Regressionsdiagnostik

Residuen, Hebelwerte, Ausreisser, VIF und was man bei Verletzungen tut.

Die vier Beispiele rechnen mit vier deterministisch konstruierten Datensätzen. Sie unterscheiden sich in genau einer Eigenschaft, damit sichtbar wird, welche Diagnose worauf anspricht.

Kurzsteckbrief

Fragestellung	Sind die Voraussetzungen eines linearen Modells erfüllt, und wenn nicht, welche und mit welchen Folgen?
Eingang	ein angepasstes Modell aus lm() beziehungsweise smf.ols()
Was die Diagnostik liefert	Residuenplots, Hebelwerte, Cook-Distanzen, studentisierte Residuen, Testverfahren für einzelne Voraussetzungen
R-Funktionen	plot(modell), lmtest::bptest(), lmtest::dwtest(), lmtest::resettest(), cooks.distance(), hatvalues(), rstudent(), car::vif()
Python-Funktionen	het_breuschpagan(), durbin_watson(), linear_reset(), modell.get_influence()
Verwandte Methoden	einfache lineare Regression, multiple lineare Regression, Modellauswahl

Wann diese Methode, wann nicht

Diagnostik ist kein eigener Arbeitsschritt am Ende, sondern gehört zu jedem angepassten Modell. Die Reihenfolge lautet: Modell anpassen, Residuen ansehen, dann Koeffizienten lesen. Wer die Koeffizienten zuerst liest, ist geneigt, die Diagnostik so lange zu deuten, bis sie zum Ergebnis passt.

Sie beantwortet nicht die Frage, ob das Modell inhaltlich sinnvoll ist. Ein Modell kann alle Voraussetzungen erfüllen und trotzdem die falsche Frage beantworten oder eine Scheinbeziehung abbilden.

Situation	Stattdessen oder zusätzlich
Zielgrösse binär	logistische Regression; dort gelten andere Voraussetzungen
Es geht um die Auswahl zwischen Modellen	Modellauswahl
Die Beobachtungen sind zeitlich geordnet	Zeitreihenverfahren, siehe Stationarität
Prädiktoren messen dasselbe	Multikollinearität, siehe multiple lineare Regression

Die fünf Fragen der Diagnostik

Frage	Voraussetzung	Womit geprüft	Folge bei Verletzung
Ist der Zusammenhang geradlinig?	Linearität	Residuen gegen Vorhersage, RESET-Test	Die Steigung ist ein Kompromiss und stimmt in keinem Bereich
Streuen die Residuen überall gleich?	Homoskedastizität	Residuenplot, Scale-Location-Plot, Breusch-Pagan	Schätzwerte bleiben brauchbar, Standardfehler und p-Werte nicht
Sind die Residuen annähernd normal?	Normalität der Residuen	QQ-Plot	betrifft Tests und Intervalle, vor allem bei kleinem n
Beeinflussen sich die Beobachtungen?	Unabhängigkeit	Durbin-Watson, Erhebungsdesign	Standardfehler zu klein, p-Werte zu optimistisch
Hängt das Ergebnis an wenigen Punkten?	keine einflussreichen Einzelfälle	Cook-Distanz, Hebelwerte, studentisierte Residuen	Die Aussage hängt an einzelnen Beobachtungen

Nur die erste und die fünfte Frage betreffen die Schätzwerte selbst. Die zweite, dritte und vierte betreffen die Unsicherheit, also Standardfehler, Konfidenzintervalle und p-Werte. Diese Unterscheidung entscheidet darüber, wie dringend eine Verletzung behandelt werden muss.

Die vier Standardplots

plot(modell) in R liefert vier Bilder. Sie beantworten der Reihe nach:

Nummer	Titel	Was er zeigt	Was gut aussieht
1	Residuals vs Fitted	Linearität	strukturlose Wolke um die Nulllinie, keine Krümmung
2	Q-Q Residuals	Normalität	Punkte nah an der Geraden, Abweichungen an den Enden sind üblich
3	Scale-Location	Homoskedastizität	waagrechte Glättungslinie, gleichmässige Streubreite
4	Residuals vs Leverage	einflussreiche Punkte	keine Punkte jenseits der Cook-Konturen

Python hat kein Gegenstück in einem Aufruf; die Plots werden einzeln gebaut.

Output lesen

Kennzahl	Was sie misst	Anhaltspunkt	Falle
Cook-Distanz	wie stark sich alle Vorhersagen ändern, wenn dieser Punkt fehlt	Werte über 1 gelten als auffällig, über 0.5 als beachtenswert	Eine einzelne grosse Cook-Distanz ist ein Befund, viele mittlere sind normal
Hebelwert (leverage)	wie extrem der x-Wert liegt	mehr als das Doppelte oder Dreifache des Mittelwerts \(p/n\)	Hoher Hebel allein schadet nicht, siehe Beispiel 4
studentisiertes Residuum	wie weit der Punkt von seiner eigenen Vorhersage abweicht, in Standardfehlern	Beträge über 3 sind auffällig	misst die y-Richtung, der Hebel die x-Richtung
Breusch-Pagan p	Hinweis auf ungleiche Streuung	grosser p-Wert ist die gute Nachricht	findet nur Streuung, die mit den Prädiktoren zusammenhängt
Durbin-Watson	Korrelation aufeinanderfolgender Residuen	Werte um 2 sind unauffällig, unter 1.5 auffällig	hängt an der Reihenfolge der Zeilen, siehe Beispiel 2
RESET p	Hinweis auf fehlende Nichtlinearität	kleiner p-Wert heisst: das Modell ist fehlspezifiziert	sagt nicht, welcher Prädiktor betroffen ist
VIF	Multikollinearität	über 5 beachten, über 10 kritisch	nur bei mehreren Prädiktoren sinnvoll

Alle diese Tests haben dasselbe Problem wie Normalitätstests: Ihre Aussagekraft hängt am Stichprobenumfang. Bei n = 5000 wird jede belanglose Abweichung signifikant, bei n = 20 bleibt jede grobe unentdeckt. Sie ergänzen die Plots, sie ersetzen sie nicht.

Interpretationsfallen

• Ein hohes R² sagt nichts über die Angemessenheit. In Beispiel 2 ist das falsche Modell dasjenige mit dem höheren R².
• Nicht jede Verletzung ist gleich schlimm. Ungleiche Streuung verzerrt die Standardfehler, aber nicht die Schätzwerte. Eine übersehene Krümmung verzerrt beides.
• Hoher Hebel ist kein Fehler. Ein Punkt weit aussen auf der x-Achse ist wertvoll, solange er auf der Linie der übrigen liegt. Gefährlich ist die Kombination aus hohem Hebel und grossem Residuum.
• Ausreisser werden geprüft, nicht entfernt. Das Entfernen ist eine inhaltliche Entscheidung mit Begründung, keine statistische Reinigung.
• Durbin-Watson misst die Zeilenreihenfolge. Bei nach x sortierten Daten zeigt er Krümmung an, nicht Autokorrelation.
• Die Residuen zählen, nicht die Rohdaten. Weder x noch y müssen normalverteilt sein.
• Eine Transformation ändert die Fragestellung. Nach dem Logarithmieren der Zielgrösse beschreibt die Steigung eine relative statt einer absoluten Änderung.

Ergebnis berichten

Die Diagnostik gehört in einen eigenen kurzen Absatz, und zwar vor den Koeffizienten:

Residuen gegen Vorhersage zeigen keine Struktur, der QQ-Plot keine nennenswerte Abweichung, der Breusch-Pagan-Test keinen Hinweis auf ungleiche Streuung (p = [p]). Die grösste Cook-Distanz beträgt [Wert]; einflussreiche Einzelfälle liegen nicht vor.

Bei einer Verletzung gehört dazu, was daraus folgte: welche Kennzahl auffällig war, welche Massnahme gewählt wurde und ob sich das Ergebnis dadurch änderte.

Abgrenzung zu verwandten Methoden

Methode	Wann diese statt der Diagnostik
Modellauswahl	die Frage lautet, welche Prädiktoren ins Modell gehören
Multiple lineare Regression	die Schätzung selbst steht im Vordergrund
Normalitätstests	Vertiefung zur Frage der Verteilungsform
Bootstrap	Intervalle ohne Verteilungsannahme, wenn die Voraussetzungen nicht zu retten sind

Beispiele

Beispiel 1: Das unauffällige Modell

Frage und Datenlage

120 Beobachtungen, ein Prädiktor, ein sauber linearer Zusammenhang mit gleichmässiger Streuung. Dieses Beispiel ist der Referenzfall: So sehen die Plots aus, wenn nichts vorliegt.

Rechnung

R

modell_a <- lm(y ~ x, data = sauber)

par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
plot(modell_a)
par(mfrow = c(1, 1))

bptest(modell_a)

    studentized Breusch-Pagan test

data:  modell_a
BP = 1.008, df = 1, p-value = 0.3154

dwtest(modell_a)

    Durbin-Watson test

data:  modell_a
DW = 2.1231, p-value = 0.7213
alternative hypothesis: true autocorrelation is greater than 0

resettest(modell_a, power = 2, type = "fitted")

    RESET test

data:  modell_a
RESET = 1.1344, df1 = 1, df2 = 117, p-value = 0.289

shapiro.test(residuals(modell_a))

    Shapiro-Wilk normality test

data:  residuals(modell_a)
W = 0.99851, p-value = 1

round(c(r2 = summary(modell_a)$r.squared,
        residualstreuung = summary(modell_a)$sigma,
        max_cook = max(cooks.distance(modell_a)),
        max_hebel = max(hatvalues(modell_a)),
        mittlerer_hebel = mean(hatvalues(modell_a)),
        max_stud_residuum = rstudent(modell_a)[which.max(abs(rstudent(modell_a)))]), 4)

                 r2    residualstreuung            max_cook           max_hebel 
             0.9023              4.0097              0.2468              0.0669 
    mittlerer_hebel max_stud_residuum.1 
             0.0167             -2.6918 

Abbildung 1: Beispiel 1: die vier Standardplots eines unauffälligen Modells.

Python

modell_a = smf.ols("y ~ x", data=sauber).fit()


def vier_plots(modell, titel):
    einfluss_obj = modell.get_influence()
    wurzel_resid = np.sqrt(np.abs(einfluss_obj.resid_studentized_internal))

    fig, achsen = plt.subplots(2, 2, figsize=(7, 5.4))
    achsen[0, 0].scatter(modell.fittedvalues, modell.resid, s=12, color="0.3")
    achsen[0, 0].axhline(0, linestyle="--", color="black")
    achsen[0, 0].set_title("Residuen gegen Vorhersage", fontsize=9)
    sm.qqplot(modell.resid, line="s", ax=achsen[0, 1])
    achsen[0, 1].set_title("QQ-Plot der Residuen", fontsize=9)
    achsen[1, 0].scatter(modell.fittedvalues, wurzel_resid, s=12, color="0.3")
    achsen[1, 0].set_title("Scale-Location", fontsize=9)
    achsen[1, 1].scatter(einfluss_obj.hat_matrix_diag,
                         einfluss_obj.resid_studentized_internal,
                         s=12, color="0.3")
    achsen[1, 1].axhline(0, linestyle="--", color="black")
    achsen[1, 1].set_xlabel("Hebelwert")
    achsen[1, 1].set_title("Residuen gegen Hebelwert", fontsize=9)
    fig.suptitle(titel, fontsize=10)
    plt.tight_layout()
    plt.show()


vier_plots(modell_a, "Beispiel 1")

Abbildung 2: Beispiel 1: dieselben Bilder in Python.

bp = het_breuschpagan(modell_a.resid, modell_a.model.exog)
print("Breusch-Pagan LM =", round(bp[0], 4), " p =", round(bp[1], 4))

Breusch-Pagan LM = 1.008  p = 0.3154

print("Durbin-Watson =", round(durbin_watson(modell_a.resid), 4))

Durbin-Watson = 2.1231

reset = linear_reset(modell_a, power=2, use_f=True)
print("RESET F =", round(reset.fvalue, 4), " p =", round(reset.pvalue, 4))

RESET F = 1.1344  p = 0.289

print("Shapiro der Residuen p =", round(stats.shapiro(modell_a.resid).pvalue, 5))

Shapiro der Residuen p = 0.99999

einfluss_a = modell_a.get_influence()
stud_a = einfluss_a.resid_studentized_external
print("R^2 =", round(modell_a.rsquared, 4),
      " Residualstreuung =", round(np.sqrt(modell_a.mse_resid), 4))

R^2 = 0.9023  Residualstreuung = 4.0097

print("max Cook =", round(einfluss_a.cooks_distance[0].max(), 4),
      " max Hebel =", round(einfluss_a.hat_matrix_diag.max(), 4),
      " mittlerer Hebel =", round(einfluss_a.hat_matrix_diag.mean(), 4))

max Cook = 0.2468  max Hebel = 0.0669  mittlerer Hebel = 0.0167

print("grösstes studentisiertes Residuum =",
      round(stud_a[np.abs(stud_a).argmax()], 4))

grösstes studentisiertes Residuum = -2.6918

Output Zeile für Zeile

Kennzahl	Wert hier	Wie er zu lesen ist
R²	0.9023	Merken für den Vergleich mit Beispiel 2.
Residualstreuung	4.01	Entspricht der eingebauten Streuung von 4.
Breusch-Pagan p	0.315	Kein Hinweis auf ungleiche Streuung. Der grosse p-Wert ist hier die gute Nachricht.
Durbin-Watson	2.123	Um 2 herum, also unauffällig.
RESET p	0.289	Kein Hinweis auf eine fehlende Nichtlinearität.
Shapiro der Residuen p	praktisch 1 (R gibt 1 aus, Python 0.99999)	Unauffällig, und derart nahe an 1, weil die Störgrösse deterministisch aus Normalverteilungsquantilen gebaut ist.
grösste Cook-Distanz	0.247	Weit unter der Schwelle von 1.
grösster Hebelwert	0.0669	Der mittlere Hebel beträgt 0.0167, also 2 geteilt durch 120. Der grösste liegt beim Vierfachen, was bei 120 Punkten am Rand des Wertebereichs normal ist.
grösstes studentisiertes Residuum	-2.69	Unter 3, also unauffällig. Bei 120 Beobachtungen erwartet man einige Werte um 2.5.

Interpretation

Alle vier Plots zeigen das erwartete Bild: strukturlose Wolke, Punkte auf der QQ-Geraden, waagrechte Scale-Location-Linie, keine Punkte ausserhalb der Cook-Konturen. Das Modell darf gelesen werden.

Wichtig ist, dieses Bild einmal gesehen zu haben. Wer nur Verletzungen kennt, neigt dazu, in jedem Residuenplot ein Muster zu erkennen.

Beispiel 2: Übersehene Krümmung

Frage und Datenlage

Derselbe Prädiktor, aber der Zusammenhang ist in Wahrheit quadratisch. Angepasst wird trotzdem eine Gerade.

Rechnung

R

modell_b <- lm(y ~ x, data = krumm)

par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
plot(krumm$x, krumm$y, pch = 16, col = "grey40", xlab = "x", ylab = "y",
     main = "Daten und Gerade")
abline(modell_b, col = "steelblue", lwd = 2)
plot(modell_b, which = 1)
par(mfrow = c(1, 1))

round(summary(modell_b)$r.squared, 4)

[1] 0.9214

bptest(modell_b)

    studentized Breusch-Pagan test

data:  modell_b
BP = 0.40991, df = 1, p-value = 0.522

dwtest(modell_b)

    Durbin-Watson test

data:  modell_b
DW = 0.52853, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is greater than 0

resettest(modell_b, power = 2, type = "fitted")

    RESET test

data:  modell_b
RESET = 382.52, df1 = 1, df2 = 117, p-value < 2.2e-16

shapiro.test(residuals(modell_b))

    Shapiro-Wilk normality test

data:  residuals(modell_b)
W = 0.90307, p-value = 2.827e-07

round(summary(modell_b)$sigma, 4)

[1] 8.2458

Abbildung 3: Beispiel 2: Gerade auf gekrümmten Daten und der Bogen im Residuenplot.

Python

modell_b = smf.ols("y ~ x", data=krumm).fit()

fig, achsen = plt.subplots(1, 2, figsize=(7, 3.2))
achsen[0].scatter(krumm["x"], krumm["y"], s=12, color="0.4")
gitter = np.linspace(krumm["x"].min(), krumm["x"].max(), 50)
achsen[0].plot(gitter, modell_b.params.iloc[0] + modell_b.params.iloc[1] * gitter,
               color="steelblue", linewidth=2)
achsen[0].set_title("Daten und Gerade", fontsize=9)
achsen[1].scatter(modell_b.fittedvalues, modell_b.resid, s=12, color="0.3")
achsen[1].axhline(0, linestyle="--", color="black")
glatt = sm.nonparametric.lowess(modell_b.resid, modell_b.fittedvalues)
achsen[1].plot(glatt[:, 0], glatt[:, 1], color="firebrick", linewidth=2)
achsen[1].set_title("Residuen gegen Vorhersage", fontsize=9)
plt.tight_layout()
plt.show()

Abbildung 4: Beispiel 2: dieselbe Gegenüberstellung in Python.

print("R^2 =", round(modell_b.rsquared, 4))

R^2 = 0.9214

bp = het_breuschpagan(modell_b.resid, modell_b.model.exog)
print("Breusch-Pagan p =", round(bp[1], 4))

Breusch-Pagan p = 0.522

print("Durbin-Watson =", round(durbin_watson(modell_b.resid), 4))

Durbin-Watson = 0.5285

print("RESET p =", linear_reset(modell_b, power=2, use_f=True).pvalue)

RESET p = 1.1155488532642536e-38

print("Shapiro der Residuen p =", stats.shapiro(modell_b.resid).pvalue)

Shapiro der Residuen p = 2.827160354309601e-07

print("Residualstreuung =", round(np.sqrt(modell_b.mse_resid), 4))

Residualstreuung = 8.2458

Output Zeile für Zeile

Kennzahl	Beispiel 1	Beispiel 2	Wie das zu lesen ist
R²	0.9023	0.9214	Das falsche Modell hat das höhere R². Wer Modelle nach R² beurteilt, wählt hier die falsche Spezifikation.
Residualstreuung	4.01	8.25	Doppelt so gross. Diese Zahl verrät das Problem, R² nicht.
Breusch-Pagan p	0.315	0.522	Unauffällig. Der Test sucht ungleiche Streuung und findet keine, denn die Streuung ist ja gleichmässig. Er ist für dieses Problem blind.
Durbin-Watson	2.123	0.529	Deutlich auffällig, und zwar aus einem Grund, der Aufmerksamkeit verdient, siehe Kasten.
RESET p	0.289	1.1e-38 (R schreibt < 2.2e-16)	Der Test, der genau für diesen Fall gebaut ist. Er prüft, ob quadrierte Vorhersagewerte noch etwas erklären, und findet massiv etwas.
Shapiro der Residuen p	0.99999	2.8e-07	Ebenfalls auffällig, aber irreführend: Die Residuen sind nicht deshalb nicht normalverteilt, weil die Störgrösse es nicht wäre, sondern weil sie den systematischen Bogen enthalten.

Warum Durbin-Watson hier anschlägt, obwohl es keine Zeitreihe ist. Der Test misst, ob aufeinanderfolgende Residuen einander ähneln, und “aufeinanderfolgend” heisst: in der Reihenfolge der Zeilen. In diesen Daten sind die Zeilen nach x sortiert. Der Bogen im Residuenplot führt dann dazu, dass benachbarte Zeilen ähnliche Residuen haben, und genau das misst Durbin-Watson.

Praktisch heisst das zweierlei. Bei nach einem Prädiktor sortierten Daten ist ein kleiner Durbin-Watson-Wert ein Hinweis auf Krümmung, nicht auf Autokorrelation. Und umgekehrt: Wer echte Autokorrelation prüfen will, muss die Zeilen in ihrer zeitlichen Reihenfolge haben, sonst misst der Test etwas anderes.

Interpretation und Ergebnissatz

Die Lehre dieses Beispiels ist die Arbeitsteilung der Diagnosen: Breusch-Pagan schweigt, RESET schlägt Alarm, Durbin-Watson schlägt aus einem anderen Grund Alarm als vermutet, und R² zeigt in die falsche Richtung. Kein einzelnes Mass deckt alle Verletzungen ab.

Der Residuenplot zeigt alles auf einen Blick: Die Residuen sind an den Rändern positiv und in der Mitte negativ, die Glättungslinie bildet einen deutlichen Bogen.

Das Vorgehen: einen quadratischen Term aufnehmen, transformieren oder ein anderes Modell wählen, und danach die Residuen erneut ansehen.

Ein rein lineares Modell beschreibt den Zusammenhang nicht angemessen. Der RESET-Test ist deutlich signifikant (F(1, 117) = 382.5, p < 0.001), und die Residuen zeigen einen ausgeprägten Bogen. Das höhere R² gegenüber dem korrekt spezifizierten Modell (0.921 gegenüber 0.902) ist kein Gütezeichen. Berichtet wird das Modell mit quadratischem Term.

Beispiel 3: Ungleiche Streuung

Frage und Datenlage

Der Zusammenhang ist geradlinig, aber die Streuung wächst mit x: Bei kleinen x-Werten liegen die Punkte eng an der Geraden, bei grossen weit auseinander. Das ist der häufigste Fall in wirtschaftlichen Daten, wo Beträge und ihre Schwankung gemeinsam wachsen.

Rechnung

R

modell_c <- lm(y ~ x, data = trichter)

par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
plot(modell_c, which = 1)
plot(modell_c, which = 3)
par(mfrow = c(1, 1))

bptest(modell_c)

    studentized Breusch-Pagan test

data:  modell_c
BP = 16.283, df = 1, p-value = 5.455e-05

resettest(modell_c, power = 2, type = "fitted")

    RESET test

data:  modell_c
RESET = 0.0068298, df1 = 1, df2 = 117, p-value = 0.9343

shapiro.test(residuals(modell_c))

    Shapiro-Wilk normality test

data:  residuals(modell_c)
W = 0.99093, p-value = 0.6198

summary(modell_c)$coefficients

            Estimate Std. Error  t value     Pr(>|t|)
(Intercept) 20.30292  1.4183213 14.31475 1.468395e-27
x            2.97770  0.1317836 22.59538 1.117751e-44

Abbildung 5: Beispiel 3: der Trichter im Residuenplot.

Python

modell_c = smf.ols("y ~ x", data=trichter).fit()

fig, achsen = plt.subplots(1, 2, figsize=(7, 3.2))
achsen[0].scatter(modell_c.fittedvalues, modell_c.resid, s=12, color="0.3")
achsen[0].axhline(0, linestyle="--", color="black")
achsen[0].set_title("Residuen gegen Vorhersage", fontsize=9)
wurzel = np.sqrt(np.abs(modell_c.get_influence().resid_studentized_internal))
achsen[1].scatter(modell_c.fittedvalues, wurzel, s=12, color="0.3")
achsen[1].set_title("Scale-Location", fontsize=9)
plt.tight_layout()
plt.show()

Abbildung 6: Beispiel 3: dieselben Bilder in Python.

bp = het_breuschpagan(modell_c.resid, modell_c.model.exog)
print("Breusch-Pagan LM =", round(bp[0], 4), " p =", bp[1])

Breusch-Pagan LM = 16.283  p = 5.455243919507552e-05

print("RESET p =", round(linear_reset(modell_c, power=2, use_f=True).pvalue, 4))

RESET p = 0.9343

print("Shapiro der Residuen p =", round(stats.shapiro(modell_c.resid).pvalue, 4))

Shapiro der Residuen p = 0.6198

# Robuste Standardfehler: derselbe Schaetzwert, andere Unsicherheit
robust = smf.ols("y ~ x", data=trichter).fit(cov_type="HC3")
print(pd.DataFrame({
    "Schätzwert": modell_c.params.round(4),
    "SE gewöhnlich": modell_c.bse.round(4),
    "SE robust (HC3)": robust.bse.round(4),
}))

           Schätzwert  SE gewöhnlich  SE robust (HC3)
Intercept     20.3029         1.4183           1.1505
x              2.9777         0.1318           0.1389

Output Zeile für Zeile

Kennzahl	Wert hier	Wie er zu lesen ist
Breusch-Pagan p	0.0000546	Auffällig. Der Test, der für diesen Fall gebaut ist, findet ihn.
RESET p	0.934	Unauffällig, korrekt: Der Zusammenhang ist ja geradlinig.
Shapiro der Residuen p	0.620	Ebenfalls unauffällig. Die Residuen sind normalverteilt, nur nicht gleich stark streuend.
Steigung	2.978	Praktisch unverändert gegenüber Beispiel 1 (3.037).
Standardfehler gewöhnlich	0.132
Standardfehler robust (HC3)	siehe Ausgabe	Die interessante Gegenüberstellung: Der Schätzwert bleibt derselbe, nur die Unsicherheit wird anders beziffert.

Der Residuenplot zeigt den namensgebenden Trichter: Die Punkte fächern nach rechts auf. Der Scale-Location-Plot macht dasselbe deutlicher, weil er die Beträge aufträgt und die Glättungslinie dann sichtbar ansteigt statt waagrecht zu verlaufen.

Interpretation und Ergebnissatz

Ungleiche Streuung ist die harmloseste der Verletzungen, weil sie die Schätzwerte nicht verzerrt. Die Steigung von 2.98 ist so brauchbar wie im sauberen Fall. Verzerrt sind nur die Standardfehler, und damit Konfidenzintervalle und p-Werte.

Drei Wege stehen offen:

Weg	Wann er passt
Robuste Standardfehler (HC3)	der Standardweg; Schätzwerte bleiben, nur die Unsicherheit wird korrekt beziffert
Zielgrösse logarithmieren	wenn die Streuung proportional zum Niveau wächst; ändert aber die Bedeutung der Steigung
Gewichtete Regression	wenn die Streuungsstruktur bekannt ist

In R braucht der erste Weg das Zusatzpaket sandwich zusammen mit lmtest::coeftest(); in Python genügt cov_type="HC3" beim Anpassen.

Der Breusch-Pagan-Test zeigt ungleiche Streuung (LM = 16.28, p < 0.001), und der Residuenplot fächert nach rechts auf. Die Schätzwerte sind davon nicht betroffen; berichtet werden robuste Standardfehler nach HC3. Linearität (RESET p = 0.93) und Normalität der Residuen (Shapiro p = 0.62) sind erfüllt.

Beispiel 4: Ein einflussreicher Punkt

Frage und Datenlage

Die sauberen Daten aus Beispiel 1, ergänzt um eine einzige Beobachtung: x = 32, also weit rechts ausserhalb des beobachteten Bereichs von rund 2 bis 18, und y = 20, also weit unterhalb der Geraden.

Rechnung

R

modell_d <- lm(y ~ x, data = einfluss)

par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
plot(einfluss$x, einfluss$y, pch = 16, col = "grey40", xlab = "x", ylab = "y",
     main = "Mit und ohne den Fall")
points(32, 20, pch = 16, col = "firebrick", cex = 1.4)
abline(modell_d, col = "firebrick", lwd = 2)
abline(modell_a, col = "steelblue", lwd = 2, lty = 2)
legend("topleft", c("mit", "ohne"), col = c("firebrick", "steelblue"),
       lty = c(1, 2), lwd = 2, bty = "n", cex = 0.8)
plot(cooks.distance(modell_d), type = "h", lwd = 2, xlab = "Beobachtung",
     ylab = "Cook-Distanz", main = "Cook-Distanzen")
abline(h = 1, lty = 2, col = "firebrick")
par(mfrow = c(1, 1))

round(c(steigung_ohne = coef(modell_a)[2], steigung_mit = coef(modell_d)[2],
        r2_ohne = summary(modell_a)$r.squared,
        r2_mit = summary(modell_d)$r.squared), 4)

steigung_ohne.x  steigung_mit.x         r2_ohne          r2_mit 
         3.0369          2.1493          0.9023          0.5413 

round(c(max_cook = max(cooks.distance(modell_d)),
        max_hebel = max(hatvalues(modell_d)),
        mittlerer_hebel = mean(hatvalues(modell_d)),
        groesstes_stud_residuum = rstudent(modell_d)[which.max(abs(rstudent(modell_d)))]), 4)

                   max_cook                   max_hebel 
                    12.4670                      0.2083 
            mittlerer_hebel groesstes_stud_residuum.121 
                     0.0165                    -21.4831 

Abbildung 7: Beispiel 4: der Punkt, die beiden Geraden und die Cook-Distanzen.

Python

modell_d = smf.ols("y ~ x", data=einfluss).fit()

fig, achsen = plt.subplots(1, 2, figsize=(7, 3.2))
achsen[0].scatter(einfluss["x"], einfluss["y"], s=12, color="0.4")
achsen[0].scatter([32], [20], s=45, color="firebrick")
gitter = np.linspace(2, 33, 50)
achsen[0].plot(gitter, modell_d.params.iloc[0] + modell_d.params.iloc[1] * gitter,
               color="firebrick", linewidth=2, label="mit")
achsen[0].plot(gitter, modell_a.params.iloc[0] + modell_a.params.iloc[1] * gitter,
               color="steelblue", linewidth=2, linestyle="--", label="ohne")
achsen[0].legend(frameon=False, fontsize=8)
achsen[0].set_title("Mit und ohne den Fall", fontsize=9)
cook = modell_d.get_influence().cooks_distance[0]
achsen[1].vlines(range(len(cook)), 0, cook, linewidth=2)
achsen[1].axhline(1, linestyle="--", color="firebrick")
achsen[1].set_xlabel("Beobachtung")
achsen[1].set_ylabel("Cook-Distanz")
achsen[1].set_title("Cook-Distanzen", fontsize=9)
plt.tight_layout()
plt.show()

Abbildung 8: Beispiel 4: dieselbe Gegenüberstellung in Python.

einfluss_obj = modell_d.get_influence()
print("Steigung ohne:", round(modell_a.params.iloc[1], 4),
      " mit:", round(modell_d.params.iloc[1], 4))

Steigung ohne: 3.0369  mit: 2.1493

print("R^2 ohne:", round(modell_a.rsquared, 4),
      " mit:", round(modell_d.rsquared, 4))

R^2 ohne: 0.9023  mit: 0.5413

print("max Cook:", round(cook.max(), 4),
      " max Hebel:", round(einfluss_obj.hat_matrix_diag.max(), 4),
      " mittlerer Hebel:", round(einfluss_obj.hat_matrix_diag.mean(), 4))

max Cook: 12.467  max Hebel: 0.2083  mittlerer Hebel: 0.0165

stud = einfluss_obj.resid_studentized_external
print("grösstes studentisiertes Residuum:",
      round(stud[np.abs(stud).argmax()], 4))

grösstes studentisiertes Residuum: -21.4831

print("Breusch-Pagan p =", het_breuschpagan(modell_d.resid, modell_d.model.exog)[1])

Breusch-Pagan p = 9.564024480342477e-07

print("RESET p =", linear_reset(modell_d, power=2, use_f=True).pvalue)

RESET p = 9.873700108167648e-28

Output Zeile für Zeile

Kennzahl	Ohne den Fall	Mit dem Fall	Wie das zu lesen ist
Steigung	3.037	2.149	Ein Punkt von 121 senkt die Steigung um knapp ein Drittel.
R²	0.9023	0.5413	Von 90 auf 54 Prozent.
Residualstreuung	4.01	8.85	Mehr als verdoppelt.
grösste Cook-Distanz	0.247	12.467	Die Schwelle liegt bei 1. Der Wert liegt eine Grössenordnung darüber, und die übrigen 120 Werte verschwinden im Diagramm optisch auf der Nulllinie.
grösster Hebelwert	0.0669	0.2083	Der mittlere Hebel beträgt 0.0165; dieser Punkt hat den 12.6-fachen Wert.
grösstes studentisiertes Residuum	-2.69	-21.48	Der Punkt liegt 21 Standardfehler unter seiner eigenen Vorhersage.
Breusch-Pagan p	0.315	0.0000010	Auffällig geworden, ohne dass sich an der Streuungsstruktur etwas geändert hätte.
RESET p	0.289	9.9e-28	Ebenfalls.

Ein einzelner Punkt lässt alle Tests anschlagen. Breusch-Pagan, RESET und Shapiro sind in Beispiel 1 unauffällig und hier alle drei hochsignifikant, obwohl sich an Streuungsstruktur, Linearität und Verteilung nichts geändert hat.

Daraus folgt eine Reihenfolge für die Diagnostik: zuerst nach einflussreichen Punkten sehen, dann die übrigen Tests lesen. Andernfalls behandelt man eine Krümmung oder eine Streuungsstruktur, die es gar nicht gibt.

Interpretation und Ergebnissatz

Hebel und Cook-Distanz messen Verschiedenes, und beide werden hier gebraucht. Der Hebel von 0.21 sagt: Der Punkt liegt weit aussen auf der x-Achse. Für sich genommen wäre das unproblematisch; ein weit aussen liegender Punkt auf der Linie der übrigen macht die Schätzung sogar genauer. Das studentisierte Residuum von -21.5 sagt: Er liegt zugleich weit von seiner Vorhersage entfernt. Erst die Kombination ergibt die Cook-Distanz von 12.5.

Was jetzt zu tun ist, ist keine statistische Frage:

Befund zur Ursache	Vorgehen
Erfassungsfehler	ausschliessen, Ausschluss und Begründung berichten
Gültig, aber ausserhalb der Fragestellung	ausschliessen, Grundgesamtheit enger fassen
Gültig, Ursache unklar	Ergebnis mit und ohne den Fall berichten

Eine von 121 Beobachtungen ist stark einflussreich (Cook-Distanz 12.5, Hebelwert 0.21 gegenüber einem mittleren Hebel von 0.017, studentisiertes Residuum -21.5). Sie senkt die Steigung von 3.04 auf 2.15 und R² von 0.90 auf 0.54. Zugleich werden Breusch-Pagan, RESET und der Shapiro-Test auf die Residuen signifikant, obwohl die übrigen 120 Beobachtungen unauffällig sind. Der Fall wurde als Messfehler ausgeschlossen; berichtet wird das Modell ohne ihn.

Verständnisfragen

Ein Modell hat R² = 0.92, ein konkurrierendes mit anderer Spezifikation R² = 0.90. Der Residuenplot des ersten zeigt einen Bogen. Welches ist besser?

Das zweite, trotz kleinerem R²

Richtig. R² misst erklärte Streuung, nicht Angemessenheit. Ein Modell mit übersehener Krümmung kann ein höheres R² haben als das korrekt spezifizierte, wie Beispiel 2 zeigt. Der Residuenplot und die Residualstreuung entscheiden.

Das erste, R² ist das Gütemass

Dann würde man systematisch die falsche Spezifikation wählen.

Nicht entscheidbar ohne p-Werte

Die p-Werte sind bei einem fehlspezifizierten Modell ohnehin nicht vertrauenswürdig.

Der Breusch-Pagan-Test ist unauffällig, der RESET-Test hochsignifikant. Was liegt vor?

Eine übersehene Nichtlinearität, keine ungleiche Streuung

Richtig. Die beiden Tests prüfen verschiedene Voraussetzungen. Breusch-Pagan sucht ungleiche Streuung, RESET sucht fehlende Terme im Modell. Dass der eine schweigt, sagt nichts über den anderen.

Ein Widerspruch, einer der beiden Tests muss falsch sein

Sie widersprechen sich nicht, sie beantworten verschiedene Fragen.

Ungleiche Streuung, die der Test übersehen hat

Möglich, aber der RESET-Test wäre dafür der falsche Anzeiger.

Eine Beobachtung hat einen Hebelwert von 0.21 bei einem mittleren Hebel von 0.017, aber ein studentisiertes Residuum von 0.3. Ist sie ein Problem?

Nein, sie liegt weit aussen und trotzdem auf der Linie der übrigen

Richtig. Hoher Hebel allein ist unproblematisch und macht die Schätzung sogar genauer. Gefährlich ist die Kombination aus hohem Hebel und grossem Residuum, und genau die misst die Cook-Distanz.

Ja, jeder Hebel über dem Dreifachen des Mittelwerts ist kritisch

Die Faustregel markiert Punkte, die man ansehen sollte, nicht solche, die man entfernen muss.

Ja, weil sie ausserhalb des Datenbereichs liegt

Ausserhalb des Bereichs liegt sie nicht, sie liegt an dessen Rand; und das ist zulässig.

Die Daten sind nach dem Prädiktor sortiert, und Durbin-Watson liefert 0.53. Was folgt daraus?

Ein Hinweis auf Krümmung, nicht auf Autokorrelation

Richtig. Durbin-Watson misst die Ähnlichkeit aufeinanderfolgender Zeilen. Bei nach x sortierten Daten führt ein Bogen im Residuenverlauf genau zu diesem Bild. Für echte Autokorrelation müssten die Zeilen zeitlich geordnet sein.

Die Beobachtungen sind voneinander abhängig

Das wäre der Schluss bei zeitlich geordneten Daten.

Der Test ist bei sortierten Daten nicht anwendbar

Anwendbar ist er, er misst nur etwas anderes als erwartet.

Verlinkte Ressourcen

• An Introduction to Statistical Learning, Kapitel 3.3
• R-Dokumentation zu plot.lm
• statsmodels: Outlier and influence measures