Der Expectation-Maximization (EM) Algorithmus ist eine generelle Technik, um Maximum-Likelihood-Schätzer zu erhalten, wenn in einem statistischen Modell nicht alle Informationen verfügbar sind, beispielsweise aufgrund fehlender Werte. Der Algorithmus kann Mittelwerte und Varianzen für fehlende Werte schätzen, wenn die Anzahl fehlender Werte definiert ist (und die Anzahl fehlenden Werte ist stets durch die beobachteten Werte der anderen Variablen (mit-)bestimmt). EM bietet damit eine analytische Lösung für die Vorhersage fehlender Werte. Bei einer stochastischen Version des Algorithmus wird die Varianz der vorhergesagten Werte durch Addition eines Fehlers erhöht. Die Fehlerwerte stammen aus einer Normalverteilung mit dem Mittelwert null und der Varianz der fehlenden Werte. Um EM anzuwenden, müssen die Variablen multivariat normalverteilt sein. In R existieren vielfältige Pakete zur Imputation fehlender Werte mit EM.
Wir starten mit einer Datenmatrix, in der einige Werte fehlen:
# Erstellen der Datenmatrix mit fehlenden Werten
d <- matrix(c(2, 3,
4, NA,
NA, 6,
6, 7),
nrow = 4, byrow = TRUE)
# Anzeigen der Datenmatrix
colnames(d) <- c("X", "Y")
d X Y
[1,] 2 3
[2,] 4 NA
[3,] NA 6
[4,] 6 7Wir berechnen die anfänglichen Schätzungen für Mittelwerte und Kovarianzmatrix basierend auf beobachteten Werten:
# Berechnung der Mittelwerte ohne NA
mu <- colMeans(d, na.rm = TRUE)
# Kovarianzmatrix basierend auf vorhandenen Werten
cov_matrix <- cov(d, use = "pairwise.complete.obs")
# Anzeigen der Anfangswerte
mu; cov_matrix X Y
4.000000 5.333333 X Y
X 4 8.000000
Y 8 4.333333Ersetzung der fehlenden Werte durch Regressionsschätzer (bedingte Mittelwerte).
Diese Regressionsschätzer lassen sich mit lm() nicht reproduzieren, weil lm() für X1 ~ X2 ausschließlich vollständige Paare (listwise) verwendet. Der hier verwendete Schätzer basiert dagegen auf Mittelwerten/Varianzen aus allen verfügbaren Einzelwerten und einer Kovarianz aus vollständigen Paaren (pairwise.complete.obs). Diese Mischung führt zu anderen Koeffizienten als ein ‘normales’ Regressionsmodell mit lm().
Wir betrachten die einfache lineare Regression von (y) auf (x):
\(\mathbb{E}(Y \mid X=x) =\hat{y} = a + b\,x\)
\(b = \frac{COV(x,y)}{Var(x)}\)
\(a = \bar{y} - b\,\bar{x}\)
Setzt man das in \(\hat{y} = a + b x\) ein, erhält man:
\(\hat{y} = (\bar{y} - b\bar{x}) + bx = \bar{y} + b(x - \bar{x})\)
\(\mathbb{E}(Y \mid X=x)\: \hat{=}\hat{y}\:=\:\bar{y}+\frac{\operatorname{Cov}(x,y)}{\operatorname{Var}(x)}\,(x-\bar{x})\)
df <- as.data.frame(d)
# Regression über Kovarianzen: Y auf X (E[Y | X])
b_Y_on_X <- cov_matrix["Y","X"] / cov_matrix["X","X"]
a_Y_on_X <- mu["Y"] - b_Y_on_X * mu["X"]
# Regression über Kovarianzen: X auf Y (E[X | Y])
b_X_on_Y <- cov_matrix["X","Y"] / cov_matrix["Y","Y"]
a_X_on_Y <- mu["X"] - b_X_on_Y * mu["Y"]
# Vorhersagen (bedingte Erwartungswerte) einsetzen
idx_X <- is.na(df$X) & !is.na(df$Y)
df$X[idx_X] <- a_X_on_Y + b_X_on_Y * df$Y[idx_X]
idx_Y <- is.na(df$Y) & !is.na(df$X)
df$Y[idx_Y] <- a_Y_on_X + b_Y_on_X * df$X[idx_Y]
# Ergebnis
df X Y
1 2.000000 3.000000
2 4.000000 5.333333
3 5.230769 6.000000
4 6.000000 7.000000Jetzt berechnen wir aus dieser vollständig gefüllten Matrix einfach neue Schätzer: Deshalb ändern sich die imputierten Werte in der nächsten Iteration typischerweise erneut.
mu_m <- colMeans(df, na.rm = TRUE)
# Kovarianzmatrix basierend auf vorhandenen Werten
cov_matrix_m <- cov(df, use = "pairwise.complete.obs")
# Anzeigen der Anfangswerte
mu_m; cov_matrix_m X Y
4.307692 5.333333 X Y
X 3.045365 2.940171
Y 2.940171 2.888889# Iteration 2 E-Schritt mit mu_m, cov_matrix_m1 (auf df1 weiterarbeiten)
miss_X <- is.na(d[, "X"]) & !is.na(d[, "Y"])
miss_Y <- is.na(d[, "Y"]) & !is.na(d[, "X"])
df2 <- df
df2$X[miss_X] <- mu_m["X"] + (cov_matrix_m["X","Y"] / cov_matrix_m["Y","Y"]) * (df2$Y[miss_X] - mu_m["Y"])
df2$Y[miss_Y] <- mu_m["Y"] + (cov_matrix_m["Y","X"] / cov_matrix_m["X","X"]) * (df2$X[miss_Y] - mu_m["X"])
df X Y
1 2.000000 3.000000
2 4.000000 5.333333
3 5.230769 6.000000
4 6.000000 7.000000# Iteration 2: M-Schritt
mu2 <- colMeans(df2)
VKM2 <- cov(df2)
mu2; VKM2 X Y
4.246548 5.259067 X Y
X 2.909811 2.910234
Y 2.910234 2.910951Nach zwei Iterationen haben wir aktualisierte Schätzwerte für die fehlenden Werte sowie für Mittelwerte und Kovarianzmatrix. Der Algorithmus kann weiter iteriert werden, bis Konvergenz erreicht ist.
library(missMethods)
d_imp <- impute_EM(d, stochastic = FALSE, maxits = 2, criterion = 0, verbose = TRUE)
d_imp X Y
[1,] 2.000000 3.000000
[2,] 4.000000 5.150989
[3,] 4.792067 6.000000
[4,] 6.000000 7.000000
attr(,"iterations")
[1] 3Etwas überrascht stellte ich fest, dass die EM-Implementation in R nicht mit meiner schrittweisen Rechnung übereinstimmt. In meiner nachträglichen Recherche sah ich, dass die EM-Schätzung im multivariat-normalen Modell die Kovarianz als Maximum-Likelihood-Schätzer bestimmt und damit eine andere Kovarianzformel verwendet als die Standardfunktioncov() (die durch \(n−1\) teilt). Außerdem werden in missMethods andere Startwerte genutzt (nicht wie bei mir die beobachteten Mittelwerte). Mit dem Paket norm ließe sich das prinzipiell explizit über Startwerte nachbilden. Darauf verzichte ich an dieser Stelle jedoch, um den Fokus des Eintrags auf die Grundidee der Iterationen zu behalten.
Schaffauer (1997) gibt ein Rechenbeispiel zum EM-Algorithmus (Kapitel 3), das ich vor einigen Jahren bereits einmal nachvollzogen habe. In R werde ich dieses Beispiel nachbauen, um eine Übereinstimmung zwischen den Formeln und der Vorgehensweise und dem R-Output zu erzeugen.
Schafer, J. L. (1997). Analysis of incomplete multivariate data. CRC press.