Read the data set from Lecture 4.
dat1 <- read.table("lecture4ex1.csv", sep=",", header=TRUE)Transform the fish length variable by 10 to make the multilevel analysis easier to converge. Next, the fish length and fish color are group mean centered.
dat1$lengthx <- dat1$length/10
dat1$avelengthx <- ave(dat1$lengthx, dat1$groupid)
dat1$difflengthx <- dat1$lengthx - dat1$avelengthx
dat1$avegold <- ave(dat1$goldcolor, dat1$groupid)
dat1$diffgold <- dat1$goldcolor - dat1$avegoldLet’s start using build-up strategy to examine the relationship between predictors and criterion. Start with the null model (Model 0).
library(lme4)
library(r2mlm)
out1b0 <- lmer(consume ~ 1 + (1|groupid), data=dat1, REML=FALSE)
r2mlm(out1b0, bargraph=FALSE)## $Decompositions
## total within between
## fixed, within 0.0000000 0 NA
## fixed, between 0.0000000 NA 0
## slope variation 0.0000000 0 NA
## mean variation 0.4939929 NA 1
## sigma2 0.5060071 1 NA
##
## $R2s
## total within between
## f1 0.0000000 0 NA
## f2 0.0000000 NA 0
## v 0.0000000 0 NA
## m 0.4939929 NA 1
## f 0.0000000 NA NA
## fv 0.0000000 0 NA
## fvm 0.4939929 NA NANext, add level-1 predictors: difflengthx (Model 1) and
diffgold (Model 2).
out1b1 <- lmer(consume ~ 1 + difflengthx + (1|groupid), data=dat1, REML=FALSE)
anova(out1b0, out1b1)## Data: dat1
## Models:
## out1b0: consume ~ 1 + (1 | groupid)
## out1b1: consume ~ 1 + difflengthx + (1 | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1b0 3 5120.1 5134.2 -2557.0 5114.1
## out1b1 4 4169.1 4187.8 -2080.5 4161.1 953.01 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2mlm(out1b1, bargraph=FALSE)## $Decompositions
## total within between
## fixed, within 0.3196216 0.7128107 NA
## fixed, between 0.0000000 NA 0
## slope variation 0.0000000 0.0000000 NA
## mean variation 0.5516039 NA 1
## sigma2 0.1287746 0.2871893 NA
##
## $R2s
## total within between
## f1 0.3196216 0.7128107 NA
## f2 0.0000000 NA 0
## v 0.0000000 0.0000000 NA
## m 0.5516039 NA 1
## f 0.3196216 NA NA
## fv 0.3196216 0.7128107 NA
## fvm 0.8712254 NA NAout1b2 <- lmer(consume ~ 1 + diffgold + (1|groupid), data=dat1, REML=FALSE)
anova(out1b0, out1b2)## Data: dat1
## Models:
## out1b0: consume ~ 1 + (1 | groupid)
## out1b2: consume ~ 1 + diffgold + (1 | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1b0 3 5120.1 5134.2 -2557.0 5114.1
## out1b2 4 5118.7 5137.5 -2555.3 5110.7 3.4003 1 0.06518 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2mlm(out1b2, bargraph=FALSE)## $Decompositions
## total within between
## fixed, within 0.002118835 0.004190369 NA
## fixed, between 0.000000000 NA 0
## slope variation 0.000000000 0.000000000 NA
## mean variation 0.494356033 NA 1
## sigma2 0.503525132 0.995809631 NA
##
## $R2s
## total within between
## f1 0.002118835 0.004190369 NA
## f2 0.000000000 NA 0
## v 0.000000000 0.000000000 NA
## m 0.494356033 NA 1
## f 0.002118835 NA NA
## fv 0.002118835 0.004190369 NA
## fvm 0.496474868 NA NAOnly fish length was significant. Thus Model 1 is retained. The interaction effect between level-1 predictors is investigated. The full model is the model with interaction (Model 4) and the restricted model is the model with both level-1 predictors without interaction.
out1b3 <- lmer(consume ~ 1 + difflengthx + diffgold + (1|groupid), data=dat1, REML=FALSE)
out1b4 <- lmer(consume ~ 1 + difflengthx + diffgold + difflengthx:diffgold + (1|groupid), data=dat1, REML=FALSE)
anova(out1b3, out1b4)## Data: dat1
## Models:
## out1b3: consume ~ 1 + difflengthx + diffgold + (1 | groupid)
## out1b4: consume ~ 1 + difflengthx + diffgold + difflengthx:diffgold + (1 | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1b3 5 4165.3 4188.8 -2077.7 4155.3
## out1b4 6 4166.6 4194.7 -2077.3 4154.6 0.7832 1 0.3762The interaction was not significant. Thus, Model 1 is still retained.
Next, check the effect of level-2 predictors: volume (Model
5), nfish (Model 6), avelengthx (Model 7),
avegold (Model 8).
out1b5 <- lmer(consume ~ 1 + difflengthx + volume + (1|groupid), data=dat1, REML=FALSE)
anova(out1b1, out1b5)## Data: dat1
## Models:
## out1b1: consume ~ 1 + difflengthx + (1 | groupid)
## out1b5: consume ~ 1 + difflengthx + volume + (1 | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1b1 4 4169.1 4187.8 -2080.5 4161.1
## out1b5 5 4145.3 4168.8 -2067.7 4135.3 25.755 1 3.877e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1out1b6 <- lmer(consume ~ 1 + difflengthx + nfish + (1|groupid), data=dat1, REML=FALSE)
anova(out1b1, out1b6)## Data: dat1
## Models:
## out1b1: consume ~ 1 + difflengthx + (1 | groupid)
## out1b6: consume ~ 1 + difflengthx + nfish + (1 | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1b1 4 4169.1 4187.8 -2080.5 4161.1
## out1b6 5 4170.2 4193.7 -2080.1 4160.2 0.8566 1 0.3547out1b7 <- lmer(consume ~ 1 + difflengthx + avelengthx + (1|groupid), data=dat1, REML=FALSE)
anova(out1b1, out1b7)## Data: dat1
## Models:
## out1b1: consume ~ 1 + difflengthx + (1 | groupid)
## out1b7: consume ~ 1 + difflengthx + avelengthx + (1 | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1b1 4 4169.1 4187.8 -2080.5 4161.1
## out1b7 5 4127.3 4150.8 -2058.7 4117.3 43.739 1 3.752e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1out1b8 <- lmer(consume ~ 1 + difflengthx + avegold + (1|groupid), data=dat1, REML=FALSE)
anova(out1b1, out1b8)## Data: dat1
## Models:
## out1b1: consume ~ 1 + difflengthx + (1 | groupid)
## out1b8: consume ~ 1 + difflengthx + avegold + (1 | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1b1 4 4169.1 4187.8 -2080.5 4161.1
## out1b8 5 4170.0 4193.5 -2080.0 4160.0 1.0117 1 0.3145The effects of avelengthx and volume were
significant. Model 9 is then created to include only both significant
effects.
out1b9 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + (1|groupid), data=dat1, REML=FALSE)Next, the interactions between level 2 predictors are investigated. Four level-2 predictors can make six two-way interactions.
avelengthx*volume (Model 10) was not significant. The
restricted model is Model 9.out1b10 <- lmer(consume ~ 1 + difflengthx + avelengthx*volume + (1|groupid), data=dat1, REML=FALSE)
anova(out1b9, out1b10)## Data: dat1
## Models:
## out1b9: consume ~ 1 + difflengthx + avelengthx + volume + (1 | groupid)
## out1b10: consume ~ 1 + difflengthx + avelengthx * volume + (1 | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1b9 6 4116.4 4144.6 -2052.2 4104.4
## out1b10 7 4117.7 4150.6 -2051.9 4103.7 0.7174 1 0.397avelengthx*nfish (Model 12) was not significant. The
restricted model is Model 9 with nfish (Model 11).volume*nfish (Model 13) was not significant. The
restricted model is Model 9 with nfish (Model 11).out1b11 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + nfish + (1|groupid), data=dat1, REML=FALSE)
out1b12 <- lmer(consume ~ 1 + difflengthx + volume + avelengthx*nfish + (1|groupid), data=dat1, REML=FALSE)
out1b13 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume*nfish + (1|groupid), data=dat1, REML=FALSE)
anova(out1b11, out1b12)## Data: dat1
## Models:
## out1b11: consume ~ 1 + difflengthx + avelengthx + volume + nfish + (1 | groupid)
## out1b12: consume ~ 1 + difflengthx + volume + avelengthx * nfish + (1 | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1b11 7 4093.9 4126.8 -2040.0 4079.9
## out1b12 8 4092.8 4130.4 -2038.4 4076.8 3.0826 1 0.07913 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(out1b11, out1b13)## Data: dat1
## Models:
## out1b11: consume ~ 1 + difflengthx + avelengthx + volume + nfish + (1 | groupid)
## out1b13: consume ~ 1 + difflengthx + avelengthx + volume * nfish + (1 | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1b11 7 4093.9 4126.8 -2040.0 4079.9
## out1b13 8 4095.9 4133.5 -2039.9 4079.9 0.0348 1 0.852avelengthx*avegold (Model 15) was not significant. The
restricted model is Model 9 with avegold (Model 14).volume*avegold (Model 16) was not significant. The
restricted model is Model 9 with avegold (Model 14).out1b14 <- lmer(consume ~ 1 + difflengthx + volume + avelengthx + avegold + (1|groupid), data=dat1, REML=FALSE)
out1b15 <- lmer(consume ~ 1 + difflengthx + volume + avelengthx*avegold + (1|groupid), data=dat1, REML=FALSE)
out1b16 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume*avegold + (1|groupid), data=dat1, REML=FALSE)
anova(out1b14, out1b15)## Data: dat1
## Models:
## out1b14: consume ~ 1 + difflengthx + volume + avelengthx + avegold + (1 | groupid)
## out1b15: consume ~ 1 + difflengthx + volume + avelengthx * avegold + (1 | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1b14 7 4116.0 4148.9 -2051.0 4102.0
## out1b15 8 4115.8 4153.4 -2049.9 4099.8 2.1313 1 0.1443anova(out1b14, out1b16)## Data: dat1
## Models:
## out1b14: consume ~ 1 + difflengthx + volume + avelengthx + avegold + (1 | groupid)
## out1b16: consume ~ 1 + difflengthx + avelengthx + volume * avegold + (1 | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1b14 7 4116.0 4148.9 -2051.0 4102.0
## out1b16 8 4115.7 4153.3 -2049.8 4099.7 2.2924 1 0.13avegold*nfish (Model 18) was not significant. The
restricted model is Model 9 with avegold and
nfish (Model 17).out1b17 <- lmer(consume ~ 1 + difflengthx + diffgold + volume + avelengthx + avegold + nfish + (1|groupid), data=dat1, REML=FALSE)
out1b18 <- lmer(consume ~ 1 + difflengthx + diffgold + volume + avelengthx + avegold*nfish + (1|groupid), data=dat1, REML=FALSE)
anova(out1b17, out1b18)## Data: dat1
## Models:
## out1b17: consume ~ 1 + difflengthx + diffgold + volume + avelengthx + avegold + nfish + (1 | groupid)
## out1b18: consume ~ 1 + difflengthx + diffgold + volume + avelengthx + avegold * nfish + (1 | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1b17 9 4088.3 4130.6 -2035.2 4070.3
## out1b18 10 4090.3 4137.3 -2035.1 4070.3 0.0559 1 0.8131All two-way interactions between level-2 predictors are not
significant. Thus, Model 9 is still retained. Next, the random effects
of two level-1 predictors are investigated. First, the difference slopes
of difflengthx (Model 19) is compared with Model 9.
out1b19 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1b9, out1b19)## Data: dat1
## Models:
## out1b9: consume ~ 1 + difflengthx + avelengthx + volume + (1 | groupid)
## out1b19: consume ~ 1 + difflengthx + avelengthx + volume + (1 + difflengthx | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1b9 6 4116.4 4144.6 -2052.2 4104.4
## out1b19 8 4103.4 4141.0 -2043.7 4087.4 17.02 2 0.0002015 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Second, the difference slopes of diffgold (Model 21) is
compared with Model 9 with the fixed effect of diffgold
(Model 20).
out1b20 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + volume + (1|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
out1b21 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + volume + (1 + goldcolor|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1b20, out1b21)## Data: dat1
## Models:
## out1b20: consume ~ 1 + difflengthx + diffgold + avelengthx + volume + (1 | groupid)
## out1b21: consume ~ 1 + difflengthx + diffgold + avelengthx + volume + (1 + goldcolor | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1b20 7 4112.7 4145.6 -2049.3 4098.7
## out1b21 9 4116.4 4158.7 -2049.2 4098.4 0.2953 2 0.8628The random slope of difflengthx was the only significant
effect. Thus, Model 19 is retained.
Next, the cross-level interactions were investigated. There are four
level-2 predictors that can predict the random slopes of
difflengthx.
difflengthx:avelengthx (Model 22) was not significant.
The restricted model is Model 19.difflengthx:volume (Model 23) was significant. The
restricted model is Model 19.out1b22 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + difflengthx:avelengthx + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1b19, out1b22)## Data: dat1
## Models:
## out1b19: consume ~ 1 + difflengthx + avelengthx + volume + (1 + difflengthx | groupid)
## out1b22: consume ~ 1 + difflengthx + avelengthx + volume + difflengthx:avelengthx + (1 + difflengthx | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1b19 8 4103.4 4141.0 -2043.7 4087.4
## out1b22 9 4104.4 4146.6 -2043.2 4086.4 1.0546 1 0.3044out1b23 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + difflengthx:volume + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1b19, out1b23)## Data: dat1
## Models:
## out1b19: consume ~ 1 + difflengthx + avelengthx + volume + (1 + difflengthx | groupid)
## out1b23: consume ~ 1 + difflengthx + avelengthx + volume + difflengthx:volume + (1 + difflengthx | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1b19 8 4103.4 4141.0 -2043.7 4087.4
## out1b23 9 4100.1 4142.4 -2041.0 4082.1 5.3026 1 0.02129 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1difflengthx:nfish (Model 25) was not significant. The
restricted model is Model 19 with nfish (Model 24).out1b24 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + nfish + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
out1b25 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + nfish + difflengthx:nfish + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1b24, out1b25)## Data: dat1
## Models:
## out1b24: consume ~ 1 + difflengthx + avelengthx + volume + nfish + (1 + difflengthx | groupid)
## out1b25: consume ~ 1 + difflengthx + avelengthx + volume + nfish + difflengthx:nfish + (1 + difflengthx | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1b24 9 4079.2 4121.5 -2030.6 4061.2
## out1b25 10 4078.2 4125.2 -2029.1 4058.2 3.0472 1 0.08088 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1difflengthx:avegold (Model 27) was not significant. The
restricted model is Model 19 with avegold (Model 26).out1b26 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + avegold + (1 + difflengthx|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
out1b27 <- lmer(consume ~ 1 + difflengthx + avelengthx + volume + avegold + difflengthx:avegold + (1 + difflengthx|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1b26, out1b27)## Data: dat1
## Models:
## out1b26: consume ~ 1 + difflengthx + avelengthx + volume + avegold + (1 + difflengthx | groupid)
## out1b27: consume ~ 1 + difflengthx + avelengthx + volume + avegold + difflengthx:avegold + (1 + difflengthx | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1b26 9 4100.7 4143.0 -2041.4 4082.7
## out1b27 10 4099.6 4146.6 -2039.8 4079.6 3.1209 1 0.07729 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1difflengthx:volume was the only significant cross-level
interaction. Therefore, Model 23 is the final model from the build-up
strategy.
The starting model (Model 0) for the tear-down strategy has all level-1 (fixed and random effects) and level-2 (fixed effects) predictors as well as all cross-level interactions (fixed effects). Note that the interactions among level-1 predictors and interactions among level-2 predictors are not included to simplify the example.
out1t0 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))First, the cross-level interactions will be investigated. The significant interactions will be retained and the nonsignificant interactions will be dropped. Let’s check each interaction as follows:
difflengthx:avelengthx was not significant. However,
the model that drop this interaction was not convergent. The problem is
fixed by multiplying difflengthx by 10. Then, the model
with the interaction (Model 0a) and without the interaction (Model 1a)
was compared.dat1a <- dat1
dat1a$difflengthx <- dat1a$difflengthx * 10
out1t1a <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold|groupid), data=dat1a, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
out1t0a <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold|groupid), data=dat1a, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1t1a, out1t0a)## Data: dat1a
## Models:
## out1t1a: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold | groupid)
## out1t0a: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1t1a 21 4077.9 4176.6 -2018 4035.9
## out1t0a 22 4079.9 4183.3 -2018 4035.9 0.0127 1 0.9104difflengthx:avegold (Model 2) was not significant when
compared with Model 0.out1t2 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1t2, out1t0)## Data: dat1
## Models:
## out1t2: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold | groupid)
## out1t0: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1t2 21 4081.5 4180.1 -2019.7 4039.5
## out1t0 22 4079.9 4183.3 -2018.0 4035.9 3.5519 1 0.05948 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1difflengthx:volume (Model 3) was not significant when
compared with Model 0.out1t3 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1t3, out1t0)## Data: dat1
## Models:
## out1t3: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold | groupid)
## out1t0: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1t3 21 4080.5 4179.2 -2019.3 4038.5
## out1t0 22 4079.9 4183.3 -2018.0 4035.9 2.6189 1 0.1056difflengthx:nfish (Model 4) was not significant when
compared with Model 0.out1t4 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1t4, out1t0)## Data: dat1
## Models:
## out1t4: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold | groupid)
## out1t0: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1t4 21 4077.9 4176.6 -2018 4035.9
## out1t0 22 4079.9 4183.3 -2018 4035.9 0.0022 1 0.9627diffgold:avelengthx (Model 5) was not significant when
compared with Model 0.out1t5 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1t5, out1t0)## Data: dat1
## Models:
## out1t5: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold | groupid)
## out1t0: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1t5 21 4077.9 4176.6 -2018 4035.9
## out1t0 22 4079.9 4183.3 -2018 4035.9 0.0012 1 0.9723diffgold:avegold (Model 6) was not significant when
compared with Model 0.out1t6 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1t6, out1t0)## Data: dat1
## Models:
## out1t6: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold | groupid)
## out1t0: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1t6 21 4077.9 4176.6 -2018 4035.9
## out1t0 22 4079.9 4183.3 -2018 4035.9 2e-04 1 0.99diffgold:volume (Model 7) was not significant when
compared with Model 0.out1t7 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:nfish + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1t7, out1t0)## Data: dat1
## Models:
## out1t7: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:nfish + (1 + difflengthx + diffgold | groupid)
## out1t0: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1t7 21 4078.2 4176.9 -2018.1 4036.2
## out1t0 22 4079.9 4183.3 -2018.0 4035.9 0.306 1 0.5801diffgold:nfish (Model 8) was not significant when
compared with Model 0.out1t8 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out1t8, out1t0)## Data: dat1
## Models:
## out1t8: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + (1 + difflengthx + diffgold | groupid)
## out1t0: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + difflengthx:avelengthx + difflengthx:avegold + difflengthx:volume + difflengthx:nfish + diffgold:avelengthx + diffgold:avegold + diffgold:volume + diffgold:nfish + (1 + difflengthx + diffgold | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1t8 21 4079.1 4177.8 -2018.5 4037.1
## out1t0 22 4079.9 4183.3 -2018.0 4035.9 1.1763 1 0.2781All cross-level interactions were not significant so they are dropped from the model as Model 9.
out1t9 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + (1 + difflengthx + diffgold|groupid), data=dat1, REML=FALSE)Next, the random effects of both level-1 predictors were investigated.
out1t10 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + (1 + diffgold|groupid), data=dat1, REML=FALSE)
anova(out1t10, out1t9)## Data: dat1
## Models:
## out1t10: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + (1 + diffgold | groupid)
## out1t9: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + (1 + difflengthx + diffgold | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1t10 11 4091.4 4143.1 -2034.7 4069.4
## out1t9 14 4075.6 4141.4 -2023.8 4047.6 21.766 3 7.298e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1out1t11 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1t11, out1t9)## Data: dat1
## Models:
## out1t11: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + (1 + difflengthx | groupid)
## out1t9: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + (1 + difflengthx + diffgold | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1t11 11 4070.9 4122.6 -2024.5 4048.9
## out1t9 14 4075.6 4141.4 -2023.8 4047.6 1.2725 3 0.7357The random effect of difflengthx was significant but the
random effect of diffgold was not significant. Therefore,
Model 11 was retained.
Next, all remaining fixed effects were tested except
difflengthx (because the random effect was
significant).
diffgold (Model 12) was significant so the
effect will be retained.out1t12 <- lmer(consume ~ 1 + difflengthx + avelengthx + avegold + volume + nfish + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1t12, out1t11)## Data: dat1
## Models:
## out1t12: consume ~ 1 + difflengthx + avelengthx + avegold + volume + nfish + (1 + difflengthx | groupid)
## out1t11: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + (1 + difflengthx | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1t12 10 4074.3 4121.3 -2027.2 4054.3
## out1t11 11 4070.9 4122.6 -2024.5 4048.9 5.4279 1 0.01982 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1avelengthx (Model 13) was significant so the
effect will be retained.out1t13 <- lmer(consume ~ 1 + difflengthx + diffgold + avegold + volume + nfish + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1t13, out1t11)## Data: dat1
## Models:
## out1t13: consume ~ 1 + difflengthx + diffgold + avegold + volume + nfish + (1 + difflengthx | groupid)
## out1t11: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + (1 + difflengthx | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1t13 10 4104.1 4151.1 -2042.1 4084.1
## out1t11 11 4070.9 4122.6 -2024.5 4048.9 35.228 1 2.932e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1avegold (Model 14) was significant so the
effect will be retained.out1t14 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + volume + nfish + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1t14, out1t11)## Data: dat1
## Models:
## out1t14: consume ~ 1 + difflengthx + diffgold + avelengthx + volume + nfish + (1 + difflengthx | groupid)
## out1t11: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + (1 + difflengthx | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1t14 10 4075.8 4122.8 -2027.9 4055.8
## out1t11 11 4070.9 4122.6 -2024.5 4048.9 6.8729 1 0.008751 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1volume (Model 15) was significant so the
effect will be retained.out1t15 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + nfish + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1t15, out1t11)## Data: dat1
## Models:
## out1t15: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + nfish + (1 + difflengthx | groupid)
## out1t11: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + (1 + difflengthx | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1t15 10 4105.3 4152.3 -2042.7 4085.3
## out1t11 11 4070.9 4122.6 -2024.5 4048.9 36.423 1 1.588e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1nfish (Model 16) was significant so the effect
will be retained.out1t16 <- lmer(consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + (1 + difflengthx|groupid), data=dat1, REML=FALSE)
anova(out1t16, out1t11)## Data: dat1
## Models:
## out1t16: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + (1 + difflengthx | groupid)
## out1t11: consume ~ 1 + difflengthx + diffgold + avelengthx + avegold + volume + nfish + (1 + difflengthx | groupid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out1t16 10 4097.3 4144.3 -2038.7 4077.3
## out1t11 11 4070.9 4122.6 -2024.5 4048.9 28.42 1 9.764e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1All fixed effects were significant. Thus, the final model is Model 11.
Read the new data set which is the original interview dataset in Lecture 4 with additional covariates.
dat2 <- read.table("lecture9ex1.csv", sep=",", header=TRUE)Based on the research hypothesis, the interaction of interviewers’ and interviewees’ sex, the random effect of the interviewee’s sex, and all fixed effects of all covariates are included in the model (Model 0).
library(lme4)
out2m0 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex|erid), data=dat2, REML=FALSE)
summary(out2m0)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex +
## erexp + eesex:ersex + (1 + eesex | erid)
## Data: dat2
##
## AIC BIC logLik deviance df.resid
## 60259.1 60345.7 -30117.6 60235.1 9988
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9117 -0.6471 -0.0151 0.6477 3.3815
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## erid (Intercept) 18.402 4.290
## eesex 1.842 1.357 0.02
## Residual 18.585 4.311
## Number of obs: 10000, groups: erid, 1000
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 73.668030 0.221986 331.859
## eesex -0.530397 0.136242 -3.893
## eesalary -0.008172 0.054475 -0.150
## eeworkexp 0.337878 0.014556 23.212
## I(eeiq - 100) 0.028475 0.003065 9.291
## ersex 0.221768 0.297455 0.746
## erexp 0.011689 0.011099 1.053
## eesex:ersex 0.957125 0.192642 4.968
##
## Correlation of Fixed Effects:
## (Intr) eesex eeslry ewrkxp I(-100 ersex erexp
## eesex -0.237
## eesalary 0.151 -0.012
## eeworkexp -0.286 0.006 -0.534
## I(eeiq-100) 0.024 -0.006 -0.059 -0.050
## ersex -0.671 0.178 0.000 0.002 -0.007
## erexp -0.141 -0.019 -0.006 -0.006 -0.004 0.000
## eesex:ersex 0.172 -0.707 0.012 -0.011 0.007 -0.251 0.000First, all covariates are checked whether they are needed in the
model. There are four fixed effects that can be dropped (excluding
eesex and ersex):
eesalary can be dropped. First, check
whether eesalary has a random effect (Model 1). The random
effect was significant so the random effect is retained.out2m1 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex + eesalary|erid), data=dat2, REML=FALSE)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00225309 (tol = 0.002, component 1)anova(out2m1, out2m0)## Data: dat2
## Models:
## out2m0: score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex | erid)
## out2m1: score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex + eesalary | erid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out2m0 12 60259 60346 -30118 60235
## out2m1 15 60257 60365 -30113 60227 8.351 3 0.03929 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1eeworkexp can be dropped. First, check
whether eeworkexp has a random effect (Model 2). The random
effect was not significant. Next, dropping eeworkexp (Model
3) from the model was significant so this fixed effect will be
retained.out2m2 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex + eeworkexp|erid), data=dat2, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out2m2, out2m0)## Data: dat2
## Models:
## out2m0: score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex | erid)
## out2m2: score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex + eeworkexp | erid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out2m0 12 60259 60346 -30118 60235
## out2m2 15 60259 60367 -30114 60229 6.6109 3 0.08539 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1out2m3 <- lmer(score ~ 1 + eesex + eesalary + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex|erid), data=dat2, REML=FALSE)
anova(out2m3, out2m0)## Data: dat2
## Models:
## out2m3: score ~ 1 + eesex + eesalary + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex | erid)
## out2m0: score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex | erid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out2m3 11 60781 60860 -30379 60759
## out2m0 12 60259 60346 -30118 60235 523.5 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1eeiq can be dropped. First, check
whether eeworkexp has a random effect (Model 4). After some
transformations to make the model convergent, the random effect was not
significant. Next, dropping eeiq (Model 5) from the model
was significant so this fixed effect will be retained.out2m0a <- lmer(score ~ 1 + eesex + eesalary + I(eeworkexp/10) + I(eeiq/10 - 10) + ersex + erexp + eesex:ersex + (1 + eesex|erid), data=dat2, REML=FALSE)
out2m4a <- lmer(score ~ 1 + eesex + eesalary + I(eeworkexp/10) + I(eeiq/10 - 10) + ersex + erexp + eesex:ersex + (1 + eesex + I(eeiq/10 - 10)|erid), data=dat2, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out2m4a, out2m0a)## Data: dat2
## Models:
## out2m0a: score ~ 1 + eesex + eesalary + I(eeworkexp/10) + I(eeiq/10 - 10) + ersex + erexp + eesex:ersex + (1 + eesex | erid)
## out2m4a: score ~ 1 + eesex + eesalary + I(eeworkexp/10) + I(eeiq/10 - 10) + ersex + erexp + eesex:ersex + (1 + eesex + I(eeiq/10 - 10) | erid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out2m0a 12 60259 60346 -30118 60235
## out2m4a 15 60264 60373 -30117 60234 0.6676 3 0.8808out2m5 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + ersex + erexp + eesex:ersex + (1 + eesex|erid), data=dat2, REML=FALSE)
anova(out2m5, out2m0)## Data: dat2
## Models:
## out2m5: score ~ 1 + eesex + eesalary + eeworkexp + ersex + erexp + eesex:ersex + (1 + eesex | erid)
## out2m0: score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex | erid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out2m5 11 60343 60422 -30161 60321
## out2m0 12 60259 60346 -30118 60235 85.881 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1erexp is the level-2 predictor so the random effect is
not tested. Dropping erexp (Model 6) was not significant so
this fixed effect will be dropped.out2m6 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + eesex:ersex + (1 + eesex|erid), data=dat2, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
anova(out2m6, out2m0)## Data: dat2
## Models:
## out2m6: score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + eesex:ersex + (1 + eesex | erid)
## out2m0: score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex | erid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out2m6 11 60258 60338 -30118 60236
## out2m0 12 60259 60346 -30118 60235 1.109 1 0.2923From the results, as shown in Model 7, the random effect of
eesalary will be added and the fixed effect of
erexp is dropped.
out2m7 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + eesex:ersex + (1 + eesex + eesalary|erid), data=dat2, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))Next, the cross-level interaction will be examined. There are three combination that can be added.
eesalary:ersex. Model 8 was not
significantly better than Model 7. Thus, this interaction will be
dropped.out2m8 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + eesex:ersex + eesalary:ersex + (1 + eesex + eesalary|erid), data=dat2, REML=FALSE)
anova(out2m8, out2m7)## Data: dat2
## Models:
## out2m7: score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + eesex:ersex + (1 + eesex + eesalary | erid)
## out2m8: score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + eesex:ersex + eesalary:ersex + (1 + eesex + eesalary | erid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out2m7 14 60256 60357 -30114 60228
## out2m8 15 60258 60366 -30114 60228 0.3538 1 0.5519eesalary:erexp (and erexp)
in the model. Model 10 is compared with Model 9, which is Model 7 with
erexp. The models were not significantly different so
eesalary:erexp is dropped.out2m9 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex + eesalary|erid), data=dat2, REML=FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
out2m10 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + eesalary:erexp + (1 + eesex + eesalary|erid), data=dat2, REML=FALSE)
anova(out2m10, out2m9)## Data: dat2
## Models:
## out2m9: score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex + eesalary | erid)
## out2m10: score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + eesalary:erexp + (1 + eesex + eesalary | erid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out2m9 15 60257 60365 -30113 60227
## out2m10 16 60258 60373 -30113 60226 0.6726 1 0.4121eesex:erexp (and erexp) in
the model. Model 11 was not significantly better than Model 9, which is
Model 7 with erexp. Then, eesalary:erexp is
dropped.out2m11 <- lmer(score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + eesex:erexp + (1 + eesex + eesalary|erid), data=dat2, REML=FALSE)
anova(out2m11, out2m9)## Data: dat2
## Models:
## out2m9: score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + (1 + eesex + eesalary | erid)
## out2m11: score ~ 1 + eesex + eesalary + eeworkexp + I(eeiq - 100) + ersex + erexp + eesex:ersex + eesex:erexp + (1 + eesex + eesalary | erid)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## out2m9 15 60257 60365 -30113 60227
## out2m11 16 60258 60374 -30113 60226 0.5213 1 0.4703All cross-level interactions were not significant. Therefore, Model 7 is the final model.