neighborhood <- read.dta("http://www.stata-press.com/data/mlmus3/neighborhood.dta")
head(neighborhood)
##   neighid schid  attain   p7vrq  p7read  dadocc dadunemp daded momed male
## 1     675     0  0.7350  21.972  12.134   2.316        0     0     0    1
## 2     647     0  0.2641  -7.028 -12.866  16.196        0     0     1    0
## 3     650     0 -1.3276 -11.028 -31.866 -23.454        1     0     0    1
## 4     650     0  0.7350   3.972   3.134   2.316        0     0     0    1
## 5     648     0 -0.1325  -2.028   0.134  -3.454        0     0     0    0
## 6     648     0  0.5610  -5.028  -0.866  -3.454        0     0     0    0
##   deprive dummy
## 1  -0.182     1
## 2   0.206     1
## 3   0.534     1
## 4   0.534     1
## 5   0.188     1
## 6   0.188     1
  1. Fit a random-intercept model with attain as the response variable and without any covariates by ML using xtmixed. What are the estimated variance components between and within neighborhoods? Obtain the estimated intraclass correlation.
model1 <- lmer(attain~(1|neighid), data=neighborhood, REML = FALSE)
summary(model1)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: attain ~ (1 | neighid)
##    Data: neighborhood
## 
##      AIC      BIC   logLik deviance df.resid 
##   6422.0   6439.2  -3208.0   6416.0     2307 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.33164 -0.65532  0.01513  0.58177  2.96174 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  neighid  (Intercept) 0.2015   0.4489  
##  Residual             0.8044   0.8969  
## Number of obs: 2310, groups:  neighid, 524
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.08202    0.02844   2.885

Calculate ICC

RanEffects<-as.data.frame(VarCorr(model1))
ICC<- RanEffects[1,4]/(RanEffects[1,4]+RanEffects[2,4])
print(ICC)
## [1] 0.2003543
  1. Include the covariate deprive in the model and interpret the estimates. Discuss the changes in the estimated standard deviations of the random intercept and level-1 residual.
model2 <- lmer(attain~ deprive+(1|neighid), data=neighborhood, REML = FALSE)
summary(model2)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: attain ~ deprive + (1 | neighid)
##    Data: neighborhood
## 
##      AIC      BIC   logLik deviance df.resid 
##   6273.6   6296.6  -3132.8   6265.6     2306 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2846 -0.6513  0.0081  0.5848  3.4957 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  neighid  (Intercept) 0.08579  0.2929  
##  Residual             0.81220  0.9012  
## Number of obs: 2310, groups:  neighid, 524
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.09995    0.02367   4.222
## deprive     -0.52077    0.03824 -13.619
## 
## Correlation of Fixed Effects:
##         (Intr)
## deprive -0.045
  1. Include the student-level covariates and interpret the estimates. Also comment on how the estimated standard deviations have changed.
model3 <- lmer(attain~ deprive+p7vrq+p7read+dadocc+dadunemp+daded+momed+male+(1|neighid), data=neighborhood, REML=FALSE)
# model4 <- lmer(attain~ p7vrq+p7read+dadocc+dadunemp+daded+momed+male+(deprive|neighid), data=neighborhood, REML=FALSE)
# model5 <- lmer(attain~ p7vrq+p7read+dadocc+dadunemp+daded+momed+male+(1|neighid)+(1|deprive), data=neighborhood, REML=FALSE)

summary(model3)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: attain ~ deprive + p7vrq + p7read + dadocc + dadunemp + daded +  
##     momed + male + (1 | neighid)
##    Data: neighborhood
## 
##      AIC      BIC   logLik deviance df.resid 
##   4797.0   4860.2  -2387.5   4775.0     2299 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9354 -0.6413 -0.0348  0.5706  3.5381 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  neighid  (Intercept) 0.006228 0.07892 
##  Residual             0.456631 0.67574 
## Number of obs: 2310, groups:  neighid, 524
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  0.080731   0.022959   3.516
## deprive     -0.148123   0.025331  -5.848
## p7vrq        0.027763   0.002256  12.306
## p7read       0.026065   0.001746  14.928
## dadocc       0.008239   0.001363   6.044
## dadunemp    -0.114896   0.046830  -2.453
## daded        0.140966   0.040810   3.454
## momed        0.062405   0.037454   1.666
## male        -0.055381   0.028434  -1.948
## 
## Correlation of Fixed Effects:
##          (Intr) depriv p7vrq  p7read dadocc dadnmp daded  momed 
## deprive  -0.026                                                 
## p7vrq    -0.108  0.057                                          
## p7read    0.097  0.083 -0.767                                   
## dadocc    0.100  0.155 -0.056 -0.083                            
## dadunemp -0.228 -0.110  0.038  0.017  0.123                     
## daded    -0.214  0.038  0.002 -0.063 -0.212 -0.009              
## momed    -0.240  0.023 -0.017 -0.023 -0.065  0.002 -0.419       
## male     -0.596  0.006  0.085 -0.050  0.011  0.020 -0.005 -0.013
# summary(model4)
# summary(model5)
  1. Obtain the overall coefficient of determination R2 for the model in step 3.
#Joop Hox - P. 69 Also same as sjstats r2(model3, model1)
(VarCorr(model1)$neighid[1]-VarCorr(model3)$neighid[1])/VarCorr(model1)$neighid[1] #R-squared (tau-00):
## [1] 0.9690982
#R Squared
(attr(VarCorr(model1),"sc")^2-attr(VarCorr(model3),"sc")^2)/attr(VarCorr(model1),"sc")^2
## [1] 0.4323131
#Same as http://web.pdx.edu/~newsomj/mlrclass/ho_r2.pdf 
(VarCorr(model3)$neighid[1]+attr(VarCorr(model3),"sc")^2)/
  (VarCorr(model1)$neighid[1]+attr(VarCorr(model1),"sc")^2)
## [1] 0.4601397