In martscht/PsyBSc7: Übungen für die Praktika im Sommersemester 2020

library(learnr)
library(fontawesome)
#source('setup.R')
library(psych)

Modelloptimierung

Bei der Regressionsanalyse hat die Modelloptimierung zum Ziel, ein Regresionsmodell zu verbessern, das heißt, möglichst viel Varianz der abhängigen Variable zu erklären. Modelloptimierung bedeutet, ein Modell zu verbessern, durch:

• Aufnehmen zusätzlicher, bedeutsamer Prädiktoren
• Ausschließen von Prädiktoren, die nicht zur Varianzaufklärung beitragen

Ziel:

• sparsames Modell, in dem jeder enthaltene Prädiktor einen Beitrag zur Varianzaufklärung des Kriteriums leistet
• kein wichtiger Prädiktor vergessen

Dafür kann man beispielsweise das Inkrement, also den Zuwachs an erklärter Kriteriumsvarianz durch Hinzunahme eines Prädiktors, sowie das Dekrement, also die Verringerung an erklärter Varianz durch Ausschluss eines Prädiktors betrachten.

Methoden:

• Forward-Selection: Hinzunahme des Prädiktors mit dem größten, signifikanten Inkrement (bis keine zusätzliche Varianzaufklärung durch weitere Hinzunahme erreicht werden kann)
• Backward-Selection: Ausschluss von Prädiktoren mit dem kleinsten Dekrement, ohne die erklärte Kriteriumsvarianz signifikant zu reduzieren
• Stepwise Regression: Wie bei der Forward-Selection wird immer der Prädiktor mit dem größten Inkrement hinzugenommen, und zusätzlich in jedem Schritt geprüft, ob ein anderer Prädiktor wieder gelöscht werden kann (nicht-signifikantes Dekrement). Problematisch, wird nicht empfohlen

Übungs-Datensatz

Die Modelloptimierung wird am gleichen Datensatz demonstriert, der auch in der Sitzung zu Regression I verwendet wurde. Eine Stichprobe von 100 Schülerinnen und Schülern hat einen Lese- und einen Mathematiktest bearbeitet, und zusätzlich dazu einen allgemeinen Intelligenztest absolviert. Im Datensatz enthalten ist zudem das Geschlecht (Variable: female, 0 = m, 1 = w). Die abhängige Variable ist die Matheleistung, die durch die anderen Variablen im Datensatz vorhergesagt werden soll.

Über das Paket PsyBSc7 laden Sie den Datensatz wie gehabt wie folgt:

# Datensatz laden
data(Schulleistungen, package = "PsyBSc7")

Der Vollständigkeit halber und zur Erinnerung: Wenn Sie mit einem Datensatz arbeiten wollen, der lokal auf ihrem Rechner gespeichert ist, könnten sie diesen wie folgt laden:

# Datensatz laden
Schulleistungen <- readRDS("Schulleistungen.rda")

Inkrement und Dekrement

Testen eines Inkrements

Vergleich eines eingeschränkten Modells $M_c$ mit weniger Prädiktoren gegen ein uneingeschränktes Modell $M_u$ mit zusätzlichen Prädiktoren. Beispiel: Inkrement durch Hinzunahme von Intelligenz. Führt die Hinzunahme von Intelligenz zu einem signifikanten Zuwachs an erklärter Varianz, wenn Lesekompetenz und Geschlecht bereits im Vorhersagemodell sind?

m.c <- lm(math ~ reading + female, data = Schulleistungen)      # constrained
m.u <- lm(math ~ reading + female + IQ, data = Schulleistungen) # unconstrained
summary(m.c)
summary(m.u)

• Determinationskoeffizient ohne IQ ("constrained"): $R_c^2=r sprintf("%.2f",summary(m.c)$r.squared)$
• Determinationskoeffizient mit IQ ("unconstrained"): $R_u^2=r sprintf("%.2f",summary(m.u)$r.squared)$

Die Differenz aus beiden Determinationskoeffizienten können wir wie folgt berechnen:

# Inkrement = Differenz in R2 aus restringiertem Modell 2 minus R2 aus unrestringiertem Modell 1
summary(m.u)$r.squared - summary(m.c)$r.squared

Wir sehen also, dass das Modell mit IQ mehr Varianz erklärt als das ohne IQ. Die Hinzunahme des IQ führt zu einem Zuwachs von r round((summary(m.u)$r.squared - summary(m.c)$r.squared)*100)% erklärter Varianz ($\Delta R^2=r round(summary(m.u)$r.squared - summary(m.c)$r.squared, 3)$). Dieses Inkrement soll auf Signifikanz getestet werden. Der Modellvergleich kann mit der anova-Funktion vorgenommen werden. Ergänzen Sie den Befehl, um die beiden zuvor erstellten Modelle zu vergleichen.

# Modellvergleich mit der anova-Funktion
anova(...)

Das Inkrement des IQs ist auf einem Alpha-Fehlerniveau von 0.05 signifikant von null verschieden ($p=r sprintf("%.3f",anova(m.c, m.u)[[6]][2])$).

Zum Vergleich finden Sie hier die Berechnung des F-Tests nach Vorlesungs-Folie 11 aus der zweiten Vorlesung zu Regression:

R2.u <- summary(m.u)$r.squared
R2.c <- summary(m.c)$r.squared
df.diff <- summary(m.u)$df[1] - summary(m.c)$df[1]
df.u <- summary(m.u)$df[2]
F.diff <- ((R2.u - R2.c) / df.diff) /
  ((1 - R2.u) / df.u)
p.diff <- 1-pf(F.diff, df.diff, df.u)
F.diff
p.diff

Testen eines Dekrements

Das Dekrement ist der Unterschied im $R^2$ zwischen dem restringierten Modell $M_c$ und dem unrestringierten Modell $M_u$. Die Testung eines Dekrements erfolgt analog dem Inkrement: das eingeschränkte Modell $M_c$ mit weniger Prädiktoren wird mit dem uneingeschränkten Modell $M_u$ mit mehr Prädiktoren verglichen. Es soll nun geprüft werden, ob das Weglassen des Geschlechts aus dem Modell zu einem signifikanten Rückgang der erklärten Varianz führt. Stellen Sie nach dem Modell von oben das uneingeschränkte und das eingeschränkte Modell auf, berechnen Sie die Differenz in der Varianzaufklärung und testen Sie, ob dieser Unterschied signifikant ist.

m.u <- lm(math ~ ..., data = Schulleistungen) # unconstrained
m.c <- lm(math ~ ..., data = Schulleistungen) # constrained

summary(m.u)$r.squared - summary(m.c)$r.squared
# Modellvergleich mit der anova-Funktion
anova(m.c, m.u)

Der Ausschluss des Geschlechts führt zu einer Verringerung von r round((summary(m.u)$r.squared - summary(m.c)$r.squared)*100)% erklärter Varianz ($\Delta R^2=r round(summary(m.u)$r.squared - summary(m.c)$r.squared, 3)$). Dieser Unterschied ist nicht signifikant von null verschieden ($p=r sprintf("%.3f",anova(m.c, m.u)[[6]][2])$)

Schrittweise Selektion von Prädiktoren

Eine "theoriefreie" schrittweise Auswahl von Prädiktoren kann in R mit der step-Funktion erfolgen. Diese macht, anders als unser zuvor demonstriertes Vorgehen, nicht von Partialkorrelationen und Inkrementen Gebrauch, sondern vom sogenannten Informationskriterium AIC (Akaike Information Criterion). Dieses basiert auf der Likelihood eines geschätzten Modells $L(\hat{\theta})$ und der Anzahl der Modellparametern $p$:

$AIC=-2L(\hat{\theta}) + 2p$

Die Likelihood bezeichnet ein Maß für die Plausibilität/Wahrscheinlichkeit eines Modells, unter Berücksichtigung der gegebenen (empirisch erhobenen) Daten. Anders ausgedrückt: Die Likelihood eines Modells, gegeben die empirischen Daten. Um das beste Modell zu finden, kann man die Likelihood verschiedener Modelle vergleichen.

Für lineare Regressionsmodelle lässt sich der AIC wie folgt darstellen:

$AIC_{\sigma}=n \cdot log(\sigma_e^2) + 2p$

Der AIC ist hier eine Funktion der Stichprobengröße $n$, der Residualvarianz $\sigma_e^2$ und der Anzahl der Parameter (= Regressionskoeffizienten) $p=m+1$. Es wird hier ersichtlich, dass der AIC von der Varianz der abhängigen Variablen abhängt, da diese wiederum die Residualvarianz beeinflusst.

Der AIC ist ein sogenanntes inverses Maß, das bedeutet, dass Modelle mit einem kleineren AIC besser sind als Modelle mit einem größeren AIC. Der AIC wird durch den Term $n \cdot log(\sigma_e^2)$ kleiner, wenn die Residualvarianz kleiner wird, also mehr Varianz erklärt wird. Durch den "Strafterm" $2p$ wird der AIC größer, wenn das Modell mehr Prädiktoren enthält. Es soll also ein Modell gefunden werden, das mit möglichst wenigen Prädiktoren möglichst viel Varianz erklärt (Sparsamkeitsprinzip).

Die Schrittweise Selektion kann "vorwärts", "rückwärts", oder in beide Richtungen erfolgen. Die Standardeinstellung der step-Funktion ist die, dass ein Modell mit allen möglichen Prädiktoren als Ausgangspunkt genommen wird. Es wird dann der Prädiktor ausgeschlossen, der die größte Reduktion des AIC erlaubt, dann der nächste usw. In jedem Schritt wird auch wieder geprüft, ob Prädiktoren, die nicht im Modell sind, bei Aufnahme wieder zu einer Reduktion des AIC führen würden. Das Verfahren stoppt, wenn:

1. nur noch Prädiktoren im Modell sind, deren Ausschluss zu einer Erhöhung des AIC führen würden und
2. nur Prädiktoren "übrig" sind, deren Einschluss den AIC nicht verbessern würde.

Einfaches Beispiel: Optimierung des Modells für Mathematikleistung, Start mit allen drei möglichen Prädiktoren:

# Modell mit allen Prädiktoren
m <- lm(math ~ reading + female + IQ, data = Schulleistungen)
# Optimierung
summary(step(m))

Es ist zu sehen dass es im Ausgangsmodell nur eine Möglichkeit gibt, das Modell zu verbessern, nämlich den Ausschluss von Lesekompetenz (reading) (AIC von 889.88 auf 887.88). Danach gibt es keine Möglichkeit zur Verbesserung mehr, beide verbleibenden Prädiktoren würden bei einem Ausschluss zu einer Verschlechterung des AIC führen (female auf 888.25 und IQ auf 952.08). Damit sind Geschlecht und IQ die Prädiktoren für das optimierte Modell.

An der Ausgabe für das "finale" Modell am Schluss ist zu sehen, dass der Effekt von Geschlecht im finalen Modell hier nicht signifikant ist. Auch oben haben wir gesehen, dass unter Betrachtung des Dekrements dieser Prädiktor wegfallen würde.

Sparsamkeit wird beim AIC im "Strafterm" $2p$ nicht so hoch gewichtet wie bei anderen Informationskriterien. In der Funktion step kann man über die Veränderung des Parameters k steuern, wie streng die Prädiktorauswahl vorgenommen wird. Wenn man hier $k = log(n)$ angibt, wird statt des AIC das sogenannte Bayessche Informationskriterium BIC (Bayesian Information Criterion) verwendet.

$BIC=-2L(\hat{\theta}) + log(n)\cdot p$

Vorsicht, in der Ausgabe der step-Funktion steht immer AIC, auch wenn dies nur mit der Standardeinstellung von $k=2$ tatsächlich dem AIC entspricht!

# Optimierung mit BIC
summary(step(m, k=log(nrow(Schulleistungen))))

Bei der Verwendung des stengeren Kriteriums wird auch Geschlecht aus dem Modell entfernt, es verbleibt nur der IQ im finalen Modell.

Weitere Möglichkeiten

Wie immer gibt es in R viele weitere Wege, zum selben Ziel zu kommen. Eine Vielzahl von Funktionen für die schrittweise Regression bietet z.B. das Paket olsrr, Im Rahmen des Praktikums verwenden wir soweit möglich die Basis-Funktionen von R und beschränken uns daher bei den Aufgaben für schrittweise Analysen auf die step-Funktion.

Das Paket olsrr beinhaltet verschiedene Funktionen, die für die Regressionsanalyse nützlich sind, u.a. auch Funktionen, die die schrittweise Auswahl von Prädiktoren auf Basis verschiedener Kriterien und nach verschiedenen Methoden (vorwärts, rückwärts, etc.) ermöglichen. Finden Sie hier mehr Informationen dazu. Die Funktion ols_step_both_p beinhaltet die Auswahl auf Basis der Signifikanz des Inkrements oder Dekrements und führt in jedem Schritt Tests für Einschluss und Ausschluss durch.
Der Input ist ein Regressionsmodell, das mit der bekannten Funktion lm erstellt wurde. Über die zusätzlichen Argumente kann gesteuert werden, wie streng bei Aufnahme und Ausschluss getestet wird. Über das Argument details können Sie den gesamten Verlauf der schrittweisen Selektion (nicht nur das finale Ergebnis) anzeigen lassen.

# install.packages("olsrr")
library(olsrr)
# pent = p enter, p-Wert zur Aufnahme ins Modell
# prem = p remove, p-Wert zum Ausschluss aus dem Modell
ols_step_both_p(m, pent = .05, prem = .10, details = TRUE)

Aufgaben zur eigenständigen Bearbeitung

1. Datensatz

Der Beispieldatensatz enthält Daten zu Lesekompetenz aus der deutschen Stichprobe der PISA-Erhebung in Deutschland 2009. Gegenüber dem letzten Termin sind in dieser Version noch mehr Variablen enthalten, insbesondere über die Eltern und Besitztümer im Elternhaus.

Im Einzelnen enthält der Datensatz die folgenden Variablen:

• Klassenstufe (Grade)
• Alter in Jahren (Age)
• Geschlecht (Female, 0=m / 1=w)
• Lesekompetenz (Reading)
• Lesefreude (JoyRead)
• Lernzeit in Minuten für Deutsch (LearnMins)
• Sozialstatus (HISEI, "highest international socio-economic index of occupational status")
• Fragebogen-Score für kulturelle Besitztümer zu Hause (CultPoss, z. B. klassische Literatur, Kunstwerke)
• Anzahl Bücher zu Hause (Books, eigentlich eine ordinale Ratingskala)
• Anzahl Fernseher zu Hause (TVs)
• Anzahl Computer zu Hause (Computers)
• Anzahl Autos zu Hause (Cars)
• Migrationshintergrund (MigHintergrund, 0=beide Eltern in D geboren, 1=min. 1 Elternteil im Ausland geboren)
• Bildungsabschluss des Vaters (FatherEdu, International Standard Classification of Education)
• Bildungsabschluss der Mutter (MotherEdu, International Standard Classification of Education)

Laden Sie den Datensatz und verschaffen Sie sich einen Eindruck über die Daten. Installieren Sie außerdem das Paket "lm.beta".

#install.packages("lm.beta")
library(lm.beta) # erforderlich für standardiserte Gewichte
# Datensatz laden 
data("PISA2009", package = "PsyBSc7")

2. Modellkonstruktion mit individuellen Merkmalen der Schüler/innen

Untersucht werden sollen verschiedene Modelle zur Vorhersage der Lesekompetenz anhand der individuellen Schülermerkmale Alter, Klassenstufe, Geschlecht und Migrationshintergrund.

2 a) Alter und Klassenstufe

Berechnen Sie drei separate Regressionsmodelle mit (1) Alter, (2) Klassenstufe und (3) Alter und Klassenstufe als Prädiktoren für Lesekompetenz. Berechnen Sie die Unterschiede in der erklärten Varianz und deren Signifikanz für die Vergleiche der Modelle 1 vs. 3 und 2 vs. 3. Interpretieren Sie die Ergebnisse. Für welches Modell würden Sie sich entscheiden? Bitte begründen Sie Ihre Antwort. Beziehen Sie sich dabei auf alles, was Sie in beiden bisherigen Sitzungen zur Regression gelernt haben.

#Drei Modelle
m1 <- lm(Reading ~ Age, data = PISA2009)
summary(m1)

m2 <- lm(Reading ~ Grade, data = PISA2009)
summary(m2)

m3 <- lm(Reading ~ Age + Grade, data = PISA2009)
summary(m3)

# Vergleich m1 vs. m2
summary(m3)$r.squared - summary(m1)$r.squared # Inkrement
anova(m1, m3)                                 # Modellvergleich

summary(m3)$r.squared - summary(m2)$r.squared # Inkrement
anova(m2, m3)                                 # Modellvergleich

Ergebnisinterpretation 2 a)

• Alter allein erklärt einen Varianzanteil von $R^2=r sprintf("%.3f",summary(m1)$r.squared)$, dieser Effekt ist nicht signifikant.
• Klassenstufe allein erklärt einen Varianzanteil von $R^2=r sprintf("%.2f",summary(m2)$r.squared)$, dieser Effekt ist signifikant.
• Beide Prädiktoren gemeinsam erklären einen signifikant von null verschiedenen Varianzanteil von $R_c^2=r sprintf("%.2f",summary(m3)$r.squared)$.
• Das Inkrement von Alter bei hinzunahme in das Modell mit Klassenstufe ($\Delta R^2=r sprintf("%.2f",summary(m3)$r.squared-summary(m2)$r.squared)$) ist signifikant ($p=r sprintf("%.3f",anova(m2, m3)[[6]][2])$), obwohl Alter für sich genommen kein signifikanter Prädiktor für Lesekompetenz ist.
• Es handelt sich bei Alter um eine Suppressorvariable. Sie ist mit der abhängigen Variablen kaum korreliert ($r_{Age,Reading}=r sprintf("%.2f",cor(PISA2009$Age, PISA2009$Reading))$), korreliert aber hoch mit der Klassenstufe ($r_{Age,Grade}=r sprintf("%.2f",cor(PISA2009$Age, PISA2009$Grade))$).
• Inhaltliche Interpretation der Modelle:
  • Bei konstantem Alter können Schüler/innen in höheren Klassen besser Lesen ($\hat{\beta}{Grade}=r sprintf("%.2f", m3$coefficients[3])$). Dieser Effekt ist etwas stärker als ohne die Kontrolle von Alter im Modell 2 ($\hat{\beta}{Grade}=r sprintf("%.2f", m2$coefficients[2])$).
  • Bei konstanter Klassenstufe zeigen ältere Schüler/innen eine niedrigere Lesekompetenz ($\hat{\beta}{Age}=r sprintf("%.2f", m3$coefficients[2])$). Dies geht vermutlich darauf zurück, dass diese Schüler/innen Klassen wiederholt haben und generell leistungsschwächer sind. Ohne die Kontrolle der Klassenstufe in Modell 1 ist dieser "Effekt" nicht sichtbar ($\hat{\beta}{Age}=r sprintf("%.2f", m1$coefficients[2])$).

# o.g. Korrelationen von Alter mit Lesekompetenz und Klassenstufe
cor.test(PISA2009$Age, PISA2009$Reading)
cor.test(PISA2009$Age, PISA2009$Grade)

2 b) Hinzunahme von Geschlecht in Migrationshintergrund

Ergänzen Sie Modell 3 (mit Alter und Klassenstufe) um die beiden Prädiktoren Geschlecht und Migrationshintergrund. Wie viel zusätzliche Varianz erklären diese beiden Prädiktoren? Interpretieren Sie auch hier die Ergebnisse und begründen, für welches Modell Sie sich entscheiden würden.

m4 <- lm(Reading ~ Age + Grade + Female + MigHintergrund, data = PISA2009)
summary(m4)

# Vergleich m3 vs. m4
summary(m4)$r.squared - summary(m3)$r.squared # Inkrement
anova(m3, m4)                                 # Modellvergleich

Ergebnisinterpretation 2 b)

• Geschlecht und Migrationshintergrund erklären beide gemeinsam einen zusätzlichen Varianzanteil von ($\Delta R^2=r sprintf("%.2f",summary(m4)$r.squared-summary(m3)$r.squared)$). Dieser Zuwachs ist signifikant ($p=r sprintf("%.3f",anova(m3, m4)[[6]][2])$).
• Im Modell mit allen vier Prädiktoren sind die Regressionsgewichte der beiden zusätzlichen Prädiktoren auf dem 5%-Niveau trotzdem nicht signifikant.

3. Schrittweise Analyse mit Eigenschaften des Elternhauses

Es soll untersucht werden, welche Charakteristika der Eltern und Besitztümer des Elternhauses die besten Prädiktoren für die Lesekompetenz der getesteten Schüler/innen sind. Aus den folgenden Variablen soll ein optimales Vorhersagemodell gebildet werden:

• Sozialstatus (HISEI)
• Bildungsabschluss des Vaters (FatherEdu)
• Bildungsabschluss der Mutter (MotherEdu)
• Kulturelle Besitztümer (CultPoss)
• Anzahl Bücher zu Hause (Books)
• Anzahl Fernseher zu Hause (TVs)
• Anzahl Computer zu Hause (Computers)
• Anzahl Autos zu Hause (Cars)

3 a) Führen Sie eine schrittweise Regressionsanalyse mit dem AIC als Auswahlkriterium durch.

Geben Sie zunächst an, welche Prädiktoren im Modell mit allen Prädiktoren signifikant sind. Führen Sie dann eine schrittweise Regressionsanalyse durch. Welche Prädiktoren werden hier ausgewählt?

# Modell mit allen o.g. Prädiktoren
test <- lm(Reading ~ HISEI + FatherEdu + MotherEdu + 
             CultPoss + Books + TVs + Computers + Cars,
           data = PISA2009)

summary(lm.beta(test))

step(test)

m.final <- lm(Reading ~ HISEI + MotherEdu + Books, data = PISA2009)

3 a) Ergebnisinterpretation

• Das Modell mit allen Prädiktoren erklärt einen Varianzanteil von $R^2=r round(summary(test)$r.squared,2)$.
• Auf dem 5%-Niveau ist jedoch nur der Effekt des Sozialstatus signifikant, die meisten Prädiktoren haben keine signifikanten Regressionsgewichte.
• Nach schrittweiser Selektion bleiben nur Sozialstatus, Bildungsabschluss der Mutter und Anzahl Bücher zu Hause im Modell. Alle drei Regressionsgewichte sind auf dem 5%-Niveau signifikant von null verschieden.
• Die drei Prädiktoren erklären gemeinsam einen Varianzanteil von $R^2=r round(summary(m.final)$r.squared,2)$, also nur unerheblich weniger als das Modell mit allen Prädiktoren.

3 b) Ermitteln Sie auch standardiserte Regressionsgewichte für das finale Modell. Wie stark sind die Effekte der ausgewählten Prädiktoren im Vergleich?

# Zusammenfassung einschließlich standardisierter Koeffizienten für das finale Modell
summary(lm.beta(m.final))

3 b) Ergebnisinterpretation

• Gemessen am standardisierten Regressionsgewicht ist der Effekt des Sozialstatus am größten.

3 c) Führen Sie zum Vergleich die schrittweise Regressionsanalyse nun auf Basis von Inkrement und Dekrement (mit dem Package olsrr) durch und vergleichen Sie den Output mit der Lösung aus Aufgabe 3a).

# Modell mit allen o.g. Prädiktoren
test <- lm(Reading ~ HISEI + FatherEdu + MotherEdu + 
             CultPoss + Books + TVs + Computers + Cars,
           data = PISA2009)

ols_step_both_p(test, details = TRUE)

3 c) Ergebnisinterpretation

• Während in Aufgabe 3b) die Selektion mit einem vollen Modell gestartet ist, werden hier Prädiktoren nacheinander hinzugefügt.
• Das final ausgewählte Modell ist in beiden Aufgaben das gleiche.

martscht/PsyBSc7 documentation built on Sept. 1, 2020, 10:50 p.m.