### R Textbook Examples Multilevel Analysis Techniques and Applications by Joop Hox Chapter 4: Some Important Methodological and Statistical Issues

library(foreign)
library(lme4)


Table 4.1 on page 57 using the popular dataset.

popdata<-read.dta("http://www.ats.ucla.edu/stat/stata/examples/mlm_ma_hox/popular.dta")
Part 1
m4.1.1 <- lmer(popular~sex + (1|school), popdata)
summary(m4.1.1)

Linear mixed model fit by REML
Formula: popular ~ sex + (1 | school)
Data: popdata
AIC  BIC logLik deviance REMLdev
4501 4523  -2246     4485    4493
Random effects:
Groups   Name        Variance Std.Dev.
school   (Intercept) 0.86222  0.92856
Residual             0.45992  0.67817
Number of obs: 2000, groups: school, 100

Fixed effects:
Estimate Std. Error t value
(Intercept)  4.89722    0.09529   51.39
sexgirl      0.84370    0.03096   27.25

Correlation of Fixed Effects:
(Intr)
sexgirl -0.158
Part 2
attach(popdata)
sex01 <- 1*(sex=="girl")

csex <- sex01 - mean(sex01)
m4.1.2 <- lmer(popular~csex + (1|school))
summary(m4.1.2)

Linear mixed model fit by REML
Formula: popular ~ csex + (1 | school)
AIC  BIC logLik deviance REMLdev
4501 4523  -2246     4485    4493
Random effects:
Groups   Name        Variance Std.Dev.
school   (Intercept) 0.86222  0.92856
Residual             0.45992  0.67817
Number of obs: 2000, groups: school, 100

Fixed effects:
Estimate Std. Error t value
(Intercept)  5.30810    0.09409   56.41
csex         0.84370    0.03096   27.25

Correlation of Fixed Effects:
(Intr)
csex 0.000

Part 3
m4.1.3 <- lmer(popular~sex + (sex|school))
summary(m4.1.3)

Linear mixed model fit by REML
Formula: popular ~ sex + (sex | school)
AIC  BIC logLik deviance REMLdev
4348 4382  -2168     4330    4336
Random effects:
Groups   Name        Variance Std.Dev. Corr
school   (Intercept) 0.94027  0.96968
sexgirl     0.27252  0.52203  -0.279
Residual             0.39244  0.62645
Number of obs: 2000, groups: school, 100

Fixed effects:
Estimate Std. Error t value
(Intercept)  4.89009    0.09902   49.39
sexgirl      0.84311    0.05964   14.14

Correlation of Fixed Effects:
(Intr)
sexgirl -0.307

Part 4
m4.1.4 <- lmer(popular~csex + (csex|school))
summary(m4.1.4)

Linear mixed model fit by REML
Formula: popular ~ csex + (csex | school)
AIC  BIC logLik deviance REMLdev
4348 4382  -2168     4330    4336
Random effects:
Groups   Name        Variance Std.Dev. Corr
school   (Intercept) 0.86758  0.93144
csex        0.27252  0.52204  -0.017
Residual             0.39244  0.62645
Number of obs: 2000, groups: school, 100

Fixed effects:
Estimate Std. Error t value
(Intercept)  5.30068    0.09424   56.24
csex         0.84311    0.05963   14.14

Correlation of Fixed Effects:
(Intr)
csex -0.015
Table 4.2 on page 60, continuing to use the popular dataset.
Part 1
m4.2.1 <- lmer(popular~sex + texp + (1 + sex|school))
summary(m4.2.1)

Linear mixed model fit by REML
Formula: popular ~ sex + texp + (1 + sex | school)
AIC  BIC logLik deviance REMLdev
4290 4329  -2138     4261    4276
Random effects:
Groups   Name        Variance Std.Dev. Corr
school   (Intercept) 0.41158  0.64155
sexgirl     0.27329  0.52278  0.062
Residual             0.39248  0.62648
Number of obs: 2000, groups: school, 100

Fixed effects:
Estimate Std. Error t value
(Intercept)  3.34001    0.16079   20.77
sexgirl      0.84315    0.05969   14.13
texp         0.10835    0.01022   10.61

Correlation of Fixed Effects:
(Intr) sexgrl
sexgirl -0.020
texp    -0.908  0.000

Part 2
m4.2.2 <- lmer(popular~sex + texp + sex*texp + (1 + sex|school))
summary(m4.2.2)

Linear mixed model fit by REML
Formula: popular ~ sex + texp + sex * texp + (1 + sex | school)
AIC  BIC logLik deviance REMLdev
4284 4329  -2134     4246    4268
Random effects:
Groups   Name        Variance Std.Dev. Corr
school   (Intercept) 0.41198  0.64186
sexgirl     0.22641  0.47582  0.077
Residual             0.39241  0.62643
Number of obs: 2000, groups: school, 100

Fixed effects:
Estimate Std. Error t value
(Intercept)   3.313521   0.161017  20.579
sexgirl       1.329594   0.133052   9.993
texp          0.110235   0.010232  10.773
sexgirl:texp -0.034035   0.008457  -4.024

Correlation of Fixed Effects:
(Intr) sexgrl texp
sexgirl     -0.046
texp        -0.909  0.042
sexgirl:txp  0.042 -0.908 -0.046

Part 3
ctexp <- texp - mean(texp)
m4.2.3 <- lmer(popular~csex + ctexp + (1 + csex|school))
summary(m4.2.3)

Linear mixed model fit by REML
Formula: popular ~ csex + ctexp + (1 + csex | school)
AIC  BIC logLik deviance REMLdev
4290 4329  -2138     4261    4276
Random effects:
Groups   Name        Variance Std.Dev. Corr
school   (Intercept) 0.49674  0.70480
csex        0.27330  0.52278  0.418
Residual             0.39248  0.62648
Number of obs: 2000, groups: school, 100

Fixed effects:
Estimate Std. Error t value
(Intercept)  5.29601    0.07193   73.63
csex         0.84315    0.05969   14.13
ctexp        0.10835    0.01022   10.61

Correlation of Fixed Effects:
(Intr) csex
csex   0.359
ctexp -0.005  0.000

Part 4
m4.2.4 <- lmer(popular~csex + ctexp + csex*ctexp + (1 + csex|school))
summary(m4.2.4)

Linear mixed model fit by REML
Formula: popular ~ csex + ctexp + csex * ctexp + (1 + csex | school)
AIC  BIC logLik deviance REMLdev
4284 4329  -2134     4246    4268
Random effects:
Groups   Name        Variance Std.Dev. Corr
school   (Intercept) 0.48849  0.69892
csex        0.22641  0.47583  0.402
Residual             0.39241  0.62643
Number of obs: 2000, groups: school, 100

Fixed effects:
Estimate Std. Error t value
(Intercept)  5.296908   0.071346   74.24
csex         0.844154   0.055616   15.18
ctexp        0.093660   0.010851    8.63
csex:ctexp  -0.034035   0.008457   -4.02

Correlation of Fixed Effects:
(Intr) csex   ctexp
csex        0.338
ctexp      -0.005 -0.002
csex:ctexp -0.002 -0.004  0.336

Figure 4.3 on page 61.
m4.3 <- lmer(popular~sex + texp + sex*texp + (1 + sex|school))
sextexp = texp * sex01
linp <- fixef(m4.3)[1] + fixef(m4.3)[2]*sex01 + fixef(m4.3)[3]*texp + fixef(m4.3)[4]*sex01*texp

plot(x = texp[sex01==1], y = linp[sex01==1], type = "l",
xlab = "teacher experience in years", ylab = "predicted values",
xlim = c(0,30), ylim = c(3.5, 7), col = "red")
points(x = texp[sex01==0], y = linp[sex01==0], type = "l", col = "blue")
legend("bottomright", c("girls", "boys"), col = c("red", "blue"), lty = c(1,1))

Table 4.3 on page 63, continuing to use the popular dataset.
Part 1
m4.3.1 <- lmer(popular~ (1 |school))
summary(m4.3.1)

Linear mixed model fit by REML
Formula: popular ~ (1 | school)
AIC  BIC logLik deviance REMLdev
5122 5138  -2558     5113    5116
Random effects:
Groups   Name        Variance Std.Dev.
school   (Intercept) 0.87981  0.93798
Residual             0.63868  0.79917
Number of obs: 2000, groups: school, 100

Fixed effects:
Estimate Std. Error t value
(Intercept)   5.3076     0.0955   55.58

Part 2
m4.3.2 <- lmer(popular~ sex + (1 |school))
summary(m4.3.2)

Linear mixed model fit by REML
Formula: popular ~ sex + (1 | school)
AIC  BIC logLik deviance REMLdev
4501 4523  -2246     4485    4493
Random effects:
Groups   Name        Variance Std.Dev.
school   (Intercept) 0.86222  0.92856
Residual             0.45992  0.67817
Number of obs: 2000, groups: school, 100

Fixed effects:
Estimate Std. Error t value
(Intercept)  4.89722    0.09529   51.39
sexgirl      0.84370    0.03096   27.25

Correlation of Fixed Effects:
(Intr)
sexgirl -0.158

Part 3
m4.3.3 <- lmer(popular~ sex + texp + (1 |school))
summary(m4.3.3)

Linear mixed model fit by REML
Formula: popular ~ sex + texp + (1 | school)
AIC  BIC logLik deviance REMLdev
4454 4482  -2222     4429    4444
Random effects:
Groups   Name        Variance Std.Dev.
school   (Intercept) 0.48596  0.69711
Residual             0.45992  0.67818
Number of obs: 2000, groups: school, 100

Fixed effects:
Estimate Std. Error t value
(Intercept)  3.56068    0.17147   20.77
sexgirl      0.84467    0.03095   27.29
texp         0.09345    0.01085    8.61

Correlation of Fixed Effects:
(Intr) sexgrl
sexgirl -0.088
texp    -0.905  0.000

Part 4
m4.3.4 <- lmer(popular~ sex + texp + (1 + sex|school))
summary(m4.3.4)

Linear mixed model fit by REML
Formula: popular ~ sex + texp + (1 + sex | school)
AIC  BIC logLik deviance REMLdev
4290 4329  -2138     4261    4276
Random effects:
Groups   Name        Variance Std.Dev. Corr
school   (Intercept) 0.41158  0.64155
sexgirl     0.27329  0.52278  0.062
Residual             0.39248  0.62648
Number of obs: 2000, groups: school, 100

Fixed effects:
Estimate Std. Error t value
(Intercept)  3.34001    0.16079   20.77
sexgirl      0.84315    0.05969   14.13
texp         0.10835    0.01022   10.61

Correlation of Fixed Effects:
(Intr) sexgrl
sexgirl -0.020
texp    -0.908  0.000

Part 5
m4.3.5 <- lmer(popular~ sex + texp + sex*texp + (1 + sex|school))
summary(m4.3.5)

Linear mixed model fit by REML
Formula: popular ~ sex + texp + sex * texp + (1 + sex | school)
AIC  BIC logLik deviance REMLdev
4284 4329  -2134     4246    4268
Random effects:
Groups   Name        Variance Std.Dev. Corr
school   (Intercept) 0.41198  0.64186
sexgirl     0.22641  0.47582  0.077
Residual             0.39241  0.62643
Number of obs: 2000, groups: school, 100

Fixed effects:
Estimate Std. Error t value
(Intercept)   3.313521   0.161017  20.579
sexgirl       1.329594   0.133052   9.993
texp          0.110235   0.010232  10.773
sexgirl:texp -0.034035   0.008457  -4.024

Correlation of Fixed Effects:
(Intr) sexgrl texp
sexgirl     -0.046
texp        -0.909  0.042
sexgirl:txp  0.042 -0.908 -0.046

Table 4.4 on page 70--we have skipped this table for now.

The content of this web site should not be construed as an endorsement of any particular web site, book, or software product by the University of California.