Chybějící data
Karel Rečka
08. 12. 2022
1 Balíčky
library(mice)
library(tidyverse)
library(naniar)
library(simputation)
library(VIM)2 Základní funkce
Nejdříve si vytvoříme cvičný vektor
x <- c(NA, NaN, Inf, ".", "missing")
x># [1] NA "NaN" "Inf" "." "missing"Základní funkce z base R pro identifikaci chybějících hodnot je
is.na(). Vrací logické vektor stejné délky jako vektor
vstupní, který obsahuje hodnoty TRUE, je-li daný prvek
NA, jinak FALSE. Balíček
naniar má ale vlastní užitečné funkce, které
vracejí jen jednu logickou hodnotu:
any_na(): je aspoň jedna hodnotaNA? Funkce vrací jen jednu logickou hodnotu.all_na(x)neboall_miss(): Jsou všechny hodnotyNA?all_complete(): Jsou všechny hodnoty kompletní (ani jedna neníNA)?
# Které hodnoty jsou NA?
is.na(x)># [1] TRUE FALSE FALSE FALSE FALSE# Je aspoň jedna hodnota NA?
any_na(x)># [1] TRUE# Jsou všechny hodnoty NA?
all_na(x)># [1] FALSEall_miss(x)># [1] FALSE# Jsou všechny hodnoty kompletní?
all_complete(x)># [1] FALSE3 Import dat
Použitý dataset pochází z výkumu, který zkoumal souběžnou validitu nového testu kognitivních schopností (LEM) u romských a neromských dětí
df <- haven::read_sav("data/lem.sav") %>%
mutate(
across(c(sch_type, school, class, sex, ethnic,
starts_with("edu_")), as_factor),
b_total = b1 + b2
)Zde je přehled použitého datasetu:
glimpse(df)># Rows: 380
># Columns: 22
># $ id <chr> "0001", "0002", "0003", "0004", "0005", "0006", "0007", "000…
># $ school <fct> 1, 1, 6, 6, 12, 12, 12, 12, 3, 3, 3, 1, 4, 4, 11, 12, 12, 1,…
># $ class <fct> 1, 2, 13, 13, 25, 26, 26, 26, 7, 7, 7, 2, 11, 11, 24, 26, 26…
># $ sch_type <fct> základní, základní, základní, základní, základní, základní, …
># $ sex <fct> chlapec, chlapec, dívka, chlapec, chlapec, chlapec, dívka, c…
># $ ethnic <fct> nerom, nerom, nerom, nerom, nerom, nerom, nerom, nerom, nero…
># $ raven <dbl> NA, NA, NA, NA, 25, 28, 29, NA, NA, NA, NA, NA, NA, NA, NA, …
># $ b1 <dbl> NA, NA, NA, NA, 25, 20, 24, NA, 18, 20, 18, NA, NA, NA, NA, …
># $ b2 <dbl> 16, 19, 24, 26, 24, 20, 24, 20, 14, 13, 10, 12, 24, 23, 16, …
># $ lem_total <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
># $ lem1 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
># $ lem2 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
># $ lem3 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
># $ lem4 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
># $ lem5 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
># $ lem6 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
># $ reading <dbl> 20, 26, 32, 39, 29, 31, 35, 30, 23, 30, 12, NA, NA, NA, NA, …
># $ maths <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
># $ edu_f <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
># $ edu_m <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
># $ edu_a <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
># $ b_total <dbl> NA, NA, NA, NA, 49, 40, 48, NA, 32, 33, 28, NA, NA, NA, NA, …summary(df)># id school class sch_type sex
># Length:380 4 : 65 13 : 24 základní:327 dívka :193
># Class :character 6 : 58 26 : 24 zvláštní: 53 chlapec:187
># Mode :character 7 : 41 2 : 21
># 12 : 41 9 : 21
># 1 : 39 12 : 21
># 3 : 32 6 : 20
># (Other):104 (Other):249
># ethnic raven b1 b2 lem_total
># nerom:294 Min. : 6.00 Min. : 8.00 Min. : 3.00 Min. : 36
># Rom : 86 1st Qu.:21.00 1st Qu.:21.00 1st Qu.:15.00 1st Qu.:144
># Median :28.00 Median :23.00 Median :21.00 Median :170
># Mean :25.84 Mean :22.19 Mean :19.42 Mean :163
># 3rd Qu.:31.00 3rd Qu.:24.00 3rd Qu.:25.00 3rd Qu.:185
># Max. :34.00 Max. :25.00 Max. :29.00 Max. :224
># NA's :61 NA's :51 NA's :44 NA's :51
># lem1 lem2 lem3 lem4
># Min. :13.00 Min. :10.00 Min. : 1.00 Min. : 1.00
># 1st Qu.:35.00 1st Qu.:26.00 1st Qu.:11.00 1st Qu.:15.00
># Median :43.00 Median :31.00 Median :15.00 Median :21.00
># Mean :41.28 Mean :30.27 Mean :14.78 Mean :18.49
># 3rd Qu.:49.00 3rd Qu.:35.00 3rd Qu.:18.00 3rd Qu.:23.00
># Max. :54.00 Max. :42.00 Max. :28.00 Max. :24.00
># NA's :53 NA's :54 NA's :59 NA's :66
># lem5 lem6 reading maths
># Min. : 9.00 Min. : 1.00 Min. : 2.00 Min. : 2.00
># 1st Qu.:36.00 1st Qu.:16.00 1st Qu.:22.00 1st Qu.:23.00
># Median :40.00 Median :22.00 Median :29.00 Median :31.00
># Mean :38.98 Mean :20.39 Mean :27.65 Mean :29.43
># 3rd Qu.:44.00 3rd Qu.:26.00 3rd Qu.:35.00 3rd Qu.:36.00
># Max. :51.00 Max. :30.00 Max. :39.00 Max. :46.00
># NA's :56 NA's :55 NA's :46 NA's :78
># edu_f edu_m edu_a b_total
># základní :69 základní :103 1 :76 Min. :17.00
># vyučen :40 vyučen : 15 2 :33 1st Qu.:36.00
># střední bez maturity:55 střední bez maturity: 43 3 :63 Median :44.00
># střední s maturitou :73 střední s maturitou :106 4 :66 Mean :41.82
># nadstavba : 1 nadstavba : 0 5 :50 3rd Qu.:49.00
># vysokoškolské :79 vysokoškolské : 67 6 :48 Max. :53.00
># NA's :63 NA's : 46 NA's:44 NA's :66# Funkce overview pro extrakci názvů sloupců, jejich labelů i labelů hodnot
# (pokud jsou # k dispozici)
overview <- function(data) {
nms <- colnames(data)
atr <- data %>%
map(attributes)
get_var_lab <- function(x) {
if (is.null(x$label)) {
""
} else {
x$label
}
}
var_label <- atr %>% map_chr(get_var_lab)
get_value_lab <- function(x) {
if (is.null(x$labels)) {
""
} else {
str_c(x$labels, " ", names(x$labels), collapse = ". ")
}
}
value_labels <- atr %>% map_chr(get_value_lab)
tibble(
var_names = nms,
var_labels = var_label,
value_labels = value_labels
)
}
overview(df) %>%
print(n = Inf)># # A tibble: 22 × 3
># var_names var_labels value_labels
># <chr> <chr> <chr>
># 1 id "" ""
># 2 school "" ""
># 3 class "" ""
># 4 sch_type "typ školy" ""
># 5 sex "pohlaví" ""
># 6 ethnic "etnicita" ""
># 7 raven "Raven" ""
># 8 b1 "Boehm1" ""
># 9 b2 "Boehm2" ""
># 10 lem_total "LEM Celkový skór" ""
># 11 lem1 "LEM třídění" ""
># 12 lem2 "LEM asociace rozpoznávání" ""
># 13 lem3 "LEM asociace pojmenování" ""
># 14 lem4 "LEM číselné řady" ""
># 15 lem5 "LEM opakování slabik" ""
># 16 lem6 "LEM analogie" ""
># 17 reading "čtení (školní test 0-39 bodů)" ""
># 18 maths "počítání (školní test 0-50 bodů)" ""
># 19 edu_f "Father's education" ""
># 20 edu_m "Mother's education" ""
># 21 edu_a "" ""
># 22 b_total "Boehm1" ""Další funkce z balíčku nanir
n_miss()zjistí celkový počet chybějících hodnot napříč proměnnými.- funkce
prop_miss()aprop_complete()zjistí podíl chybějících a nechybějících hodnot
n_miss(df)># [1] 893prop_miss(df)># [1] 0.1068182prop_complete(df)># [1] 0.89318184 Tabulky
Funkce miss_var_summary() zjistí počet a podíl
chybějících hodnot v jednotlivých proměnných a funkce
miss_var_table() ukáže, kolika proměnným chybí kolik
hodnot.
# Počet a podíl chybějících hodnot v jednotlivých proměnných
miss_var_summary(df)># # A tibble: 22 × 3
># variable n_miss pct_miss
># <chr> <int> <dbl>
># 1 maths 78 20.5
># 2 lem4 66 17.4
># 3 b_total 66 17.4
># 4 edu_f 63 16.6
># 5 raven 61 16.1
># 6 lem3 59 15.5
># 7 lem5 56 14.7
># 8 lem6 55 14.5
># 9 lem2 54 14.2
># 10 lem1 53 13.9
># # … with 12 more rows# Kolika proměnným chybí kolik hodnot
miss_var_table(df)># # A tibble: 13 × 3
># n_miss_in_var n_vars pct_vars
># <int> <int> <dbl>
># 1 0 6 27.3
># 2 44 2 9.09
># 3 46 2 9.09
># 4 51 2 9.09
># 5 53 1 4.55
># 6 54 1 4.55
># 7 55 1 4.55
># 8 56 1 4.55
># 9 59 1 4.55
># 10 61 1 4.55
># 11 63 1 4.55
># 12 66 2 9.09
># 13 78 1 4.55Otázka: Ve kterých proměnných chybí nejvíce hodnot?
Pomocí group_by můžeme srovnat různé skupiny, např. Romy a neromy, z hlediska podílu chybějících hodnot v jednotlivých proměnných:
df %>%
group_by(ethnic) %>% # Rozdělit podle ethnic
miss_var_summary() %>% # Vypočíst počet a podíl chybějících hodnot
pivot_wider(names_from = ethnic, values_from = c(n_miss, pct_miss)) %>% # Převést na wider formát
mutate(diff = pct_miss_nerom - pct_miss_Rom) %>% # Vypočíst rozdíl v podílu NA
arrange(-abs(diff)) # seřadit podle vypočteného rozdílu># # A tibble: 21 × 6
># variable n_miss_nerom n_miss_Rom pct_miss_nerom pct_miss_Rom diff
># <chr> <int> <int> <dbl> <dbl> <dbl>
># 1 edu_f 42 21 14.3 24.4 -10.1
># 2 lem5 50 6 17.0 6.98 10.0
># 3 lem_total 46 5 15.6 5.81 9.83
># 4 lem6 49 6 16.7 6.98 9.69
># 5 lem3 52 7 17.7 8.14 9.55
># 6 lem2 48 6 16.3 6.98 9.35
># 7 lem1 47 6 16.0 6.98 9.01
># 8 maths 65 13 22.1 15.1 6.99
># 9 b1 44 7 15.0 8.14 6.83
># 10 raven 51 10 17.3 11.6 5.72
># # … with 11 more rowsOtázka: Ve kterých se obě skupiny nejvíce liší, co se týče podílu chybějících hodnot? Jaké vás napadá vysvětlení?
Funkce miss_case_summary() ukáže počet NAs pro každého
respondenta a funkce miss_case_table() kolika respondentům
chybí kolik hodnot
# počet NAs pro každého respondenta
miss_case_summary(df)># # A tibble: 380 × 3
># case n_miss pct_miss
># <int> <int> <dbl>
># 1 12 15 68.2
># 2 13 15 68.2
># 3 14 15 68.2
># 4 15 15 68.2
># 5 16 15 68.2
># 6 17 15 68.2
># 7 1 14 63.6
># 8 2 14 63.6
># 9 3 14 63.6
># 10 4 14 63.6
># # … with 370 more rows# kolika respondentům chybí kolik hodnot
miss_case_table(df)># # A tibble: 16 × 3
># n_miss_in_case n_cases pct_cases
># <int> <int> <dbl>
># 1 0 200 52.6
># 2 1 56 14.7
># 3 2 5 1.32
># 4 3 24 6.32
># 5 4 12 3.16
># 6 5 6 1.58
># 7 6 19 5
># 8 7 29 7.63
># 9 8 2 0.526
># 10 9 1 0.263
># 11 10 3 0.789
># 12 11 5 1.32
># 13 12 3 0.789
># 14 13 2 0.526
># 15 14 7 1.84
># 16 15 6 1.58# Případy s největším počtem NA
df %>%
slice(12:17)># # A tibble: 6 × 22
># id school class sch_type sex ethnic raven b1 b2 lem_t…¹ lem1 lem2
># <chr> <fct> <fct> <fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
># 1 0012 1 2 základní chla… nerom NA NA 12 NA NA NA
># 2 0020 4 11 základní dívka nerom NA NA 24 NA NA NA
># 3 0021 4 11 základní chla… nerom NA NA 23 NA NA NA
># 4 0022 11 24 zvláštní dívka nerom NA NA 16 NA NA NA
># 5 0025 12 26 základní dívka nerom NA NA 25 NA NA NA
># 6 0027 12 26 základní chla… nerom NA NA 25 NA NA NA
># # … with 10 more variables: lem3 <dbl>, lem4 <dbl>, lem5 <dbl>, lem6 <dbl>,
># # reading <dbl>, maths <dbl>, edu_f <fct>, edu_m <fct>, edu_a <fct>,
># # b_total <dbl>, and abbreviated variable name ¹lem_totalOtázka: V čem se shodují žáci s nejvtším počtem chybějících hodnot (žáci na řádcích 12 až 17)
Funkcí add_n_miss() můžete přidat do datasetu sloupec s
počtem chybějících hodnot v daném řádku. Defaultně počítá pro každý
řádek počet chybějících hodnot ve všech sloupcích, ale můžeme
specifikovat konkrétní sloupce stejnými způsoby jako ve funkci
select().
df <- df %>%
add_n_miss()
df %>%
select(n_miss_all, everything()) %>%
arrange(desc(n_miss_all))># # A tibble: 380 × 23
># n_miss_all id school class sch_type sex ethnic raven b1 b2 lem_t…¹
># <int> <chr> <fct> <fct> <fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl>
># 1 15 0012 1 2 základní chla… nerom NA NA 12 NA
># 2 15 0020 4 11 základní dívka nerom NA NA 24 NA
># 3 15 0021 4 11 základní chla… nerom NA NA 23 NA
># 4 15 0022 11 24 zvláštní dívka nerom NA NA 16 NA
># 5 15 0025 12 26 základní dívka nerom NA NA 25 NA
># 6 15 0027 12 26 základní chla… nerom NA NA 25 NA
># 7 14 0001 1 1 základní chla… nerom NA NA 16 NA
># 8 14 0002 1 2 základní chla… nerom NA NA 19 NA
># 9 14 0003 6 13 základní dívka nerom NA NA 24 NA
># 10 14 0004 6 13 základní chla… nerom NA NA 26 NA
># # … with 370 more rows, 12 more variables: lem1 <dbl>, lem2 <dbl>, lem3 <dbl>,
># # lem4 <dbl>, lem5 <dbl>, lem6 <dbl>, reading <dbl>, maths <dbl>,
># # edu_f <fct>, edu_m <fct>, edu_a <fct>, b_total <dbl>, and abbreviated
># # variable name ¹lem_total5 Grafy
Prozkoumat chybějící hodnoty můžeme i graficky. Funkce
gg_miss_var() zobrazí počet nebo podíl chybějících hodnot v
jendotlivých proměnných. Pomocí argumentu facet můžeme
specifikovat, zda se má graf rozdělit do panelů podle úrovní nějaké
proměnné a argumentem show_pct specifikujeme, zda chceme
zobrazit relativní četnosti (TRUE), nebo absolutní
(FALSE)
gg_miss_var(df)># Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
># use `guide = "none"` instead.
gg_miss_var(df, facet = ethnic, show_pct = TRUE)># Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
># use `guide = "none"` instead.
gg_miss_var(df, facet = sch_type, show_pct = TRUE)># Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
># use `guide = "none"` instead.
gg_miss_var(df, facet = sex, show_pct = TRUE)># Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
># use `guide = "none"` instead.
Otázka: V jakých proměnných se z hlediska podílu chybějících hodnot nejvíce liší Romové vs neromové a žáci ze zvláštní vs základní školy?
Graf jednotlivých chybějících hodnot pro jednotlivé případy můžeme
získat pomocí funkce vis_miss(). Argumentem
cluster = TRUE specifikujeme, že chceme případy seřadit
pomocí hierarchické clusterové analýzy dle vzorce chybějících .
vis_miss(df, cluster = TRUE) +
theme(axis.text.x = element_text(angle = 90)) # natočit labely, aby se nepřekrývaly># Warning: `gather_()` was deprecated in tidyr 1.2.0.
># Please use `gather()` instead.
Nejčetnější vzorce chybějících hodnot můžeme získat pomocí funkce
gg_miss_upset().
gg_miss_upset(df)
gg_miss_upset(df,
nsets = 10, # maximální počet proměnných
nintersects = 10) # maximální počet vzorců
Otázka: Hodnoty kterých proměnných mají tendenci chybět společně?
Srovnat více skupin v podílu chybějících hodnot umožňuje také heatmapa.
gg_miss_fct(df, fct = ethnic)
gg_miss_fct(df, fct = school)
gg_miss_fct(df, fct = sch_type)
Otázka: Která škola school (kvůli
anonymizaci známe jen její numerické id) se od ostatních nejvíce
odlišuje z hlediska podílu chybějících hodnot a o které proměnné se
jedná?
Scatterploty s chybějícími hodnotami můžeme získat pomocí funkce
geom_miss_point(). Chybějící hodnoty budou umístěny na
okraje grafu mimo pozorované rozmezí hodnot. Tak můžeme prozkoumat, zda
chybění v nějaké spojité proměnné závisí na hodnotě v jiné spojité
proměnné.
ggplot(df, aes(raven, lem_total)) +
geom_miss_point(position = position_jitter(width = 0.5))
ggplot(df, aes(maths, lem_total)) +
geom_miss_point(position = position_jitter(width = 0.5))
Otázka: Vypadá to, že by chybění v
lem_total záviselo na hodnotách raven nebo
maths nebo naopak chybění v raven nebo
maths záviselo na hodnotě lem_total?
6 Nabular data
Tzv. shadow matrix obsahuje stejný počet řádků a sloupců
jako původní dataset, ale místo původních hodnot obsahují pouze
informace o tom, zda je nějaká hodnota chybějící (NA), nebo
není (!NA). Nabular data pak vzniknout připojením
shadow matrix k původnímu datasetu jako nové sloupce.
# Pouze shadow matrix
as_shadow(df) ># # A tibble: 380 × 23
># id_NA school_NA class_NA sch_typ…¹ sex_NA ethni…² raven…³ b1_NA b2_NA lem_t…⁴
># <fct> <fct> <fct> <fct> <fct> <fct> <fct> <fct> <fct> <fct>
># 1 !NA !NA !NA !NA !NA !NA NA NA !NA NA
># 2 !NA !NA !NA !NA !NA !NA NA NA !NA NA
># 3 !NA !NA !NA !NA !NA !NA NA NA !NA NA
># 4 !NA !NA !NA !NA !NA !NA NA NA !NA NA
># 5 !NA !NA !NA !NA !NA !NA !NA !NA !NA NA
># 6 !NA !NA !NA !NA !NA !NA !NA !NA !NA NA
># 7 !NA !NA !NA !NA !NA !NA !NA !NA !NA NA
># 8 !NA !NA !NA !NA !NA !NA NA NA !NA NA
># 9 !NA !NA !NA !NA !NA !NA NA !NA !NA NA
># 10 !NA !NA !NA !NA !NA !NA NA !NA !NA NA
># # … with 370 more rows, 13 more variables: lem1_NA <fct>, lem2_NA <fct>,
># # lem3_NA <fct>, lem4_NA <fct>, lem5_NA <fct>, lem6_NA <fct>,
># # reading_NA <fct>, maths_NA <fct>, edu_f_NA <fct>, edu_m_NA <fct>,
># # edu_a_NA <fct>, b_total_NA <fct>, n_miss_all_NA <fct>, and abbreviated
># # variable names ¹sch_type_NA, ²ethnic_NA, ³raven_NA, ⁴lem_total_NA# argumentem (only_miss = TRUE) vyřadíme ze shadow matrix proměnné bez chybějících hodnot
df_nabular <- bind_shadow(df, only_miss = TRUE)
glimpse(df_nabular)># Rows: 380
># Columns: 39
># $ id <chr> "0001", "0002", "0003", "0004", "0005", "0006", "0007", "…
># $ school <fct> 1, 1, 6, 6, 12, 12, 12, 12, 3, 3, 3, 1, 4, 4, 11, 12, 12,…
># $ class <fct> 1, 2, 13, 13, 25, 26, 26, 26, 7, 7, 7, 2, 11, 11, 24, 26,…
># $ sch_type <fct> základní, základní, základní, základní, základní, základn…
># $ sex <fct> chlapec, chlapec, dívka, chlapec, chlapec, chlapec, dívka…
># $ ethnic <fct> nerom, nerom, nerom, nerom, nerom, nerom, nerom, nerom, n…
># $ raven <dbl> NA, NA, NA, NA, 25, 28, 29, NA, NA, NA, NA, NA, NA, NA, N…
># $ b1 <dbl> NA, NA, NA, NA, 25, 20, 24, NA, 18, 20, 18, NA, NA, NA, N…
># $ b2 <dbl> 16, 19, 24, 26, 24, 20, 24, 20, 14, 13, 10, 12, 24, 23, 1…
># $ lem_total <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ lem1 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ lem2 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ lem3 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ lem4 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ lem5 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ lem6 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ reading <dbl> 20, 26, 32, 39, 29, 31, 35, 30, 23, 30, 12, NA, NA, NA, N…
># $ maths <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ edu_f <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ edu_m <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ edu_a <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ b_total <dbl> NA, NA, NA, NA, 49, 40, 48, NA, 32, 33, 28, NA, NA, NA, N…
># $ n_miss_all <int> 14, 14, 14, 14, 11, 11, 11, 14, 12, 12, 12, 15, 15, 15, 1…
># $ raven_NA <fct> NA, NA, NA, NA, !NA, !NA, !NA, NA, NA, NA, NA, NA, NA, NA…
># $ b1_NA <fct> NA, NA, NA, NA, !NA, !NA, !NA, NA, !NA, !NA, !NA, NA, NA,…
># $ b2_NA <fct> !NA, !NA, !NA, !NA, !NA, !NA, !NA, !NA, !NA, !NA, !NA, !N…
># $ lem_total_NA <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ lem1_NA <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ lem2_NA <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ lem3_NA <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ lem4_NA <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ lem5_NA <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ lem6_NA <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ reading_NA <fct> !NA, !NA, !NA, !NA, !NA, !NA, !NA, !NA, !NA, !NA, !NA, NA…
># $ maths_NA <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ edu_f_NA <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ edu_m_NA <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ edu_a_NA <fct> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
># $ b_total_NA <fct> NA, NA, NA, NA, !NA, !NA, !NA, NA, !NA, !NA, !NA, NA, NA,…Shadow matrix pak můžeme využít pro porovnání skupin vymezených podle chybějících dat v nějaké proměnné, abychom zjistili, zda chybění nějak systematicky nesouvisí s jinou proměnnou, jejíž hodnoty známe.
my_summarise <- function(data, group_vars, summary_vars) {
data %>%
group_by(across({{group_vars}})) %>%
summarise(across({{ summary_vars }}, ~sum(!is.na(.x)), .names = "n_{.col}"),
across({{ summary_vars }}, mean, na.rm = TRUE, .names = "M_{.col}"),
across({{ summary_vars }}, sd, na.rm = TRUE, .names = "SD_{.col}"))
}
my_summarise(df_nabular,
group_vars = raven_NA,
summary_vars = lem_total)># # A tibble: 2 × 4
># raven_NA n_lem_total M_lem_total SD_lem_total
># <fct> <int> <dbl> <dbl>
># 1 !NA 286 162. 32.3
># 2 NA 43 168. 27.6my_summarise(df_nabular,
group_vars = maths_NA,
summary_vars = lem_total)># # A tibble: 2 × 4
># maths_NA n_lem_total M_lem_total SD_lem_total
># <fct> <int> <dbl> <dbl>
># 1 !NA 274 161. 32.2
># 2 NA 55 171. 28.4my_summarise(df_nabular,
group_vars = edu_f_NA,
summary_vars = lem_total)># # A tibble: 2 × 4
># edu_f_NA n_lem_total M_lem_total SD_lem_total
># <fct> <int> <dbl> <dbl>
># 1 !NA 288 165. 31.1
># 2 NA 41 149. 33.3Otázka: Na chybějí ve které proměnné závisí průměrný
skór lem_total nejsilněji? Jak to může souviset s
následujícícm vztahem?
df_nabular %>%
count(sch_type, edu_f_NA) %>%
group_by(sch_type) %>%
mutate(p = n/sum(n))># # A tibble: 4 × 4
># # Groups: sch_type [2]
># sch_type edu_f_NA n p
># <fct> <fct> <int> <dbl>
># 1 základní !NA 278 0.850
># 2 základní NA 49 0.150
># 3 zvláštní !NA 39 0.736
># 4 zvláštní NA 14 0.264Samozřejmě shadow matrix můžeme využít i k tvorbě grafů. Zde je např.
graf hustoty pravděpodobnosti lem_total podle chybění v
reading (barevně) a podle ethnic (fazety)
df_nabular %>%
ggplot(aes(lem_total, color = reading_NA)) +
geom_density() +
facet_wrap(~ethnic)># Warning: Removed 51 rows containing non-finite values (stat_density).
Otázka: Liší se nějak výrazně rozdělení
lem_total v závislosti na chybějí v proměnné
reading u Romů nebo neromů?
7 Imputace dat
V analýze se zaměříme primárně na to, jaký efekt má škola (sch_type) na zvládnutí základních prvků gramotnosti - zaměříme se na čtení (reading).
Čtení bylo testováno školním testem. Předpokládáme tedy, že v základní škole se dítě naučí více než ve zvláštní. V duchu projektu budeme zvažovat, že efekt typu školy (tj. jak dobře dokáže škola základy žáka naučit) se liší podle etnicity (ethnic).
I když bychom si to úplně nepřáli, je možné, že roli bude hrát také pohlaví žáka (sex). Konečně, musíme zohlednit také intelekt (lem), protože schopnost čtení je jím samozřejmě ovlivněna.
Jednoduchá forma imputace dat, která se hodí pro exploraci vzorců
chybějících hodnot, je doplnění chybějícíh hodnot tak, aby se doplněné
hodnoty pohybovaly mimo rozmezí chybějících dat a šlo je snadno
rozpoznat. To uděláme pomocí funkce impute_below_all() a
poté pmocí add_label_shadow() doplníme sloupec
any_missing jako jednoduchý indikátor toho, zda se na danm
řádku vyskytují nějaké chybějí hodnoty.
df_imputed <- df_nabular %>% # připojení shadow matrix
impute_below_all() %>% # imputace chybějících hodnot mimo rozmezí
add_label_shadow() # přidat proměnnou "any_missing"Pak můžeme např. prozkoumat, jestli chybějí v reading
závisí hodnotě lem_total nebo jestli chybění v
lem_total závisí na chybění v reading a také
zda respondenti s aspoň jednou chybějící hodnotou se liší od ostatních
třeba nějak specifickou kombinací hodnot v reading a
lem_total.
df_imputed %>%
ggplot(aes(lem_total, reading, color = any_missing)) +
geom_point()
Otázka: Lze v grafu výše pozorovat něco z toho, co bylo popsáno v minulém odstavci?
7.1 Jednodušší metody imputace
7.1.1 Imputace průměru
Nejjednoduší metodou imputace je jako chybějící hodnoty imputovat průměr dané proměnné (vypočtený na základě nechybějících hodnot). Pro kategorické proměnné je imputace modu obdobou imputace průměru u numerických proměnných.
Takto bychom mohli chybějící hodnoty v numerických proměnných nahradit průměry:
df_imp_mean <- df %>%
bind_shadow() %>%
add_label_shadow() %>%
mutate(
across(where(is.numeric), impute_mean)
)Zásadní problém je ten, že imputované hodnoty jsou vlastně konstantou, což samozřejmě může zásadně ovlivnit univariančí i multivariační rozdělení našich dat (viz níže)
ggplot(df_imp_mean, aes(lem_total, reading, color = any_missing)) +
geom_point()
ggplot(df_imp_mean, aes(reading, fill = reading_NA)) +
geom_histogram()># `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Kvůli nižší variabilitě se sníží vztahy s ostatními proměnnými.
x <- df %>%
select(where(is.numeric)) %>%
cor(use = "pairwise") %>%
round(digits = 2) %>%
.[, c("reading", "lem_total")]
y <- df_imp_mean %>%
select(where(is.numeric)) %>%
cor() %>%
round(digits = 2)%>%
.[, c("reading", "lem_total")]
cbind(x, y)># reading lem_total reading lem_total
># raven 0.68 0.64 0.64 0.60
># b1 0.71 0.60 0.68 0.51
># b2 0.80 0.68 0.77 0.61
># lem_total 0.68 1.00 0.61 1.00
># lem1 0.58 0.71 0.52 0.70
># lem2 0.36 0.62 0.32 0.62
># lem3 0.45 0.69 0.39 0.66
># lem4 0.63 0.78 0.53 0.71
># lem5 0.55 0.74 0.49 0.71
># lem6 0.69 0.80 0.60 0.77
># reading 1.00 0.68 1.00 0.61
># maths 0.71 0.58 0.66 0.51
># b_total 0.81 0.70 0.77 0.62
># n_miss_all -0.08 -0.07 -0.06 -0.037.1.2 Hot-deck imputace
Dalším jednoduchým typem imputace je hot-deck. Principem je jednoduše
nahradit chybějící hodnoty poslední nechybějící hodnotou. Pro doplnění
“validnějších” hodnot má smysl stanovit si tzv. “doménovou” proměnnou
nebo proměnné, aby se chybějící hodnoty doplňovaly pouze v rámci téže
“domény”. Pokud si jako domain_var stanovíme typ školy a
etnikum, budou se případům s chybějícími hodnotami imputovat chybějící
hodnoty pouze od osob na tomtéž typu školy a téhož etnika. Dále si
můžeme vybrat proměnnou, podle které se má před imputací dataset
seřadit. Měla by to být proměnná, která co nejtěsněji souvisí s
proměnnými, jejichž chybějící hodnoty chceme imputovat, protože pak
imputujeme plauzibilnější hodnoty. Kdybychom např. měřili tělesnou výšku
a hmotnost a měli bychom chybějící hodnoty v hmotnosti a za ně chtěli
imputovat validní hodnoty, dávalo by smysl stanovit si jako “doménu”
pohlaví, aby se ženám s chybějícími hodnotami imputovaly pouze hodnoty
od jiných žen a mužům pouze hodnoty od jiných mužů, a navíc si seřadit
data podle výšky, aby se případům s chybějícími hodnotami v tělesné
hmotnosti imputovaly hodnoty tělesné hmotnosti od případů s podobnou
tělesnou výškou.
Níže nahrazujeme chybějící hodnoty v proměnných
lem_total a reading pomocí hot-deck imputace,
přičemž typ školy a etnikum stanovujeme jako doménové proměnné a
seřazujeme dataset podle zscore (protože každý žák vyplnit
aspoň jeden kognitivní subtest, můžeme všechny subtesty převést na
z-skóry, vypočíst z nich průměr a podle toho dataset seřadit).
df <- df %>%
mutate(zscore = select(., raven, b1, b2, lem1:lem6) %>%
mutate(across(.fns = scale)) %>%
rowMeans(na.rm = TRUE))
df_imp_hd <- df %>%
VIM::hotdeck(variable = c("lem_total", "reading"), # do kterých proměnných imputovat data
domain_var = c("sch_type", "ethnic"), # vymezení "domén"
ord_var = "zscore") # seřadit podleK samotné imputaci byla použita funkce marginplot() z
balíčku VIM, která kromě imputace samotné přidá nové sloupce, v našem
případě lem_total_imp a reading_imp s
logickými hodnotami indikujícími, zda daná hodnota byla byla imputována
(protože chyběla), nebo nikoli (protože nechyběla).
Pak se opět můžeme podívat na vztah obou proměnných po imputaci.
Můžeme nyní zkusit použít funkci marginplot() z balíčku
VIM.
# Draw a margin plot
df_imp_hd %>%
select(lem_total, reading, ends_with("_imp")) %>%
VIM::marginplot(delimiter = "_imp")
Otázka: Ovlivnila v tomto případě hotdeck imputace zásadně podobu univariačního a bivariačního rozdělení obou proměnných?
7.1.3 kNN imputace
Další jednodušší metodou imputace je kNN (k-Nearest Neighbours)
imputace. Nejprve si vybereme soubor proměnných, na základě nich se
vypočte matice vzdáleností pro jednotlivé případy (tzv. Gower Distance).
Čím si jsou nějaké případy podobnější z hlediska vybranných proměnných,
tím menší “vzdálenost” mezi nimi je. Pak si stanovíme počet neighbour
čili počet nejpodobnějších případů. Chybějící hodnoty určitého případu
se pak budou doplňovat na základě jeho nejbližších “sousedů” s validními
hodnotami.
Níže vypočteme matici vzdáleností na základě typu školy, etnika a skórů
v Ravenovi a Boehmově testu, pro každý případ s chybějícími hodnotami v
proměnných lem_total nebo reading se najde pět
nejbližších sousedů a místo chybějících hodnot se doplní vážený průměr z
těchto pěti nejbližších sousedů (vážený vzdáleností tak, aby větší váhu
měly co nejpodobnější případy).
df_imp_knn <- df %>%
VIM::kNN(variable = c("lem_total", "reading"),
dist_var = c("sch_type", "ethnic",
"raven", "b_total"),
k = 5, # počet neighbors
numFun = weighted.mean,
weightDist = TRUE)S tímto nastavením by metoda kNN vedla k vytvoření artefaktru: až
příliš mnoho hodnot imputovaných do reading se soustřeďuje
okolo mediánu.
df_imp_knn %>%
select(lem_total, reading, ends_with("_imp")) %>%
marginplot(delimiter = "_imp")
7.1.4 Lineárně-regresní imputace
Jako poslední “nerefinovanou” metodu imputace si představíme imputaci pomocí lineární regrese.
Nejprve si napíšeme funkci pro vápočet průměrné absolutní změny hodnot pro posouzení stability imputace.
mapc <- function(a, b) {
mean(abs(b - a) / a, na.rm = TRUE)
}Pak si vytvoříme dva zatím prázdné vektory pro zaznamenání průměrné absolutní změny:
diff_lem_total <- double(10)
diff_reading <- double(10)Dále si vytvoříme vektory pro zaznamenání toho, které z původních
hodnot v lem_total a reading byly
chybějící:
na_lem_total <- is.na(df$lem_total)
na_reading <- is.na(df$reading)Aby se nám regresní model odhadl se všemi případy, musíme prozatimně zaplnit chybějící hodnoty pomocí nějaké jednoduché metody, např. imputací průměru. impute_s # Inicializace chybějících hodnot pomocí prosté hotdeck imputece
df_imp_lm <- df %>%
mutate(lem_total_imp = is.na(lem_total),
reading_imp = is.na(reading),
across(c(lem_total, reading), impute_mean))Pak provedena pomocí for loop iterativní imputaci (celkem 10krát, ale mohli bychom i vícekrát). Každá iterace proběhne v těchto krocích:
- Uloží data z předchozí iterace.
- Resetuje chybějící hodnoty v proměnné
reading(ty, které byly původněNA, budou znovuNA). - Do proměnné
readingimputuje predikované hodnoty s lineární regrese s pohlavím, typem školy, etnicitou (plus interakcí mezi typem školy a etnicitou) alem_totaljako prediktory. - Resetuje chybějící hodnoty v proměnné
lem_total(ty, které byly původněNA, budou znovuNA). - Do proměnné
lem_totalimputuje predikované hodnoty s lineární regrese s pohlavím, typem školy, etnicitou (plus interakcí mezi typem školy a etnicitou) areadingjako prediktory. - Pro každou iteraci uloží průměrnou absolutní změnu hodnot proměnných
lem_totalareading.
Proč vlastně používáme iterativní odhad? Model pro imputaci
reading obsahuje jako prediktor lem_total a
pokud se změní hodnoty v lem_total, změní se i koeficienty
modelu a taktéž predikované hodnoty. Podobně Model pro imputaci
reading obsahuje jako prediktor reading a
pokud se změní hodnoty v reading změní se i koeficienty
modelu a taktéž predikované hodnoty určené pro imputaci. Postupně by ale
měla být změna čím dále menší a imputované hodnoty i koeficienty modelu
by se měly ustálit.
for (i in 1:10) {
# Uložení předchozí iterace
prev_iter <- df_imp_lm
# Nastavit na NA hodnoty, které byly původně NA
df_imp_lm$reading[na_reading] <- NA
# Imputovat reading pomocí regrese (predikované hodnoty)
df_imp_lm <- df_imp_lm %>%
impute_lm(reading ~ sex + sch_type*ethnic + lem_total)
# Nastavit na NA hodnoty, které byly původně NA
df_imp_lm$lem_total[na_lem_total] <- NA
# Imputovat lem_total pomocí regrese (predikované hodnoty)
df_imp_lm <- df_imp_lm %>%
impute_lm(lem_total ~ sex + sch_type*ethnic + reading)
# Výpočet MAPC
diff_reading[i] <- mapc(prev_iter$reading, df_imp_lm$reading)
diff_lem_total[i] <- mapc(prev_iter$lem_total, df_imp_lm$lem_total)
}Jak vidíme níže, 10 iterací s rezervou stačilo, aby se imputované hodnoty ustálily.
plot(diff_reading, type = "b")
plot(diff_lem_total, type = "b")
Opět se můžeme podívat na univariační i bivariační rozdělení imputovaných proměnných.
df_imp_lm %>%
as.data.frame() %>%
select(lem_total, reading, ends_with("_imp")) %>%
marginplot(delimiter = "_imp")
Otázka: Co je na imputovaných hodnotách v grafu výše podezřelého?
7.2 Vícenásobné imputace pomocí balíčku mice
Slabinou všech předchozích metod, které jsme si ukázali, je to, že zanedbávají nejistotu imputovaných hodnot, ale pracují s nimi, jako bych se jednalo o skutečně pozorované hodnoty a my měli k dispozici kompletní data. To obvykle vede k podhodnocení standardních chyb, protože ty počítají pouze s výběrovou variabilitou, nikoli imputační variabilitou – a tu je třeba zohlednit.
Jedním z řešení je vícenásobná imputace a odhad této “imputační” variability na základě toho, jak se různí více imputovaných datasetů.
Ukážeme si jen jeden z možných příkladů toho, jak lze provést vícenásobnou imputaci pomocí balíčku mice.
Nejprve si definujeme modely sloužící k imutaci hodnot. Každá proměnná, která obsahuje chybějící data, může mít svůj vlastní imputační model. Doporučení pro specifikaci imputačních modelů jsou následující:
- Zahrňme všechny proměnné které budou součastní finálního modelu, jejž chceme fitovat a to včetně závislé proměnné (ta bude také použita jako predikor chybějících dat v ostatních proměnných) a odůvodněných interakcí. Jinak mohou být odhady finálního modelu zkreslené (směrem k nule, tedy podhodnocené)
- Zahrňme všechny proměnné, které souvisejí s absencí odpovědí v cílových proměnných. Můžeme tak učinit na základě teoretické úvahy a na empirickém základě. To znamená, že bychom měli zahrnout proměnné, jejichž distribuce se výrazně liší mezi skupinami s/bez chybějících dat. Obvykle to můžeme zjistit, když blíže prozkoumáme korelace různých proměnných s indikátorem chybějících hodnot (!NA vs NA) v cílovových (imputovaných) proměnných.
- Zahrňme proměnné vysvětlující významné podíl rozptylu v proměnných, které plánujem imputovat. Tím snižujeme variabilitu imputovaných hodnot I v tomto případě můžeme prozkoumat, jak korelují ostatní proměnné s cílovými (imputovanými) proměnnými.
- Vyřaďte ty proměnné z kroku 2 a 3, které mají příliš chybějících hodnot v podskupině nekompletních případů (tj. případů, které mají nějakou chybějící hodnotu v proměnné cílového modelu). Jednoduchým ukazatelem je podíl chybějících hodnot v této skupině.
Dejme tomu, že chceme imputovat hodnoty pěti numerických proměnných a imputační modely specifikujeme takto:
form_list <- c(
raven = raven ~ sex + sch_type*ethnic + b_total + reading + maths + lem_total,
b_total = b_total ~ sex + sch_type*ethnic + raven + reading + maths + lem_total,
reading = reading ~ sex + sch_type*ethnic + raven + b_total + maths + lem_total,
maths = maths ~ sex + sch_type*ethnic + raven + b_total + reading + lem_total,
lem_total = lem_total ~ sex + sch_type*ethnic + raven + b_total + reading + maths
)Když jsme si definovali imputační modely, měli bychom si zvolit metodu imputace, která závisí mj. na povaze dat. Protože máme numerická data, je jednou z vhodných metod predictive mean matchning (PMM). Tato metoda vychází opět z lineární regrese. Princip této metody je následujícíc: I pro případy, které mají v závislé proměnné chybějící hodnoty, lze vypočíst predikovanou hodnotu a najít případy s co nejpodobnější predikovanou hodnotou, které ale zároveň mají k dispozici i hodnotu pozorovanou (tj. jejichž hodnota závislé proměnné není missing). Dejme tomu, že idenfitikujeme vždy pět takovýchto nejpodobnějších případů (donors), jeden z nich náhodně vybereme a jeho pozorovanou hodnotu použijeme k imputaci. To je celý princip PMM.
Kódem níže pak generujeme deset imputovaných datasetů (proto jde o vícenásobnou imputaci)
df_imp_pmm <- mice(df_nabular,
method = "pmm", # linear regression, predicted values
formulas = form_list, # Specifikace imputačních modelů.
auxiliary = FALSE, # Nepřidávat jako prediktory žádné další proměnné
donors = 5, # počet potenciálních donors
m = 10, # Počet multiple imputations
maxit = 10, # Počet iterací pro "ustálení" imputovaných hodnot
seed = 123) # Iniciální hodnota generátoru náhodných číselVýsledný objekt df_imp_pmm je třídy mids
(multiple imputed data sets) a obsahuje jak původní data, tak 10
imputovaných datasetů. Funkcí complete() můžeme extrahovat
jednotlivé imputované datasety anebo i všechny z nich, a to v různých
formátech (podívejte se nápovědu ?mice::complete).
Tatko sloučíme všechny imputované datasety do jednoho v longer
formátu. Sloupec .imp udává číslo imputace (1 až 10).
df_imp_pmm_long <- df_imp_pmm %>%
complete(action = "long") %>%
as_tibble()
df_imp_pmm_long># # A tibble: 3,800 × 41
># .imp .id id school class sch_t…¹ sex ethnic raven b1 b2 lem_t…²
># <int> <int> <chr> <fct> <fct> <fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl>
># 1 1 1 0001 1 1 základ… chla… nerom 29 NA 16 168
># 2 1 2 0002 1 2 základ… chla… nerom 23 NA 19 174
># 3 1 3 0003 6 13 základ… dívka nerom 21 NA 24 185
># 4 1 4 0004 6 13 základ… chla… nerom 31 NA 26 204
># 5 1 5 0005 12 25 základ… chla… nerom 25 25 24 183
># 6 1 6 0006 12 26 základ… chla… nerom 28 20 20 163
># 7 1 7 0007 12 26 základ… dívka nerom 29 24 24 193
># 8 1 8 0008 12 26 základ… chla… nerom 33 NA 20 176
># 9 1 9 0009 3 7 základ… dívka nerom 20 18 14 128
># 10 1 10 0010 3 7 základ… dívka nerom 27 20 13 123
># # … with 3,790 more rows, 29 more variables: lem1 <dbl>, lem2 <dbl>,
># # lem3 <dbl>, lem4 <dbl>, lem5 <dbl>, lem6 <dbl>, reading <dbl>, maths <dbl>,
># # edu_f <fct>, edu_m <fct>, edu_a <fct>, b_total <dbl>, n_miss_all <int>,
># # raven_NA <fct>, b1_NA <fct>, b2_NA <fct>, lem_total_NA <fct>,
># # lem1_NA <fct>, lem2_NA <fct>, lem3_NA <fct>, lem4_NA <fct>, lem5_NA <fct>,
># # lem6_NA <fct>, reading_NA <fct>, maths_NA <fct>, edu_f_NA <fct>,
># # edu_m_NA <fct>, edu_a_NA <fct>, b_total_NA <fct>, and abbreviated …Můžeme si dále vytvoři sloupce indikující, zda v
lem_total nebo reading byla původně chybějící
hodnota.
df_imp_pmm_long <- df_imp_pmm_long %>%
mutate(lem_or_reading_na = lem_total_NA == "NA" | reading_NA == "NA")A opět se podívat na jejich bivariační rozdělení:
df_imp_pmm_long %>%
ggplot(aes(lem_total, reading, color = lem_or_reading_na)) +
geom_point() +
facet_wrap(~.imp, ncol = 2) +
theme(legend.position = "top")
Otázka: Poznali byste nyní (kdyby nebylo barevného odlišení), které hodnoty byly imputovány a které nikoli?
Také bychom měli ověřit, že se imputované hodnoty “ustálily”. Můžeme to ověřit nepřímo např. tak, že se podíváme, zda průměry a směrodatné odchylky imputovaných proměnných jsou stabilní – přesněji řečeno, že hlavně v posledních iteracích už jen náhodně kolísají okolo nějaké pomyslné střední hodnoty a nevykazují žádný (vzestupný/sestupný, divergentní/konvergentní) trend.
plot(df_imp_pmm)
Otázka: Ustálily se impoutované hodnoty?
Nakonec si můžeme ukázat, jak imputované datasety použít k odhadu
lineárně regresního modelu. V modelu 1 jsou prediktory úrovně čtení
pohlaví, typ školy, etnikum a interakce etnika a typu školy, v modelu 2
jsme přidali jako prediktor lem_total (měřítko obecné
inteligence). K odhadu použijeme funkci with(). Prvním
argumentem je list s imputovanými datasety, druhým argumentem je funkce
k odhadu modelu společně s definicí regresní rovnice.
# Model 1
fit_pmm1 <- with(
df_imp_pmm,
lm(reading ~ sex + sch_type*ethnic)
)
# Model 2
fit_pmm2 <- with(
df_imp_pmm,
lm(reading ~ sex + sch_type*ethnic + lem_total)
)Pomocí funkce pool() pak dostaneme souhrn s
výsledky.
pool(fit_pmm1)># Class: mipo m = 10
># term m estimate ubar b t
># 1 (Intercept) 10 31.3373219 0.2574552 0.005612977 0.2636294
># 2 sexchlapec 10 0.2003536 0.4279400 0.050033607 0.4829770
># 3 sch_typezvláštní 10 -15.7693330 1.9904730 0.301935869 2.3226024
># 4 ethnicRom 10 -9.6502099 0.8862580 0.067586638 0.9606033
># 5 sch_typezvláštní:ethnicRom 10 9.3304598 4.0349942 0.423086154 4.5003890
># dfcom df riv lambda fmi
># 1 375 356.3679 0.02398194 0.02342028 0.02885526
># 2 375 223.7909 0.12860907 0.11395360 0.12176738
># 3 375 185.1766 0.16685957 0.14299884 0.15210732
># 4 375 280.0115 0.08388675 0.07739438 0.08391430
># 5 375 239.3324 0.11533964 0.10341212 0.11081177pool(fit_pmm2)># Class: mipo m = 10
># term m estimate ubar b
># 1 (Intercept) 10 12.0925731 4.5876130194 8.985722e-01
># 2 sexchlapec 10 -0.1694591 0.3511863435 6.706865e-02
># 3 sch_typezvláštní 10 -9.8855321 2.0354009205 5.916185e-01
># 4 ethnicRom 10 -6.7492211 0.8234698952 1.250496e-01
># 5 lem_total 10 0.1111303 0.0001459347 3.045974e-05
># 6 sch_typezvláštní:ethnicRom 10 6.9162139 3.3654144455 4.467235e-01
># t dfcom df riv lambda fmi
># 1 5.5760424269 374 147.95982 0.2154561 0.1772636 0.1881637
># 2 0.4249618551 374 151.48082 0.2100751 0.1736050 0.1843040
># 3 2.6861813222 374 99.31339 0.3197308 0.2422697 0.2570817
># 4 0.9610244187 374 184.72598 0.1670426 0.1431332 0.1522621
># 5 0.0001794404 374 139.29179 0.2295938 0.1867233 0.1981544
># 6 3.8568103084 374 204.73981 0.1460135 0.1274099 0.1358107Ten obsahuje tyto sloupce:
term= prediktor;m= počet imputací;estimate= odhad koeficientu modelu;ubar= průměrný vnitroimputační rozptyl (rozptylem se zde myslí error variance čili čtverec standardní chyby);b= meziimputační rozptyl;t= celkový rozptyl souhrného odhadu (ubar + b) čili celkový čtverec standardní chyby;dfcom= reziduální stupně volnosti, pokud bychom měli kompletní data;df= reziduální stupně volnosti (pro testy regresních koeficientů);riv= odhad relativního nárůstu chybového rozptylu v důsledku chybějících dat;lambda= Odhad podílu chybového rozptylu způsobeného chybějícími daty;fmi= Rovněž odhad podílu chybového rozptylu způsobeného chybějícími daty, ale s korekcí podhodnocení (které je dáno tím, že namáme nekonečně mnoho imputovaných datasetů)
Odhady standardních chyb, testové statistiky a p-hodnot pak
dostaneme, pokud output funkce pool() pošleme do funkce
summary().
summary(pool(fit_pmm1), conf.int = TRUE, conf.level = .95)># term estimate std.error statistic df
># 1 (Intercept) 31.3373219 0.5134486 61.033030 356.3679
># 2 sexchlapec 0.2003536 0.6949655 0.288293 223.7909
># 3 sch_typezvláštní -15.7693330 1.5240087 -10.347273 185.1766
># 4 ethnicRom -9.6502099 0.9801037 -9.846111 280.0115
># 5 sch_typezvláštní:ethnicRom 9.3304598 2.1214120 4.398231 239.3324
># p.value 2.5 % 97.5 %
># 1 0.000000e+00 30.327552 32.347092
># 2 7.733893e-01 -1.169160 1.569867
># 3 0.000000e+00 -18.775985 -12.762681
># 4 0.000000e+00 -11.579517 -7.720903
># 5 1.643099e-05 5.151436 13.509483summary(pool(fit_pmm2), conf.int = TRUE, conf.level = .95)># term estimate std.error statistic df
># 1 (Intercept) 12.0925731 2.36136453 5.1210108 147.95982
># 2 sexchlapec -0.1694591 0.65189098 -0.2599501 151.48082
># 3 sch_typezvláštní -9.8855321 1.63895739 -6.0315980 99.31339
># 4 ethnicRom -6.7492211 0.98031853 -6.8847225 184.72598
># 5 lem_total 0.1111303 0.01339554 8.2960656 139.29179
># 6 sch_typezvláštní:ethnicRom 6.9162139 1.96387635 3.5217156 204.73981
># p.value 2.5 % 97.5 %
># 1 9.323636e-07 7.42621712 16.7589291
># 2 7.952553e-01 -1.45743159 1.1185134
># 3 2.793431e-08 -13.13745215 -6.6336120
># 4 8.760104e-11 -8.68328090 -4.8151612
># 5 8.304468e-14 0.08464538 0.1376151
># 6 5.285581e-04 3.04419922 10.7882286Pomocí funkce pool.r.squared() dostaneme odhad podílu
vysvětleného rozptylu.
pool.r.squared(fit_pmm1)># est lo 95 hi 95 fmi
># R^2 0.4736259 0.3906602 0.551194 0.1834618pool.r.squared(fit_pmm1, adjusted = TRUE)># est lo 95 hi 95 fmi
># adj R^2 0.4680082 0.384669 0.5461029 0.1853549pool.r.squared(fit_pmm2)># est lo 95 hi 95 fmi
># R^2 0.5714234 0.4930143 0.6418039 0.2399352pool.r.squared(fit_pmm2, adjusted = TRUE)># est lo 95 hi 95 fmi
># adj R^2 0.5656909 0.4866271 0.636799 0.2418318A srovnat fit obou modelů můžeme pomocí funkce D3(),
která provádí likelihood-ratio test.
mice::D3(fit_pmm2, fit_pmm1)># test statistic df1 df2 dfcom p.value riv
># 1 ~~ 2 56.45248 1 50.69628 374 8.728029e-10 0.378293Nakonec se můžeme podívat na to, zda by se naše závěry lišily, pokud bychom použili listwise deletion.
fit_listwise <- lm(reading ~ sex + sch_type*ethnic + lem_total, data = df)
summary(fit_listwise)>#
># Call:
># lm(formula = reading ~ sex + sch_type * ethnic + lem_total, data = df)
>#
># Residuals:
># Min 1Q Median 3Q Max
># -15.7432 -3.7130 0.4689 3.8737 17.3758
>#
># Coefficients:
># Estimate Std. Error t value Pr(>|t|)
># (Intercept) 12.59338 2.41548 5.214 3.56e-07 ***
># sexchlapec -0.21552 0.67494 -0.319 0.749719
># sch_typezvláštní -9.62594 1.57781 -6.101 3.43e-09 ***
># ethnicRom -7.21977 1.00240 -7.203 5.29e-12 ***
># lem_total 0.10969 0.01363 8.049 2.27e-14 ***
># sch_typezvláštní:ethnicRom 6.90882 2.02359 3.414 0.000733 ***
># ---
># Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>#
># Residual standard error: 5.73 on 285 degrees of freedom
># (89 observations deleted due to missingness)
># Multiple R-squared: 0.5864, Adjusted R-squared: 0.5792
># F-statistic: 80.82 on 5 and 285 DF, p-value: < 2.2e-16summary(pool(fit_pmm2))># term estimate std.error statistic df
># 1 (Intercept) 12.0925731 2.36136453 5.1210108 147.95982
># 2 sexchlapec -0.1694591 0.65189098 -0.2599501 151.48082
># 3 sch_typezvláštní -9.8855321 1.63895739 -6.0315980 99.31339
># 4 ethnicRom -6.7492211 0.98031853 -6.8847225 184.72598
># 5 lem_total 0.1111303 0.01339554 8.2960656 139.29179
># 6 sch_typezvláštní:ethnicRom 6.9162139 1.96387635 3.5217156 204.73981
># p.value
># 1 9.323636e-07
># 2 7.952553e-01
># 3 2.793431e-08
># 4 8.760104e-11
># 5 8.304468e-14
># 6 5.285581e-04Vypadá to, že v tomto konkrétním případě nikoli, ale obecně tomu tak být nemusí.