Kapitel 40 Fehlende Daten
Unkorrigierte, unfertige Version 24 December, 2022
Sie lernen in diesem Kapitel folgendes:
- Missing Completely At Random (MCAR)
- Missing At Random
- Missing Not At Random
To Do:
- Warum man fehlende Daten nicht einfach mit dem Mittelwert der bestehenden Daten ersetzen soll*
- Warum Last Observation Carried Forward oft nicht die beste Methode ist
- Die Vorteile der Multiple Imputation und deren Voraussetzungen
Sie finden hier einen vertieften Text über fehlende Werte.
Praktisch kein Datensatz ist komplett, meistens fehlen Daten.Daten können aus verschiedenen Gründen fehlen. Im besten Falle ist es bloss purer Zufall, es gibt kein spezieller Grund, warum ein Wert fehlt, respektive der Grund hat nichts mit den Werten anderer Variablen oder mit den Eigenschaften der Studienteilnehmenden zu tun. In diesem Fall werden die Analysen nicht durch diese fehlenden Werte verfälscht, einzig die statistische Präzision sinkt. Im Fachjargon nennt man dieses Muster Missing Completely At Random (MCAR).
Manchmal gibt es Gründe, die in den Eigenschaften der Studienteilnehmenden liegen und wir können die fehlenden Werte durch diese Eigenschaften schätzen. Die fehlenden Werte sind also nicht komplett zufällig, aber innerhalb der durch die Eigenschaften definierten Untergruppen sind die fehlen die Werte zufällig. Wir nenenn diese Situation im Fachjargon Missing At Random. Wenn wir die fehlenden Werte mit den zugrundeliegenden Variablen schätzen können, so werden die Analysen nicht verfälscht sein.
Wenn es Gründe für die fehlenden Werte gibt, die wir nicht kennen, oder die leider sogar mit dem wahren aber unbekannten (weil fehlenden Wert) zu tun haben, so sind die Werte Missing Not At Random. Beispiel: Wenn wir eine Studie über Müdigkeit durchführen, kann es sein, dass müdere Patient:innen mehr fehlende Werte im Müdigkeitsfragebogen aufweisen, als weniger müde Patient:innen. Dadurch können wir leider die fehlenden Werte nicht gut schätzen und die Analysen werden durch die fehlenden Werte verfälscht (Bias).
Wir werden verschiedene Beispiele benutzen um die Situationen Missing Completely At Random (MCAR), Missing At Random (MAR), and Missing Not At Random zu illustrieren.
40.0.1 Wir benutzen folgende Pakete in diesem Kapitel:
knitr::opts_chunk$set(echo = TRUE)
r <- getOption("repos")
r["CRAN"] <- "https://stat.ethz.ch/CRAN/"
options(repos = r)
list.of.packages <- c("patchwork","bookdown","rmarkdown" ,"knitr","rio", "psych","janitor",
"tidyverse","jtools","summarytools", "qgraph", "gtsummary" , "viridis", "wesanderson", "missMethods", "ggpubr", "ggrepel", "naniar", "finalfit", "missMethods", "rpart", "rpart.plot")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
library(summarytools)
library(psych)
library(janitor)
library(sjlabelled)
library(tidyverse)
library(gtsummary)
library(viridis)
library(wesanderson)
library(missMethods)
library(ggpubr)
library(ggrepel)
library(naniar)
library(finalfit)
library(rpart)
library(rpart.plot)
library(mice)
library(knitr)
library(patchwork)40.1 Missing Completely at Random (MCAR)
Zur Illustration benutzen wir Daten einer im Journal Plos One veröffentlichten Studie. Die Daten dieser Studie sind komplett, aber wir werden zur Illustration zufällig Werte löschen.
Wir laden die Daten direkt vom Internet:
df1<-rio::import("https://doi.org/10.1371/journal.pone.0262238.s008", format="xlsx")
df1<-df1 %>%
rename(WalkingDistance_m_6min=Distance30) %>%
select(Number, ActiveSmoking, Age, Sex, COPDduration,FEV1, WalkingDistance_m_6min ) %>%
mutate(Sex=factor(Sex, levels=c(1,2), labels=c("Men", "Women")))40.1.1 Komplette Daten
Sie sehen hier eine Zusammenfassung der kompletten Daten (enthalten keine fehlenden Daten)
options(width = 300)
summary(df1)## Number ActiveSmoking Age Sex COPDduration FEV1 WalkingDistance_m_6min
## Min. : 1.00 Min. : 0.0 Min. :53.00 Men :46 Min. : 1.00 Min. :0.440 Min. :214.7
## 1st Qu.:14.25 1st Qu.: 0.0 1st Qu.:64.00 Women: 4 1st Qu.: 7.25 1st Qu.:1.015 1st Qu.:304.8
## Median :27.50 Median : 0.0 Median :68.50 Median :16.00 Median :1.540 Median :371.4
## Mean :28.10 Mean :160.1 Mean :69.08 Mean :25.40 Mean :1.496 Mean :359.8
## 3rd Qu.:41.75 3rd Qu.: 1.0 3rd Qu.:76.00 3rd Qu.:44.25 3rd Qu.:1.805 3rd Qu.:408.0
## Max. :56.00 Max. :999.0 Max. :80.00 Max. :72.00 Max. :2.530 Max. :555.3Wir sehen mit dem Summary Befehl, das keine fehlende Werte angegeben werden (da es keine hat).
40.1.2 Streudiagramm FEV1 - Six Minute Walking Distance (keine fehlenden Werte)
Wir zeichnen ein Streudiagramm der Distanz im sechs Minuten Gehtest über die FEV1 (Der FEV1-Wert (Sekundenkapazität) wird zur Beurteilung der Lungenfunktion benutzt. Kann weniger Luft in einer Sekunge ausgeatmet werden, deutet dies auf eine Verengung der Bronchien hin). Die Korrelation (Pearson’s r) schreiben wir direkt in die Graphik hinein. Der Datensatz enthält viel mehr Männer als Frauen, deshalb ist die Schätzung der Korrelation für die Frauen sehr unsicher, was wir am breiten 95% Konfidenzband (grauer Bereich um die Regressionslinie) erkennen.
names(df1)## [1] "Number" "ActiveSmoking" "Age" "Sex" "COPDduration" "FEV1" "WalkingDistance_m_6min"correlation_men<-cor.ci(df1[df1$Sex=="Men",c("FEV1", "WalkingDistance_m_6min")], plot=FALSE)
correlation_women<-cor.ci(df1[df1$Sex=="Women",c("FEV1", "WalkingDistance_m_6min")], plot=FALSE)
cor_ci_men<-paste0(round(correlation_men$rho[2,1],2), " (",round(correlation_men$ci[2],2), " to ", round(correlation_men$ci[4],2),")")
cor_ci_women<-paste0(round(correlation_women$rho[2,1],2), " (",round(correlation_women$ci[2],2), " to ", round(correlation_women$ci[4],2),")")
p1<-ggplot(data=df1, aes(FEV1, WalkingDistance_m_6min, group=Sex, colour=Sex))+
geom_point()+
# stat_cor(method = "pearson", label.x = 1, label.y = 100)+ # This line does overplot both correlations
geom_smooth(method="lm")+
annotate("text", x = 1.8, y = 365,label = paste("r = " , cor_ci_men))+
annotate("text", x = 0.85, y = 255,label = paste("r = " ,cor_ci_women))
p1
Abbildung 40.1: Scatterplot with regression lines seperately for men and women.
40.1.3 Korrelationsmatrix (Daten enthalten keine fehlenden Werte)
options(width = 300)
round(cor(df1[df1$Sex=="Women",c(-1, -4)], use="all.obs"),2)## ActiveSmoking Age COPDduration FEV1 WalkingDistance_m_6min
## ActiveSmoking 1.00 -0.11 0.62 -0.65 -0.79
## Age -0.11 1.00 -0.75 -0.02 -0.34
## COPDduration 0.62 -0.75 1.00 -0.62 -0.36
## FEV1 -0.65 -0.02 -0.62 1.00 0.88
## WalkingDistance_m_6min -0.79 -0.34 -0.36 0.88 1.00round(cor(df1[df1$Sex=="Men",c(-1, -4)], use="all.obs"),2)## ActiveSmoking Age COPDduration FEV1 WalkingDistance_m_6min
## ActiveSmoking 1.00 -0.01 -0.25 -0.15 -0.01
## Age -0.01 1.00 0.11 -0.21 -0.55
## COPDduration -0.25 0.11 1.00 -0.04 -0.12
## FEV1 -0.15 -0.21 -0.04 1.00 0.28
## WalkingDistance_m_6min -0.01 -0.55 -0.12 0.28 1.0040.2 Zufälliges Löschen von Daten (d.h. Missing Completely At Random)
Mit dem Paket missMethods können wir fehlende Werte produzieren, die Missing Completely At Random sind. Paket missMethods, siehe hier für ein Tutorial (hier klicken).
Wir löschen 30% der Werte in den Variablen WalkingDistance_m_6min, FEV1, und Age:
set.seed(12345)
df1_mcar <- delete_MCAR(df1, 0.15, c("WalkingDistance_m_6min", "FEV1", "Age"))40.2.1 Korrelationsmatrix (Daten enthalten fehlende Werte)
R kann natürlich die Korrelation nur mit den nicht-fehlenden Werten berechnen. Wir müssen R sagen, wie die Fälle ausgewählt werden:
- pairwise.complete.obs: Hier benutzen wir alle Beobachtungen, die im Paar zweier Variablen keine fehlenden Werte haben.
- complete.obs: hier werden nur Fälle benutzt, die in allen Variablen keine fehlenden Werte haben (wird auch Listwise Deletion genannt)
Korrelationen für die Frauen:
round(cor(df1_mcar[df1_mcar$Sex=="Women",c(-1, -4)], use="pairwise.complete.obs"),2)## ActiveSmoking Age COPDduration FEV1 WalkingDistance_m_6min
## ActiveSmoking 1.00 -0.11 0.62 -0.65 -0.79
## Age -0.11 1.00 -0.75 -0.02 -0.34
## COPDduration 0.62 -0.75 1.00 -0.62 -0.36
## FEV1 -0.65 -0.02 -0.62 1.00 0.88
## WalkingDistance_m_6min -0.79 -0.34 -0.36 0.88 1.00Korrelationen für die Männer
round(cor(df1_mcar[df1_mcar$Sex=="Men",c(-1, -4)], use="pairwise.complete.obs"),2)## ActiveSmoking Age COPDduration FEV1 WalkingDistance_m_6min
## ActiveSmoking 1.00 0.05 -0.25 -0.06 0.16
## Age 0.05 1.00 0.17 -0.32 -0.54
## COPDduration -0.25 0.17 1.00 -0.11 -0.03
## FEV1 -0.06 -0.32 -0.11 1.00 0.08
## WalkingDistance_m_6min 0.16 -0.54 -0.03 0.08 1.00Wir sehen, dass die Korrelationen nicht mehr die gleichen sind (für die Variablen mit fehlenden Werten). Aber der Unterschied ist nur zufällig, und wir dürfen die Daten mit den fehlenden Werten verwenden. Es gibt keine Verzerrung durch die fehlenden Werte.
40.2.2 Scatterplot (Missing Completely At Random)
And we plot again the same scatter plot as above:
names(df1_mcar)## [1] "Number" "ActiveSmoking" "Age" "Sex" "COPDduration" "FEV1" "WalkingDistance_m_6min"correlation_men<-cor.ci(df1_mcar[df1_mcar$Sex=="Men",c("FEV1", "WalkingDistance_m_6min")], plot=FALSE)
correlation_women<-cor.ci(df1_mcar[df1_mcar$Sex=="Women",c("FEV1", "WalkingDistance_m_6min")], plot=FALSE)
cor_ci_men<-paste0(round(correlation_men$rho[2,1],2), " (",round(correlation_men$ci[2],2), " to ", round(correlation_men$ci[4],2),")")
cor_ci_women<-paste0(round(correlation_women$rho[2,1],2), " (",round(correlation_women$ci[2],2), " to ", round(correlation_women$ci[4],2),")")
p2<-ggplot(data=df1_mcar, aes(FEV1, WalkingDistance_m_6min, group=Sex, colour=Sex))+
geom_point()+
# stat_cor(method = "pearson", label.x = 1, label.y = 100)+ # This line does overplot both correlations
geom_smooth(method="lm")+
annotate("text", x = 1.8, y = 365,label = paste("r = " , cor_ci_men))+
annotate("text", x = 0.85, y = 255,label = paste("r = " ,cor_ci_women))
combined<-p1 + p2 +
plot_layout(widths = c(4, 4), guides = "collect") &
theme(
legend.position = "top", legend.title = element_blank(),
legend.margin = margin(5, 0, 5, 0), legend.box.margin = margin(0, 0, -10, 0),
legend.justification = "right",
plot.margin = margin(10, 10, 10, 10)
)
combined
Abbildung 40.2: Die Korrelationen mit dem Datensatz mit fehlenden Werten sind nicht genau die gleichen wie die ohne fehlende Werte. Der Unterschied ist jedoch nur ein zufälliger Unterschied und das Ergebnis ist nicht verzerrt.
Da der Wert der Variablen oder die Werte anderer Variablen - unabhängig davon, ob wir sie gemessen haben oder nicht - keinen Einfluss darauf haben, ob unsere Daten fehlen, können wir sagen, dass die Ergebnisse immer noch repräsentativ für die Grundgesamtheit sind, d. h. wir haben keine Verzerrung (d.h. keinen Bias). Aber wir verlieren an statistischer Präzision.
Wir können dies testen:
40.3 Testen ob fehlende Werte Missing Completely At Random oder nicht
Statistisch können wir nur testen, ob die fehlenden Daten Missing Completely At Random sind oder nicht. Dazu benutzen wir den Little’s test hier klicken für den Originalartikel.
Im naniar-Paket finden wir den Befehl mcar_test. Ist der P-Wert dieses Tests klein, so können wir nicht ausschliessen, dass die fehlenden Daten nicht Missing Completely At Random sind.In unserem Beispiel ist der P-Wert sehr hoch, deswegen haben wir keinen Verdacht, dass die Daten nicht Missing Completely At Random” sind. Was wir in unserem Beispiel ja schon wissen, da wir die Daten so manipuliert haben, dass die fehlenden Werte ganz sicher Missing Completely At Random* sind.
Sie finden hier mehr Informationen zum Test.
littleltest<-naniar::mcar_test(df1_mcar)Der nächste Code-chunk ist nur zum Erstellen von Textvarianten, die wir je nach P-Wert weiter unten im Text benutzen.
if (littleltest$p.value < 0.05) {
interpretation<-"there is evidence that the data are not Missing Completely At Random, therefore, we assume that they are either Missing At Random or Missing Not At Random."
} else {
interpretation<-"there is no evidence that the data are not Missing Completely At Random."
}
if (littleltest$p.value < 0.001) {
pvalue="p < 0.001"
} else {
pvalue=round(littleltest$p.value, 3)
}We see that the p-value is: 0.356, therefore we can conclude that there is no evidence that the data are not Missing Completely At Random.
40.3.1 Beispiel mit mehr Daten
Das obige Beispiel ist vielleicht etwas verwirrend: Mit den Missing Completely At Random Daten erhielten wir für die Männer einen anderen Korrelationskoeffizienten (d.h. r = 0.28 versus r = 0.08). Dieser Unterschied ist jedoch kein Bias, sondern nur ein zufälliger Fehler. Um das noch besser zu illustrieren, erstellen wir einen grösseren Datensatz aus den oben benutzten Daten. Wir vervierfachen einfach unseren Datensatz.
df1<-bind_rows(df1, df1, df1, df1)Wir löschen 30% der Werte in den Variablen WalkingDistance_m_6min, FEV1, und Age:
set.seed(12345)
df1_mcar <- delete_MCAR(df1, 0.15, c("WalkingDistance_m_6min", "FEV1", "Age"))In der nächsten Abbildung die Korrelation zwischen FEV1 und Sechs-Minuten Gehtest im kompletten Datensatz und in dem mit Missing Completely At Random. Da wir hier ein viel grösseres Datenset haben, sehen wir praktische keinen Untreschied in der Korrelation zwischen den Daten ohne fehlende Werte und den daten mit MCAR fehlenden Werten, obschon die 30% der Daten fehlen. Dieses Beispiel hier ist praktisch identisch mit dem obigen Beispiel, nur mit dem Unterschied, dass wir die obigen Daten vier Mal kopiert haben. Sind die Daten wirklich MCAR, verursachen die fehlenden Werte keine Bias, sondern nur einen Verlust der statistischen Präzision - was zu zufällig unterschiedlcihen Resultaten führt.
correlation_men<-cor.ci(df1[df1$Sex=="Men",c("FEV1", "WalkingDistance_m_6min")], plot=FALSE)
correlation_women<-cor.ci(df1[df1$Sex=="Women",c("FEV1", "WalkingDistance_m_6min")], plot=FALSE)
cor_ci_men<-paste0(round(correlation_men$rho[2,1],2), " (",round(correlation_men$ci[2],2), " to ", round(correlation_men$ci[4],2),")")
cor_ci_women<-paste0(round(correlation_women$rho[2,1],2), " (",round(correlation_women$ci[2],2), " to ", round(correlation_women$ci[4],2),")")
p1<-ggplot(data=df1, aes(FEV1, WalkingDistance_m_6min, group=Sex, colour=Sex))+
geom_point()+
# stat_cor(method = "pearson", label.x = 1, label.y = 100)+ # This line does overplot both correlations
geom_smooth(method="lm")+
annotate("text", x = 1.8, y = 365,label = paste("r = " , cor_ci_men))+
annotate("text", x = 0.85, y = 255,label = paste("r = " ,cor_ci_women))+
labs(title="Keine fehlenden Werte")
correlation_men<-cor.ci(df1_mcar[df1_mcar$Sex=="Men",c("FEV1", "WalkingDistance_m_6min")], plot=FALSE)
correlation_women<-cor.ci(df1_mcar[df1_mcar$Sex=="Women",c("FEV1", "WalkingDistance_m_6min")], plot=FALSE)
cor_ci_men<-paste0(round(correlation_men$rho[2,1],2), " (",round(correlation_men$ci[2],2), " to ", round(correlation_men$ci[4],2),")")
cor_ci_women<-paste0(round(correlation_women$rho[2,1],2), " (",round(correlation_women$ci[2],2), " to ", round(correlation_women$ci[4],2),")")
p2<-ggplot(data=df1_mcar, aes(FEV1, WalkingDistance_m_6min, group=Sex, colour=Sex))+
geom_point()+
# stat_cor(method = "pearson", label.x = 1, label.y = 100)+ # This line does overplot both correlations
geom_smooth(method="lm")+
annotate("text", x = 1.8, y = 365,label = paste("r = " , cor_ci_men))+
annotate("text", x = 0.85, y = 255,label = paste("r = " ,cor_ci_women))+
labs(title="30% MCAR fehlende Werte")
combined<-p1 + p2 +
plot_layout(widths = c(4, 4), guides = "collect") &
theme(
legend.position = "top", legend.title = element_blank(),
legend.margin = margin(5, 0, 5, 0), legend.box.margin = margin(0, 0, -10, 0),
legend.justification = "right",
plot.margin = margin(10, 10, 10, 10)
)
combined
Abbildung 40.3: Streudiagramme mit und ohne fehlende Werte. Da wir hier ein grösseres Datenset haben, ist der zufällige Fehler hier deutlich kleiner als im obigen Beispiel. Wir sehen praktisch keinen Unterschied zwischen den Korrelationskoeffizienten mit und ohne fehlende Werte.
40.4 Missing At Random
40.4.1 Daten erstellen für dieses MAR Beispiel
Für dieses Beispiel extrahieren wir Daten aus der Graphik eines Buches - das Sie übrigens unbedingt lesen sollten.
- Multiple Imputation and its Application, Seite 13.
Wir benutzen Webplotdigitizer um die Daten zu extrahieren. Hier ein Link zu einem Video mit einer Erkärung, wie das gemacht werden kann.
Im Buch werden nur die Mittelwerte der Untergruppen angegeben und nicht die Standardabweichungen. Mit Mittelwerten und Standardabweichungen hätten wir sonst die Daten direkt mmit rnorm simulieren können.
Das spezielle an diesem Beispiel ist folgendes: Wir kennen in diesem Beispiel die fehlenden Werte - was wir ja in der realen Welt praktisch nie kennen. Im Beispiel werden die Löhne für zwei Arbeitsstellen angegeben.Wie könnten wir in der realen Welt an die Werte der fehlenden Daten gelangen? Wir könnten zu den Studienteilnehmenden zurückgehen, die die Löhne nicht angegeben haben, und sie fragen, ob sie nicht uns die Löhne trotzdem nennen würden. Vielleicht würden Sie es tun, wenn Sie wüssten, dass wir ihre Löhne nur für diese Übung benutzen würden. Aber eben: Es gab ja wohl Gründe, warum sie den Lohn nicht angegeben haben, und deshalb werden sie wohl auch beim direkten Nachfragen den Lohn nicht preisgeben.
data<-read.table(text="Count,Job_A_Observed,Job_A_Missing,Job_B_Observed,Job_B_Missing
1,48.79310345,48.79310345,20.28735632,22.57471264
2,49.94252874,50.86206897,21.66666667,21.72413793
3,50.86206897,52.70114943,22.12643678,22.72413793
4,52.01149425,52.70114943,22.35632184,27.87356322
5,52.70114943,53.62068966,22.81609195,27.87356322
6,53.3908046,54.54022989,22.81609195,29.86206897
7,54.31034483,55.91954023,23.27586207,33.62068966
8,54.77011494,56.37931034,23.73563218,34.77011494
9,55.22988506,56.37931034,24.65517241,36.91954023
10,55.91954023,57.29885057,24.65517241,38.90804598
11,56.37931034,57.29885057,25.34482759,33.28735632
12,56.83908046,70.2183908,26.49425287,NA
13,57.29885057,60.74712644,19.49425287,NA
14,57.29885057,61.66666667,22.72413793,NA
15,58.44827586,61.66666667,21.95402299,NA
16,58.67816092,71.35632184,20.18390805,NA
17,59.59770115,63.50574713,22.18390805,NA
18,60.05747126,65.34482759,21.64367816,NA
19,60.97701149,65.8045977,27.64367816,NA
20,60.97701149,66.72413793,28.33333333,NA
21,61.66666667,70.64367816,28.56321839,NA
22,62.12643678,67.64367816,29.25287356,NA
23,63.04597701,63.27586207,21.25287356,NA
24,63.73563218,62.5862069,30.17241379,NA
25,64.1954023,61.20689655,30.17241379,NA
26,65.11494253,57.29885057,30.86206897,NA
27,66.03448276,56.14942529,31.55172414,NA
28,66.72413793,60.22988506,32.01149425,NA
29,67.87356322,60.74712644,32.70114943,NA
30,69.94252874,62.12643678,33.62068966,NA
31,70.86206897,57.90068966,28.31034483,NA
32,72.01149425,59.22988506,34.31034483,NA
33,63.50574713,NA,35.22988506,NA
34,62.5862069,NA,35.91954023,NA
35,61.66666667,NA,36.6091954,NA
36,60.51724138,NA,37.06896552,NA
37,59.82758621,NA,38.2183908,NA
38,58.67816092,NA,38.67816092,NA
39,57.98850575,NA,38.67816092,NA
40,57.06896552,NA,39.36781609,NA
41,59.13793103,NA,25.51724138,NA
42,60.51724138,NA,35.51724138,NA
43,61.43678161,NA,40.97701149,NA
44,62.35632184,NA,36.66666667,NA
45,63.50574713,NA,35.45977011,NA
46,63.50574713,NA,34.77011494,NA
47,61.66666667,NA,34.54022989,NA
48,60.74712644,NA,34.08045977,NA
49,59.59770115,NA,33.3908046,NA
50,58.67816092,NA,33.16091954,NA
51,57.75862069,NA,33.16091954,NA
52,56.14942529,NA,32.47126437,NA
53,55.22988506,NA,32.24137931,NA
54,54.08045977,NA,31.55172414,NA
55,53.3908046,NA,31.32183908,NA
56,57.52873563,NA,30.63218391,NA
57,62.5862069,NA,30.40229885,NA
58,64.65517241,NA,30.17241379,NA
59,65.34482759,NA,29.71264368,NA
60,66.26436782,NA,28.79310345,NA
61,66.95402299,NA,28.33333333,NA
62,63.96551724,NA,27.64367816,NA
63,63.04597701,NA,27.4137931,NA
64,61.66666667,NA,30.63218391,NA
65,60.51724138,NA,31.55172414,NA
66,59.36781609,NA,31.7816092,NA
67,58.67816092,NA,32.70114943,NA
68,57.29885057,NA,32.93103448,NA
69,NA,NA,33.3908046,NA
70,NA,NA,31.62068966,NA
71,NA,NA,32.08045977,NA
72,NA,NA,24.31034483,NA
73,NA,NA,35,NA
74,NA,NA,26.22988506,NA
75,NA,NA,31.22988506,NA
76,NA,NA,35.68965517,NA
77,NA,NA,35,NA
78,NA,NA,30.08045977,NA
79,NA,NA,30.16091954,NA
80,NA,NA,32.70114943,NA
81,NA,NA,32.24137931,NA
82,NA,NA,31.55172414,NA
83,NA,NA,28.09195402,NA
84,NA,NA,30.63218391,NA
85,NA,NA,22.25287356,NA
86,NA,NA,29.02298851,NA
87,NA,NA,30.40229885,NA
88,NA,NA,31.32183908,NA
89,NA,NA,32.47126437,NA",
header=TRUE, sep=",")40.4.2 Umorganisieren der Daten von Wide zu Long
Wir müssen die Daten reorganisieren - dieser Prozess wird oft als Reshape der Daten bezeichnet, oder in r nennt man das Pivotieren der Daten.
head(data)## Count Job_A_Observed Job_A_Missing Job_B_Observed Job_B_Missing
## 1 1 48.79310 48.79310 20.28736 22.57471
## 2 2 49.94253 50.86207 21.66667 21.72414
## 3 3 50.86207 52.70115 22.12644 22.72414
## 4 4 52.01149 52.70115 22.35632 27.87356
## 5 5 52.70115 53.62069 22.81609 27.87356
## 6 6 53.39080 54.54023 22.81609 29.86207summary(data)## Count Job_A_Observed Job_A_Missing Job_B_Observed Job_B_Missing
## Min. : 1 Min. :48.79 Min. :48.79 Min. :19.49 Min. :21.72
## 1st Qu.:23 1st Qu.:57.24 1st Qu.:56.32 1st Qu.:26.49 1st Qu.:25.30
## Median :45 Median :60.29 Median :60.49 Median :30.63 Median :29.86
## Mean :45 Mean :60.08 Mean :60.08 Mean :30.00 Mean :30.01
## 3rd Qu.:67 3rd Qu.:63.16 3rd Qu.:63.33 3rd Qu.:33.16 3rd Qu.:34.20
## Max. :89 Max. :72.01 Max. :71.36 Max. :40.98 Max. :38.91
## NA's :21 NA's :57 NA's :78Keine Panik, falls Sie diesen pivot_longer nicht gleich verstehen, der Autor dieses Kapitels benötigt jeweils mehrere Stunden um diese paar Zeilen zu verstehen und an einen anderen Datensatz anzupassen. 😄. Vielleicht hilft ein Video: Click here.
Sie finden auch ein Kapitel zu diesen pivot Befehlen.
data_long<-data %>%
pivot_longer(cols=-Count,
names_pattern = "[Job]_(.*)_(.*)",
names_to=c( "Job", "Missingness"),
values_to="Salary_CHF") %>%
filter(!is.na(Salary_CHF))The long data:
head(data_long)## # A tibble: 6 × 4
## Count Job Missingness Salary_CHF
## <int> <chr> <chr> <dbl>
## 1 1 A Observed 48.8
## 2 1 A Missing 48.8
## 3 1 B Observed 20.3
## 4 1 B Missing 22.6
## 5 2 A Observed 49.9
## 6 2 A Missing 50.940.4.3 Beschreibende Statistik
Wir sehen die Löhne separat für die Daten ohne fehlende Werte und für die fehlenden Werte (Achtung: diese fehlenden Daten haben wir ja in Wahrheit praktisch nie, für diese Übung kennen wir sie) 🤕.
summary_table<-data_long %>%
group_by(Missingness, Job) %>%
summarise(mean_Salary=mean(Salary_CHF, na.rm=TRUE),
sd_Salary=sd(Salary_CHF, na.rm=TRUE),
n_nonmissing=sum(!is.na(Salary_CHF)),
n=n())## `summarise()` has grouped output by 'Missingness'. You can override using the `.groups` argument.kable(summary_table)| Missingness | Job | mean_Salary | sd_Salary | n_nonmissing | n |
|---|---|---|---|---|---|
| Missing | A | 60.08023 | 5.751603 | 32 | 32 |
| Missing | B | 30.01254 | 5.996058 | 11 | 11 |
| Observed | A | 60.07776 | 4.851583 | 68 | 68 |
| Observed | B | 29.99948 | 5.050440 | 89 | 89 |
Wir können die Daten auh mit dem Befehl tbl_summary des Paketes gtsummary analysieren. Wir sehen die Anzahl Daten, die eigentlich fehlen (Missing). Die Mittelwerte der Löhne beider Jobs (A und B) finden wir in der untersten Zeile. Diese Werte entsprechen den Löhnen, wenn es keine fehlenden Werte gäbe.
tbl_summary(data_long, by=Job,include=c(Missingness, Salary_CHF), statistic=list(Salary_CHF~"{mean}"), digits=Salary_CHF ~ 1)| Characteristic | A, N = 1001 | B, N = 1001 |
|---|---|---|
| Missingness | ||
| Missing | 32 (32%) | 11 (11%) |
| Observed | 68 (68%) | 89 (89%) |
| Salary_CHF | 60.1 | 30.0 |
| 1 n (%); Mean | ||
40.4.4 Differenz in den Löhnen zwischen beobachteten und “wahren” Lohn
Mit “wahrem” Lohn meinen wir hier den Lohn, wenn wir alle Daten haben, d.h. auch die Werte für die fehlenden Daten.
Der wahre Lohn über beide Arbeitsstellen (A und B) ist: 45.04Der Lohn für beide Arbeitsstellen, wenn wir nur die beobachteten Werte analysieren: 43.03
Der Unterschied ist darauf zurückzuführen, dass sich der Mittelwert des Gehalts der Teilnehmer mit fehlenden Gehaltsangaben von denen unterscheidet, die Gehaltsangaben gemacht haben: 52.39
Dies bedeutet, dass der Gesamtmittelwert für das Gehalt nach unten verzerrt ist (d. h. zu klein im Vergleich zum wahren Mittelwert, den wir nur hätten, wenn wir auch die fehlenden Werte analysieren könnten).
Unsere Daten fehlen also nicht völlig zufällig. Aber sind sie wenigstens “Missing At Random”?40.5 Missing At Random bedingt auf den Job
- Innerhalb beiden Berufen ist die Verteilung der fehlenden und beobachteten Gehälter ähnlich.
means<-aggregate(Salary_CHF~Job+Missingness, data_long, mean)
ggplot(data_long, aes(factor(Job), Salary_CHF))+
geom_point(data=data_long[data_long$Missingness=="Observed",], position=position_nudge(x=-0.02), aes(colour=Missingness))+
geom_point(data=data_long[data_long$Missingness=="Missing",], position=position_nudge(x=0.02), aes(colour=Missingness))+
theme_classic()+
theme(
legend.position=c(0.5, 0.5))+
xlab("Job")+
ylab("Salary (in 1000 CHF)")+
geom_text(data=means[means$Missingness=="Observed",], aes(label=paste0("Mean\n Observed: \n ",round(Salary_CHF)), y= Salary_CHF), nudge_x=-0.15)+
geom_text(data=means[means$Missingness=="Missing",], aes(label=paste0("Mean\n Missing: \n", round(Salary_CHF)), y= Salary_CHF), nudge_x=0.15)+
annotate("text", x=1.8, y=70,label=
"Over both Jobs, the probability of missing values \n increases with the salary. Within Job A and within \n Job B, the probability of missing values does not \n depend on the salary. Therefore, we say that \n values are Missing At Random dependent on the job." )+
annotate("text", x=1, y=15, label="Observed: 68 \n Missing: 32")+
annotate("text", x=2, y=15, label="Observed: 89 \n Missing: 11")+
annotate("text", x=1, y=30, label="In Job A, the distribution of \n observed and unobserved \n salaries is similar")+
annotate("text", x=2, y=55, label="In Job B, the distribution of \n observed and unobserved \n salaries is similar")+
# scale_color_viridis(discrete = TRUE, option="turbo")+
scale_color_manual(values = wes_palette("BottleRocket1", n = 2))
Abbildung 11.1: Wir sehen, dass innerhalb der beiden Stellen die fehlenden Werte nicht vom Gehalt abhängt. In Job A haben diejenigen mit höheren Gehältern keine höhere Wahrscheinlichkeit für fehlende Werte. Das Gleiche gilt für die Personen in Job B. Insgesamt (d.h. über beide Arbeitsstellen) ist die -Wahrscheinlichkeit für fehlende Werte jedoch mit dem Gehalt verbunden. Insgesamt hängt die Wahrscheinlichkeit eines fehlenden Wertes seinem wahren Wert (d. h. je höher das Gehalt, desto höher die Wahrscheinlichkeit eines fehlenden Wertes - mit anderen Worten: Personen mit höheren Gehältern geben ihr Gehalt häufiger nicht an). Daher sagen wir, dass die Werte in Abhängigkeit von der Stelle zufällig fehlen (MAR bedingt auf die Stelle)
40.5.1 Das Problem
In der Realität können wir nicht prüfen, ob die Verteilung der fehlenden Gehälter im Vergleich zu den nicht fehlenden Gehältern ähnlich ist. Denn in der Realität fehlen die fehlenden Gehälter (d. h. wir wissen nicht, wie sie sind, d. h. wir haben keine Informationen über die roten Punkte in der obigen Abbildung).
40.6 Beschreibende Statistik
Nachfolgend nur der Code (klicken Sie auf die Schaltfläche “Anzeigen”, um den Code zu sehen) zur Berechnung der Mittelwerte für die verschiedenen Untergruppen.
cat("Mean Job A: ", mean(data_long$Salary_CHF[data_long$Job=="A"]))## Mean Job A: 60.07855cat("Mean Job B: ", mean(data_long$Salary_CHF[data_long$Job=="B"]))## Mean Job B: 30.00092cat("Mean Job A, only observed: ", mean(data_long$Salary_CHF[data_long$Job=="A"& data_long$Missingness=="Observed"]))## Mean Job A, only observed: 60.07776cat("Mean Job B, only observed: ", mean(data_long$Salary_CHF[data_long$Job=="B"& data_long$Missingness=="Observed"]))## Mean Job B, only observed: 29.99948cat("Mean overall, only observed: ", mean(data_long$Salary_CHF[data_long$Missingness=="Observed"]))## Mean overall, only observed: 43.02702cat("Mean overall, all data (including missing) ", mean(data_long$Salary_CHF))## Mean overall, all data (including missing) 45.0397340.7 Beweis
Der Beweis, dass die Fehlenden Gehälter von den Gehältern über beide Stellen, aber nicht von den Gehältern innerhalb der Stellen abhängig sind:
Auch hier müssen wir bedenken, dass dieser “Beweis” in der realen Welt nicht möglich ist, weil die Daten fehlen, die wir benötigen, um die Annahme zu testen, dass unsere Daten in Abhängigkeit von der Stelle fehlen. 😠 Wir können dies nur in dieser Übung tun, weil wir die Gehälter derjenigen kennen, die fehlende Werte beim Gehalt haben. In der realen Welt könnte dies erreicht werden, indem wir Personen mit fehlenden Gehaltsangaben kontaktieren und sie fragen, ob sie aus akademischen oder didaktischen Gründen bereit wären, Informationen über ihr Gehalt zu geben.
ggplot(data_long, aes(x=Salary_CHF))+
geom_density(aes(fill=Missingness), alpha=0.4)
Abbildung 11.2: We see that there are more missing values in people with higher salaries
ggplot(data_long, aes(x=Salary_CHF))+
geom_density(aes(fill=Missingness), alpha=0.4)+
facet_grid(~Job)
Abbildung 11.3: If we analyse both jobs seperately, we see that the salaries of the missing and the observed are similarly distributed
40.7.1 Logistische Regression auf fehlende Werte
Wir können nun eine logistische Regression durchführen, um zu ermitteln, ob das Gehalt mit der Wahrscheinlichkeit fehlender Werte zusammenhängt. Wenn das Gehalt einen Zusammenhang mit der Wahrscheinlichkeit fehlender Werte hat, dann sagen wir, dass die Missingness der Variable vom wahren (d.h. unbeobachteten) Wert der Variable abhängt. In diesem Fall sehen wir, dass - wenn beide Stellen zusammen analysiert werden - die Fehlendheit des Gehalts vom Gehalt abhängt. Beachten Sie: In einer realen Situation würden wir dies nicht wissen, da wir den Wert der fehlenden Gehälter nicht kennen würden. Wir könnten dies herausfinden, wenn wir in einer Studie Personen mit fehlenden Gehaltswerten kontaktieren und sie nach ihrem Gehalt fragen würden. Oder wir fragen einige 👽, die wissen das alles.
Wir erstellen eine Variable, die 1 ist, wenn der Wert fehlt, und 0, wenn der Wert nicht fehlt. Dann führen wir eine logistische Regression mit dieser Variable als abhängige Variable und dem Gehalt als unabhängige Variable durch, um zu sehen, ob die Wahrscheinlichkeit, dass Werte fehlen, mit steigendem Gehalt zunimmt.
In der folgenden Ausgabe sehen wir, dass die Wahrscheinlichkeit fehlender Werte mit steigendem Gehalt zunimmt. Die Wahrscheinlichkeit, fehlende Werte zu haben, steigt um 4%, wenn man 1’000 CHF mehr verdient (Odds Ratio von 1.04 pro Einheit Gehaltserhöhung, eine Einheit ist hier 1’000 CHF).
data_long<-data_long %>%
mutate(missing=case_when(
Missingness=="Missing"~1,
Missingness=="Observed"~0))
logreg_missing_overall<-glm(missing~Salary_CHF, data=data_long, family=binomial)
summary(logreg_missing_overall)##
## Call:
## glm(formula = missing ~ Salary_CHF, family = binomial, data = data_long)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0370 -0.8186 -0.5243 -0.4343 2.2003
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.18659 0.62884 -5.067 4.03e-07 ***
## Salary_CHF 0.03954 0.01193 3.314 0.000919 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 208.20 on 199 degrees of freedom
## Residual deviance: 196.11 on 198 degrees of freedom
## AIC: 200.11
##
## Number of Fisher Scoring iterations: 4jtools::summ(logreg_missing_overall, confint=TRUE, exp=TRUE)| Observations | 200 |
| Dependent variable | missing |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(1) | 12.10 |
| Pseudo-R² (Cragg-Uhler) | 0.09 |
| Pseudo-R² (McFadden) | 0.06 |
| AIC | 200.11 |
| BIC | 206.70 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 0.04 | 0.01 | 0.14 | -5.07 | 0.00 |
| Salary_CHF | 1.04 | 1.02 | 1.06 | 3.31 | 0.00 |
| Standard errors: MLE |
40.7.2 Abbildung MAR
Anhand eines Diagramms wollen wir zeigen, dass die fehlenden Werte vom Gehalt abhängig sind, wenn wir die Informationen über die Art der Tätigkeit nicht verwenden.
Zunächst stellen wir das Gehalt derjenigen die fehlen separat von denen die nicht fehlen dar. Wichitg: in der realen Welt kennen wir natürlich das Gehalt derjenigen, die denn Lohn nicht angegeben haben, nicht.
p1<-ggplot(data_long, aes(Salary_CHF, missing))+
geom_point(alpha=0.3)
# Here the probability over all, we see that the probabilty for having missing values in salaryincreases with salary
p2<-p1+
geom_smooth(method=glm, method.args=list(family="binomial"))+
labs(
title="Logistic Regression Model",
x="Salary 1000 CHF",
y="Probability of having Missing Values in Salary_CHF"
)
print(p2)## `geom_smooth()` using formula = 'y ~ x'
Abbildung 22.4: The probability of missing values in salaryincreases with increasing salary. People with higher values do more often miss to provide information on salary.
40.7.3 Job-Informationen mitberücksichtigen
Wir sehen hier unten, dass es keinen Zusammenhang zwischen dem Gehalt und der Wahrscheinlichkeit fehlender Gehaltsdaten gibt, wenn man die Stelle (Job A oder B) berücksichtigt. Lassen Sie sich nicht von den verschiedenen Methoden verwirren, die im folgenden Code-Block verwendet werden. Im ersten Modell gehen wir davon aus, dass es keine Wechselwirkung zwischen der Stelle und dem Gehalt gibt, was die fehlenden Gehaltsangaben betrifft. Dies ist in der Tat nicht die beste Art, dies zu analysieren. Wir sollten eine Interaktion einbeziehen (oder die Assoziation getrennt für beide Stellen analysieren). Im zweiten und dritten Modell haben wir die Assoziation getrennt für diejenigen mit Job A und Job B analysiert. Im vierten Modell haben wir eine Wechselwirkung zwischen der Stelle und dem Gehalt modelliert, um den Zusammenhang mit der Wahrscheinlichkeit fehlender Gehaltsangaben zu ermitteln. Der fünfte Ansatz zeigt nur, wie wir die Modelle zwei und drei mit dem Purrr-Ansatz durchführen können (dies ist wirklich nur für die sehr interessierten Leser:innen).
logreg_missing_per_Job<-glm(missing~Salary_CHF+Job, data=data_long, family=binomial)
summary(logreg_missing_per_Job)##
## Call:
## glm(formula = missing ~ Salary_CHF + Job, family = binomial,
## data = data_long)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.8792 -0.8780 -0.4829 -0.4824 2.1019
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.7670618 2.1074489 -0.364 0.716
## Salary_CHF 0.0002212 0.0348951 0.006 0.995
## JobB -1.3303162 1.1178086 -1.190 0.234
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 208.20 on 199 degrees of freedom
## Residual deviance: 194.68 on 197 degrees of freedom
## AIC: 200.68
##
## Number of Fisher Scoring iterations: 4data_A<-data_long %>%
filter(Job=="A")
logreg_missing_Job_A<-glm(missing~Salary_CHF, data=data_A, family=binomial)
summary(logreg_missing_Job_A)##
## Call:
## glm(formula = missing ~ Salary_CHF, family = binomial, data = data_A)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.8787 -0.8783 -0.8781 1.5094 1.5101
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.595e-01 2.534e+00 -0.300 0.764
## Salary_CHF 9.508e-05 4.202e-02 0.002 0.998
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 125.37 on 99 degrees of freedom
## Residual deviance: 125.37 on 98 degrees of freedom
## AIC: 129.37
##
## Number of Fisher Scoring iterations: 4data_B<-data_long %>%
filter(Job=="B")
logreg_missing_Job_B<-glm(missing~Salary_CHF, data=data_B, family=binomial)
summary(logreg_missing_Job_B)##
## Call:
## glm(formula = missing ~ Salary_CHF, family = binomial, data = data_B)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.4840 -0.4831 -0.4828 -0.4820 2.1028
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.1057881 1.9069362 -1.104 0.269
## Salary_CHF 0.0005015 0.0626423 0.008 0.994
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 69.303 on 99 degrees of freedom
## Residual deviance: 69.303 on 98 degrees of freedom
## AIC: 73.303
##
## Number of Fisher Scoring iterations: 4logreg_missing_per_Job_interaction<-glm(missing~Salary_CHF*Job, data=data_long, family=binomial)
summary(logreg_missing_per_Job)##
## Call:
## glm(formula = missing ~ Salary_CHF + Job, family = binomial,
## data = data_long)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.8792 -0.8780 -0.4829 -0.4824 2.1019
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.7670618 2.1074489 -0.364 0.716
## Salary_CHF 0.0002212 0.0348951 0.006 0.995
## JobB -1.3303162 1.1178086 -1.190 0.234
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 208.20 on 199 degrees of freedom
## Residual deviance: 194.68 on 197 degrees of freedom
## AIC: 200.68
##
## Number of Fisher Scoring iterations: 4# with purrr --> this is the same as doing both anaylsis seperatly (i.e. per Job)
library(purrr)
data_long %>%
split(.$Job) %>%
map(~glm(missing~Salary_CHF, family=binomial, data=.)) %>%
map(summary)## $A
##
## Call:
## glm(formula = missing ~ Salary_CHF, family = binomial, data = .)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.8787 -0.8783 -0.8781 1.5094 1.5101
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.595e-01 2.534e+00 -0.300 0.764
## Salary_CHF 9.508e-05 4.202e-02 0.002 0.998
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 125.37 on 99 degrees of freedom
## Residual deviance: 125.37 on 98 degrees of freedom
## AIC: 129.37
##
## Number of Fisher Scoring iterations: 4
##
##
## $B
##
## Call:
## glm(formula = missing ~ Salary_CHF, family = binomial, data = .)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.4840 -0.4831 -0.4828 -0.4820 2.1028
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.1057881 1.9069362 -1.104 0.269
## Salary_CHF 0.0005015 0.0626423 0.008 0.994
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 69.303 on 99 degrees of freedom
## Residual deviance: 69.303 on 98 degrees of freedom
## AIC: 73.303
##
## Number of Fisher Scoring iterations: 4In der nächsten Ausgabe des Modells mit der Interaktion sehen wir, dass es keine Interaktion und keine Assoziation des Gehalts mit dem Fehlen des Gehalts gibt, wenn wir die Informationen über den Arbeitsplatz verwenden.
jtools::summ(logreg_missing_per_Job_interaction, confint=TRUE, exp=TRUE)| Observations | 200 |
| Dependent variable | missing |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(3) | 13.53 |
| Pseudo-R² (Cragg-Uhler) | 0.10 |
| Pseudo-R² (McFadden) | 0.06 |
| AIC | 202.68 |
| BIC | 215.87 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 0.47 | 0.00 | 67.12 | -0.30 | 0.76 |
| Salary_CHF | 1.00 | 0.92 | 1.09 | 0.00 | 1.00 |
| JobB | 0.26 | 0.00 | 130.19 | -0.42 | 0.67 |
| Salary_CHF:JobB | 1.00 | 0.86 | 1.16 | 0.01 | 1.00 |
| Standard errors: MLE |
Die Assoziation des Gehalts mit der Wahrscheinlichkeit fehlender Werte im Gehalt für diejenigen mit Job A:
jtools::summ(logreg_missing_Job_A, confint=TRUE, exp=TRUE)| Observations | 100 |
| Dependent variable | missing |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(1) | 0.00 |
| Pseudo-R² (Cragg-Uhler) | 0.00 |
| Pseudo-R² (McFadden) | 0.00 |
| AIC | 129.37 |
| BIC | 134.58 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 0.47 | 0.00 | 67.12 | -0.30 | 0.76 |
| Salary_CHF | 1.00 | 0.92 | 1.09 | 0.00 | 1.00 |
| Standard errors: MLE |
Die Assoziation des Gehalts mit der Wahrscheinlichkeit fehlender Werte beim Gehalt für Personen mit Job B:
jtools::summ(logreg_missing_Job_B, confint=TRUE, exp=TRUE)| Observations | 100 |
| Dependent variable | missing |
| Type | Generalized linear model |
| Family | binomial |
| Link | logit |
| χ²(1) | 0.00 |
| Pseudo-R² (Cragg-Uhler) | 0.00 |
| Pseudo-R² (McFadden) | 0.00 |
| AIC | 73.30 |
| BIC | 78.51 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 0.12 | 0.00 | 5.11 | -1.10 | 0.27 |
| Salary_CHF | 1.00 | 0.88 | 1.13 | 0.01 | 0.99 |
| Standard errors: MLE |
# now we put all together, but now we had seperate curves for both jobs
ggplot(data_long, aes(Salary_CHF, missing), group=Job)+
geom_point(aes(colour=Job), alpha=0.3)+
geom_smooth(aes(colour=Job), method=glm, method.args=list(family="binomial"))+
labs(
title="Logistic Regression Model",
x="Salary in CHF",
y="Probability of having Missing Values in Salary_CHF"
)## `geom_smooth()` using formula = 'y ~ x'
Abbildung 4.8: Bei einer getrennten Analyse pro Stelle zeigt sich, dass es keinen Zusammenhang zwischen dem Gehalt und der Wahrscheinlichkeit des Fehlens von Gehaltsangaben gibt. Deshalb können wir sagen, dass die Gehaltswerte nach dem Zufallsprinzip fehlen (abhängig von der Stelle).
In einer realen Situation können wir nicht testen, ob die Daten zufällig fehlen, da wir keine Informationen über die fehlenden Werte haben. Es handelt sich also um eine Vermutung.
In der Praxis, wenn wir die fehlenden Gehälter imputieren wollen, wären unsere Imputationen falsch, wenn wir die Arbeitsplatzinformationen verwerfen würden (d.h. wenn wir den Mittelwert der Gehälter über alle Arbeitsplätze imputieren würden), wäre es wichtig, die Arbeitsplatzvariable in der Imputationsmethode zu berücksichtigen.
40.8 Einfache Imputation
Obwohl einfache Imputationen nicht die beste Lösung sind, verwenden wir sie hier, um zu zeigen, dass es wichtig ist, dass wir die relevanten Variablen berücksichtigen.
Zunächst löschen wir die fehlenden Werte in unserem Datensatz. Zur Erinnerung: In dieser didaktischen Übung hatten wir Daten zum Gehalt, die wir als fehlend deklariert hatten - aber wir hatten diese Werte immer noch - so dass wir analysieren konnten, was die Gehälter der fehlenden Gehälter waren - sorry, ein bisschen kompliziert 🤔.
data_missing<-data_long %>%
mutate(Salary_CHF=ifelse(Missingness=="Missing", NA, Salary_CHF))
summary(data_missing)## Count Job Missingness Salary_CHF missing
## Min. : 1.00 Length:200 Length:200 Min. :19.49 Min. :0.000
## 1st Qu.:13.75 Class :character Class :character 1st Qu.:30.17 1st Qu.:0.000
## Median :30.00 Mode :character Mode :character Median :35.52 Median :0.000
## Mean :34.73 Mean :43.03 Mean :0.215
## 3rd Qu.:54.00 3rd Qu.:59.14 3rd Qu.:0.000
## Max. :89.00 Max. :72.01 Max. :1.000
## NA's :43Wir sehen, dass der p-Wert sehr niedrig ist, daher können wir nicht behaupten, dass die Daten völlig zufällig fehlen.
naniar::mcar_test(data_missing[1:4])## # A tibble: 1 × 4
## statistic df p.value missing.patterns
## <dbl> <dbl> <dbl> <int>
## 1 200 3 0 2Denken Sie daran, dass es keinen Test zur Unterscheidung zwischen nicht zufällig fehlenden Werten und zufällig fehlenden Werten gibt.
Nun werden die fehlenden Werte der Gehälter mit dem Gesamtmittelwert der Gehälter imputiert (das ist falsch, aber wir tun es hier trotzdem). Erinnern Sie sich daran, dass der wahre Mittelwert des Gehalts ist 45.04.
data_missing_imputed1<-data_missing %>%
mutate(Salary_CHF=case_when(
is.na(Salary_CHF)~mean(Salary_CHF, na.rm=TRUE),
TRUE~Salary_CHF))Nach der Imputation mit dem mittleren Gesamtgehalt ergibt sich das mittlere Gesamtgehalt: 43.03. Dieser Mittelwert ist verzerrt. In diesem Beispiel wissen wir das, denn wir können den wahren Mittelwert berechnen, der lautet 45.04.
Nun unterstellen wir die mittlere Abhängigkeit des Jobs:
data_missing_imputed2<-data_missing %>%
group_by(Job) %>%
mutate(Salary_CHF=case_when(
is.na(Salary_CHF)~mean(Salary_CHF, na.rm=TRUE),
TRUE~Salary_CHF ))Nach der Imputation mit dem berufsbezogenen Mittelwert ergibt sich der Gesamtmittelwert des Gehalts: 45.04, das entspricht dem wahren Lohn.
40.9 Missing Not At Random
Wenn die Daten nicht zufällig fehlen, hängt die Wahrscheinlichkeit eines fehlenden Wertes vom wahren (zugrunde liegenden) Wert der Variablen ab, und dieser Zusammenhang bleibt bestehen, selbst wenn für alle anderen beobachteten Variablen kontrolliert wird. Das bedeutet, dass wir im Gegensatz zum Gehaltsbeispiel nicht in der Lage sind, eine Situation zu schaffen, in der die Daten “Missing At Random” sind.
Wir verwenden unsere Daten aus dem Artikel, den wir für das Missint Completely At Random-Beispiel verwendet haben und simulieren Daten, die Missing Not At Random sind.
names(df1)## [1] "Number" "ActiveSmoking" "Age" "Sex" "COPDduration" "FEV1" "WalkingDistance_m_6min"df_mnar<-missMethods::delete_MNAR_censoring(df1, 0.3, c("ActiveSmoking", "Age", "Sex","COPDduration", "FEV1", "WalkingDistance_m_6min"))Werfen wir einen Blick auf das Muster der fehlenden Werte: Wir sehen, dass es mehr fehlende Werte bei denjenigen gibt, die eine niedrigere Identifikationsnummer haben, d.h. Patienten, die früher in die Studie eingeschlossen wurden, hatten mehr fehlende Werte, und diejenigen mit einer höheren Gehstrecke im Sechs-Meter-Gehtest hatten eine höhere Wahrscheinlichkeit, fehlende Werte in anderen Variablen zu haben. Was wir natürlich nicht beobachten können, ist, ob der “wahre” Wert einer Variable, zum Beispiel der Gehstrecke, einen Einfluss darauf hatte, ob dieser Wert beobachtet wurde oder fehlte. Aber die Wahrscheinlichkeit ist hoch, dass dies der Fall sein könnte (d.h. dass die Gehstrecke bei den Patienten mit einer guten Gehstrecke fehlte), da diejenigen mit einer höheren Gehstrecke im Allgemeinen mehr fehlende Werte aufwiesen.
df_mnar %>%
add_prop_miss() %>% # Adds Column Containing Proportion Of Missing Data Values
rpart(prop_miss_all ~ ., data = ., model=TRUE) %>% # the rpart function does the regression tree
prp(type = 4, extra = 101, prefix = "Prop. Miss = ") # plots the tree
df_mnar %>%
missing_plot()
df_mnar %>%
missing_pairs(position = "fill", )
Testen wir nun, ob die Daten vollständig zufällig fehlen: Wir sehen, dass der p-Wert sehr niedrig ist, daher können wir die Daten nicht als “Missing Completely At Random” betrachten.
littleltest<-naniar::mcar_test(df_mnar[1:4])
littleltest## # A tibble: 1 × 4
## statistic df p.value missing.patterns
## <dbl> <dbl> <dbl> <int>
## 1 20.4 16 0.203 8if (littleltest$p.value < 0.05) {
interpretation<-"there is evidence that the data are not Missing Completely At Random, therefore, we assume that they are either Missing At Random or Missing Not At Random."
} else {
interpretation<-"there is no evidence that the data are not Missing Completely At Random."
}
if (littleltest$p.value < 0.001) {
pvalue="p < 0.001"
} else {
pvalue=round(littleltest$p.value, 3)
}Wir können folgende Schlussfolgerung ziehen: We see that the p-value is: 0.203, therefore we can conclude that there is no evidence that the data are not Missing Completely At Random.
rio::export(data_long, "data_long.RData")
rio::export(data_missing, "data_missing.RData")
rio::export(df_mnar, "df_mnar.RData")
rio::export(df1_mcar,"df1_mcar.RData")