knitr::opts_chunk$set(message=FALSE, warning=FALSE)
library(lme4); library(texreg); library(sjstats); library(sjPlot)
## Loading required package: Matrix
## Version: 1.37.5
## Date: 2020-06-17
## Author: Philip Leifeld (University of Essex)
##
## Consider submitting praise using the praise or praise_interactive functions.
## Please cite the JSS article in your publications -- see citation("texreg").
load('Europe2.RData')
load('imm23.RData')
• Vary the effect of parents’ education among the schools. Check if the
random slope model is better.
# null model
m0 = lmer(MATH ~ 1|SCHID, data = data)
summary(m0)
## Linear mixed model fit by REML ['lmerMod']
## Formula: MATH ~ 1 | SCHID
## Data: data
##
## REML criterion at convergence: 3798.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.68160 -0.72864 -0.01926 0.73172 2.67329
##
## Random effects:
## Groups Name Variance Std.Dev.
## SCHID (Intercept) 26.12 5.111
## Residual 81.24 9.014
## Number of obs: 519, groups: SCHID, 23
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 50.759 1.151 44.09
m1 = lmer(MATH~SES+WHITE+HOMEWORK+PARENTED+SEX+(1|SCHID), data = data)
m2 = lmer(MATH~SES+WHITE+HOMEWORK+PARENTED+SEX+(1+PARENTED|SCHID), data = data)
anova(m1, m2)
## Data: data
## Models:
## m1: MATH ~ SES + WHITE + HOMEWORK + PARENTED + SEX + (1 | SCHID)
## m2: MATH ~ SES + WHITE + HOMEWORK + PARENTED + SEX + (1 + PARENTED |
## m2: SCHID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1 8 3682.5 3716.5 -1833.2 3666.5
## m2 10 3680.3 3722.8 -1830.2 3660.3 6.1607 2 0.04594 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Random slope model for PARENTED is better - p-value < 0.5
• Check if the type of school affects the results in math.
m3 = lmer(MATH~SES+WHITE+HOMEWORK+PARENTED+SEX+PUBLIC+(1+PARENTED|SCHID), data = data)
screenreg(m3)
##
## =============================================
## Model 1
## ---------------------------------------------
## (Intercept) 37.02 ***
## (2.87)
## SES -0.22
## (0.98)
## WHITEwhite 3.37 **
## (1.09)
## HOMEWORK 2.14 ***
## (0.26)
## PARENTED 2.53 ***
## (0.61)
## SEXfemale -0.37
## (0.72)
## PUBLICpublic -0.37
## (1.56)
## ---------------------------------------------
## AIC 3674.81
## BIC 3721.58
## Log Likelihood -1826.40
## Num. obs. 519
## Num. groups: SCHID 23
## Var: SCHID (Intercept) 37.33
## Var: SCHID PARENTED 1.04
## Cov: SCHID (Intercept) PARENTED -5.83
## Var: Residual 63.52
## =============================================
## *** p < 0.001; ** p < 0.01; * p < 0.05
Type of school has no significant effect on math.
• Check if the effect of parental education differs depending on the
school type.
m4 = lmer(MATH~SES+WHITE+HOMEWORK+PARENTED+SEX+MEANSES+PARENTED*PUBLIC+(1+PARENTED|SCHID), data = data)
screenreg(m4)
##
## =============================================
## Model 1
## ---------------------------------------------
## (Intercept) 40.38 ***
## (3.54)
## SES -0.34
## (0.98)
## WHITEwhite 3.04 **
## (1.09)
## HOMEWORK 2.13 ***
## (0.26)
## PARENTED 1.56 *
## (0.76)
## SEXfemale -0.49
## (0.73)
## MEANSES 2.62
## (1.93)
## PUBLICpublic -3.96
## (3.60)
## PARENTED:PUBLICpublic 1.39
## (0.77)
## ---------------------------------------------
## AIC 3670.50
## BIC 3725.78
## Log Likelihood -1822.25
## Num. obs. 519
## Num. groups: SCHID 23
## Var: SCHID (Intercept) 26.40
## Var: SCHID PARENTED 0.63
## Cov: SCHID (Intercept) PARENTED -3.71
## Var: Residual 63.84
## =============================================
## *** p < 0.001; ** p < 0.01; * p < 0.05
Interaction effect of parental education and school type is not significant.
• Indicate which school has the lowest impact of parents’ education.
# random effects
re <- ranef(m2)
# add beta from model to HOMEWORK to obtain true beta-coef per school
re$SCHID$PARENTED = 2.53+re$SCHID$PARENTED
re$SCHID[which.min(abs(re$SCHID$PARENTED - 0)),]
## (Intercept) PARENTED
## 62821 13.9222 0.4143569
School with ID 62821 has the lowest impact of parent’s education
• Plot the effects of your final model.
# fixed effects
plot_model(m2, type = 'est')
# random effects
plot_model(m2, type = 're', sort.est = 'PARENTED')
• Plot the interaction term and comment on the results
plot_model(m4, type = "pred", terms = c("PARENTED", "PUBLIC"))
Interaction effect is not significant nevertheless seems that for public schools parental education have greater slope and increase in parental education can cause greater increase in math while intercept is slightly lower.
• Calculate R2.
performance::r2(m2)
## # R2 for Mixed Models
##
## Conditional R2: 0.415
## Marginal R2: 0.306