Lineares Modell, Moderationsanalyse (Interaktion)

Unit, Examples, etc

Rmd

Moderationsanalyse lässt sich im Kern reduzieren auf die Integration eines Interaktionsterms im linearen Modell.

Allgemeines

Moderationsanalyse - Interaktion - Wechselwirkung

Im ALM (Allgemeinen Linearen Modell) sind Interaktionsterme möglich.

Moderationsanalyse ist statistisch gesehen eine Interaktion (Wechselwirkung). Moderation ist eine inhaltliche Sicht auf Interaktion. Ein (bestimmter) Prädiktor wird als Moderator bezeichnet.

lm_moderation_concept

Prädiktor und Moderator müssen neben der Interaktion ins Modell aufgenommen werden, damit Interaktionsterm sinnvoll wird.

lm_moderation_statistical_model.png

Interaktionsterme können durch Multiplikation der entsprechenden Variablen in das Modell aufgenommen werden. Dies geschieht in den verschiedenen Paketen teilweise implizit, kann aber auch von Hand erreicht werden (Datenaufbereitung).

Aufnahme eines Interaktionsterms in ein LM führt zu unterschiedlichen Steigungen an verschiedenen Stellen. Bei einem dichotomen oder kategorialen Prädiktor führt dies zu unterschiedlichen Steigungen pro Teilgruppe, bei stetigen Prädiktoren zu einer kontinuierlichen Verschiebung der Steigung. Grafisch wird durch die Interaktion die Vorhersage-Ebene verbogen.

lm_moderation_binary.png

lm_interaction.png

Der Koeffizient eines Interaktionseffekts im LM ist interpretierbar als die Steigung des einen Prädiktors (Moderator), wenn der andere bei 0 fixiert ist. Zentrieren der Prädiktoren kann zu einer besseren Interpretierbarkeit der Koeffizienten führen (Steigung beim Mittelwert des Prädiktors, Mittlerer Effekt über den Gesamtbereich des Prädiktors).

Alle Vorausstzungstests müssen durchgeführt sein, wie üblich.

Einführendes Beispiel GPA SAT

Vergleiche gpa_sat

Szenario, Hintergrund

SAT: Scholastic Assessment Test

In den USA nationaler Zugangstest zu Universitäten

GPA: Grade Point Average

Übersetzung der american grades (A+ bis F) in Zahlenwerte (0-4, höher ist besser). GPA ist schulabhängig.

Hintergrund, Szenario, Overview

When deciding whether to admit an applicant, colleges take lots of factors, such as grades, sports, activities, leadership positions, awards, teacher recommendations, and test scores, into consideration. Using SAT scores as a basis of whether to admit a student or not has created some controversy. Among other things, people question whether the SATs are fair and whether they predict college performance.

This study examines the SAT and GPA information of 105 students who graduated from a state university with a B.S. in computer science. Using the grades and test scores from high school, can you predict a student’s college grades?

Descriptions of Variables

Variable    Description
high_GPA    High school grade point average
math_SAT    Math SAT score
verb_SAT    Verbal SAT score
comp_GPA    Computer science grade point average
univ_GPA    Overall university grade point average

Alle Bewertungen (SAT, GPAs) sind umso besser, je höher der jeweilige Wert ist.

Data

oder lokal

Beispiel GPA-SAT: Interaktion zweier stetiger Variablen: Vorhersage des universitären Gesamterfolgs aus SAT-Testwerten (Mathematik, Sprache)

Das Beispielmodell als Formel:

\[ y_i = \beta_0 + \beta_1x_{1i} + \beta_2x_{2i} + \beta_3x_{1i}x_{2i} + \epsilon_i \]

df <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/sat.txt")
head(df)

##   high_GPA math_SAT verb_SAT comp_GPA univ_GPA
## 1     3.45      643      589     3.76     3.52
## 2     2.78      558      512     2.87     2.91
## 3     2.52      583      503     2.54     2.40
## 4     3.67      685      602     3.83     3.47
## 5     3.24      592      538     3.29     3.47
## 6     2.10      562      486     2.64     2.37

# we test mathematics SAT only
r.m <- lm(univ_GPA ~ math_SAT, data=df)
summary(r.m)

## 
## Call:
## lm(formula = univ_GPA ~ math_SAT, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.78921 -0.23669  0.01962  0.23730  0.82748 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.2622924  0.3838157  -0.683    0.496    
## math_SAT     0.0055132  0.0006137   8.983 1.34e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3365 on 103 degrees of freedom
## Multiple R-squared:  0.4393, Adjusted R-squared:  0.4338 
## F-statistic: 80.69 on 1 and 103 DF,  p-value: 1.344e-14

# we check residuals
require(psych)

## Lade nötiges Paket: psych

psych::describe(r.m$residuals)

##    vars   n mean   sd median trimmed  mad   min  max range  skew kurtosis   se
## X1    1 105    0 0.33   0.02       0 0.33 -0.79 0.83  1.62 -0.18    -0.56 0.03

# now we control predictive power of verbal SAT only
r.v <- lm(univ_GPA ~ verb_SAT, data=df)
summary(r.v)

## 
## Call:
## lm(formula = univ_GPA ~ verb_SAT, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.94856 -0.19745  0.03271  0.20373  0.79492 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.4080980  0.3199775   1.275    0.205    
## verb_SAT    0.0046187  0.0005316   8.688 6.03e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3414 on 103 degrees of freedom
## Multiple R-squared:  0.4229, Adjusted R-squared:  0.4173 
## F-statistic: 75.48 on 1 and 103 DF,  p-value: 6.034e-14

# we combine mathematical and verbal SAT additively
r.mv <- lm(univ_GPA ~ math_SAT + verb_SAT, data=df)
summary(r.mv)

## 
## Call:
## lm(formula = univ_GPA ~ math_SAT + verb_SAT, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7808 -0.1844  0.0568  0.2054  0.7499 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -0.2375336  0.3750378  -0.633  0.52792   
## math_SAT     0.0032909  0.0010902   3.019  0.00321 **
## verb_SAT     0.0022718  0.0009308   2.441  0.01638 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3287 on 102 degrees of freedom
## Multiple R-squared:  0.4702, Adjusted R-squared:  0.4598 
## F-statistic: 45.27 on 2 and 102 DF,  p-value: 8.488e-15

# we add interaction
r.mv.i <- lm(univ_GPA ~ math_SAT + verb_SAT + math_SAT * verb_SAT, data=df)
# or shorter
r.mv.i <- lm(univ_GPA ~ math_SAT * verb_SAT, data=df)
summary(r.mv.i)

## 
## Call:
## lm(formula = univ_GPA ~ math_SAT * verb_SAT, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.85523 -0.18050  0.03869  0.18225  0.93201 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)   
## (Intercept)       -1.088e+01  4.015e+00  -2.711  0.00788 **
## math_SAT           2.030e-02  6.475e-03   3.135  0.00225 **
## verb_SAT           2.002e-02  6.725e-03   2.976  0.00365 **
## math_SAT:verb_SAT -2.814e-05  1.057e-05  -2.663  0.00902 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3193 on 101 degrees of freedom
## Multiple R-squared:  0.505,  Adjusted R-squared:  0.4903 
## F-statistic: 34.34 on 3 and 101 DF,  p-value: 2.196e-15

Grafische Darstellung als 3D-Plot

Mit der library(rockchalk)

df <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/sat.txt")
m.i <- lm(univ_GPA ~ math_SAT * verb_SAT, data=df)

require(rockchalk)

## Lade nötiges Paket: rockchalk

plotPlane(model = m.i, plotx1 = "math_SAT", plotx2 = "verb_SAT")

# or to emphasize the view of the interaction
plotPlane(model = m.i, plotx1 = "math_SAT", plotx2 = "verb_SAT", theta = 0, phi = 10)

plotPlane(model = m.i, plotx1 = "math_SAT", plotx2 = "verb_SAT", drawArrows = T)

Mit rgl() können wir dynamische dreidimensionale Grafiken machen, die sich interaktiv drehen lassen. Der Code unten ist nur Beispielcode, der in einem RMD nicht gerendert werden kann und deshalb ist eval=FALSE gesetzt.

df <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/sat.txt")
m.i <- lm(univ_GPA ~ math_SAT * verb_SAT, data=df)

require(rgl)
open3d()
plot3d(x=df$math_SAT, y=df$verb_SAT, z=df$univ_GPA, col="red", size=4)

grd <- expand.grid(math_SAT = sort(unique(df$math_SAT)), verb_SAT = sort(unique(df$verb_SAT)) )
grd$pred <-predict(m.i, newdata=grd)
persp3d(x=unique(grd[[1]]), y=unique(grd[[2]]), 
              z=matrix(grd[[3]],length(unique(grd[[1]])),length(unique(grd[[2]]))), add=TRUE)

Modellvergleiche

Modellvergleiche hierarchischer Modelle funktionieren wie gewohnt

df <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/sat.txt")
m.0 <- lm(univ_GPA ~ 1, data=df)
m.a <- lm(univ_GPA ~ math_SAT + verb_SAT, data=df)
m.i <- lm(univ_GPA ~ math_SAT * verb_SAT, data=df)

summary(m.0)

## 
## Call:
## lm(formula = univ_GPA ~ 1, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.0929 -0.1629  0.1171  0.2971  0.6371 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.17286    0.04364    72.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4472 on 104 degrees of freedom

summary(m.a)

## 
## Call:
## lm(formula = univ_GPA ~ math_SAT + verb_SAT, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7808 -0.1844  0.0568  0.2054  0.7499 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -0.2375336  0.3750378  -0.633  0.52792   
## math_SAT     0.0032909  0.0010902   3.019  0.00321 **
## verb_SAT     0.0022718  0.0009308   2.441  0.01638 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3287 on 102 degrees of freedom
## Multiple R-squared:  0.4702, Adjusted R-squared:  0.4598 
## F-statistic: 45.27 on 2 and 102 DF,  p-value: 8.488e-15

summary(m.i)

## 
## Call:
## lm(formula = univ_GPA ~ math_SAT * verb_SAT, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.85523 -0.18050  0.03869  0.18225  0.93201 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)   
## (Intercept)       -1.088e+01  4.015e+00  -2.711  0.00788 **
## math_SAT           2.030e-02  6.475e-03   3.135  0.00225 **
## verb_SAT           2.002e-02  6.725e-03   2.976  0.00365 **
## math_SAT:verb_SAT -2.814e-05  1.057e-05  -2.663  0.00902 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3193 on 101 degrees of freedom
## Multiple R-squared:  0.505,  Adjusted R-squared:  0.4903 
## F-statistic: 34.34 on 3 and 101 DF,  p-value: 2.196e-15

anova(m.0, m.a)

## Analysis of Variance Table
## 
## Model 1: univ_GPA ~ 1
## Model 2: univ_GPA ~ math_SAT + verb_SAT
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    104 20.798                                  
## 2    102 11.018  2    9.7797 45.267 8.488e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(m.0, m.i)

## Analysis of Variance Table
## 
## Model 1: univ_GPA ~ 1
## Model 2: univ_GPA ~ math_SAT * verb_SAT
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    104 20.798                                  
## 2    101 10.296  3    10.502 34.343 2.196e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# does the interaction term result in a better model?
anova(m.a, m.i)

## Analysis of Variance Table
## 
## Model 1: univ_GPA ~ math_SAT + verb_SAT
## Model 2: univ_GPA ~ math_SAT * verb_SAT
##   Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
## 1    102 11.018                                
## 2    101 10.296  1   0.72281 7.0908 0.009018 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# ... yes!

Vorhersage von Werten

Hier dargestellt an dem einfachen Fall von zwei Vorhersagevariablen, so dass die Verhältnisse noch in einem dreidimensionalen Raum vorstellbar sind. Wenn ein Interaktionsterm mit im Spiel ist, ist die Vorhersage abhängig von der Kombination der Vorhersagevariablen. Halten wir Vorhersagevariable fest und lassen die andere laufen, bekommen wir eine Regressionsgerade mit einer bestimmten Steigung. Halten wir diese Vorhersagevariablen auf einem anderen Wert fest und lassen die andere laufen, bekommen wir eine Regressionsgerade mit einer anderen Steigung.

df <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/sat.txt")
m.i <- lm(univ_GPA ~ math_SAT * verb_SAT, data=df)

# take a look at the relevant dat columns
head(cbind(df$univ_GPA, df$math_SAT, df$verb_SAT))

##      [,1] [,2] [,3]
## [1,] 3.52  643  589
## [2,] 2.91  558  512
## [3,] 2.40  583  503
## [4,] 3.47  685  602
## [5,] 3.47  592  538
## [6,] 2.37  562  486

# we predict subj 4 'manually'
m.i$coefficients['(Intercept)'] + 
  df$math_SAT[4] * m.i$coefficients['math_SAT'] +
  df$verb_SAT[4] * m.i$coefficients['verb_SAT'] +
  df$math_SAT[4]  * df$verb_SAT[4] * m.i$coefficients['math_SAT:verb_SAT']

## (Intercept) 
##    3.468149

# head(predict.lm(m.i))
# head(residuals.lm(m.i))
head(cbind(df$univ_GPA, predict.lm(m.i), residuals.lm(m.i)))

##   [,1]     [,2]         [,3]
## 1 3.52 3.302059  0.217941319
## 2 2.91 2.653214  0.256786385
## 3 2.40 2.768033 -0.368032764
## 4 3.47 3.468149  0.001851412
## 5 3.47 2.940844  0.529155948
## 6 2.37 2.567547 -0.197547273

Beispiel GPA-SAT: Interaktion einer stetigen mit einer dichotomen Variable: Vorhersage des universitären Gesamterfolgs aus SAT-Testwerten (Mathematik, Sprache)

Nur zur Veranschaulichung. Es wird mit einer künstlich erzeugtendichotomen Variable gender gerechnet. Hauptsächlich geht es um die Visualisierung.

df <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/sat.txt")

df$gender <- 'f'
df$gender[df$verb_SAT <= mean(df$verb_SAT)] <- 'm'
df$gender <- as.factor(df$gender)

m.i <- lm(univ_GPA ~ math_SAT * gender, data=df)
summary(m.i)

## 
## Call:
## lm(formula = univ_GPA ~ math_SAT * gender, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9232 -0.2068  0.0279  0.1980  0.9541 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)  
## (Intercept)       1.463092   0.782132   1.871   0.0643 .
## math_SAT          0.002962   0.001175   2.520   0.0133 *
## genderm          -2.594627   1.146410  -2.263   0.0258 *
## math_SAT:genderm  0.003990   0.001853   2.153   0.0337 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3293 on 101 degrees of freedom
## Multiple R-squared:  0.4735, Adjusted R-squared:  0.4579 
## F-statistic: 30.28 on 3 and 101 DF,  p-value: 4.789e-14

Visualisierung Gender

df <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/sat.txt")
df$gender <- 'f'
df$gender[df$verb_SAT <= mean(df$verb_SAT)] <- 'm'
df$gender <- as.factor(df$gender)
m.a <- lm(univ_GPA ~ math_SAT + gender, data=df)
m.i <- lm(univ_GPA ~ math_SAT * gender, data=df)

require("ggplot2")

## Lade nötiges Paket: ggplot2

## 
## Attache Paket: 'ggplot2'

## Die folgenden Objekte sind maskiert von 'package:psych':
## 
##     %+%, alpha

# raw scatterplot 
ggplot(df, aes(x=math_SAT, y=univ_GPA)) + geom_point() 

# graphical look at gender specifity. Gender must have type factor to be used adequately. This is given already
str(df$gender)

##  Factor w/ 2 levels "f","m": 2 2 2 1 2 2 2 2 2 2 ...

# plot values in different color and shape 
# ggplot(df, aes(x=math_SAT, y=univ_GPA, color=gender, shape=gender)) + geom_point() 
pp <- ggplot(df, aes(x=math_SAT, y=univ_GPA))
# include global regression line
pp + geom_point(aes(color=gender, shape=gender)) + stat_smooth(method = "lm", se=FALSE)

## `geom_smooth()` using formula 'y ~ x'

# include global regression line with confidence interval
pp + geom_point(aes(color=gender, shape=gender)) + stat_smooth(method = "lm")

## `geom_smooth()` using formula 'y ~ x'

# include group specific regression lines without differences in slope (linear model without interaction)
# for this we need the coefficients of the model to plot
cbPalette = c('red', 'turquoise')
pp + geom_point(aes(color=gender, shape=gender)) +
  geom_abline(intercept = m.a$coefficients['(Intercept)'],                               slope =  m.a$coefficients['math_SAT'], color=cbPalette[1]) + 
  geom_abline(intercept = m.a$coefficients['(Intercept)'] + m.a$coefficients['genderm'], slope =  m.a$coefficients['math_SAT'], color=cbPalette[2]) 

# add scatterplot and gender specific regression lines as estimated by model with interaction m.i
pp + geom_point(aes(color=gender, shape=gender)) + stat_smooth(aes(color=gender, shape=gender), method = "lm", se=FALSE)

## Warning: Ignoring unknown aesthetics: shape

## `geom_smooth()` using formula 'y ~ x'

# gender specific regression lines as estimated by model with interaction m.i including confidence intervals and global regression line
pp + geom_point(aes(color=gender, shape=gender)) + stat_smooth(method = "lm", se=FALSE) + stat_smooth(aes(color=gender, shape=gender), method = "lm")

## Warning: Ignoring unknown aesthetics: shape

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using formula 'y ~ x'

Ausflug: Zweidimensionale Graphische Annäherung durch farbliche Schichten

Dargestellt im Beispiel body_comparisonl

Ausflug: Äquivalenz von Interaktion mit neuer erklärender Variable, die das Produkt der verwendeten erklärenden Variablen ist

df <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/sat.txt")

# we combine mathematical and verbal SAT additively
m.mv <- lm(univ_GPA ~ math_SAT + verb_SAT, data=df)
# summary(m.mv)

# we add interaction
m.mv.i <- lm(univ_GPA ~ math_SAT + verb_SAT + math_SAT * verb_SAT, data=df)
# or shorter, as it results in the same model
m.mv.i <- lm(univ_GPA ~ math_SAT * verb_SAT, data=df)
#summary(m.mv.i)

# we create a new variable math_verb_SAT
df$math_verb_SAT <- df$math_SAT * df$verb_SAT
# an additive model of math_SAT, verb_SAT, and math_verb_SAT ist equivalent to m.mv.i
m.mv.i2 <- lm(univ_GPA ~ math_SAT + verb_SAT + math_verb_SAT, data=df)

# both summaries are the same
summary(m.mv.i)

## 
## Call:
## lm(formula = univ_GPA ~ math_SAT * verb_SAT, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.85523 -0.18050  0.03869  0.18225  0.93201 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)   
## (Intercept)       -1.088e+01  4.015e+00  -2.711  0.00788 **
## math_SAT           2.030e-02  6.475e-03   3.135  0.00225 **
## verb_SAT           2.002e-02  6.725e-03   2.976  0.00365 **
## math_SAT:verb_SAT -2.814e-05  1.057e-05  -2.663  0.00902 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3193 on 101 degrees of freedom
## Multiple R-squared:  0.505,  Adjusted R-squared:  0.4903 
## F-statistic: 34.34 on 3 and 101 DF,  p-value: 2.196e-15

summary(m.mv.i2)

## 
## Call:
## lm(formula = univ_GPA ~ math_SAT + verb_SAT + math_verb_SAT, 
##     data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.85523 -0.18050  0.03869  0.18225  0.93201 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   -1.088e+01  4.015e+00  -2.711  0.00788 **
## math_SAT       2.030e-02  6.475e-03   3.135  0.00225 **
## verb_SAT       2.002e-02  6.725e-03   2.976  0.00365 **
## math_verb_SAT -2.814e-05  1.057e-05  -2.663  0.00902 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3193 on 101 degrees of freedom
## Multiple R-squared:  0.505,  Adjusted R-squared:  0.4903 
## F-statistic: 34.34 on 3 and 101 DF,  p-value: 2.196e-15

df <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/sat.txt")
head(df)

##   high_GPA math_SAT verb_SAT comp_GPA univ_GPA
## 1     3.45      643      589     3.76     3.52
## 2     2.78      558      512     2.87     2.91
## 3     2.52      583      503     2.54     2.40
## 4     3.67      685      602     3.83     3.47
## 5     3.24      592      538     3.29     3.47
## 6     2.10      562      486     2.64     2.37

# we test mathematics SAT only
r.m <- lm(univ_GPA ~ math_SAT, data=df)
summary(r.m)

## 
## Call:
## lm(formula = univ_GPA ~ math_SAT, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.78921 -0.23669  0.01962  0.23730  0.82748 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.2622924  0.3838157  -0.683    0.496    
## math_SAT     0.0055132  0.0006137   8.983 1.34e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3365 on 103 degrees of freedom
## Multiple R-squared:  0.4393, Adjusted R-squared:  0.4338 
## F-statistic: 80.69 on 1 and 103 DF,  p-value: 1.344e-14

# we check residuals
require(psych)
psych::describe(r.m$residuals)

##    vars   n mean   sd median trimmed  mad   min  max range  skew kurtosis   se
## X1    1 105    0 0.33   0.02       0 0.33 -0.79 0.83  1.62 -0.18    -0.56 0.03

# now we control predictive power of verbal SAT only
r.v <- lm(univ_GPA ~ verb_SAT, data=df)
summary(r.v)

## 
## Call:
## lm(formula = univ_GPA ~ verb_SAT, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.94856 -0.19745  0.03271  0.20373  0.79492 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.4080980  0.3199775   1.275    0.205    
## verb_SAT    0.0046187  0.0005316   8.688 6.03e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3414 on 103 degrees of freedom
## Multiple R-squared:  0.4229, Adjusted R-squared:  0.4173 
## F-statistic: 75.48 on 1 and 103 DF,  p-value: 6.034e-14

# we combine mathematical and verbal SAT additively
r.mv <- lm(univ_GPA ~ math_SAT + verb_SAT, data=df)
summary(r.mv)

## 
## Call:
## lm(formula = univ_GPA ~ math_SAT + verb_SAT, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7808 -0.1844  0.0568  0.2054  0.7499 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -0.2375336  0.3750378  -0.633  0.52792   
## math_SAT     0.0032909  0.0010902   3.019  0.00321 **
## verb_SAT     0.0022718  0.0009308   2.441  0.01638 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3287 on 102 degrees of freedom
## Multiple R-squared:  0.4702, Adjusted R-squared:  0.4598 
## F-statistic: 45.27 on 2 and 102 DF,  p-value: 8.488e-15

# we add interaction
r.mv.i <- lm(univ_GPA ~ math_SAT + verb_SAT + math_SAT * verb_SAT, data=df)
# or shorter
r.mv.i <- lm(univ_GPA ~ math_SAT * verb_SAT, data=df)
summary(r.mv.i)

## 
## Call:
## lm(formula = univ_GPA ~ math_SAT * verb_SAT, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.85523 -0.18050  0.03869  0.18225  0.93201 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)   
## (Intercept)       -1.088e+01  4.015e+00  -2.711  0.00788 **
## math_SAT           2.030e-02  6.475e-03   3.135  0.00225 **
## verb_SAT           2.002e-02  6.725e-03   2.976  0.00365 **
## math_SAT:verb_SAT -2.814e-05  1.057e-05  -2.663  0.00902 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3193 on 101 degrees of freedom
## Multiple R-squared:  0.505,  Adjusted R-squared:  0.4903 
## F-statistic: 34.34 on 3 and 101 DF,  p-value: 2.196e-15

Shiny

Kann nur nach entsprechender Installation durchgeführt werden. multReg muss installiert sein. Hier nur zu Demonstrations- und Visualisierungszwecken.

require("shiny")

## Lade nötiges Paket: shiny

oldwd <- getwd()
# setwd("/Users/pzezula/lehre/lehre_ss_2017/shiny")
# a dynamic visualization, needs shiny and shinyRGL to be installed, only on PC PZ for demonstration purposes
##runApp("multReg")
## show it on shiny server (Andreas Cordes)
## https://elder.uni-psych.gwdg.de:3838/acordes/multReg/
setwd(oldwd)

26. B. Visualisieren des Beispiels openness_happyness_affect (vgl. Beispiele unten)

Beispiele

Aggression, Computer-Spielzeit und mangelndes Einfühlungsvermögen

Quelle: @field2012discover

html

Beispiel Big-Five Offenheit, prospektive Lebenszufriedenheit und negativer Affekt

openness_happyness_affect html

Beispiel Autowasch Daten: Beispiel aus Everitt

car_wash

Beispiel Körpervergleich

Hier gibt es einen längeren Teil zu Moderationsanalysen bzw. Interaktionen

html

Links

–>

Erklärung von Interaktionstermen an einem Beispiel von Pflanzenwachstum, Sonne und Bakterien. [https://www.theanalysisfactor.com/interpreting-interactions-in-regression/] [https://www.theanalysisfactor.com/clarifications-on-interpreting-interactions-in-regression/]

Referenzen

Field (2012, p. 312) Kapitel 7: Regression

Check

• Interaktion ist die Aufnahme einer Variable ins Modell, die durch die Multiplikation von erklärenden Variablen miteinander gebildet wird.

• Bei einem Interaktionseffekt hängt die Vorhersage der Reaktionsvariable einer Beobachtung durch eine erklärende Variable vom Wert der anderen erklärenden Variable ab

• Grafisch führt eine Interaktion zu einer Verwindung der Vorhersageebene.

• Mit 3d Abbildungen können wir diese Verwindung bei zwei erklärenden Variablen visualisieren, beispielsweise mit library(rockchalk). ggplot() konzentriert sich auf 2d Abbildungen und ist weniger geeignet.

• Für jede einzelne Stelle einer erklärenden Variable bleibt die Vorhersage aber eine Linie

• Bei mehr als 2 erklärenden Variablen bewegen wir uns im vieldimensionalen Raum, das ist etwas schwerer vorstellbar

• vorstellbar als MLM nur mit fixed, also ohne random Effects

dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/sat.txt")
# we store the multiplication of two explanatary variables in a new variable
dd$iSAT = dd$math_SAT * dd$verb_SAT
m.mv <- lm(univ_GPA ~ math_SAT + verb_SAT + math_SAT:verb_SAT, data=dd)
m.mv$coefficients == lm(univ_GPA ~ math_SAT * verb_SAT, data=dd)$coefficients

##       (Intercept)          math_SAT          verb_SAT math_SAT:verb_SAT 
##              TRUE              TRUE              TRUE              TRUE

m.mv$coefficients == lm(univ_GPA ~ math_SAT + verb_SAT + iSAT, data=dd)$coefficients

##       (Intercept)          math_SAT          verb_SAT math_SAT:verb_SAT 
##              TRUE              TRUE              TRUE              TRUE

m.mv$coefficients

##       (Intercept)          math_SAT          verb_SAT math_SAT:verb_SAT 
##     -1.088382e+01      2.030158e-02      2.001761e-02     -2.814271e-05

head(cbind(dd$univ_GPA, model.matrix(m.mv), predict(m.mv), m.mv$residuals))

##        (Intercept) math_SAT verb_SAT math_SAT:verb_SAT                      
## 1 3.52           1      643      589            378727 3.302059  0.217941319
## 2 2.91           1      558      512            285696 2.653214  0.256786385
## 3 2.40           1      583      503            293249 2.768033 -0.368032764
## 4 3.47           1      685      602            412370 3.468149  0.001851412
## 5 3.47           1      592      538            318496 2.940844  0.529155948
## 6 2.37           1      562      486            273132 2.567547 -0.197547273

# visualisation of the twisted prediction plane
rockchalk::plotPlane(model = m.mv, plotx1 = "math_SAT", plotx2 = "verb_SAT", theta = 60, phi = 20)

# estimate reaction variable univ_GPA values of new observations
predict(m.mv, data.frame(math_SAT = c(580, 630, 680), verb_SAT = c(560, 590, 710)))

##        1        2        3 
## 2.960201 3.255915 3.546451

Screencast

• Multiple Regression - check: moderation, interaction StudIP - ownCloud

Version: 18 Mai, 2022 09:02