Um auf einflussreiche Fälle zu prüfen, nutzen wir wieder die Resdiuen. Hier sind zwei Werte relevant: Leverage und Cook’s D(istance). Da Cook’s D(istance) in Teilen auf dem leverage beruht, fokussieren wir direkt Cook’s D.
{{% tabs %}} {{% tab name="Leverage" %}}
Ist die partielle Ableitung des geschätzten \(i\)-ten Wert der abhängigen Variable \(\hat{y_i}\) in Abhängigkeit des gemessenen \(i\)-ten Wert der abhängigen Variable \(y_i\)
Der Grad, in dem der \(i\)-te gemessene Wert den \(i\)-ten gefitteten Wert beeinflusst wird
“Ungewöhnlichkeit” der \(x\)’s
\(h_{ii} = \frac{\partial \hat{y_i}}{\partial y_i}\)
Fall gilt als einflussreich, falls \(h_{ii} > .2\) {{% /tab %}} {{% tab name="Cook’s D" %}}
Ist die Distanz der Veränderung, die ein geschätzter Wert innerhalb der Konfidenz-Ellipse bewegt, die die Region plausibler Werte für den Parameter annimmt
Misst den Einfluss (wie stark ein Fall das Regressionsmodell ändert)
\(D_i = \frac{ \sum_{j=1}^n (\hat{y}_j - \hat{y}_{j(i)})}{p * s^2}\), wobei \(p\) Anzahl der Parameter ist
Klassifiziert als einflussreich, wenn \(D_i > \bar{D} *3\) {{% /tab %}} {{% /tabs %}}
Zuerst können wir einfach Cook’s D plotten. Dazu nutzt du wieder die Funktion plot(), gibst das Regressionsmodell an und setzt im zweiten Argument den Wert 4:
plot(
olsModel2,
4
)
Grenzwerte:
\(D_i < 1\): keine einzige Beobachtung
\(\frac{4}{4890} \approx 0.0008179959\): viele Fälle
\(\bar{D} * 3 \approx 0.0005895857\): viele Fälle
Wir können auch direkt Cook’s D gegen das leverage plotten und Linien für die Grenzwerte einziehen:
plot(
olsModel2,
6
)
abline(
h = 3 * mean(cooks.distance(olsModel2), na.rm = TRUE),
col = "chartreuse4",
lty = 6, # Grenze Cook's D
lwd = 2
)
abline(
v = 2 * (2 + 1) / 2393,
col = "chartreuse4",
lty = 6, # Grenze leverage
lwd = 2
) Um nun zu prüfen, ob diese theoretisch potentiell einflussreichen Fälle wirklich einflussreich sind, berechnen wir das Modell erneut unter Ausschluss dieser Beobachtungen.
Dafür benötigen wir einen Datensatz, der dieselbe Länge hat wie unsere Modellschätzung:
nobs(olsModel2)## [1] 4893ccol <- c(
"stfdem",
"stfeco",
"trstlgl"
)
pss2 <- pss[complete.cases(pss[,ccol]),]
dim(pss2)## [1] 4893 14Nun können wir die Werte des leverage (hat values) und Cook’s D dem Datensatz hinzufügen.
pss2$hats <- hatvalues(olsModel2)
pss2$cooksd <- cooks.distance(olsModel2) Alles, was jetzt folgt, benötigt ein tieferes Verständnis, und insbesonder das Verständnis für Schleifen und if-Bedingungen aus dem Lernblock 2. Aber wie geschrieben, ist dieses Kapitel extra und wird nicht Teil der verpflichtenden Lab-Übung sein.
Um nun die potentiell einflussreichen Fälle auszuschließen, benötigen wir einen Dummy, der uns anzeigt, ob der Fall nach den Kriterien einflussreich ist (\(1\)) oder nicht (\(0\)): In der Schleife prüfen wir zuerst, ob der jeweilige Wert des Falles über dem Grenzwert ist. Wenn ja, bekommt für das leverage die Variable leverage den Wert 1 bzw. für Cook’s D die Variable influence.
for (i in 1:dim(pss2)[1]) {
if (pss2$hats[i] > (2*(2+1)/dim(pss2)[1])) {
pss2$leverage[i] = 1
} else {
pss2$leverage[i] = 0
}
}
for (i in 1:dim(pss2)[1]) {
if (pss2$cooksd[i] > (3*mean(pss2$cooksd))) {
pss2$influence[i] = 1
} else {
pss2$influence[i] = 0
}
}Nun können wir das Modell mit dem verkleinerten Datensatz erneut rechnen: Hierbei beschränken wir den Datensatz auf die Fälle, die jeweils in den neu geschaffenen Variablen leverage und influence den Wert 0 haben, also nicht einflussreich sind.
olsModel2Re <- lm(
stfdem ~ 1 + stfeco + trstlgl,
data = pss2[pss2$leverage==0 & pss2$influence==0,]
)Um herauszufinden, wie viele Fälle wir ausgeschlossen haben, kannst du folgende einfache Subtraktion machen:
nobs(olsModel2Re) - nobs(olsModel2)## [1] -736Es wurden also 736 Fälle ausgeschlossen. Eine ganze Menge!
Schauen wir uns das Ergebnis des Re-definierten Modells mit summary() an:
Re-definiertes Modelll
summary(olsModel2Re)##
## Call:
## lm(formula = stfdem ~ 1 + stfeco + trstlgl, data = pss2[pss2$leverage ==
## 0 & pss2$influence == 0, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.8164 -0.9224 0.1008 1.1405 4.2267
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.34644 0.10316 3.358 0.000791 ***
## stfeco 0.93705 0.01508 62.147 < 2e-16 ***
## trstlgl -0.04306 0.01494 -2.883 0.003959 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.536 on 4154 degrees of freedom
## Multiple R-squared: 0.482, Adjusted R-squared: 0.4818
## F-statistic: 1933 on 2 and 4154 DF, p-value: < 2.2e-16Doch weicht das nun ab vom alten Modell? Dazu nehmen wir einfach die Funktion modelsummary(), die wir schon kennengelernt haben. Denn in einer Tabelle können wir die Werte besser vergleichen. Wichtig ist hierbei zu prüfen, ob sich die Standardfehler bzw. die Koeffizienten ändern.
modelsummary(
list(
olsModel2,
olsModel2Re
),
stars = TRUE
)| Model 1 | Model 2 | |
|---|---|---|
| (Intercept) | 0.677*** | 0.346*** |
| (0.093) | (0.103) | |
| stfeco | 0.874*** | 0.937*** |
| (0.014) | (0.015) | |
| trstlgl | −0.042** | −0.043** |
| (0.013) | (0.015) | |
| Num.Obs. | 4893 | 4157 |
| R2 | 0.460 | 0.482 |
| R2 Adj. | 0.460 | 0.482 |
| AIC | 19265.4 | 15372.7 |
| BIC | 19291.4 | 15398.0 |
| Log.Lik. | −9628.717 | −7682.341 |
| F | 2081.151 | 1932.889 |
| RMSE | 1.73 | 1.54 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
\(\Rightarrow\) Du siehst, dass sich der Koeffizient für den intercept als auch für stfeco ändert. Der intercept wird kleiner, der Effekt von stfdem dagegen größer. Die Standardfehler (Werte in Klammern) ändern sich dagegen nicht.
Wenn wir große Änderungen erkennen, müssen die gefilterten Fälle genauer inspiziert werden. Denn wir können Fälle nicht einfach aufgrund mathematischer Bedingungen exkludieren.
Um herauszufinden, um welche Fälle es sich handelt, kann man wieder eine Schleife programmieren:
levInfluential <- c()
cookInfluential <- c()
allInfluential <- c()
# nur für leverage
for (i in 1:dim(pss2)[1]) {
if (pss2$leverage[i] == 1) {
levInfluential <- c(levInfluential, pss2$idno[i])
}
}
# nur für Cook's D
for (i in 1:dim(pss2)[1]) {
if (pss2$influence[i] == 1) {
cookInfluential <- c(cookInfluential, pss2$idno[i])
}
}
# für beides
for (i in 1:dim(pss2)[1]) {
if (pss2$leverage[i] == 1 & pss2$influence[i] == 1) {
allInfluential <- c(allInfluential, pss2$idno[i])
}
}Als Ergebnis bekommst du drei Vektoren, die jeweils die einflussreichen Fälle beinhalten:
levInfluential## [1] 10008 10018 10023 10034 10041 10081 10083 10094 10116 10125 10147 10180
## [13] 10205 10216 10220 10239 10240 10255 10259 10266 10268 10274 10278 10281
## [25] 10296 10302 10305 10316 10348 10353 10364 10370 10389 10398 10418 10423
## [37] 10432 10446 10466 10476 10491 10494 10517 10520 10530 10536 10547 10548
## [49] 10559 10568 10571 10587 10592 10612 10637 10642 10654 10662 10668 10688
## [61] 10692 10705 10706 10707 10730 10772 10773 10777 10781 10797 10812 10825
## [73] 10831 10835 10841 10846 10865 10866 10888 10898 10931 10934 10975 10989
## [85] 10995 20011 20014 20055 20058 20071 20089 20092 20103 20122 20142 20147
## [97] 20149 20158 20159 20172 20179 20208 20231 20242 20255 20256 20263 20267
## [109] 20279 20318 20329 20339 20352 20371 20373 20381 20390 20398 20400 20402
## [121] 20407 20414 20420 20459 20462 20487 20500 20514 20539 20552 20555 20556
## [133] 20563 20588 20591 20600 20632 20633 20639 20642 20648 20651 20668 20678
## [145] 20684 20705 20708 20743 20744 20745 20763 20767 20776 20777 20780 20797
## [157] 20817 20822 20839 20855 20866 20876 20882 20896 20916 20917 20938 20944
## [169] 20969 20992 30005 30022 30043 30061 30068 30081 30083 30171 30184 30195
## [181] 30201 30222 30232 30244 30251 30252 30287 30295 30299 30310 30314 30332
## [193] 30341 30348 30360 30366 30375 30379 30383 30401 30402 30412 30430 30434
## [205] 30456 30482 30493 30498 30500 30501 30523 30533 30537 30547 30558 30564
## [217] 30601 30644 30664 30669 30683 30698 30715 30718 30731 30743 30754 30764
## [229] 30771 30773 30775 30782 30788 30813 30833 30856 30859 30863 30864 30866
## [241] 30874 30880 30886 30901 30908 30913 30927 40016 40037 40101 40104 40113
## [253] 40124 40137 40160 40163 40171 40173 40202 40207 40213 40221 40230 40233
## [265] 40248 40255 40270 40291 40305 40313 40314 40317 40327 40336 40357 40365
## [277] 40367 40376 40388 40416 40429 40437 40441 40447 40457 40481 40490 40525
## [289] 40528 40530 40541 40552 40571 40599 40607 40618 40619 40632 40652 40671
## [301] 40690 40694 40709 40712 40718 40722 40723 40726 40728 40733 40739 40742
## [313] 40745 40774 40794 40803 40809 40849 40850 40855 40862 40877 40890 40892
## [325] 40921 40930 40972 40978 40982 40994 40995 50000 50009 50033 50046 50053
## [337] 50066 50089 50092 50103 50105 50114 50153 50170 50183 50192 50198 50199
## [349] 50207 50210 50212 50220 50230 50236 50244 50252 50254 50264 50265 50273
## [361] 50303 50309 50320 50324 50328 50331 50352 50382 50387 50394 50413 50414
## [373] 50415 50423 50441 50444 50452 50456 50465 50468 50469 50483 50486 50492
## [385] 50502 50505 50507 50511 50513 50541 50544 50556 50574 50576 50584 50586
## [397] 50590 50596 50617 50623 50631 50639 50641 50643 50682 50686 50691 50694
## [409] 50697 50699 50710 50720 50750 50763 50774 50782 50788 50793 50807 50826
## [421] 50833 50834 50860 50882 50902 50921 50923 50928 50933 50954 50979 50986
## [433] 50988 50998cookInfluential## [1] 10007 10008 10010 10024 10050 10054 10068 10072 10076 10088 10094 10100
## [13] 10116 10119 10125 10135 10138 10146 10147 10150 10152 10216 10248 10255
## [25] 10257 10278 10291 10316 10346 10361 10395 10398 10410 10413 10432 10434
## [37] 10450 10459 10460 10475 10483 10495 10497 10530 10544 10546 10559 10585
## [49] 10613 10617 10620 10624 10633 10638 10655 10662 10677 10679 10680 10682
## [61] 10708 10709 10720 10726 10738 10759 10771 10777 10781 10782 10787 10789
## [73] 10790 10794 10798 10807 10812 10817 10820 10821 10879 10924 10944 10952
## [85] 10989 10998 10999 20054 20058 20063 20075 20084 20103 20118 20147 20170
## [97] 20172 20191 20228 20248 20255 20257 20263 20267 20301 20328 20329 20335
## [109] 20363 20365 20372 20373 20381 20390 20398 20400 20417 20447 20462 20465
## [121] 20472 20487 20494 20512 20513 20514 20530 20535 20579 20580 20600 20604
## [133] 20623 20636 20637 20639 20659 20668 20684 20707 20732 20736 20746 20748
## [145] 20763 20775 20776 20855 20860 20863 20864 20873 20907 20909 20916 20925
## [157] 20932 20950 20974 20975 20984 20998 30014 30052 30059 30072 30097 30118
## [169] 30122 30165 30184 30188 30200 30207 30210 30218 30232 30240 30249 30251
## [181] 30253 30254 30255 30292 30304 30312 30315 30321 30322 30336 30337 30348
## [193] 30370 30376 30383 30386 30412 30420 30428 30430 30445 30474 30479 30482
## [205] 30484 30485 30486 30501 30530 30537 30546 30547 30558 30564 30613 30614
## [217] 30621 30628 30638 30683 30692 30694 30695 30697 30701 30725 30731 30754
## [229] 30763 30768 30788 30811 30812 30813 30820 30831 30856 30867 30877 30898
## [241] 30910 30914 30928 30939 30956 30958 30959 30970 40037 40050 40052 40056
## [253] 40059 40067 40073 40089 40094 40113 40116 40121 40162 40172 40173 40186
## [265] 40192 40193 40205 40212 40213 40215 40218 40220 40244 40261 40278 40284
## [277] 40291 40297 40301 40305 40320 40328 40329 40331 40336 40355 40363 40369
## [289] 40376 40377 40401 40424 40440 40443 40447 40450 40451 40452 40487 40517
## [301] 40521 40525 40559 40568 40590 40607 40618 40638 40658 40672 40690 40726
## [313] 40731 40733 40739 40745 40770 40776 40779 40796 40803 40813 40814 40825
## [325] 40834 40850 40855 40868 40895 40927 40954 40974 40976 40993 50007 50034
## [337] 50046 50072 50085 50092 50096 50128 50171 50191 50207 50210 50224 50246
## [349] 50274 50283 50290 50302 50303 50328 50340 50368 50376 50392 50413 50414
## [361] 50415 50416 50430 50452 50456 50483 50487 50513 50553 50555 50562 50575
## [373] 50587 50591 50596 50609 50621 50623 50639 50640 50648 50649 50652 50669
## [385] 50682 50691 50696 50708 50710 50720 50727 50750 50755 50762 50780 50785
## [397] 50795 50815 50826 50833 50834 50843 50845 50853 50867 50902 50931 50954
## [409] 50982 50986allInfluential## [1] 10008 10094 10116 10125 10147 10216 10255 10278 10316 10398 10432 10530
## [13] 10559 10662 10777 10781 10812 10989 20058 20103 20147 20172 20255 20263
## [25] 20267 20329 20373 20381 20390 20398 20400 20462 20487 20514 20600 20639
## [37] 20668 20684 20763 20776 20855 20916 30184 30232 30251 30348 30383 30412
## [49] 30430 30482 30501 30537 30547 30558 30564 30683 30731 30754 30788 30813
## [61] 30856 40037 40113 40173 40213 40291 40305 40336 40376 40447 40525 40607
## [73] 40618 40690 40726 40733 40739 40745 40803 40850 40855 50046 50092 50207
## [85] 50210 50303 50328 50413 50414 50415 50452 50456 50483 50513 50596 50623
## [97] 50639 50682 50691 50710 50720 50750 50826 50833 50834 50902 50954 50986Anschließend könntest du dir die Wertepaare auf den Variablen in der Regression pro Fall genauer anschauen!
Das war’s zur Regressionsdiagnostik. In den sechs Schritten kannst du also sichergehen, ob deine Regressionsrechnung nicht aufgrund einer Annahmenverletzung verzerrt ist.