Alle Beispiele rechnen mit denselben 180 Personen aus 12 Schulklassen. Die Daten entstehen in beiden Sprachen aus demselben Lehmer-Generator, R und Python sehen also exakt dieselben Werte.
| Fragestellung | Wie hängt eine Zielgrösse von Prädiktoren ab, wenn die Beobachtungen in Gruppen gebündelt sind? |
| Zielvariable | metrisch (für binäre Zielgrössen gibt es die verallgemeinerte Form) |
| Datenstruktur | Personen in Klassen, Messungen in Personen, Filialen in Regionen |
| Was die Methode liefert | feste Effekte mit Standardfehlern, Varianz der Gruppen und des Rests, ICC |
| R-Funktion | lme4::lmer(y ~ x + (1 \| gruppe)) |
| Python-Funktion | smf.mixedlm("y ~ x", data, groups=data["gruppe"]).fit() |
| Verwandte Methoden | multiple lineare Regression, gepaarter t-Test, Varianzanalyse |
Gemischte Modelle sind zuständig, sobald Beobachtungen nicht unabhängig sind, weil sie aus derselben Gruppe stammen: Schülerinnen derselben Klasse, wiederholte Messungen an derselben Person, Verkäufe derselben Filiale. Innerhalb einer Gruppe ähneln sich die Werte, und eine gewöhnliche Regression tut so, als wäre jede Zeile neue, unabhängige Information.
| Situation | Stattdessen |
|---|---|
| Keine Gruppenstruktur, alle Beobachtungen unabhängig | multiple lineare Regression |
| Nur zwei Messungen je Person, keine Prädiktoren | gepaarter t-Test |
| Sehr wenige Gruppen (unter etwa fünf) | Gruppe als fester Faktor in der Regression |
| Binäre Zielgrösse mit Gruppen | verallgemeinertes gemischtes Modell, glmer(..., family = binomial) |
| Messungen in gleichmässigem zeitlichem Abstand, Vorhersage steht im Zentrum | Zeitreihenmodelle |
Jede Klasse bekommt ihren eigenen Achsenabschnitt, aber diese Achsenabschnitte werden nicht einzeln geschätzt, sondern als Streuung um einen gemeinsamen Mittelwert modelliert:
\[y_{ij} = \beta_0 + u_{j} + \beta_1 x_{ij} + \varepsilon_{ij}, \qquad u_j \sim N(0, \sigma_u^2), \quad \varepsilon_{ij} \sim N(0, \sigma^2)\]
| Bestandteil | Bedeutung |
|---|---|
| \(\beta_0, \beta_1\) | feste Effekte: gelten für alle, werden wie in der Regression gelesen |
| \(u_j\) | zufälliger Effekt der Klasse \(j\): wie weit ihr Niveau vom Durchschnitt abweicht |
| \(\sigma_u^2\) | Varianz zwischen den Klassen |
| \(\sigma^2\) | Varianz innerhalb der Klassen, der Rest |
| \(\text{ICC} = \sigma_u^2 / (\sigma_u^2 + \sigma^2)\) | Anteil der Gesamtstreuung, der auf die Klassen entfällt |
Kann auch die Steigung je Klasse verschieden sein, kommt ein zweiter zufälliger Effekt dazu, der zufällige Steigung heisst: (stunden | klasse).
| Voraussetzung | Was sie bedeutet | Prüfung | Folge bei Verletzung | Vorgehen |
|---|---|---|---|---|
| Gruppenstruktur richtig abgebildet | Jede Ebene, auf der sich Werte ähneln, ist im Modell | Erhebungsdesign | Standardfehler zu klein | fehlende Ebene ergänzen |
| Genügend Gruppen | Die Varianz zwischen Gruppen braucht Gruppen, um geschätzt zu werden | Anzahl Gruppen | Varianz auf null geschätzt, Warnung singular fit |
ab etwa zehn Gruppen verlässlich, sonst fester Faktor |
| Linearität, konstante Restvarianz | wie in der Regression | Residuen gegen angepasste Werte | verzerrte Effekte, falsche Unsicherheit | Regressionsdiagnostik |
| Annähernd normalverteilte zufällige Effekte und Residuen | QQ-Plot der Residuen und der geschätzten \(u_j\) | bei grossen Stichproben meist harmlos | robuste Standardfehler oder Bootstrap |
R, summary(lmer(...)):
| Abschnitt | Was er sagt |
|---|---|
REML criterion at convergence |
Anpassungsmass, nur zum Vergleich von Modellen mit denselben festen Effekten |
Random effects, Variance und Std.Dev. |
\(\sigma_u^2\) je zufälligem Effekt und \(\sigma^2\) als Residual |
Number of obs, groups |
Anzahl Zeilen und Anzahl Gruppen; entscheidend ist oft die zweite Zahl |
Fixed effects, Estimate, Std. Error, t value |
wie in der Regression; lme4 gibt bewusst keine p-Werte aus |
Python, MixedLM.fit().summary():
| Zeile | Was sie sagt | Unterschied zu R |
|---|---|---|
Scale |
\(\sigma^2\), die Restvarianz | in R die Zeile Residual |
Group Var |
\(\sigma_u^2\) | in R (Intercept) unter Random effects |
Coef., Std.Err., z, P>\|z\| |
feste Effekte mit Wald-Test | R zeigt nur t value, ohne p-Wert |
Converged |
ob die Schätzung konvergiert ist | in R eine Warnung, falls nicht |
Warum lme4 keinen p-Wert zeigt: Die Freiheitsgrade eines festen Effekts sind in gemischten Modellen nicht eindeutig definiert. Python rechnet mit der Normalverteilung, das ist bei vielen Gruppen eine gute Näherung und bei wenigen zu optimistisch.
anova() automatisch, Python nicht.singular fit ist keine Formalie. Es heisst, dass eine Varianz auf null oder eine Korrelation auf plus oder minus eins geschätzt wurde. Das Modell ist dann zu komplex für die Daten.[Zielgrösse] wurde mit einem linearen gemischten Modell mit zufälligem Achsenabschnitt für [Gruppe] ausgewertet ([Anzahl] Beobachtungen in [Anzahl] Gruppen, REML). Je [Einheit Prädiktor] steigt [Zielgrösse] um [Schätzwert] (SE = [SE], 95-Prozent-Konfidenzintervall [unten] bis [oben]). Die Gruppen unterscheiden sich mit einer Standardabweichung von [σu], ICC = [ICC].
Immer genannt werden: Struktur der zufälligen Effekte, Zahl der Gruppen, Schätzverfahren, feste Effekte mit Intervall und die Varianzkomponenten.
| Methode | Wann diese statt eines gemischten Modells |
|---|---|
| Multiple lineare Regression | Beobachtungen sind unabhängig |
| Regression mit Gruppe als festem Faktor | wenige Gruppen, und nur genau diese Gruppen interessieren |
| Gepaarter t-Test | genau zwei Messungen je Person, keine weiteren Prädiktoren |
| Varianzanalyse mit Messwiederholung | klassisches Design, balanciert, keine fehlenden Werte |
| Cluster-robuste Standardfehler | nur die Standardfehler sollen korrigiert werden, die Gruppenvarianz interessiert nicht |
Jedes Beispiel steht in einem eigenen Reiter.
Hängt die Punktzahl in einem Test von den Lernstunden ab, und welcher Anteil der Unterschiede geht auf die Klasse zurück? 180 Personen aus 12 Klassen, je 15.
c(zeilen = nrow(schule), klassen = nlevels(schule$klasse)) zeilen klassen
180 12 round(c(mittel = mean(schule$punkte), sd = sd(schule$punkte)), 2)mittel sd
61.25 11.04 boxplot(punkte ~ klasse, data = schule, xlab = "Klasse", ylab = "Punkte")
import matplotlib.pyplot as plt
print("Zeilen:", len(schule), "Klassen:", schule["klasse"].nunique())Zeilen: 180 Klassen: 12print("Mittel:", round(schule["punkte"].mean(), 2), "SD:", round(schule["punkte"].std(), 2))Mittel: 61.25 SD: 11.04fig, achse = plt.subplots(figsize=(7, 3.4))
_ = achse.boxplot([schule.loc[schule["klasse"] == k, "punkte"] for k in range(1, 13)])
_ = achse.set_xlabel("Klasse")
_ = achse.set_ylabel("Punkte")
plt.show()
Die Klassen liegen sichtbar auf verschiedenem Niveau. Das ist genau die Struktur, die ein gewöhnliches Regressionsmodell übersieht.
m0 <- lmer(punkte ~ stunden + (1 | klasse), data = schule)
summary(m0)Linear mixed model fit by REML ['lmerMod']
Formula: punkte ~ stunden + (1 | klasse)
Data: schule
REML criterion at convergence: 1297.3
Scaled residuals:
Min 1Q Median 3Q Max
-2.2868 -0.7036 0.0335 0.6510 3.5559
Random effects:
Groups Name Variance Std.Dev.
klasse (Intercept) 34.85 5.903
Residual 70.23 8.380
Number of obs: 180, groups: klasse, 12
Fixed effects:
Estimate Std. Error t value
(Intercept) 50.1142 2.3653 21.187
stunden 2.2244 0.3029 7.344
Correlation of Fixed Effects:
(Intr)
stunden -0.641vk <- as.data.frame(VarCorr(m0))
round(c(sd_klasse = vk$sdcor[1], sd_rest = vk$sdcor[2],
icc = vk$vcov[1] / sum(vk$vcov)), 3)sd_klasse sd_rest icc
5.903 8.380 0.332 round(confint(m0, parm = "beta_", method = "Wald"), 3) 2.5 % 97.5 %
(Intercept) 45.478 54.750
stunden 1.631 2.818m0 = smf.mixedlm("punkte ~ stunden", schule, groups=schule["klasse"]).fit(reml=True)
print(m0.summary())
var_klasse = float(m0.cov_re.iloc[0, 0])
var_rest = m0.scale
print("sd_klasse:", round(var_klasse ** 0.5, 3), "sd_rest:", round(var_rest ** 0.5, 3),
"icc:", round(var_klasse / (var_klasse + var_rest), 3))| Ausgabe | Wert hier | Bedeutung | Wie er zu lesen ist |
|---|---|---|---|
| Achsenabschnitt | 50.114 | erwartete Punktzahl bei null Lernstunden in einer durchschnittlichen Klasse | extrapoliert, niemand lernt null Stunden; nur als Rechengrösse |
stunden |
2.224 | je zusätzlicher Stunde rund 2.2 Punkte mehr | fester Effekt, gilt über alle Klassen |
Standardfehler stunden |
0.303 | Unsicherheit der Steigung | Das Wald-Intervall reicht von 1.631 bis 2.818. |
| Varianz Klasse | 34.849 | Streuung der Klassenniveaus, als Standardabweichung 5.903 Punkte | Eine Klasse eine Standardabweichung über dem Schnitt liegt knapp 6 Punkte höher. |
| Restvarianz | 70.227 | Streuung innerhalb der Klassen, als Standardabweichung 8.380 | in Python die Zeile Scale |
| ICC | 0.332 | ein Drittel der Streuung liegt zwischen den Klassen | weit entfernt von null; die Gruppenstruktur zu ignorieren wäre ein Fehler |
Die Lernstunden wirken klar, und zugleich unterscheiden sich die Klassen erheblich. Ein Drittel der Unterschiede zwischen Personen erklärt sich schon dadurch, in welche Klasse sie gehen.
Die Punktzahl wurde mit einem linearen gemischten Modell mit zufälligem Achsenabschnitt je Klasse ausgewertet (180 Personen in 12 Klassen, REML). Je Lernstunde steigt die Punktzahl um 2.224 Punkte (SE = 0.303, 95-Prozent-Konfidenzintervall 1.631 bis 2.818). Die Klassen unterscheiden sich mit einer Standardabweichung von 5.903 Punkten, ICC = 0.332.
Hat die Klassengrösse einen Einfluss auf die Punktzahl? Die Klassengrösse ist ein Merkmal der Klasse, nicht der Person: Alle 15 Personen einer Klasse haben denselben Wert. Die Daten wurden ohne jeden Effekt der Klassengrösse erzeugt.
ols <- lm(punkte ~ stunden + groesse, data = schule)
gem <- lmer(punkte ~ stunden + groesse + (1 | klasse), data = schule)
k_ols <- summary(ols)$coefficients["groesse", ]
k_gem <- summary(gem)$coefficients["groesse", ]
round(rbind(
ols = c(schaetzung = k_ols[[1]], se = k_ols[[2]], p = k_ols[[4]]),
gemischt = c(schaetzung = k_gem[[1]], se = k_gem[[2]],
p = 2 * pnorm(-abs(k_gem[[3]])))
), 4) schaetzung se p
ols 0.3481 0.4088 0.3956
gemischt 0.3610 1.0162 0.7224ols = smf.ols("punkte ~ stunden + groesse", schule).fit()
gem = smf.mixedlm("punkte ~ stunden + groesse", schule,
groups=schule["klasse"]).fit(reml=True)
print(pd.DataFrame({
"schaetzung": [ols.params["groesse"], gem.params["groesse"]],
"se": [ols.bse["groesse"], gem.bse["groesse"]],
"p": [ols.pvalues["groesse"], gem.pvalues["groesse"]],
}, index=["ols", "gemischt"]).round(4))| Modell | Schätzung | Standardfehler | p | Wie es zu lesen ist |
|---|---|---|---|---|
| OLS | 0.348 | 0.409 | 0.3956 | behandelt 180 Zeilen als unabhängig |
| gemischt | 0.361 | 1.016 | 0.7224 | berücksichtigt, dass es nur 12 verschiedene Klassengrössen gibt |
Die Schätzung ist fast gleich, der Standardfehler aber im gemischten Modell zweieinhalbmal so gross. Hier sind beide Ergebnisse nicht signifikant, und das ist richtig, weil die Daten ohne Effekt erzeugt wurden. Der zu kleine Standardfehler von OLS bleibt aber ein Fehler: Über viele solche Datensätze meldet OLS für ein Klassenmerkmal deutlich öfter einen Scheineffekt, als das Signifikanzniveau zulässt.
Der Unterschied betrifft vor allem Merkmale der Gruppe. Für die Lernstunden, die innerhalb jeder Klasse variieren, liegen die Standardfehler von OLS und gemischtem Modell näher beieinander.
Ein Effekt der Klassengrösse liess sich nicht nachweisen (Schätzung 0.361 Punkte je Person, SE = 1.016, p = 0.72, lineares gemischtes Modell mit zufälligem Achsenabschnitt je Klasse, 12 Klassen).
Wirkt eine Lernstunde in allen Klassen gleich, oder unterscheidet sich die Steigung zwischen den Klassen? Verglichen werden zwei Modelle mit einem Likelihood-Quotienten-Test, beide mit ML geschätzt.
m_int <- lmer(punkte ~ stunden + (1 | klasse), data = schule, REML = FALSE)
m_slope <- lmer(punkte ~ stunden + (stunden | klasse), data = schule, REML = FALSE)
anova(m_int, m_slope)Data: schule
Models:
m_int: punkte ~ stunden + (1 | klasse)
m_slope: punkte ~ stunden + (stunden | klasse)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
m_int 4 1307.7 1320.5 -649.84 1299.7
m_slope 6 1311.1 1330.2 -649.53 1299.1 0.6178 2 0.7342m_int = smf.mixedlm("punkte ~ stunden", schule, groups=schule["klasse"]).fit(reml=False)
m_slope = smf.mixedlm("punkte ~ stunden", schule, groups=schule["klasse"],
re_formula="~stunden").fit(reml=False)
lr = 2 * (m_slope.llf - m_int.llf)
print("Chisq:", round(lr, 3), " Df: 2 p:", round(stats.chi2.sf(lr, 2), 4))| Ausgabe | Wert hier | Bedeutung | Wie er zu lesen ist |
|---|---|---|---|
Chisq |
0.618 | doppelte Differenz der Log-Likelihoods | wie viel besser das komplexere Modell die Daten beschreibt |
Df |
2 | zwei zusätzliche Parameter | die Varianz der Steigung und ihre Korrelation mit dem Achsenabschnitt |
Pr(>Chisq) |
0.7342 | p-Wert des Vergleichs | kein Hinweis, dass die Steigung zwischen den Klassen variiert |
AIC (nur R) |
Informationskriterium | kleiner ist besser; der Vergleich steht in der Tabelle |
Die Daten wurden mit einer leicht unterschiedlichen Steigung je Klasse erzeugt. Der Test findet sie trotzdem nicht: Zwölf Klassen reichen nicht, um eine Streuung von 0.8 Punkten je Stunde von null zu unterscheiden. Das ist Interpretationsfalle drei in der Praxis. Das einfachere Modell ist hier die vertretbare Wahl, ein Beweis für gleiche Steigungen ist das nicht.
Ein Modell mit zufälliger Steigung für die Lernstunden verbesserte die Anpassung nicht (χ²(2) = 0.62, p = 0.73); berichtet wird das Modell mit zufälligem Achsenabschnitt.
Ein Datensatz enthält 2000 Messungen von 20 Filialen. Ein Merkmal der Filiale (Verkaufsfläche) soll den Umsatz erklären. Wie viele unabhängige Beobachtungen stehen für diesen Effekt ungefähr zur Verfügung?
lme4 meldet boundary (singular) fit. Was ist die sinnvollste Reaktion?
Der ICC eines Modells beträgt zwei Prozent. Was bedeutet das?