1 Szenario

Ein große Modelagentur interessiert sich für die möglichen Faktoren, die das Einkommen von Supermodels vorhersagen. Es wurden 231 Supermodels nach ihrem Jahreseinkommen (in Mio Dollar) gefragt. Die erhobenen Prädiktoren sind Alter (in Jahren; “age” im Datensatz), die Berufsjahre (in Jahren; “years” im Datensatz) und die “Schönheit” (erfasst durch eine Umfrage; “beauty” im Datensatz).

Da man sich in der Modelagentur nicht sehr oft mit Statistik befasst, ist unklar, ob diese Maße geeignet sind, um das Einkommen von Supermodels gut zu erfassen.

Man geht davon aus, dass die Prädiktoren “Age” und “Schönheit” positiv mit dem Kriterium “Einkommen” zusammenhängen. Der Prädiktor “years” (also die Berufsjahre) sollten dagegen negativ mit dem Kriterium zusammenhängen.

2 Aufgaben

Die Daten befinden sich in der Datei “Daten_Supermodels.csv”

  1. Stellen Sie für die drei Hypothesen statistische Hypothesenpaare auf. Beziehen Sie dabei die Hypothesen auf die Regressionskoeffizienten.
  2. Berechnen Sie eine multiple Regressionsanalyse. Interpretieren Sie das Ergebnis.
  3. Berechnen Sie die vorhergesagten Werte, die Residuen, standardisierten und studentisierten Residuen. Fügen Sie diese der Datenmatrix hinzu.
  4. Prüfen Sie ob die folgenden Annahmen erfüllt sind:
  • Linearität der Zusammenhänge
  • Normalverteilung der Residuen
  • Homoskedastizität
  • keine Multikollinearität.

Interpretieren Sie das Ergebnis ihrer Analysen.

  1. Was bedeutet das Ergebnis der Prüfung der Annahmen für die Analyse der Daten?
# R-Commander laden

library(Rcmdr)
#library(RcmdrMisc)
library(lmtest)

3 Daten importieren und erste Zeilen anzeigen

Daten <- read.delim("Daten_Supermodels.csv")


head(Daten, 10) # dieser Code zeigt euch die ersten (deshalb "head") 10 Zeilen des Datensatzes.
##         salary      age    years   beauty
## 1   0.37039711 16.66699 3.148455 78.25149
## 2  53.72478935 20.34707 5.506886 68.56999
## 3   1.46015912 18.20307 5.330748 75.04376
## 4   0.02433401 15.35626 3.840088 65.14253
## 5  95.33807011 24.17183 8.532050 71.77039
## 6  14.63547789 18.26022 4.393158 78.05224
## 7   8.67333220 17.69861 4.396633 72.06817
## 8   2.64927429 17.48589 4.110735 75.32745
## 9   7.54732264 17.06954 3.515077 72.03374
## 10  1.20251467 20.07689 6.828727 71.93139

Wundern Sie sich nicht, dass die Werte so viele Nachkommastellen haben (da hatte wohl jemand ein sehr sehr genaues Messinstrument).

Ab hier können Sie jetzt das Script erweitern und selbständig weiterrechnen.

4 Hypothesen:

Es wird folgendes lineares Modell an die Daten angepasst: salary = beta0 + beta1 x age + beta2 x beauty + beta3 x years.

Für die statistischen Hypothesen über Regressionskoeffizienten ergibt sich daher:

  • H1,1: beta1 > 0, H0,1: beta1 <= 0
  • H1,2: beta2 > 0, H0,2: beta2 <= 0
  • H1,3: beta3 < 0, H0,3: beta3 >= 0

5 Analysen

5.1 Regressionsanalyse

5.1.1 Regressionsmodell erstellen

RegModel.1 <- lm(salary~age+beauty+years, data=Daten)
summary(RegModel.1)
## 
## Call:
## lm(formula = salary ~ age + beauty + years, data = Daten)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.853  -7.950  -4.197   4.605  68.085 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -60.8897    16.4966  -3.691  0.00028 ***
## age           6.2344     1.4112   4.418 1.54e-05 ***
## beauty       -0.1964     0.1524  -1.289  0.19871    
## years        -5.5612     2.1222  -2.621  0.00937 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.57 on 227 degrees of freedom
## Multiple R-squared:  0.184,  Adjusted R-squared:  0.1733 
## F-statistic: 17.07 on 3 and 227 DF,  p-value: 4.973e-10

Entscheidung über die statistischen Hypothesen:

  • H1,1 angenommen (beta 1 ist positiv)
  • H0,2 beibehalten (beta2 ist nicht signifikant geworden)
  • H1,3 angenommen (beta3 ist negativ) - aber Achtung, der Zusammenhang sieht im Graphen eigentlich positiv aus. Daher jetzt mal aufmerksam die Annahmen prüfen :) ).

Reminder: der t-test im Output erfolgt ungerichtet. Für gerichtete Hypothesen, die auch in der vorhergesagten Richtung bestätigt werden konnten, darf der p-Wert also noch halbiert werden.

Koeffizienten standardisieren

library(QuantPsyc)

lm.beta(RegModel.1)
##         age      beauty       years 
##  0.94214234 -0.08299604 -0.54779846

6 Prüfung der Voraussetzungen

6.1 Berechnen der vorhergesagten Werte, der Residuen, standardisierten und studentisierten Residuen. Werte zu den Daten hinzufügen.

Daten<- within(Daten, {
  fitted.RegModel.1 <- fitted(RegModel.1)
  residuals.RegModel.1 <- residuals(RegModel.1)
  rstudent.RegModel.1 <- rstudent(RegModel.1)
  hatvalues.RegModel.1 <- hatvalues(RegModel.1)
  cooks.distance.RegModel.1 <- cooks.distance(RegModel.1)
  obsNumber <- 1:nrow(Daten) 
})

(Hier wurden zusätzlich auch die hat values und cook’s distance hinzugefügt.)

Was noch fehlt und nicht über “Füge Regressionsstatistiken” hinzugefügt werden kann, sind die standardisierten Residuen. Wir haben aber bereits die Residuen hinzugefügt und wir können diese jetzt einfach selbst standardisieren und als neue Spalte ergänzen.

Das geht in R-Commander über Datenmanagement.

Daten <- local({
  .Z <- scale(Daten[,c("residuals.RegModel.1")])
  within(Daten, {
    Z.residuals.RegModel.1 <- .Z[,1] 
  })
})

6.2 Prüfen der Annahmen

6.2.1 Linearität der Zusammenhänge

Als erstes: Streudiagrammmatrix:

scatterplotMatrix(~age+beauty+salary+years, regLine=TRUE, smooth=list(span=0.5, spread=FALSE), diagonal=list(method="density"), data=Daten)

Als zweites schauen wir uns einen Plot an, der die vorhergesagten und die tatsächlichen Werte anzeigt:

scatterplot(salary~fitted.RegModel.1, regLine=TRUE, smooth=list(span=0.5, spread=FALSE), boxplots=FALSE, data=Daten)

als statistischen Test für die Prüfung der Annahme gibt es den “Reset-Test”:

Logik: Regression der quadrierten Residuen auf Potenzen der vorhergesagten Werte des Kriteriums

library(zoo, pos=20)
library(lmtest, pos=20)
resettest(salary ~ age + beauty + years, power=2:3, type="regressor", 
  data=Daten)
## 
##  RESET test
## 
## data:  salary ~ age + beauty + years
## RESET = 2.2335, df1 = 6, df2 = 221, p-value = 0.04106

Test ist signifikant, d.h., Zusammenhänge scheinen nicht gut durch lineares Modell beschreibbar zu sein

6.2.2 Normalverteilung der Residuen

grafische Analysen:

  1. Histogramm der Residuen
  2. Q-Q Plot
# Histogramm
with(Daten, Hist(Z.residuals.RegModel.1, scale="frequency", breaks="Sturges", col="darkgray"))

Den ersten QQ-Plot erhält man im Commander unter Grafiken -> QQ-Diagramm. In der Liste die Variable für die Residuen auswählen.

Im Video haben wir das Vorgehen über Modelle -> Grafiken -> Residuen-Quantil-Vergleichsgrafik gezeigt. Hier erhält man auch einen QQ-Plot, allerdings werden die studentisierten Residuen mit den Quantilen der t-Verteilung verglichen. Das kann man auch machen. Der dadurch erstellte Code ist der zweite im nächsten Chunk.

# Q-Q Plot über Grafiken...
with(Daten, qqPlot(residuals.RegModel.1, dist="norm", id=list(method="y", n=2, labels=rownames(Daten))))

## [1] 135   5
# Q-Q Plot über Modelle...
qqPlot(RegModel.1, simulate=TRUE, id=list(method="y", n=2))

## [1]   5 135

Statistisch können wir die Normalverteilung der Residuen mit dem “Shapiro-Wilk” Test prüfen:

shapiro.test(Daten$residuals.RegModel.1)
## 
##  Shapiro-Wilk normality test
## 
## data:  Daten$residuals.RegModel.1
## W = 0.84745, p-value = 2.462e-14

Der Test ist signifikant; die H0 (“die Residuen sind normalverteilt”) wird also abgelehnt.

6.2.3 Homoskedastizität

Homoskedastizität = Gleichheit der Varianzen des Fehlers und damit des Kriteriums bei allen Ausprägungen der Prädiktoren.

Heißt: der Fehler ist unabhängig von den Prädiktoren und hat eine konstante Streuung.

Bei Verletzung dieser Annahme werden die Standardfehler der Koeffizienten nicht korrekt geschätzt.

Für eine grafische Analyse bietet es sich an, die Modellvorhersagen gegen den Fehler in einem Streudiagramm abzutragen:

scatterplot(residuals.RegModel.1~fitted.RegModel.1, regLine=TRUE, smooth=FALSE, boxplots=FALSE, data=Daten)

Man sieht hier, dass die vorhergesagten Werte positiv mit den Residuen korrelieren. Je höher die Prädiktorwerte, umso höher der Fehler.

Der dazugehörige statistische Test ist der Breusch-Pagan Test.

Grundidee der Testung: Regression der (quadrierten) Residuen auf die Prädiktoren. In anderen Worten: Regressionsanalyse, die schaut, ob die Prädiktoren die Residuen des ursprünglichen Modells vorhersagen. (Wenn Homoskedastizität vorliegt, dann sollte das nicht so sein.)

bptest(salary ~ age + beauty + years, varformula = ~ fitted.values(RegModel.1), studentize=FALSE, data=Daten)
## 
##  Breusch-Pagan test
## 
## data:  salary ~ age + beauty + years
## BP = 97.131, df = 1, p-value < 2.2e-16

Der Test ist signifikant, d.h., wir nehmen die H1 (“es liegt Heteroskedastizität vor”) an.

6.2.4 keine Multikollinearität

Dazu sollten wir uns zunächst eine Streudiagramm-Matrix ansehen.

scatterplotMatrix(~age+beauty+salary+years, regLine=TRUE, smooth=list(span=0.5, spread=FALSE), diagonal=list(method="density"), data=Daten)

Statistisch sollte man sich die VIFs (Variance Inflation Factors) ansehen.

Logik: Jeder Prädiktor im Modell wird einmal als Kriterium genommen und es wird analysiert, ob dessen Varianz durch die anderen (übrigen) Prädiktoren erklärt werden kann (wenn ja, dann korrlieren die Prädiktoren). (s. Folie 12)

Mittelwert der VIFs sollte ungefähr 1 sein und jeder einzelne VIF sollte < 10 sein.

vif(RegModel.1)
##       age    beauty     years 
## 12.652841  1.153364 12.156757

Mittelwert der VIFs:

mean(vif(RegModel.1))
## [1] 8.65432

Wir können auch checken, wie groß die Korrelation zwischen “age” und “years” ist.

cor.test(Daten$age, Daten$years)
## 
##  Pearson's product-moment correlation
## 
## data:  Daten$age and Daten$years
## t = 48.592, df = 229, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9417555 0.9649314
## sample estimates:
##       cor 
## 0.9547716

6.2.5 Ausreißer und zu einflussreiche Messwerte

Hier sollte man Cook’s distances bestimmen, da dies standardisierte Werte sind.

Grundidee hinter der Analyse: Regression jedes vorhergesagten Werts auf alle beobachteten Werte. Die beobachteten Werte werden sozusagen zu Prädiktoren der vorhergesagten Werte. Wenn einzelne Werte einen zu hohen Einfluss haben, dann kriegen sie ein deutlich größeres Regressionsgewicht als die übrigen Werte.

Orientierung: Werte > 1 sollte man überprüfen. Werte > 4 haben einen zu großen Einfluss.

Grafische Analyse über Histogramm:

with(Daten, Hist(cooks.distance.RegModel.1, scale="frequency", breaks="Sturges", col="darkgray"))

Alle Werte sind < 1. Es scheint also keine modellverzerrenden Ausreißerwerte zu geben.

Andere Möglichkeit der Prüfung: Bestimmung des Maximalwerts von Cook’s D

max(Daten$cooks.distance.RegModel.1)
## [1] 0.226934
Daten$cooks.distance.RegModel.1 >= 4 
##   [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [37] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [49] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [85] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [97] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [109] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [121] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [133] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [145] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [157] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [169] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [181] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [193] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [205] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [217] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [229] FALSE FALSE FALSE

7 Neue Analyse ohne den Prädiktor “age”

Da die beiden Prädiktoren age und years stark korrelieren, bietet es sich an, einen der beiden aus dem Modell zu entfernen. Das ist eine inhaltliche Entscheidung (siehe Erläuterungen am Ende des Markdowns).

7.1 Reg-Modell anpassen

RegModel.1 <- lm(salary~years+beauty, data=Daten)
summary(RegModel.1)
## 
## Call:
## lm(formula = salary ~ years + beauty, data = Daten)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -21.293  -8.537  -4.464   4.905  70.646 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6.02001   11.28918  -0.533    0.594    
## years        3.40763    0.64263   5.303  2.7e-07 ***
## beauty       0.02282    0.14978   0.152    0.879    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.15 on 228 degrees of freedom
## Multiple R-squared:  0.1139, Adjusted R-squared:  0.1061 
## F-statistic: 14.65 on 2 and 228 DF,  p-value: 1.033e-06

7.2 Berechnen der vorhergesagten Werte, der Residuen, standardisierten und studentisierten Residuen. Werte zu den Daten hinzufügen.

Daten<- within(Daten, {
  fitted.RegModel.1 <- fitted(RegModel.1)
  residuals.RegModel.1 <- residuals(RegModel.1)
  rstudent.RegModel.1 <- rstudent(RegModel.1)
  hatvalues.RegModel.1 <- hatvalues(RegModel.1)
  cooks.distance.RegModel.1 <- cooks.distance(RegModel.1)
  obsNumber <- 1:nrow(Daten) 
})

Was noch fehlt und nicht über “Füge Regressionsstatistiken” hinzugefügt werden kann, sind die standardisierten Residuen. Wir haben aber bereits die Residuen hinzugefügt und wir können diese jetzt einfach selbst standardisieren und als neue Spalte ergänzen.

Das geht in R-Commander über Datenmanagement.

Daten <- local({
  .Z <- scale(Daten[,c("residuals.RegModel.1")])
  within(Daten, {
    Z.residuals.RegModel.1 <- .Z[,1] 
  })
})

7.3 Prüfen der Annahmen

7.3.1 Linearität der Zusammenhänge

Als erstes: Streudiagrammmatrix:

scatterplotMatrix(~years+beauty+salary, regLine=TRUE, smooth=list(span=0.5, spread=FALSE), diagonal=list(method="density"), data=Daten)

Als zweites schauen wir uns einen Plot an, der die vorhergesagten und die tatsächlichen Werte anzeigt:

scatterplot(salary~fitted.RegModel.1, regLine=TRUE, smooth=list(span=0.5, spread=FALSE), boxplots=FALSE, data=Daten)

als statistischen Test für die Prüfung der Annahme gibt es den “Reset-Test”:

Logik: Regression der quadrierten Residuen auf Potenzen der vorhergesagten Werte des Kriteriums

resettest(salary ~ beauty + years, power=2:3, type="regressor", 
  data=Daten)
## 
##  RESET test
## 
## data:  salary ~ beauty + years
## RESET = 2.1959, df1 = 4, df2 = 224, p-value = 0.07035

Test ist (knapp) nicht signifikant, d.h., Zusammenhänge scheinen durch lineares Modell beschreibbar zu sein

7.3.2 Normalverteilung der Residuen

grafische Analysen:

  1. Histogramm der Residuen
  2. Q-Q Plot
# Histogramm
with(Daten, Hist(Z.residuals.RegModel.1, scale="frequency", breaks="Sturges", col="darkgray"))

# Q-Q Plot
with(Daten, qqPlot(Z.residuals.RegModel.1, dist="norm", id=list(method="y", n=2, labels=rownames(Daten))))

## [1]   5 135

Statistisch können wir die Normalverteilung der Residuen mit dem “Shapiro-Wilk” Test prüfen:

shapiro.test(Daten$residuals.RegModel.1)
## 
##  Shapiro-Wilk normality test
## 
## data:  Daten$residuals.RegModel.1
## W = 0.82446, p-value = 1.932e-15

Der Test ist signifikant; die H0 (“die Residuen sind normalverteilt”) wird also abgelehnt.

7.3.3 Homoskedastizität

Homoskedastizität = Gleichheit der Varianzen des Fehlers und damit des Kriteriums bei allen Ausprägungen der Prädiktoren.

Heißt: der Fehler ist unabhängig von den Prädiktoren und hat eine konstante Streuung.

Bei Verletzung dieser Annahme werden die Standardfehler der Koeffizienten nicht korrekt geschätzt.

Für eine grafische Analyse bietet es sich an, die Modellvorhersagen gegen den Fehler in einem Streudiagramm abzutragen:

scatterplot(residuals.RegModel.1~fitted.RegModel.1, regLine=TRUE, smooth=FALSE, boxplots=FALSE, data=Daten)

Man sieht hier, dass die vorhergesagten Werte positiv mit den Residuen korrelieren. Je höher die Prädiktorwerte, umso höher der Fehler.

Der dazugehörige statistische Test ist der Breusch-Pagan Test.

Grundidee der Testung: Regression der (quadrierten) Residuen auf die Prädiktoren. In anderen Worten: Regressionsanalyse, die schaut, ob die Prädiktoren die Residuen des ursprünglichen Modells vorhesagen. (Wenn Homoskedastizität vorliegt, dann sollte das nicht so sein.)

bptest(salary ~ age + beauty + years, varformula = ~ fitted.values(RegModel.1), studentize=FALSE, data=Daten)
## 
##  Breusch-Pagan test
## 
## data:  salary ~ age + beauty + years
## BP = 74.344, df = 1, p-value < 2.2e-16

Der Test ist signifikant, d.h., wir nehmen die H1 (“es liegt Heteroskedastizität vor”) an.

7.3.4 keine Multikollinearität

Dazu sollten wir uns zunächst eine Streudiagramm-Matrix ansehen.

scatterplotMatrix(~age+beauty+salary+years, regLine=TRUE, smooth=list(span=0.5, spread=FALSE), diagonal=list(method="density"), data=Daten)

Statistisch sollte man sich die VIFs (Variance Inflation Factors) ansehen.

Logik: Jeder Prädiktor im Modell wird einmal als Kriterium genommen und es wird analysiert, ob dessen Varianz durch die anderen (übrigen) Prädiktoren erklärt werden kann (wenn ja, dann korrlieren die Prädiktoren). (s. Folie 12)

Mittelwert der VIFs sollte ungefähr 1 sein und jeder einzelne VIF sollte < 10 sein.

vif(RegModel.1)
##    years   beauty 
## 1.031029 1.031029

Mittelwert der VIFs:

mean(vif(RegModel.1))
## [1] 1.031029

Wir können auch checken, wie groß die Korrelation zwischen “years” und “beauty” ist (den verbleibenden zwei Prädiktoren).

cor.test(Daten$beauty, Daten$years)
## 
##  Pearson's product-moment correlation
## 
## data:  Daten$beauty and Daten$years
## t = 2.6656, df = 229, p-value = 0.008232
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.04541768 0.29592953
## sample estimates:
##       cor 
## 0.1734784

7.3.5 Ausreißer und zu einflussreiche Messwerte

Hier sollte man Cook’s distances bestimmen, da dies standardisierte Werte sind.

Grundidee hinter der Analyse: Regression jedes vorhergesagten Werts auf alle beobachteten Werte. Die beobachteten Werte werden sozusagen zu Prädiktoren der vorhergesagten Werte. Wenn einzelne Werte einen zu hohen Einfluss haben, dann kriegen sie ein deutlich größeres Regressionsgewicht als die übrigen Werte.

Orientierung: Werte > 1 sollte man überprüfen. Werte > 4 haben einen zu großen Einfluss.

Grafische Analyse über Histogramm:

with(Daten, Hist(cooks.distance.RegModel.1, scale="frequency", breaks="Sturges", col="darkgray"))

Alle Werte sind < 1. Es scheint also keine modellverzerrenden Ausreißerwerte zu geben.

Andere Möglichkeit der Prüfung: Bestimmung des Maximalwerts von Cook’s D

max(Daten$cooks.distance.RegModel.1)
## [1] 0.2844381
Daten$cooks.distance.RegModel.1 >= 4 
##   [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [37] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [49] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [85] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [97] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [109] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [121] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [133] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [145] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [157] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [169] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [181] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [193] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [205] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [217] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [229] FALSE FALSE FALSE

8 Fazit

In den Daten liegen mehrere Eigenschaften vor, die die Anpassung eines linearen Modells erschweren.

Zunächst (im ersten Modell) waren die Prädiktoren “age” und “years” zu hoch korreliert (je älter die Person, desto länger im Geschäft). Dadurch werden die Regressionskoeffizienten unzuverlässig geschätzt (siehe VL). Im Output zeigt sich das zum Beispiel daran, dass “years” angeblich einen negativen Zusammenhang mit salary aufweist, während sich in der Streudiagrammmatrix jedoch eher ein positiver Zusammenhang zeigt. Zudem erschwert es die inhaltliche Interpretation, da die beiden Variablen stark konfundiert sind: Ist das Einkommen zu einem bestimmten Zeitpunkt nun aufgrund des Alters oder aufgrund der Jahre im Geschäft höher?

Aus inhaltlichen Erwägungen (vermutlich sind es eher die Jahre im Geschäft, die durch Erfahrung, Kontakte etc. zu einem höheren Einkommen führen als das reine Alter), haben wir uns anschließend dazu entschlossen, den Prädiktor “age” aus dem Modell zu entfernen. Das Problem der Multikollinearität wird damit gelöst (auf statistischer Ebene - inhaltlich können wir nur mutmaßen, dass years tatsächlich der wichtigere Prädiktor ist).

Zudem scheinen die Zusammenhänge zwischen den Prädiktoren und salary (im Originalmodell) zumindest teilweise nicht linear zu verlaufen. Bei höheren Werten der Prädiktoren scheint das Einkommen plötzlich stärker als nur linear anzusteigen. Das sehen wir im Streudiagramm des Kriteriums in Abhängigkeit von den vorhergesagten Werten, es deutet sich aber auch bereits in der Streudiagramm-Matrix an.

Der RESET-Test auf Linearität bestätigt diesen Eindruck.

Wenn die Zusammenhänge zwischen Prädiktoren und Kriterium nicht linear sind, ist zu erwarten, dass sich das auch auf die Residuen niederschlägt. Und genau das sehen wir, wenn wir die Annahmen der Normalverteilung der Residuen und der Homoskedastizität prüfen (graphisch und statistisch). Die Residuen sind nicht normalverteilt, es gibt zu viele positive Residuen. Das liegt daran, dass das Modell die tatsächlichen Kriteriumswerte in einem bestimmten Bereich systematisch unterschätzt (Kriteriumswerte zu oft über der Regressionsgeraden -> viele positive Residuen, die nicht durch genug negative ausgeglichen werden).

Die Annahme der Homoskedastizität ist ebenfalls verletzt: Im niedrigen Bereich der vorhergesagten Werte macht das Modell einen viel kleineren Fehler als im höheren Bereich. Das ergibt sich ebenfalls aus dem bereits beschriebenen Grund.

Wir können also festhalten, dass ein lineares Modell nicht gut geeignet ist, um den Zusammenhang zwischen Jahren im Geschäft, beauty und Einkommen zu erklären. Stattdessen sollte man ein Modell wählen, dass den plötzlich viel stärkeren Anstieg des Einkommens ab einer bestimmten Anzahl von Geschäftsjahren einfangen kann. Hierfür könnte beispielsweise eine Exponentialfunktion geeignet sein. Bei dieser Gelegenheit könnte man auch beauty aus dem Modell entfernen, da es sich ja bisher nicht als signifikanter Prädiktor erwies (zur Übung kann bezüglich des Prädiktors beauty natürlich auch nochmal ein Vergleich der bisherigen linearen Modelle berechnet werden, siehe letzte Woche).

9 Exkurs: Warum so viele separate Annahmen prüfen, wenn ein einziges Problem in den Daten (der nicht-lineare Zusammenhang) sich hier in der Verletzung von gleich drei Annahmen spiegelt?

Manche Annahmen sind relativ(!) unabhängig von den anderen (zb Multikollinearität oder ob einzelne Ausreißer vorliegen, die sich nicht groß auf das Muster der Residuen niederschlagen), andere hängen oft zusammen (zb Linearität und NV der Residuen sowie Homoskedastizität). Wenn bei letzteren eine der Annahmen verletzt ist, sind häufig auch die anderen verletzt - aber nicht immer.

Man kann sich zb ein Modell vorstellen, bei dem alle Residuen entweder +5 oder -5 sind. Die Hälfte der tatsächlichen Daten liegt also 5 Punkte über der Regressionsgeraden und die andere Hälfte liegt 5 Punkte darunter. Ob ein Wert +5 oder -5 ist, kann dabei völlig unabhängig von den Ausprägungen der Prädiktoren sein. In diesem Fall wäre Homoskedastizität gegeben, weil die Streuung der Residuen unabhängig von den Ausprägungen der Prädiktoren ist. Allerdings wären die Residuen nicht normalverteilt, sondern bimodal (die Hälfte ist +5, die andere Hälfte -5). Inhaltlich gäbe es dann aus irgendwelchen Gründen zwei Subgruppen von Leuten, wobei der Faktor, der diese Gruppen unterscheidet, offenbar noch nicht identifiziert und ins Modell aufgenommen wurde (denn dann wäre dieser Unterschied ja schon erklärt und würde nicht als Fehler auftauchen). Die Linearität der Zusammenhänge könnte formal trotzdem gegeben sein, weil die Linie, die genau zwischen diesen beiden Subgruppen liegt, die Daten immer noch besser beschreibt als ein Modell mit quadratischen oder kubischen Komponenten. In diesem Fall würde man also nur anhand der Prüfung der Normalverteilung der Residuen merken, dass etwas nicht stimmt.

Weiterhin könnten in einem Modell die Residuen auch normalverteilt sein, aber keine konstante Streuung aufweisen. Das könnte passieren, wenn die vielen Werte, die nur wenig von der Modellvorhersage abweichen, alle bei den niedrigen Werten der Prädiktoren liegen. Die wenigen Werte mit großen Abweichungen könnten aber alle bei den hohen Bereichen der Prädiktoren liegen. Wenn man nur die Normalverteilung der Residuen prüfen würde (die uns zeigt, grob gesagt, wie viele große und kleine Residuen es gibt, unabhängig von ihrer “Lage”), würde man das nicht feststellen können.