Lineare Regression
8. Oktober 2002
1 Lineare erkl¨ arende Variablen
Wir nehmen zuerst an, dass der Erwartungswert des Preises eines Hauses in Gainesville linear von der Anzahl der Badezimmer, bath, und vom Status newabh¨angt. Die Anzahl der H¨auser in den einzelnen Kategorien erh¨alt man mittels
attach(houses); table(new, bath) bath
new 1 2 3 0 10 51 4 1 0 26 2
Interessanterweise ist kein neues Haus mit nur einem Badezimmer in den Daten enthalten. Die mittleren Preise in den Kategorien liefert
tapply(price, list(factor(new), factor(bath)), mean)
1 2 3
0 40.32 91.94 177.25 1 NA 113.69 249.70
So ist beispielsweise der mittlere Preis eines neuen Hauses mit zwei Badezimmern gleich 113.690$.
Passen wir jetzt ein lineares Regressionsmodell in den beiden erkl¨arenden Variablenbathundnewdiesen Daten an, d.h. betrachten wir das Modell
E(pricei) =β1+βbbathi+βnnewi, i= 1, . . . , n= 93, (1) so liefert die Sch¨atzung
var.mod <- lm(price ~ bath + new); summary(var.mod) Coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept) -47.113 14.869 -3.169 0.00209 **
bath 71.480 7.553 9.463 3.77e-15 ***
new 22.454 6.793 3.305 0.00136 **
---
Residual standard error: 29.55 on 90 degrees of freedom Multiple R-Squared: 0.5625, Adjusted R-squared: 0.5528 F-statistic: 57.85 on 2 and 90 DF, p-value: 1.11e-016
Der Hauspreis steigt f¨ur jedes weitere Badezimmer um jeweils 71.480$. Dieser Preiszuwachs ist signifikant (p < 0.001). Generell kosten neue H¨auser um 22.454$ mehr als gebrauchte (p = 0.001). Der Intercept (Modellkonstante) steht f¨ur den theoretischen Preis eines gebrauchten Hauses ohne Bad. Die gesch¨atzte Standardabweichung der Preise ist 29.550$. Das Modell erkl¨art nur 56% der Preisvariabilit¨at. Die globale Hypothese, dass keine der beiden Variablen etwas vom Hauspreis erkl¨art, kann mit dem overall F-Test aufH0:βb=βn= 0 verworfen werden (p <0.001).
Um Sch¨atzungen (Vorhersagen, predictions) f¨ur die zu erwartenden Hauspreise zu erhalten, definiert man vorerst neue (erkl¨arende) Daten, f¨ur die dann das obige Modell (1) gerechnet wird.
p.data <- data.frame(bath = c(1, 2, 3, 1, 2, 3), new = c(0, 0, 0, 1, 1, 1)); p.data bath new
1 1 0
2 2 0
3 3 0
4 1 1
5 2 1
6 3 1
var.pred <- predict(var.mod, p.data); matrix(var.pred, byrow=T, 2, 3) [,1] [,2] [,3]
[1,] 24.367 95.847 167.328 [2,] 46.821 118.301 189.782
So kostet ein neues Haus mit zwei Badezimmer unter dem Modell (1) gesch¨atzte 118.301$. Es ist auch m¨oglich, eine mittlere Preissch¨atzung f¨ur ein neues Haus mit einem Bad anzugeben (46.821$). F¨ur diese Kombination waren keine Daten verf¨ugbar. Diese Tabelle zeigt auch den linearen Einfluss der Anzahl der Badezimmer auf den mittleren Hauspreis. Egal, ob ein Haus gebraucht oder neu ist, steigt der mittlere Preis f¨ur jedes weitere Bad um 71.480$. Die Tabelle mit den mittleren Hauspreisen (ohne Modellannah- men) zeigte aber eine Verteuerung, die von einem Bad auf zwei B¨ader deutlich geringer ausfiel (etwa 50.000$ f¨ur ein gebrauchtes Haus) als der Preiszuwachs von zwei aud drei Badezimmer (etwa 86.000$ bei einem gebrauchten Haus, und sogar 136.000$ bei neuen H¨ausern). Dies ist ein Hinweis darauf, dassbath nicht als lineare erkl¨arende Gr¨oße mit nur einem Parameter in das Modell eingehen sollte, sondern als Faktor, dessen verschiedene Stufen jeweils durch eigene Parameter charakterisiert sind.
2 Kategoriale erkl¨ arende Faktoren
Es wird nun angenommen, dass bathals dreistufiger, kategorialer Faktor und newals zweistufiger Indi- kator in das Regressionsmodell eingehen. Wir betrachten daher das Modell
E(pricei) =β1+βbathi=2+βbathi=3+βnewi=1. (2) Standardm¨assig wird in R die Dummy- oderTreatment-Kodierung f¨ur Faktoren verwendet. Die Erwar- tungswerte werden hierbei folgenderweise aufgebaut:
bath
new 1 2 3
0 β1 β1+βbath=2 β1+βbath=3
1 β1+βnew=1 β1+βnew=1+βbath=2 β1+βnew=1+βbath=3
Z.B. beschreibt der Parameter βbath=2 die Differenz E(price|bath=2, new)−E(price|bath=1, new) f¨ur eine beliebige Stufe von new. Dieser Parameter wird im R-Output factor(bath)2 genannt. Als Mo- dellsch¨atzung folgt
fac.mod <- lm(price ~ factor(bath) + factor(new)); summary(fac.mod) Coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept) 40.320 8.888 4.536 1.78e-05 ***
factor(bath)2 50.391 9.701 5.194 1.29e-06 ***
factor(bath)3 152.615 14.677 10.398 < 2e-16 ***
factor(new)1 25.394 6.525 3.892 0.000192 ***
---
Residual standard error: 28.11 on 89 degrees of freedom Multiple R-Squared: 0.6085, Adjusted R-squared: 0.5953 F-statistic: 46.11 on 3 and 89 DF, p-value: 0
Die Preisdifferenz zwischen H¨auser mit einem und mit zwei Badezimmer ist damit 50.391$. Wir m¨ussen weitere 152.615$ - 50.391$ = 102.224$ zahlen, falls wir an einem Haus mit drei Badezimmer interessiert sind. Unter dem Modell (1) waren beide Differenzen gleich 71.480$. Neue H¨auser kosten jetzt generell um 25.394$ mehr als gebrauchte. Dieser Unterschied h¨angt dabei nicht von der Anzahl der Badezimmer ab.
Die gesch¨atzten Erwartungswerte sind
fac.pred <- predict(fac.mod, p.data); matrix(fac.pred, byrow=T, 2, 3) [,1] [,2] [,3]
[1,] 40.320 90.711 192.935 [2,] 65.714 116.105 218.330
Diese passen viel besser zu den beobachteten Mittelwerten als die Prognosen unter dem Modell (1). Da wir jetzt eine nichtlineare Struktur inbatherlauben, ben¨otigt das Modell (2) um einen Parameter mehr und hat eine etwas bessere Prognosekraft von R2 = 0.61. Auch die gesch¨atzte Standardabweichung von priceist jetzt mit 28.110$ etwas geringer als zuvor.
2.1 Modellieren ohne Intercept
Betrachtet man das Modell (2) ohne Intercept, so ergibt dies
E(pricei) =βbathi=1+βbathi=2+βbathi=3+βnewi=1. (3) Definierte in (2) der Interceptβ1 den mittleren Preis eines gebrauchten Hauses mit einem Bad, so ¨uber- nimmt dies jetzt der Parameterβbathi=1. Die beiden anderen Bad-Parameterβbathi=2 undβbathi=3 sind hier keine Differenzen sondern modellieren den Preis direkt, d.h. man hat jetzt
bath
new 1 2 3
0 βbath=1 βbath=2 βbath=3
1 βnew=1+βbath=1 βnew=1+βbath=2 βnew=1+βbath=3
fac.noInt.mod <- lm(price ~ -1 + factor(bath) + factor(new)); summary(fac.noInt.mod) Coefficients:
Estimate Std. Error t value Pr(>|t|) factor(bath)1 40.320 8.888 4.536 1.78e-05 ***
factor(bath)2 90.711 3.888 23.332 < 2e-16 ***
factor(bath)3 192.935 11.679 16.520 < 2e-16 ***
factor(new)1 25.394 6.525 3.892 0.000192 ***
---
Residual standard error: 28.11 on 89 degrees of freedom Multiple R-Squared: 0.9361, Adjusted R-squared: 0.9333 F-statistic: 326.1 on 4 and 89 DF, p-value: 0
Die Prognosen sind ident denen von (2).R2 und die F-Statistik werden aber unter (3) anders gerechnet.
Ein interceptfreies Modell ist ¨aquivalent mit der Annahme, dass der globale Erwartungswert µder Re- sponsevariableny gleich Null ist. Unter dieser, hier sicherlich falschen, Annahme resultiert die ANOVA- Zerlegung
Xn i=1
(yi−0)2= Xn i=1
(yi−µˆi+ ˆµi−0)2= Xn i=1
(yi−µˆi)2+ Xn i=1
ˆ µ2i. Der gemischte TermP
i(yi−µˆi)ˆµi= [(I−H)y]0Hy=y0Hy−y0Hy= 0 verschwindet wiederum. W¨ahrend SSE(ˆµ) =P
i(yi−µˆi)2der ¨ublichen Definition entspricht, sind also die Quadratsummen SSR(ˆµ) =P
iµˆ2i und SST anders definiert. R2 beschreibt jetzt den Anteil der vom Modell (ohne Intercept) erkl¨arten Variablit¨at, wobei hier die Responsevariabilit¨at unter der Annahme µ = 0 gesch¨atzt ist. Aus diesem Grund ist auchR2 nicht mehr mit dem unter Modell (2) vergleichbar.
2.2 Interaktionen
Unter dem Modell (2) wurde angenommen, dass sich die Zellenmittelwerte additiv bilden. Dabei wurde f¨ur jede erkl¨arende, kategoriale Gr¨oße im Prediktor eine Referenzkategorie definiert und durch den Intercept parametrisiert. Im Modell (2) war dies die Zelle mit bath = 1 und new= 0. Unabh¨angig von der new Kategorie wurden die beiden Abweichungen vonbathzur Refernzstufe (ein Badezimmer) parametrisiert.
Dies ben¨otigte daher zwei Parameter; einen weniger als der Faktor Stufen hatte. Auch f¨ur den anderen Faktornewist man gleich vorgegangen. Die eine Abweichung von der Referenzstufe (gebrauchtes Haus) wurde durch einen Parameter beschrieben. Derartige Modell werden Modelle ohne Interaktionen genannt.
Sie beinhalten nur Haupteffekte. Manchmal ist es aber notwendig, zus¨atzlich zu diesen Haupteffekten auch Interaktionseffekte im Modell zu ber¨ucksichtigen. Dies erlaubt dann zell-spezifische Abweichungen g¨anzlich frei zu modellieren. F¨ur unserem Fall folgt dann
E(pricei) =β1+βbathi=2+βbathi=3+βnewi=1+βbathi=2,newi=1+βbathi=3,newi=1. (4) Dieses Modell kann jetzt den Erwartungswert in jeder Zelle frei beschreiben und es gilt
bath
new 1 2 3
0 β1 β1+βbath=2 β1+βbath=3
1 β1+βnew=1 β1+βnew=1+βbath=2+βbath=2,new=1 β1+βnew=1+βbath=3+βbath=3,new=1
F¨ur die sechs Zellen stehen jetzt sechs Parameter zur Verf¨ugung. Solche Modell nennt man auch volle oder maximale Modelle. Angepasst werden sie durch
fac.IA.mod <- lm(price ~ factor(bath) * factor(new)); summary(fac.IA.mod) Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 40.320 8.734 4.616 1.32e-05 ***
factor(bath)2 51.621 9.552 5.404 5.50e-07 ***
factor(bath)3 136.930 16.340 8.380 7.73e-13 ***
factor(new)1 72.450 23.920 3.029 0.00322 **
factor(bath)2:factor(new)1 -50.699 24.828 -2.042 0.04415 * ---
Residual standard error: 27.62 on 88 degrees of freedom Multiple R-Squared: 0.6262, Adjusted R-squared: 0.6092 F-statistic: 36.86 on 4 and 88 DF, p-value: 0
Da in einer Zelle keine Beobachtungen vorhanden sind, kann einer der beiden Interaktions-Parameter (hier der Parameter factor(bath)3:factor(new)1) auch nicht gesch¨atzt werden. Bemerke, dass trotzdem jeder beobachtete Zellenmittelwert exakt durch dieses Modell reproduziert werden kann.
fac.IA.pred <- predict(fac.IA.mod, p.data); matrix(fac.IA.pred, byrow=T, 2, 3) [,1] [,2] [,3]
[1,] 40.32 91.94 177.25 [2,] 112.77 113.69 249.70
Weiters liefert das Modell (4) auch eine Prognose f¨ur den Hauspreis in der nichtbeobachteten Zelle (112.770$). Da die Interaktion hier nur aus einem Parameter besteht, ist es m¨oglich, die Notwendigkeit dieser Interaktion ¨uber den Wald-Test zu best¨atigen (p= 0.044). Der Wert vonR2 steigt von 0.61 unter dem Modell (2) auf 0.63. Die gesch¨atzte Standardabweichung reduziert sich dabei geringf¨ugig von 28.110$
auf 27.620$.
3 Erkl¨ arende lineare Variablen und Faktoren
Weitaus realistischer wird das Modell, wenn wir zus¨atzlich nochareaals erkl¨arende Variable einsetzen.
Zuerst wird in jeder der sechs Zellen dieselbearea-Steigung angenommen. Die Modelle f¨ur die einzelnen Zellen unterscheiden sich dann nur in den Intercepts und wir erhalten die sechs parallele Geraden
E(pricei) =β1+βbathi=2+βbathi=3+βnewi=1+βarea∗areai. (5) Bemerke, dass auch die sechs Intercepts hier nicht g¨anzlich frei w¨ahlbar sind sondern mittels nur vier Parameter (Haupteffekte) aufgebaut sind. Als Sch¨atzung resultiert
vf.mod <- lm(price ~ area + factor(bath) + factor(new)); summary(vf.mod) Coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept) -17.133 6.360 -2.694 0.00846 **
area 59.847 4.202 14.243 < 2e-16 ***
factor(bath)2 10.489 6.054 1.733 0.08668 . factor(bath)3 49.434 10.881 4.543 1.75e-05 ***
factor(new)1 20.144 3.628 5.552 2.96e-07 ***
---
Residual standard error: 15.55 on 88 degrees of freedom Multiple R-Squared: 0.8816, Adjusted R-squared: 0.8762 F-statistic: 163.7 on 4 and 88 DF, p-value: 0
plot(area, fitted(vf.mod), xlim=c(0,4), ylim=c(0,310))
abline(-17.133, 59.847); abline(-17.133+10.489, 59.847) abline(-17.133+49.434, 59.847); abline(-17.133+20.144, 59.847) abline(-17.133+10.489+20.144, 59.847); abline(-17.133+49.434+20.144, 59.847)
0 1 2 3 4
050100150200250300
area
fitted(vf.mod)
Der eine Steigungsparameter ist signifikant von Null verschieden (p <0.001). Die Interceptunterschiede sind auch signifikant (p-Werte alle kleiner als 0.09). Die Modellanpassung steigt sogar auf R2 = 0.88.
Es stellt sich aber die Frage, ob nicht vielleicht status-spezifische Anstiege oder f¨ur jede Stufe vonbath unterschiedliche Steigungsparameter noch bessere Ergebnisse liefern w¨urden.
Unterschiedliche Steigungsparameter f¨ur gebrauchte und neue H¨auser erlaubt das Modell E(pricei) =β1+βbathi=2+βbathi=3+βnewi=1+¡
βarea+βarea, newi=1
¢∗areai. (6) Dieses Modell erh¨alt man durch Hinzunahme der Interaktion zwischenareaundfactor(new). Diese Art der Parametrisierung liefert die beiden Steigungenβarea f¨ur gebrauchte H¨auser und βarea+βarea, new=1
bei neuen. Die Gr¨oßeβarea, new=1gibt daher gerade den Unterschied in den beiden Steigungen an.
vf.mod.nslope <- lm(price ~ factor(bath) + area*factor(new)); summary(vf.mod.nslope) Coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept) -12.509 6.196 -2.019 0.04659 * factor(bath)2 14.466 5.870 2.464 0.01570 * factor(bath)3 47.967 10.333 4.642 1.21e-05 ***
area 55.031 4.249 12.953 < 2e-16 ***
factor(new)1 -25.290 14.272 -1.772 0.07991 . area:factor(new)1 25.610 7.808 3.280 0.00149 **
---
Residual standard error: 14.75 on 87 degrees of freedom Multiple R-Squared: 0.8946, Adjusted R-squared: 0.8885 F-statistic: 147.7 on 5 and 87 DF, p-value: 0
Nat¨urlich kann auf die Notwendigkeit von zwei Steigungen auch getestet werden. Im obigen Output ergibt daf¨ur der Wald-Test einen p-Wert von 0.0015 f¨ur die Hypothese H0 :βarea, new=1 = 0. ¨Aquivalent dazu liefert die ANOVA-Zerlegung dasselbe Ergebnis
anova(vf.mod.nslope) Analysis of Variance Table Response: price
Df Sum Sq Mean Sq F value Pr(>F) factor(bath) 2 97328 48664 223.623 < 2.2e-16 ***
area 1 53553 53553 246.090 < 2.2e-16 ***
factor(new) 1 7451 7451 34.239 8.405e-08 ***
area:factor(new) 1 2341 2341 10.759 0.001493 **
Residuals 87 18933 218
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
H0 kann also verworfen werden. Vergr¨oßert man area um eine Einheit, so steigt der Preis f¨ur ein ge- brauchtes Haus um 55.031$ w¨ahrend er bei einem neuen Haus um sogar 55.031$ + 25.610$ = 80.641$
ansteigt. Diese beiden Zuw¨achse sind dabei von der Anzahl der Badezimmer unabh¨angig. ¨Uberraschen- derweise scheinen unter dem Modell (6) generell neue H¨auser 25.290$ weniger als gebrauchte zu kosten.
Da aber f¨ur zwei Badezimmer die beiden Regressionsgeraden
E(price|bath=2, new=0) =−12.509 + 14.466 + 55.031∗area = 1.957 + 55.031∗area
E(price|bath=2, new=1) =−12.509−25.290 + 14.466 + (55.031 + 25.610)∗area =−23.333 + 80.641∗area resultieren, sind neue H¨auser dann teurer als gebrauchte, falls 25.29<25.61∗area, d.h. falls 0.99<area gilt (was in den Daten immer der Fall ist).
plot(area, fitted(vf.mod.nslope), xlim=c(0,4), ylim=c(0,310))
abline(-12.509, 55.031); abline(-12.509+14.466, 55.031)
abline(-12.509+47.967, 55.031); abline(-12.509-25.290, 55.031+25.610) abline(-12.509+14.466-25.290, 55.031+25.610
abline(-12.509+47.967-25.290, 55.031+25.610)
0 1 2 3 4
050100150200250300
area
fitted(vf.mod.nslope)
Alternativ dazu kann man untersuchen, ob H¨auser mit unterschiedlich vielen Badezimmern unterschied- liche Abh¨angigkeiten vonarea aufweisen, d.h. ob das Modell
E(pricei) =β1+βbathi=2+βbathi=3+βnewi=1+¡
βarea+βarea, bathi=2+βarea, bathi=3
¢∗areai (7) sinnvoll ist. Das Sch¨atzresultat lautet
vf.mod.bslope <- lm(price ~ factor(new) + area*factor(bath)); summary(vf.mod.bslope) Coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept) -1.674 17.186 -0.097 0.9227 factor(new)1 19.801 3.375 5.866 8.1e-08 ***
area 43.743 17.257 2.535 0.0131 *
factor(bath)2 6.190 18.777 0.330 0.7424 factor(bath)3 -53.195 30.519 -1.743 0.0849 . area:factor(bath)2 9.435 17.825 0.529 0.5979 area:factor(bath)3 48.272 19.488 2.477 0.0152 * ---
Residual standard error: 14.46 on 86 degrees of freedom Multiple R-Squared: 0.8999, Adjusted R-squared: 0.8929 F-statistic: 128.9 on 6 and 86 DF, p-value: 0
DerR2 Wert ist bereits 90%, die gesch¨atzte Standardabweichung der Hauspreise unter dem Modell (7) gar nur noch 14.460$. Die Steigung bei drei Badezimmern scheint sich von der bei einem Bad signifikant zu unterscheiden (p= 0.015). Bez¨uglich der Steigung bei zwei Badezimmern und der bei einem Bad kann kein signifikanter Unterschied festgestellt werden. Da der Faktorbathdrei Stufen hat, kann man nur mit der ANOVA-Zerlegung auf die Notwendigkeit von drei unterschiedlichen Steigungen testen.
anova(vf.mod.bslope) Analysis of Variance Table Response: price
Df Sum Sq Mean Sq F value Pr(>F) factor(new) 1 22833 22833 109.247 < 2.2e-16 ***
area 1 129539 129539 619.783 < 2.2e-16 ***
factor(bath) 2 5960 2980 14.257 4.492e-06 ***
area:factor(bath) 2 3299 1650 7.893 0.0007126 ***
Residuals 86 17975 209
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Der p-Wert ist mit 0.0007 derart klein, sodass die beiden Parameter βarea, bath=2 und βarea, bath=3 von Null verschieden zu sein scheinen. Dies ist haupts¨achlich in der großen Steigung beibath = 3begr¨undet.
Wir haben also die sechs Modelle
E(price|bath=1, new=0) =−1.67 + 43.74∗area, E(price|bath=1, new=1) = 18.13 + 43.74∗area, E(price|bath=2, new=0) = 4.52 + 53.18∗area, E(price|bath=2, new=1) = 24.32 + 53.18∗area, E(price|bath=3, new=0) =−54.87 + 92.02∗area, E(price|bath=3, new=1) =−35.07 + 92.02∗area plot (area, fitted(vf.mod.bslope), xlim=c(0,4), ylim=c(0,310))
plot (area, fitted(vf.mod.bslope), xlim=c(0,4), ylim=c(0,310))
abline( -1.674, 43.743)
abline( -1.674+19.801, 43.743)
abline( -1.674+ 6.190, 43.743+ 9.435) abline( -1.674+19.801+ 6.190, 43.743+ 9.435) abline( -1.674-53.195, 43.743+48.272) abline( -1.674+19.801-53.195, 43.743+48.272)
0 1 2 3 4
050100150200250300
area
fitted(vf.mod.bslope)
Somit scheinen einerseits die beiden status-spezifischen Steigungen aus dem Modell (6), sowie auch die drei Steigungen je Badezimmeranzahl im Modell (7) notwendig zu sein. Nat¨urlich kann man auch beides simultan ber¨ucksichtigen. Das dazu passende Modell ist dann
E(pricei) = β1+βbathi=2+βbathi=3+βnewi=1
+¡
βarea+βarea, bathi=2+βarea, bathi=3+βarea, newi=1
¢∗areai (8) und f¨uhrt zu den Sch¨atzungen
vf.mod.nb.slope <- lm(price ~ (factor(new) + factor(bath)) * area) summary(vf.mod.nb.slope)
Coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept) -1.674 16.419 -0.102 0.91906 factor(new)1 -19.887 13.463 -1.477 0.14332 factor(bath)2 13.117 18.084 0.725 0.47023 factor(bath)3 -43.317 29.338 -1.476 0.14351
area 43.743 16.487 2.653 0.00952 **
factor(new)1:area 22.394 7.375 3.036 0.00318 **
factor(bath)2:area 5.574 17.077 0.326 0.74491 factor(bath)3:area 41.433 18.754 2.209 0.02985 * ---
Residual standard error: 13.81 on 85 degrees of freedom Multiple R-Squared: 0.9097, Adjusted R-squared: 0.9023 F-statistic: 122.4 on 7 and 85 DF, p-value: 0
anova(vf.mod.nb.slope) Analysis of Variance Table Response: price
Df Sum Sq Mean Sq F value Pr(>F) factor(new) 1 22833 22833 119.6890 < 2.2e-16 ***
factor(bath) 2 86460 43230 226.6033 < 2.2e-16 ***
area 1 49039 49039 257.0541 < 2.2e-16 ***
factor(new):area 1 2341 2341 12.2732 0.0007354 ***
factor(bath):area 2 2717 1358 7.1207 0.0013831 **
Residuals 85 16216 191
DasR2 steigt weiter auf 91% und die gesch¨atzte Standardabweichung von Preis sinkt auf 13.810$. Die factor(bath):areaZeile der obigen ANOVA-Zerlegung testet auf die Notwendigkeit von Badezimmer- spezifischen Steigungen (p = 0.001), gegeben das Modell enth¨alt bereits beide Haupteffekte und zwei status-spezifische Steigungen. Sind aber Badezimmer-spezifische Steigungen notwendig, falls das Modell bereits status-spezifische enth¨alt? Dazu braucht man nur die Reihenfolge der beiden Effekte in der Pre- diktordefinition (lm-Aufruf) ¨andern und man erh¨alt
anova(lm(price ~ (factor(bath) + factor(new)) * area)) Analysis of Variance Table
Response: price
Df Sum Sq Mean Sq F value Pr(>F) factor(bath) 2 97328 48664 255.0882 < 2.2e-16 ***
factor(new) 1 11965 11965 62.7193 8.216e-12 ***
area 1 49039 49039 257.0541 < 2.2e-16 ***
factor(bath):area 2 3299 1650 8.6474 0.0003815 ***
factor(new):area 1 1759 1759 9.2198 0.0031780 **
Residuals 85 16216 191
also wiederum Signifikanz (p= 0.003).
Noch immer sind die sechs Intercepts sowie die sechs Steigungen im Modell (8) nicht g¨anzlich frei w¨ahlbar.
Dies erm¨oglicht das maximale Modell
E(pricei) = β1+βbathi=2+βbathi=3+βnewi=1+βbathi=2,newi=1+βbathi=3,newi=1
+ ¡
βarea+βarea, bathi=2+βarea, bathi=3+βarea, newi=1 (9) +βarea, bathi=2,newi=1+βarea, bathi=3,newi=1
¢∗areai
Modell (9) erlaubt f¨ur die sechs Intercepts in den Zellen genau sechs Parameter wie auch weitere sechs Parameter f¨ur die sechs Steigungen. Angepasst wird es einfach mittels
vf.mod.max <- lm(price ~ factor(new) * factor(bath) * area) summary(vf.mod.max)
Coefficients: (2 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.674 16.082 -0.104 0.91737
factor(new)1 -203.102 81.432 -2.494 0.01461 *
factor(bath)2 10.799 17.862 0.605 0.54711
factor(bath)3 -25.578 29.695 -0.861 0.39153
area 43.743 16.148 2.709 0.00820 **
factor(new)1:factor(bath)2 196.430 83.527 2.352 0.02106 * factor(new)1:area 84.983 27.786 3.058 0.00299 **
factor(bath)2:area 6.966 16.795 0.415 0.67938 factor(bath)3:area 34.836 18.602 1.873 0.06463 . factor(new)1:factor(bath)2:area -70.332 29.823 -2.358 0.02071 * ---
Residual standard error: 13.53 on 83 degrees of freedom Multiple R-Squared: 0.9154, Adjusted R-squared: 0.9063 F-statistic: 99.83 on 9 and 83 DF, p-value: 0
Es resultiert das mit den vorliegenden Daten maximal m¨ogliche R2= 0.91, sowie 13.530$ als gesch¨atzte Preisstreuung. Im Vergleich zum Modell (8) von zuvor haben wir hier zwei zus¨atzliche Intercept-Parameter und zwei weitere Steigungs-Parameter, von denen jeweils einer nicht gesch¨atzt werden kann.
anova(vf.mod.max)
Analysis of Variance Table Response: price
Df Sum Sq Mean Sq F value Pr(>F) factor(new) 1 22833 22833 124.7710 < 2.2e-16 ***
factor(bath) 2 86460 43230 236.2248 < 2.2e-16 ***
area 1 49039 49039 267.9685 < 2.2e-16 ***
factor(new):factor(bath) 1 1518 1518 8.2964 0.005053 **
factor(new):area 1 875 875 4.7804 0.031601 *
factor(bath):area 2 2674 1337 7.3055 0.001195 **
factor(new):factor(bath):area 1 1018 1018 5.5616 0.020710 *
Residuals 83 15189 183
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Interpretation:
• Das Modell erlaubt 12 Parameter. Jedoch sind zwei davon wegen linearer Abh¨angigkeiten in der Designmatrix (X singul¨ar) nicht sch¨atzbar.
• Gegeben das Modell enth¨alt bereits alle drei Haupteffekte (new, bath, area) ben¨otigen wir f¨ur die Intercepts auch den einen sch¨atzbaren Interaktionsparameter zwischennewundbath(p= 0.0050).
• Gegeben das Modell erlaubt bereits sechs frei w¨ahlbare Intercepts sind auch die beiden status- spezifischen Steigungen notwendig (p= 0.0316).
• Zus¨atzlich sind auch weitere zwei Steigungsparameter f¨ur das unterschiedliche Verhalten in den Badezimmer-Kategorien signifikant (p= 0.0012).
• Der eine sch¨atzbare Steigungs-Interaktionsparameter ist dar¨uberhinaus auch notwendig (p= 0.0271).
• Somit hat das finale Modell jeweils f¨unf Intercept- und f¨unf Steigungs-Parameter.
plot (area, fitted(vf.mod.max), xlim=c(0,4), ylim=c(0,310)) abline(-1.674, 43.743) # new=0 & bath=1 abline(-1.674+10.799, 43.743+ 6.966) # new=0 & bath=2 abline(-1.674-25.578, 43.743+34.836) # new=0 & bath=3
abline(-1.674-203.102, 43.743+84.983) # new=1 & bath=1 abline(-1.674-203.102+10.799+196.430, 43.743+84.983+6.966-70.332) # new=1 & bath=2 abline(-1.674-203.102-25.578, 43.743+84.983+34.836) # new=1 & bath=3
0 1 2 3 4
050100150200250300
area
fitted(vf.mod.max)
Die Variable bedhaben wir noch nicht im Modell verwendet. Solltefactor(bed) als Haupteffekt sowie in Interaktion mit den ¨ubrigen Effekten im Modell vertreten sein. Dies kann sehr einfach durch den ANOVA-Vergleich der entsprechenden Modelle getestet werden.
anova(vf.mod.max, lm(price ~ factor(new) * factor(bath) * factor (bed) * area)) Analysis of Variance Table
Model 1: price ~ factor(new) + factor(bath) + area + factor(new):factor(bath) + factor(new):area + factor(bath):area + factor(new):factor(bath):area
Model 2: price ~ factor(new) + factor(bath) + factor(bed) + area + factor(new):factor(bath) + factor(new):factor(bed) + factor(new):area + factor(bath):factor(bed) +
factor(bath):area + factor(bed):area + factor(new):factor(bath):factor(bed) + factor(new):factor(bath):area + factor(new):factor(bed):area +
factor(bath):factor(bed):area + factor(new):factor(bath):factor(bed):area
Res.Df RSS Df Sum of Sq F Pr(>F) 1 83 15189
2 71 12930 12 2259 1.0337 0.4286
Es scheint nicht n¨otig zu sein, unser gefundenes Modell (9) dermassen zu erweitern (p= 0.43). Alternativ k¨onntebedals lineare Z¨ahlvariable vielleicht ins Modell eingehen.
anova(vf.mod.max, lm(price ~ factor(new) * factor(bath) * bed * area)) Analysis of Variance Table
Model 1: price ~ factor(new) + factor(bath) + area + factor(new):factor(bath) + factor(new):area + factor(bath):area + factor(new):factor(bath):area
Model 2: price ~ factor(new) + factor(bath) + bed + area + factor(new):factor(bath) + factor(new):bed + factor(new):area + factor(bath):bed + factor(bath):area + bed:area + factor(new):factor(bath):bed + factor(new):factor(bath):area +
factor(new):bed:area + factor(bath):bed:area + factor(new):factor(bath):bed:area Res.Df RSS Df Sum of Sq F Pr(>F)
1 83 15189
2 76 13062 7 2127 1.77 0.11
Auch dies ist nur aufp= 0.11 signifikant. Dar¨uberhinaus zeigt auch der Residuenplot keine auffallenden Muster in Abh¨angigkeit vonbed. Nur acht dern = 93 standardisierten Residuen liegen ausserhalb des zentralen 95% Bereichs der Standard-Normalverteilung.
r <- rstandard(vf.mod.max); large.r <- as.numeric(rownames(houses)[abs(r) > 2]) plot(bed, r); abline(0, 0)
abline( qnorm(0.975), 0, lty=2); abline( qnorm(0.025), 0, lty=2) text(bed[large.r], r[large.r], large.r, pos=4, cex=0.8)
sum(abs(rstandard(vf.mod.max)) > qnorm(0.975)) [1] 8
1 2 3 4 5
−3−2−1012
bed
r
4
7
3041 65
86 93
Wir wollen noch einige weitere diagnostische Verfahren und Masse auf das Modell (9) anwenden, um dadurch Hinweise auf Besonderheiten in der Sch¨atzung zu bekommen. Zuerst betrachten wir die Diago- nalelemente der HatmatrixH und pr¨ufen, f¨ur welche Datenhii>2p/ngilt (high-leverage Punkte).
h <- lm.influence(vf.mod.max)$hat; is.inf <- as.numeric(rownames(houses)[h > 20/93]) cbind(houses, h)[is.inf, ]
price area bed bath new h 4 137.0 2.40 3 3 0 0.2691 5 309.4 3.30 4 3 1 1.0000 6 17.5 0.40 1 1 0 0.5469 7 19.6 1.28 3 1 0 0.2459 22 72.5 1.28 3 1 0 0.2459 84 142.0 2.12 3 3 0 0.3585 89 149.0 3.05 4 2 0 0.2534 90 150.0 2.04 3 3 0 0.3974 92 190.0 2.57 4 3 1 1.0000 93 280.0 3.85 4 3 0 0.9750
Erstaunlicherweise sind dies alle sechs H¨auser mit jeweils drei Badezimmern. Die drei Beobachtungen 5, 92 und 93 scheinen dabei besonders einflussreich zu sein. Ein sehr kleines, gebrauchtes Haus (Beobachtung 6) mit nur einem Bad und einem Schlafzimmer ist auch darunter.
index <- 1:n; plot(index, h, ylab="h(ii)"); abline(20/93, 0, lty=2) text(index[is.inf], h[is.inf], is.inf, pos=4, cex=0.8)
0 20 40 60 80
0.00.20.40.60.81.0
index
h(ii)
4 5
6
7 22
84
89 90
92 93
Der sechste sch¨atzbare Parameter im Modell steht f¨ur den Interaktionsterm im Intercept w¨ahrend der zehnte Parameter die Interaktion in der Steigung beschreibt. Speziell f¨ur diese beide Parameter werden
die ¨Anderungen in deren Sch¨atzungen berechnet, die durch Weglassen deri-ten Beobachtung entstehen.
Belsley, Kuh and Welsch (1980) definieren diese Gr¨oßen komponentenweise durch DFBETASj(i)= βˆj−βˆj(i)
ˆ σ(i)
q
(X0X)−1jj
, j= 1, . . . , p.
Dieses Mass wird groß genannt, falls DFBETASj(i)>2/√ ngilt.
dfbetas5 <- dfbetas(vf.mod.max)[ , 5]; n <- length(price)
is.inf5 <- as.numeric(rownames(houses)[abs(dfbetas5) > 2/sqrt(n)]) dfbetas10 <- dfbetas(vf.mod.max)[ ,10]
is.inf10 <- as.numeric(rownames(houses)[abs(dfbetas10) > 2/sqrt(n)]) plot(index, dfbetas5, ylab="DFBETAS(new:bath)")
abline(0, 0); abline(-2/sqrt(n), 0, lty=2); abline(2/sqrt(n), 0, lty=2) text(index[is.inf5], dfbetas5[is.inf5], is.inf5, pos=4, cex=0.8)
plot(index, dfbetas10, ylim=c(-0.5, 0.5), ylab="DFBETAS(new:bath:area)") abline(0, 0); abline(-2/sqrt(n), 0, lty=2); abline(2/sqrt(n), 0, lty=2) text(index[is.inf10], dfbetas10[is.inf10], is.inf10, pos=4, cex=0.8)
0 20 40 60 80
−1.0−0.50.00.5
index
DFBETAS(new:bath)
7 16
22
0 20 40 60 80
−0.4−0.20.00.20.4
index
DFBETAS(new:bath:area)
90 91
Die Beobachtung 5 hat bei beiden Parametern einen fehlenden Wert f¨ur das DFBETAS. Dies ist darin begr¨undet, dass ohne dieser Beobachtung die beiden Parameter gar nicht sch¨atzbar sind. Den st¨arksten Einfluss auf die Interaktionssch¨atzung im Intercept hat die Beobachtung 7, w¨ahrend die Wirkung auf die Steigungsinteraktion bei den Beobachtungen 92 und 93 mit Werten von−22.64 und 3.776 (nicht mehr dargestellt) am gr¨oßten ist.
Weiters hat die Beobachtung 7 den gr¨oßten Einfluss auf die Sch¨atzung der Standardabweichung der Response. Ohne diese einzelne Zeile k¨onnte ˆσ= 13.530 sogar auf ˆσ(7)= 12.870 reduziert werden mit
ˆ
σ(i)2 = 1 n−1−p
X
k6=i
(yk−x0kβˆ(i))2. Ahnliches gilt f¨ur die DFFITS Werte, wobei¨
DFFITSi=µˆi−µˆi(i)
ˆ σ(i)√
hii
. Auffallende Datenzeilen sind hierbei jene mit DFFITSi>2p
p/n. Wiederum resultiert ein fehlender Wert beim Weglassen der Beobachtung 5. Die Beobachtungen 92 und 93 sind mit Werten von−36.43 und 14.14
extrem und werden deshalb nicht dargestellt. Wiederum fallen dabei die Datens¨atze 4, 7, 90, 92, und 93 besonders auf.
dfsigma <- lm.influence(vf.mod.max)$sigma
is.infsigma <- as.numeric(rownames(houses)[abs(dfsigma) < 13.53]) plot(index, dfsigma)
abline(13.53, 0)
text(index[is.infsigma], dfsigma[is.infsigma], is.infsigma, pos=4, cex=0.8) dffits <- dffits(vf.mod.max)
is.inffits <- as.numeric(rownames(houses)[abs(dffits) > 2*sqrt(10/n)]) plot(index, dffits, ylim=c(-2, 2), ylab="DFFITS")
abline(0, 0)
abline(-2*sqrt(10/n), 0, lty=2); abline(2*sqrt(10/n), 0, lty=2)
text(index[is.inffits], dffits[is.inffits], is.inffits, pos=4, cex=0.8)
0 20 40 60 80
13.013.213.413.6
index
dfsigma
4
7
30 41
65 70
86 93
0 20 40 60 80
−2−1012
index
DFFITS
4
7 22
65 89
90
91
3.1 Zusammenfassen von Kategorien
All diese diagnostischen Methoden weisen darauf hin, dass H¨auser mit drei Badezimmer sehr problema- tisch im Modell sind und einen enormen Einfluss auf die Sch¨atzungen aufweisen. Die Datens¨atze 5 und 92 beschreiben die einzigen beiden neuen H¨auser mit drei Bedezimmern. Man sollte sich daher ¨uberlegen, anstelle des dreistufugen Faktors factor(bath)nur einen zweistufigen zu verwenden in dem die zweite und dritte Stufe zusammengefasst sind zur Stufemehr als ein Badezimmer. Dazu wird eine neue Variable bath12definiert mit
bath12 <- pmin(bath,2) table(bath12, new)
new bath12 0 1
1 10 0 2 55 28
Da eine Zelle gar nicht besetzt ist, kann keine Steigunsinteraktionen gesch¨atzt werden. Wir betrachten daher das zu Modell (7) entsprechende Modell ohne Interaktion zwischenbath12undnew, d.h.
E(pricei) =β1+βbath12i=2+βnewi=1+¡
βarea+βarea, bath12i=2+βarea, newi=1¢
∗areai. (10) Dieses entspricht dem maximalen Modell, beinhaltet jedoch davon nur die sch¨atzbaren Parameter. Als Sch¨atzung ergibt sich daf¨ur
summary(lm(price ~ (factor(new) + factor(bath12)) * area)) Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.67 19.22 -0.09 0.93083
factor(new)1 -37.44 15.34 -2.44 0.01668 * factor(bath12)2 -9.32 20.68 -0.45 0.65348
area 43.74 19.30 2.27 0.02593 *
factor(new)1:area 31.94 8.41 3.80 0.00027 ***
factor(bath12)2:area 20.32 19.78 1.03 0.30707 ---
Residual standard error: 16.2 on 87 degrees of freedom Multiple R-Squared: 0.873, Adjusted R-squared: 0.866 F-statistic: 120 on 5 and 87 DF, p-value: 0
Im obigensummarysieht man, dass der Wald-Test darauf hinweist, dass der neue Faktor bath12sowohl als Haupteffekt als auch in der Interaktion irrelevant zu sein scheint. Dies kann man testen mittels anova(lm(price ~ (factor(new) + factor(bath12)) * area), lm(price ~ factor(new) * area)) Analysis of Variance Table
Model 1: price ~ factor(new) + factor(bath12) + area + factor(new):area + factor(bath12):area
Model 2: price ~ factor(new) + area + factor(new):area Res.Df RSS Df Sum of Sq F Pr(>F)
1 87 22753
2 89 23794 -2 -1041 1.99 0.14
Wie bereits vermutet kann aufbath12im Modell verzichtet werden (p= 0.14). Somit folgt als einfaches Alternativmodell
vf.final <- lm(price ~ factor(new) * area) summary(vf.final)
Coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept) -16.60 6.21 -2.67 0.00894 **
factor(new)1 -31.83 14.82 -2.15 0.03445 *
area 66.60 3.69 18.03 < 2e-16 ***
factor(new)1:area 29.39 8.19 3.59 0.00055 ***
---
Residual standard error: 16.4 on 89 degrees of freedom Multiple R-Squared: 0.868, Adjusted R-squared: 0.863 F-statistic: 194 on 3 and 89 DF, p-value: 0
anova(vf.final)
Analysis of Variance Table Response: price
Df Sum Sq Mean Sq F value Pr(>F) factor(new) 1 22833 22833 85.4 1.2e-14 ***
area 1 129539 129539 484.5 < 2e-16 ***
factor(new):area 1 3439 3439 12.9 0.00055 ***
Residuals 89 23794 267
Dieses simple Modell hat mit nur vier Parameter sogar ein R2 = 0.87 und die gesch¨atzte Standard- abweichung der Preise betr¨agt 16.400$. Alle enthaltenen Effekte k¨onnen vollst¨andig gesch¨atzt werden und scheinen relevant zu sein. Es beschreibt durch zwei frei w¨ahlbare Geraden (gebraucht und neu) die mittleren Preise der H¨auser.
plot (area[new==0], price[new==0], xlim=c(0,4), xlab="area", ylim=c(0,310), ylab="fitted(vf.final)", pch="*", col=1) points(area[new==1], price[new==1], col=2)
abline(-16.60, 66.60) # new=0
abline(-16.60-31.83, 66.60+29.39, col=2) # new=1 extremes <- c(5, 7, 36, 41, 65, 89, 92, 93)
text(area[extremes], price[extremes], extremes, pos=4, cex=0.8, col=(new[extremes]+1)) h <- lm.influence(vf.final)$hat
is.inf <- as.numeric(rownames(houses)[h > 8/93]) plot(index, h, ylab="h(ii)"); abline(8/93, 0, lty=2) text(index[is.inf], h[is.inf], is.inf, pos=4, cex=0.8) r <- rstandard(vf.final)
plot(index, r); abline(0, 0)
abline( qnorm(0.975), 0, lty=2); abline( qnorm(0.025), 0, lty=2) is.large.r <- as.numeric(rownames(houses)[abs(r) > qnorm(0.975)]) text(index[is.large.r], r[is.large.r], is.large.r, pos=4, cex=0.8) rstud <- rstudent(vf.final)
plot(index, rstud); abline(0, 0)
abline( qnorm(0.975), 0, lty=2); abline( qnorm(0.025), 0, lty=2) is.large.r <- as.numeric(rownames(houses)[abs(rstud) > qnorm(0.975)]) text(index[is.large.r], r[is.large.r], is.large.r, pos=4, cex=0.8) dfbetas3 <- dfbetas(vf.final)[ , 3]
is.inf3 <- as.numeric(rownames(houses)[abs(dfbetas3) > 2/sqrt(n)]) dfbetas4 <- dfbetas(vf.final)[ , 4]
is.inf4 <- as.numeric(rownames(houses)[abs(dfbetas4) > 2/sqrt(n)]) plot(index, dfbetas3, ylab="DFBETAS(area)")
abline(0, 0); abline(-2/sqrt(n), 0, lty=2); abline(2/sqrt(n), 0, lty=2) text(index[is.inf3], dfbetas3[is.inf3], is.inf3, pos=4, cex=0.8)
plot(index, dfbetas4, ylim=c(-1, 3.5), ylab="DFBETAS(new:area)")
abline(0, 0); abline(-2/sqrt(n), 0, lty=2); abline(2/sqrt(n), 0, lty=2) text(index[is.inf4], dfbetas4[is.inf4], is.inf4, pos=4, cex=0.8)
dfsigma <- lm.influence(vf.final)$sigma
is.infsigma <- as.numeric(rownames(houses)[abs(dfsigma) < 16.2]) plot(index, dfsigma)
abline(16.4, 0)
text(index[is.infsigma], dfsigma[is.infsigma], is.infsigma, pos=4, cex=0.8) dffits <- dffits(vf.final)
is.inffits <- as.numeric(rownames(houses)[abs(dffits) > 2*sqrt(4/n)]) plot(index, dffits, ylim=c(-1, 4), ylab="DFFITS")
abline(0, 0)
abline(-2*sqrt(4/n), 0, lty=2); abline(2*sqrt(4/n), 0, lty=2)
text(index[is.inffits], dffits[is.inffits], is.inffits, pos=4, cex=0.8)
**
*
*
* * *
***
*
** *** ****** **** **
***
** * * *
** * ** *
***** * ** ****
** * *
** ** *
*
*
0 1 2 3 4
050100150200250300
area
fitted(vf.final)
5
7
36 41
65
89 92
93 5
36
65 92
5
7
36 41
65
89 92
93
0 20 40 60 80
0.00.10.20.30.40.5
index
h(ii)
5
6
36 89
92 93
0 20 40 60 80
−3−2−10123
index
r
5
7
41
65 89
93
0 20 40 60 80
−3−2−101234
index
rstud
5
7
41
65 89
93
Wiederum sind die Beobachtungen 5 und 93 diejenigen mit den gr¨oßten Werten von hii, welche jedoch weitaus kleiner sind als zuvor.
Die standardisierten Residuen sind alle sehr ¨ahnlich den studentisierten Residuen. Die beiden extrem teuren H¨auser 5 und 93 werden etwas durch das Modell untersch¨atzt und generieren daher positive Residuen. Das gebrauchte Haus Nummer 7 ist ¨außerst billig mit 19.600$ relativ zur Fl¨ache 1.28.
0 20 40 60 80
−0.50.00.51.01.5
index
DFBETAS(area) 7
41 76
89 93
0 20 40 60 80
−10123
index
DFBETAS(new:area)
5
36
65
89
93
0 20 40 60 80
15.415.615.816.016.216.4
index
dfsigma
5 7
30
41
65 70
89 90
93
0 20 40 60 80
−101234
index
DFFITS
5
7
36
65
89 93
Haus Nummer 5 ist neu und hat daher nur einen geringen Einfluss aus die Steigung aber einen großen auf die Steigungsdifferenz. Im Gegensatz dazu ist 93 ein altes Haus mit großem Einfluss auf die Steigung und kleinem auf die Steigungsdifferenz. Weiters f¨allt hier noch das gebrauchte Haus 89 auf, das auf einem großen Grundst¨uck steht (3.05) und daf¨ur relativ billig ist (149.000$).
Die H¨auser 5, 7, und 93 leisten auch einen großen Beitrag zur Sch¨atzung der Standardabweichung. Ohne diese Daten k¨onnte jeweils diese Sch¨azung um 1.145$, 828$, und 522$ reduziert werden.
Bez¨uglich der Erwartungswerte werden die Preise der H¨auser 5 und 93 viel h¨oher eingesch¨atzt wenn die beiden Beobachtungen zur Sch¨atzung verwendet werden, als wenn sie daf¨ur weggelassen werden. Beide ziehen also deren Regressionsgeraden jeweils stark nach oben. Im Gegensatz dazu hat das Haus 89 eine negative Tendenz auf den mittleren Hauspreis.
3.2 Alternative Faktorstufen-Kodierung
Bis jetzt haben wir uns mit Faktoren besch¨aftigt, die in der Treatment-Kodierung in das Modell einge- gangen sind. Die Parameter haben dabei jeweils die Differenz des Einflusses einer Stufe im Vergleich zum Einfluss einer Referenzstufe (Basisstufe) dieses Faktors beschrieben. Als Basis wird standardm¨assig die alphabetisch erste Stufe verwendet. F¨urbathist dies die Stufe1. Die Basis kann einfach ge¨andert werden.
M¨ochte manbath = 3als Basis haben, d.h. Unterschiede der Stufenbath = 1undbath = 2zur Stufe bath = 3 als Parameter im Modell haben, so verwende man C(factor(bath), base=3) anstelle von factor(bath)in der Modelldefinition.
Erinnern wir uns an das Modell (2) mit den beiden Faktoren als Haupteffekte summary(lm(price ~ C(factor(bath), base=1) + factor(new)))
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 40.320 8.888 4.536 1.78e-05 ***
C(factor(bath), base = 1)2 50.391 9.701 5.194 1.29e-06 ***
C(factor(bath), base = 1)3 152.615 14.677 10.398 < 2e-16 ***
factor(new)1 25.394 6.525 3.892 0.000192 ***
Diese Parametrisierung erlaubt es, die Abweichung zwischen bath=2 oder bath=3 von der Basisstufe bath=1(dem Intercept Parameter) zu verifizieren. So unterscheidet sich ein Haus mit zwei Badezimmern von einem Haus mit nur einem Bad preislich um 50.391$ (p < 0.001). Weiters sind H¨auser mit drei Badezimmern im Mittel um 152.615$ teurer als H¨auser mit nur einem Bad (p <0.001).
Nichts kann jedoch dar¨uber ausgesagt oder getestet werden, wie sichbath=2 vonbath=3 unterscheidet.
Dies liegt aber nur an der verwendeten Kodierung.
Die numerische Kodierung der Stufen erh¨alt man durch contrasts(factor(bath))
2 3 1 0 0 2 1 0 3 0 1
contrasts(factor(new)) 1
0 0 1 1
Die Zeilen beschreiben hierbei die einzelnen Kategorien des Faktors. Die Spalten beinhalten die zur Beschreibung aller Kategorien notwendigen Dummyvariablen zu diesem Faktor. Als Basisstufe versteht man jene Zeile, die nur Nullen hat. Zusammen mit dem Intercept bedeutet dies f¨ur ein Haus mit nur einem Bad, dass die entsprechende Zeile der DesignmatrixX in der Dummy2 und Dummy3 Spalte jeweils eine Null aufweist. Der Intercept beschreibt daher den Einfluß bei nur einem Badezimmer. Hat das Haus 3 Badezimmer, so ist eine Null in der Dummy2 Spalte und eine Eins in der Dummy3 Spalte. Als Basis vonnewwirdnew=0verwendet.
Alternativ kann nun das Modell mit contrasts(C(factor(bath), base=3))
1 2 1 1 0 2 0 1 3 0 0
contrasts(C(factor(new), base=1)) 2
0 0 1 1
gesch¨atzt werden. Daf¨ur resultiert
summary(lm(price ~ C(factor(bath), base=3) + C(factor(new), base=1))) Coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept) 192.94 11.68 16.52 < 2e-16 ***
C(factor(bath), base = 3)1 -152.62 14.68 -10.40 < 2e-16 ***
C(factor(bath), base = 3)2 -102.22 11.91 -8.58 2.8e-13 ***
C(factor(new), base = 1)2 25.39 6.53 3.89 0.00019 ***
Hier wird ersichtlich, dass sich bath=2 von der Basis bath=3 um−102.224$ unterscheidet (p <0.001) oder dass gebrauchte H¨auser um um 25.390$ teurer sind als neue.
Falls wir nicht an Unterschieden zwischen den Kategorien sondern an den Mittelwerten der Kategorien interessiert sind, so sollte man den Intercept aus dem Modell nehmen. Jedoch wird dadurch oft die ANOVA-Zerlegung nicht interpretierbar. Es folgt daf¨ur hier
summary(lm(price ~ factor(bath) + factor(new) - 1)) Coefficients:
Estimate Std. Error t value Pr(>|t|) factor(bath)1 40.320 8.888 4.536 1.78e-05 ***
factor(bath)2 90.711 3.888 23.332 < 2e-16 ***
factor(bath)3 192.935 11.679 16.520 < 2e-16 ***
factor(new)1 25.394 6.525 3.892 0.000192 ***
und E(price|new=0, bath=3) = 192.935$ oder E(price|new=1, bath=3) = 192.935 + 25.394 = 218.329$
k¨onnen (fast) direkt abgelesen werden.
Anstelle der Treatment-Kodierung kann auch die Summen-Kodierung eingesetzt werden. Am einfachsten schaltet man generell darauf um, indem die entsprechende Option neu gesetzt wird, d.h.
options(contrasts=c("contr.treatment", "contr.poly")) # default treatment-coding options(contrasts=c("contr.sum", "contr.poly")) # switch to sum-coding Dies f¨uhrt zu zwei andere Dummy Variablen mit
contrasts(factor(bath)) [,1] [,2]
1 1 0
2 0 1
3 -1 -1
Bemerke, dass hierbei alle Spaltensummen Null sind. Nat¨urlich hat dieses generelle Umschalten auch f¨ur den zweiten Faktornewdie Summen-Kodierung zur Folge
contrasts(factor(new)) [,1]
0 1
1 -1
Zusammen ergibt dies
summary(lm(price ~ factor(bath) + factor(new))) Coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept) 120.686 5.273 22.889 < 2e-16 ***
factor(bath)1 -67.669 7.281 -9.294 9.33e-15 ***
factor(bath)2 -17.278 5.341 -3.235 0.001707 **
factor(new)1 -12.697 3.263 -3.892 0.000192 ***
F¨ur ein Haus mit drei Badezimmer resultiert E(price|new=0, bath=3) = 120.686 + 67.669 + 17.278− 12.697 = 192.936$ oder E(price|new=1, bath=3) = 120.686 + 67.669 + 17.278 + 12.697 = 218.330$.
