library(haven)
popula <- read_dta("https://stats.idre.ucla.edu/stat/stata/examples/mlm_ma_hox/popular.dta")
head(popula)library(nlme)
m4_lme <- lme(popular ~ sex*texp, data = popula, random = ~ sex|school, method = "ML")
summary(m4_lme)Linear mixed-effects model fit by maximum likelihood
Data: popula
AIC BIC logLik
4261.85 4306.657 -2122.925
Random effects:
Formula: ~sex | school
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.6347377 (Intr)
sex 0.4692521 0.08
Residual 0.6264320
Fixed effects: popular ~ sex * texp
Value Std.Error DF t-value p-value
(Intercept) 3.313651 0.15954654 1898 20.769180 0e+00
sex 1.329479 0.13183479 1898 10.084432 0e+00
texp 0.110229 0.01013882 98 10.872007 0e+00
sex:texp -0.034025 0.00837995 1898 -4.060303 1e-04
Correlation:
(Intr) sex texp
sex -0.046
texp -0.909 0.042
sex:texp 0.042 -0.908 -0.046
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.93387434 -0.64733763 0.02294381 0.53272602 3.49150352
Number of Observations: 2000
Number of Groups: 100 library(magrittr)
library(dplyr)
popula %<>% mutate(ctexp = texp - mean(texp))
m4a_lme <- lme(popular ~ sex*ctexp, data = popula, random = ~ sex|school, method = "ML")
summary(m4a_lme)Linear mixed-effects model fit by maximum likelihood
Data: popula
AIC BIC logLik
4261.85 4306.657 -2122.925
Random effects:
Formula: ~sex | school
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.6347377 (Intr)
sex 0.4692521 0.08
Residual 0.6264320
Fixed effects: popular ~ sex * ctexp
Value Std.Error DF t-value p-value
(Intercept) 4.885851 0.06660875 1898 73.35149 0e+00
sex 0.844178 0.05510824 1898 15.31855 0e+00
ctexp 0.110229 0.01013882 98 10.87201 0e+00
sex:ctexp -0.034025 0.00837995 1898 -4.06030 1e-04
Correlation:
(Intr) sex ctexp
sex -0.044
ctexp -0.006 0.000
sex:ctexp 0.000 -0.004 -0.046
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.93387434 -0.64733763 0.02294381 0.53272602 3.49150352
Number of Observations: 2000
Number of Groups: 100 library(texreg)
htmlreg(list(m4_lme, m4a_lme))| Model 1 | Model 2 | ||
|---|---|---|---|
| (Intercept) | 3.31*** | 4.89*** | |
| (0.16) | (0.07) | ||
| sex | 1.33*** | 0.84*** | |
| (0.13) | (0.06) | ||
| texp | 0.11*** | ||
| (0.01) | |||
| sex:texp | -0.03*** | ||
| (0.01) | |||
| ctexp | 0.11*** | ||
| (0.01) | |||
| sex:ctexp | -0.03*** | ||
| (0.01) | |||
| AIC | 4261.85 | 4261.85 | |
| BIC | 4306.66 | 4306.66 | |
| Log Likelihood | -2122.92 | -2122.92 | |
| Num. obs. | 2000 | 2000 | |
| Num. groups | 100 | 100 | |
| p < 0.001, p < 0.01, p < 0.05 | |||
The first model on slide 31 of the lecture demonstrates a cross-level interaction between sex and teachers experience to see the effect on popularity. The second models refits the model by centering the covariates on slide 29 and uses the centered covariates in the second model. Centering the independent variables (teachers experience in this case) redefines the 0 point for the predictor variable to be the mean of teacher’s experience. This shifts the scale over but retains the unit. The effect of the slope between teachers experience and popularity does not’t change at all but the interpretation of the intercept does. The intercept of the second model is the mean popularity when all the independent variables equal 0. Since the second model has turned “0” into the average score through centering,the intercept is now recognized as the value of popularity when teacher’s experience is set to the average.