Vorab:

Auch hier laeuft das Skript nicht in Gaenze automatisch durch, sondern die einzelnen Schritte (“chunks”) sollten Schritt fuer Schritt nachgearbeitet werden.

Kap. 25.6.9: Multiple lineare Regression

… wieder anhand des Datensatzes “Rollet 1888”. Ziel: Koerperhoehenschaetzung anhand des Oberarmknochens (H1), des Oberschenkelknochens (F1) und des Schienbeins (T1). Nach allen Erfahrungen (gerne selbst versuchen) ist die Hinzunahme der Unterarmknochen (R1) nicht sinnvoll. Zur Erlaeuterung siehe im Buch S. 371 ff.

> setwd("D:/RData")
> #
> Rollet1888 <- readSPSS("D:/RData/5_Rollet1888.sav", rownames=FALSE, 
+   stringsAsFactors=TRUE, tolower=FALSE)
> #
> RegMult <- lm(Khoehe~F1+H1+T1, data=Rollet1888)
> summary(RegMult)

Call:
lm(formula = Khoehe ~ F1 + H1 + T1, data = Rollet1888)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.9825 -2.1371 -0.0433  2.0831  9.5935 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 50.82918    4.97472  10.218  < 2e-16 ***
F1           0.07255    0.03286   2.208  0.02966 *  
H1           0.16066    0.04259   3.772  0.00028 ***
T1           0.08123    0.03761   2.160  0.03327 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.359 on 96 degrees of freedom
Multiple R-squared:  0.8412,    Adjusted R-squared:  0.8362 
F-statistic: 169.5 on 3 and 96 DF,  p-value: < 2.2e-16

Die ueblichen weiteren Pruefungen. Beachte: “Warte auf Bestaetigung des Seitenwechsels…” bezieht sich auf das Plot-Fenster des R-Commanders, wo man links oben - satt des ueblichen “Datei - History - Rsize” ein “Weiter” liest und dieses anklicken muss, um die Bildfolge sehen zu koennen.

> Rollet1888<- within(Rollet1888, {
+   fitted.RegMult <- fitted(RegMult)
+   residuals.RegMult <- residuals(RegMult)
+   rstudent.RegMult <- rstudent(RegMult)
+   hatvalues.RegMult <- hatvalues(RegMult)
+   cooks.distance.RegMult <- cooks.distance(RegMult)
+   obsNumber <- 1:nrow(Rollet1888) 
+ })
> #
> influencePlot(RegMult, id=list(method="noteworthy", n=5))

      StudRes        Hat       CookD
5   2.1676070 0.03921355 0.046162880
15 -0.3794232 0.12160947 0.005027562
18  0.7298767 0.09999314 0.014869025
27 -1.5557797 0.11553206 0.077889288
40  2.9975922 0.01663036 0.035072626
43  2.0397053 0.05978384 0.064027116
51 -2.8208490 0.03612197 0.069512608
59 -2.3676793 0.01222808 0.016555230
63 -1.9147769 0.04868974 0.045645050
73  0.6285525 0.09123702 0.009979043
77  1.0879650 0.10349886 0.034097654
85 -2.3299243 0.02507206 0.033362236
> #
> influenceIndexPlot(RegMult, id=list(method="y", n=5), vars=c("Cook", 
+   "Studentized", "Bonf", "hat"))

> #
> plot(RegMult)

> #
> outlierTest(RegMult)
No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:
   rstudent unadjusted p-value Bonferroni p
40 2.997592          0.0034725      0.34725

Pruefen, ob die Residuen signifikant von einer Normalverteilung abweichen:

> normalityTest(~residuals.RegMult, test="shapiro.test", data=Rollet1888)

    Shapiro-Wilk normality test

data:  residuals.RegMult
W = 0.99363, p-value = 0.9229

Pruefen auf Autokorrelation der Eingangsvariablen (S. 373):

> dwtest(Khoehe ~ F1 + H1 + T1, alternative="greater", data=Rollet1888)

    Durbin-Watson test

data:  Khoehe ~ F1 + H1 + T1
DW = 1.9559, p-value = 0.3755
alternative hypothesis: true autocorrelation is greater than 0

Pruefen auf Homoskedaszitaet resp. Hetereoskedaszitaet, erst die Eingangsvariablen, dann die vorhergesagten Werte (S. 373-374):

> bptest(Khoehe ~ F1 + H1 + T1, studentize=FALSE, data=Rollet1888)

    Breusch-Pagan test

data:  Khoehe ~ F1 + H1 + T1
BP = 0.70375, df = 3, p-value = 0.8723
> #
> bptest(Khoehe ~ F1 + H1 + T1, varformula = ~ fitted.values(RegMult), 
+   studentize=FALSE, data=Rollet1888)

    Breusch-Pagan test

data:  Khoehe ~ F1 + H1 + T1
BP = 0.033353, df = 1, p-value = 0.8551

Pruefen auf Multikollinearitaet (p. 374):

> vif(RegMult)
      F1       H1       T1 
8.298703 8.260075 7.658881 
> round(cov2cor(vcov(RegMult)), 3) # Correlations of parameter estimates
            (Intercept)     F1     H1     T1
(Intercept)       1.000 -0.268  0.028 -0.116
F1               -0.268  1.000 -0.518 -0.460
H1                0.028 -0.518  1.000 -0.456
T1               -0.116 -0.460 -0.456  1.000

Ermitteln der 95% Bootstrap-Intervalle (S. 375) - braucht etwas Rechen-Zeit:

> .bs.samples <- Boot(RegMult, R=2000, method="case")
> plotBoot(.bs.samples)

> confint(.bs.samples, level=0.95, type="bca")
Bootstrap bca confidence intervals

                   2.5 %     97.5 %
(Intercept) 41.464282661 59.0175064
F1           0.008494248  0.1369580
H1           0.080340207  0.2499102
T1           0.005603232  0.1507187
> remove(.bs.samples)

finis