This document reproduces examples in lecture and illustrates how to use various functions that I wrote. These functions do the following:
icc: Computes interclass correlation for 2-level modelsbic.hlm(model,cluster.id): Computes BIC statistic that balences between and withinrobust(model.fit, response, idvar, dftype): Computes Sandwich standard errors and reports model based and sandwich bases statistics.contrast(model,L) where L is c x p matrix of contrast. This gives Wald and F-tests results hlmRsq(dataset, model, id): Computes R_1^2 and R_2^2All of these functions are in a single file called “All_functions.txt”. This can be set as a source file. For more information on these functions see the lectures notes and notes in “All_functions.txt”.
For some of these functions, I implicitly assume that the data were such that rows from same group or cluster are grouped together in a block of the data frame. (For large data sets this might reduce computation time needed to fit models to data).
setwd("D:\\Dropbox\\edps587\\lectures\\6 inference1")
source("All_functions.txt")library(lme4)## Loading required package: Matrixlibrary(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steplibrary(texreg)## 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").library(optimx)
library(xtable) # for Latex output of data frame# Read in data files
hsb1 <- read.table("HSB1.txt", header=TRUE)
hsb2 <- read.table("HSB2.txt", header=TRUE)
# Merge into one file
hsb.tmp <- merge(hsb1,hsb2, by=c('id'))
# Center ses
hsb.tmp$cses <- hsb.tmp$ses - hsb.tmp$meanses
# Make sure data are ordered/sorted by id.
hsb <- hsb.tmp[order(hsb.tmp$id),]
hsb$id[1:50]## [1] 1224 1224 1224 1224 1224 1224 1224 1224 1224 1224 1224 1224 1224 1224 1224
## [16] 1224 1224 1224 1224 1224 1224 1224 1224 1224 1224 1224 1224 1224 1224 1224
## [31] 1224 1224 1224 1224 1224 1224 1224 1224 1224 1224 1224 1224 1224 1224 1224
## [46] 1224 1224 1288 1288 1288One last step for this example: to avoid “warning different re-scale”
hsb$zsize <- (hsb$size)/sd(hsb$size)We will work with the following set of models
summary(model1 <- lmer(mathach ~ 1 + cses + female + minority + meanses + zsize + sector + (1 |id), data=hsb, REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + female + minority + meanses + zsize + sector +
## (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46307.3 46369.3 -23144.7 46289.3 7176
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2641 -0.7225 0.0410 0.7599 2.9191
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.61 1.269
## Residual 35.89 5.991
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.4000 0.3576 172.2298 34.672 < 2e-16 ***
## cses 1.9127 0.1084 7076.1669 17.644 < 2e-16 ***
## female -1.2402 0.1586 5173.8657 -7.821 6.29e-15 ***
## minority -2.8871 0.2009 3610.0411 -14.373 < 2e-16 ***
## meanses 3.9634 0.3324 170.3640 11.924 < 2e-16 ***
## zsize 0.4154 0.1355 159.0317 3.065 0.00256 **
## sector 2.1179 0.2979 153.0838 7.109 4.16e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses female minrty meanss zsize
## cses -0.016
## female -0.245 0.050
## minority -0.041 0.153 0.016
## meanses 0.122 0.041 0.048 0.256
## zsize -0.837 -0.013 0.019 -0.089 -0.056
## sector -0.653 -0.023 -0.006 -0.148 -0.360 0.437summary(model2 <- lmer(mathach ~ 1 + cses + female + minority + meanses + zsize + sector + cses*zsize + (1 |id), data=hsb, REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + female + minority + meanses + zsize + sector +
## cses * zsize + (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46301.6 46370.4 -23140.8 46281.6 7175
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3048 -0.7227 0.0349 0.7618 2.9425
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.61 1.269
## Residual 35.85 5.987
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.3967 0.3576 172.2227 34.667 < 2e-16 ***
## cses 1.3867 0.2176 7027.2675 6.373 1.97e-10 ***
## female -1.2422 0.1585 5176.8006 -7.838 5.53e-15 ***
## minority -2.8687 0.2009 3608.6721 -14.280 < 2e-16 ***
## meanses 3.9711 0.3324 170.3710 11.948 < 2e-16 ***
## zsize 0.4161 0.1355 159.0309 3.070 0.00252 **
## sector 2.1139 0.2979 153.0911 7.097 4.46e-11 ***
## cses:zsize 0.2909 0.1044 7023.0769 2.787 0.00533 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses female minrty meanss zsize sector
## cses -0.005
## female -0.245 0.029
## minority -0.041 0.048 0.015
## meanses 0.122 0.013 0.048 0.256
## zsize -0.837 -0.008 0.019 -0.089 -0.056
## sector -0.653 -0.007 -0.006 -0.148 -0.360 0.437
## cses:zsize -0.003 -0.867 -0.005 0.033 0.008 0.002 -0.005summary(model3 <- lmer(mathach ~ 1 + cses + female + minority + meanses + zsize + sector
+ cses*zsize + cses*sector + female*meanses + female*zsize + female*sector
+ minority*meanses + minority*zsize + minority*sector
+ (1 + cses|id), data=hsb, REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + female + minority + meanses + zsize + sector +
## cses * zsize + cses * sector + female * meanses + female *
## zsize + female * sector + minority * meanses + minority *
## zsize + minority * sector + (1 + cses | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46268.7 46399.4 -23115.4 46230.7 7166
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1998 -0.7180 0.0320 0.7585 3.0068
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 1.6336 1.2781
## cses 0.1567 0.3959 0.08
## Residual 35.5190 5.9598
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.01798 0.43268 318.08997 27.776 < 2e-16 ***
## cses 2.26181 0.31214 151.43779 7.246 2.04e-11 ***
## female -0.35821 0.42171 6761.19516 -0.849 0.395678
## minority -4.08707 0.61827 3074.11945 -6.610 4.50e-11 ***
## meanses 4.23444 0.46291 404.78301 9.147 < 2e-16 ***
## zsize 0.69323 0.17582 359.44252 3.943 9.68e-05 ***
## sector 1.73886 0.38075 308.07750 4.567 7.16e-06 ***
## cses:zsize 0.03215 0.12291 157.93699 0.262 0.794000
## cses:sector -0.99850 0.25375 151.90533 -3.935 0.000126 ***
## female:meanses -0.05433 0.44329 2820.95011 -0.123 0.902460
## female:zsize -0.41608 0.16767 6582.42896 -2.481 0.013109 *
## female:sector -0.23677 0.39242 2664.68437 -0.603 0.546322
## minority:meanses -0.80456 0.49179 2167.15439 -1.636 0.101991
## minority:zsize 0.06961 0.22129 2386.42217 0.315 0.753109
## minority:sector 2.06971 0.49486 3273.08184 4.182 2.96e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itsummary(model4 <- lmer(mathach ~ 1 + female + minority + cses + meanses + zsize + sector
+ cses*zsize + cses*sector + female*meanses + female*zsize + female*sector
+ minority*meanses + minority*zsize + minority*sector
+ (1 + female + minority + cses |id), data=hsb, REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + female + minority + cses + meanses + zsize + sector +
## cses * zsize + cses * sector + female * meanses + female *
## zsize + female * sector + minority * meanses + minority *
## zsize + minority * sector + (1 + female + minority + cses | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46272.8 46451.7 -23110.4 46220.8 7159
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2037 -0.7149 0.0334 0.7609 2.9991
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 2.2268 1.4923
## female 0.7022 0.8380 -0.76
## minority 0.8756 0.9358 -0.12 0.24
## cses 0.1446 0.3803 0.26 -0.32 -0.35
## Residual 35.3183 5.9429
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.02587 0.46923 139.24621 25.629 < 2e-16 ***
## female -0.34012 0.46377 123.25042 -0.733 0.464723
## minority -4.20661 0.67149 157.13987 -6.265 3.43e-09 ***
## cses 2.28237 0.31118 148.15982 7.335 1.35e-11 ***
## meanses 4.22063 0.50040 181.84030 8.435 1.01e-14 ***
## zsize 0.67978 0.18960 156.33377 3.585 0.000449 ***
## sector 1.73592 0.41738 134.29183 4.159 5.67e-05 ***
## cses:zsize 0.01961 0.12264 154.94416 0.160 0.873183
## cses:sector -1.00330 0.25288 147.66573 -3.967 0.000113 ***
## female:meanses -0.03191 0.48383 165.52237 -0.066 0.947498
## female:zsize -0.42547 0.18382 132.57553 -2.315 0.022174 *
## female:sector -0.30052 0.42847 142.90967 -0.701 0.484214
## minority:meanses -0.77917 0.53918 120.09086 -1.445 0.151030
## minority:zsize 0.11075 0.24035 141.49886 0.461 0.645657
## minority:sector 2.11884 0.54302 132.65245 3.902 0.000151 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itWe could relax convergence criteria on the gradient to make it look converged; for example, change tol = .002 to tol=.007.
Other things we can do is to Check.conv.singular; that is, check whether hat(T) (estimated matrix of \(\tau\)s) is a “good” covariance matrix. Bad covrariance matrices are non-singular and those where some parameters are on the boundary of the feasible space. For example, random effects variances equal to 0 or correlations between random effects equal to +/- 1.0)”
Other things that you can try:
o control= lmerControl(optimizer = “nloptwrap”,restart_edge = FALSE, boundary.tol = 0) o control= lmerControl(optimizer = “Nelder-Mead”) o check to make sure model is not misspecified o Try package brms and use Bayesian estimation
summary(model4.exta <- lmer(mathach ~ 1 + female + minority + cses + meanses + zsize + sector
+ cses*zsize + cses*sector + female*meanses + female*zsize + female*sector
+ minority*meanses + minority*zsize + minority*sector
+ (1 + female + minority + cses |id), data=hsb, REML=FALSE,
control= lmerControl(
check.conv.singular= .makeCC(action = "message", tol = 1e-4),
check.conv.hess = .makeCC(action = "warning", tol = 1e-6),
check.conv.grad = .makeCC("warning", tol = .002, relTol = NULL))))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + female + minority + cses + meanses + zsize + sector +
## cses * zsize + cses * sector + female * meanses + female *
## zsize + female * sector + minority * meanses + minority *
## zsize + minority * sector + (1 + female + minority + cses | id)
## Data: hsb
## Control: lmerControl(check.conv.singular = .makeCC(action = "message",
## tol = 1e-04), check.conv.hess = .makeCC(action = "warning",
## tol = 1e-06), check.conv.grad = .makeCC("warning", tol = 0.002,
## relTol = NULL))
##
## AIC BIC logLik deviance df.resid
## 46272.8 46451.7 -23110.4 46220.8 7159
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2037 -0.7149 0.0334 0.7609 2.9991
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 2.2268 1.4923
## female 0.7022 0.8380 -0.76
## minority 0.8756 0.9358 -0.12 0.24
## cses 0.1446 0.3803 0.26 -0.32 -0.35
## Residual 35.3183 5.9429
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.02587 0.46923 139.24621 25.629 < 2e-16 ***
## female -0.34012 0.46377 123.25042 -0.733 0.464723
## minority -4.20661 0.67149 157.13987 -6.265 3.43e-09 ***
## cses 2.28237 0.31118 148.15982 7.335 1.35e-11 ***
## meanses 4.22063 0.50040 181.84030 8.435 1.01e-14 ***
## zsize 0.67978 0.18960 156.33377 3.585 0.000449 ***
## sector 1.73592 0.41738 134.29183 4.159 5.67e-05 ***
## cses:zsize 0.01961 0.12264 154.94416 0.160 0.873183
## cses:sector -1.00330 0.25288 147.66573 -3.967 0.000113 ***
## female:meanses -0.03191 0.48383 165.52237 -0.066 0.947498
## female:zsize -0.42547 0.18382 132.57553 -2.315 0.022174 *
## female:sector -0.30052 0.42847 142.90967 -0.701 0.484214
## minority:meanses -0.77917 0.53918 120.09086 -1.445 0.151030
## minority:zsize 0.11075 0.24035 141.49886 0.461 0.645657
## minority:sector 2.11884 0.54302 132.65245 3.902 0.000151 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itTo me this model looks fine, but changing criteria to make things look good if not good practice.
relgrad <- with(model4@optinfo$derivs, solve(Hessian, gradient))
max(abs(relgrad)) ## [1] 0.0001473563isSingular(model4)## [1] FALSEThe wald statistic for a parameter, for example \(H_{o}: \gamma_{kl}=\gamma_{kl}^{*}\) where \(\gamma_{kl}^{*}\) is a fix values (usually 0) equals
\[ z = \frac{\gamma_{ik} - \gamma_{kl}^{*}}{\hat{se}} \sim N(0,1)\] But usually, \[ X^2 = z^2 \sim \chi^2_1 \] which is what is computed below.
# Put summary into object
s4 <- summary(model4)
# Get the cofficients from the summary
s4 <- as.data.frame(s4[10])
# Give them names
names(s4)<- c("Estimate","StdError","df t","t","Pr(>|t|)")
# Compute wald statistics for each
s4$df.Wald <- rep(1,nrow(s4))
s4$Wald4 <- s4$t**2
# Nicer output
options(scipen = 999)
print(s4, type = "html", digits=3)## Estimate StdError df t t
## (Intercept) 12.0259 0.469 139 25.6288
## female -0.3401 0.464 123 -0.7334
## minority -4.2066 0.671 157 -6.2646
## cses 2.2824 0.311 148 7.3346
## meanses 4.2206 0.500 182 8.4345
## zsize 0.6798 0.190 156 3.5854
## sector 1.7359 0.417 134 4.1591
## cses:zsize 0.0196 0.123 155 0.1599
## cses:sector -1.0033 0.253 148 -3.9674
## female:meanses -0.0319 0.484 166 -0.0659
## female:zsize -0.4255 0.184 133 -2.3146
## female:sector -0.3005 0.428 143 -0.7014
## minority:meanses -0.7792 0.539 120 -1.4451
## minority:zsize 0.1107 0.240 141 0.4608
## minority:sector 2.1188 0.543 133 3.9019
## Pr(>|t|)
## (Intercept) 0.00000000000000000000000000000000000000000000000000000142
## female 0.46472260894266570474542277224827557802200317382812500000
## minority 0.00000000342875657637011156848860782275778547045774757862
## cses 0.00000000001348180100413990077212050833210810196760576218
## meanses 0.00000000000001007720518222340721570598476297675460955361
## zsize 0.00044942501349583621750305439945805119350552558898925781
## sector 0.00005665999029530202705726149581266781751764938235282898
## cses:zsize 0.87318306263939660105677376122912392020225524902343750000
## cses:sector 0.00011296674286750000627300993816248819712200202047824860
## female:meanses 0.94749751697782957915450197106110863387584686279296875000
## female:zsize 0.02217352941763535076336033569077699212357401847839355469
## female:sector 0.48421402319891826415698687924304977059364318847656250000
## minority:meanses 0.15102953653566428271481925094121834263205528259277343750
## minority:zsize 0.64565698529181547726807366416323930025100708007812500000
## minority:sector 0.00015104602842428818177712701587012134041287936270236969
## df.Wald Wald4
## (Intercept) 1 656.83521
## female 1 0.53784
## minority 1 39.24523
## cses 1 53.79627
## meanses 1 71.14103
## zsize 1 12.85503
## sector 1 17.29775
## cses:zsize 1 0.02556
## cses:sector 1 15.74044
## female:meanses 1 0.00435
## female:zsize 1 5.35723
## female:sector 1 0.49192
## minority:meanses 1 2.08836
## minority:zsize 1 0.21233
## minority:sector 1 15.22499options(scipen = 0)
# For Carolyn's lecture notes #
# xtable(s4,digits=3)s4$upper <- s4$Estimate - qnorm(.025)*s4$StdError
s4$lower <- s4$Estimate - qnorm(.975)*s4$StdError
options(scipen = 999) # turn of scientic notation
round(s4[,8:9],digits=2)## upper lower
## (Intercept) 12.95 11.11
## female 0.57 -1.25
## minority -2.89 -5.52
## cses 2.89 1.67
## meanses 5.20 3.24
## zsize 1.05 0.31
## sector 2.55 0.92
## cses:zsize 0.26 -0.22
## cses:sector -0.51 -1.50
## female:meanses 0.92 -0.98
## female:zsize -0.07 -0.79
## female:sector 0.54 -1.14
## minority:meanses 0.28 -1.84
## minority:zsize 0.58 -0.36
## minority:sector 3.18 1.05options(scipen = 0) # turn it back one
# For LaTex
# xtable(s4[,8:9],digits=2)summary(model5 <- lmer(mathach ~ 1 + female + minority + cses + meanses + zsize + sector
+ cses*zsize + cses*sector + female*meanses + female*zsize + female*sector
+ minority*meanses + minority*zsize + minority*sector
+ (1 + female + cses |id), data=hsb, REML=FALSE) )## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + female + minority + cses + meanses + zsize + sector +
## cses * zsize + cses * sector + female * meanses + female *
## zsize + female * sector + minority * meanses + minority *
## zsize + minority * sector + (1 + female + cses | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46268.5 46419.8 -23112.2 46224.5 7163
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1742 -0.7150 0.0297 0.7602 2.9993
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 2.2491 1.4997
## female 0.6621 0.8137 -0.69
## cses 0.1436 0.3789 0.22 -0.49
## Residual 35.4175 5.9513
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.00065 0.47017 166.45639 25.524 < 2e-16 ***
## female -0.33177 0.46182 123.22654 -0.718 0.473866
## minority -4.13739 0.61408 2714.13003 -6.738 1.96e-11 ***
## cses 2.27126 0.31098 149.76016 7.304 1.54e-11 ***
## meanses 4.20243 0.49905 206.80067 8.421 6.23e-15 ***
## zsize 0.69898 0.18989 188.84073 3.681 0.000303 ***
## sector 1.74109 0.41675 152.12986 4.178 4.94e-05 ***
## cses:zsize 0.02701 0.12246 156.09767 0.221 0.825712
## cses:sector -0.99632 0.25269 149.68479 -3.943 0.000123 ***
## female:meanses -0.02628 0.47965 167.01772 -0.055 0.956372
## female:zsize -0.42767 0.18305 132.50884 -2.336 0.020976 *
## female:sector -0.28671 0.42591 143.17928 -0.673 0.501926
## minority:meanses -0.83179 0.48772 1890.00909 -1.705 0.088273 .
## minority:zsize 0.09301 0.21930 2043.69480 0.424 0.671511
## minority:sector 2.10585 0.49283 3073.13947 4.273 1.99e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itSummarize what we’ve fit so far
screenreg(list(model1,model2,model3,model4,model5),
custom.model.names=c("Model 1", "Model 2","Model 3","Model 4", "Model 5") )##
## =======================================================================================================
## Model 1 Model 2 Model 3 Model 4 Model 5
## -------------------------------------------------------------------------------------------------------
## (Intercept) 12.40 *** 12.40 *** 12.02 *** 12.03 *** 12.00 ***
## (0.36) (0.36) (0.43) (0.47) (0.47)
## cses 1.91 *** 1.39 *** 2.26 *** 2.28 *** 2.27 ***
## (0.11) (0.22) (0.31) (0.31) (0.31)
## female -1.24 *** -1.24 *** -0.36 -0.34 -0.33
## (0.16) (0.16) (0.42) (0.46) (0.46)
## minority -2.89 *** -2.87 *** -4.09 *** -4.21 *** -4.14 ***
## (0.20) (0.20) (0.62) (0.67) (0.61)
## meanses 3.96 *** 3.97 *** 4.23 *** 4.22 *** 4.20 ***
## (0.33) (0.33) (0.46) (0.50) (0.50)
## zsize 0.42 ** 0.42 ** 0.69 *** 0.68 *** 0.70 ***
## (0.14) (0.14) (0.18) (0.19) (0.19)
## sector 2.12 *** 2.11 *** 1.74 *** 1.74 *** 1.74 ***
## (0.30) (0.30) (0.38) (0.42) (0.42)
## cses:zsize 0.29 ** 0.03 0.02 0.03
## (0.10) (0.12) (0.12) (0.12)
## cses:sector -1.00 *** -1.00 *** -1.00 ***
## (0.25) (0.25) (0.25)
## female:meanses -0.05 -0.03 -0.03
## (0.44) (0.48) (0.48)
## female:zsize -0.42 * -0.43 * -0.43 *
## (0.17) (0.18) (0.18)
## female:sector -0.24 -0.30 -0.29
## (0.39) (0.43) (0.43)
## minority:meanses -0.80 -0.78 -0.83
## (0.49) (0.54) (0.49)
## minority:zsize 0.07 0.11 0.09
## (0.22) (0.24) (0.22)
## minority:sector 2.07 *** 2.12 *** 2.11 ***
## (0.49) (0.54) (0.49)
## -------------------------------------------------------------------------------------------------------
## AIC 46307.34 46301.58 46268.70 46272.84 46268.49
## BIC 46369.26 46370.37 46399.42 46451.72 46419.84
## Log Likelihood -23144.67 -23140.79 -23115.35 -23110.42 -23112.24
## Num. obs. 7185 7185 7185 7185 7185
## Num. groups: id 160 160 160 160 160
## Var: id (Intercept) 1.61 1.61 1.63 2.23 2.25
## Var: Residual 35.89 35.85 35.52 35.32 35.42
## Var: id cses 0.16 0.14 0.14
## Cov: id (Intercept) cses 0.04 0.15 0.13
## Var: id female 0.70 0.66
## Var: id minority 0.88
## Cov: id (Intercept) female -0.95 -0.84
## Cov: id (Intercept) minority -0.17
## Cov: id female minority 0.19
## Cov: id female cses -0.10 -0.15
## Cov: id minority cses -0.13
## =======================================================================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05# Easy way to get likelihood ratio test for nested models
anova(model1,model2)## Data: hsb
## Models:
## model1: mathach ~ 1 + cses + female + minority + meanses + zsize + sector +
## model1: (1 | id)
## model2: mathach ~ 1 + cses + female + minority + meanses + zsize + sector +
## model2: cses * zsize + (1 | id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model1 9 46307 46369 -23145 46289
## model2 10 46302 46370 -23141 46282 7.7638 1 0.00533 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This is a test whether the interaction is significant or not.
This will use the function “contrast’ from the”All_functions.txt” source file.
The model we will work with here
cmodel <- lmer(mathach ~ 1 + cses + female + minority + meanses + zsize + sector
+ cses*zsize + cses*sector + female*meanses + female*zsize
+ female*sector + minority*meanses + minority*zsize
+ minority*sector
+ (1 + cses + female | id), data=hsb, REML=FALSE)
summary(cmodel)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + female + minority + meanses + zsize + sector +
## cses * zsize + cses * sector + female * meanses + female *
## zsize + female * sector + minority * meanses + minority *
## zsize + minority * sector + (1 + cses + female | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46268.5 46419.8 -23112.2 46224.5 7163
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1742 -0.7150 0.0297 0.7602 2.9993
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 2.2497 1.4999
## cses 0.1437 0.3791 0.22
## female 0.6628 0.8141 -0.69 -0.49
## Residual 35.4173 5.9512
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.00063 0.47021 166.42528 25.522 < 2e-16 ***
## cses 2.27127 0.31099 149.77952 7.303 1.54e-11 ***
## female -0.33175 0.46186 123.20631 -0.718 0.473932
## minority -4.13740 0.61408 2714.18953 -6.738 1.96e-11 ***
## meanses 4.20244 0.49909 206.75627 8.420 6.26e-15 ***
## zsize 0.69899 0.18991 188.80267 3.681 0.000303 ***
## sector 1.74107 0.41679 152.10249 4.177 4.95e-05 ***
## cses:zsize 0.02701 0.12246 156.11671 0.221 0.825740
## cses:sector -0.99632 0.25270 149.70315 -3.943 0.000123 ***
## female:meanses -0.02626 0.47969 166.98281 -0.055 0.956409
## female:zsize -0.42768 0.18307 132.48531 -2.336 0.020984 *
## female:sector -0.28671 0.42594 143.15476 -0.673 0.501952
## minority:meanses -0.83175 0.48772 1889.99079 -1.705 0.088288 .
## minority:zsize 0.09301 0.21930 2043.69457 0.424 0.671509
## minority:sector 2.10585 0.49283 3073.21670 4.273 1.99e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itYou need to create a matrix where the number of columns equals the number of fixed effect. A fast way to get number of coefficents is
length(fixef(cmodel)) ## [1] 15And if you want the names either check output summary table of fixed effects, or
rownames(summary(cmodel)$coefficients)## [1] "(Intercept)" "cses" "female" "minority"
## [5] "meanses" "zsize" "sector" "cses:zsize"
## [9] "cses:sector" "female:meanses" "female:zsize" "female:sector"
## [13] "minority:meanses" "minority:zsize" "minority:sector"The number of rows equals the number of tests you want to perform (simultaneously).
We will test the \(\gamma_{kl}\) for the interaction between minority and sector. The contrast matrix L for this
L1 <- matrix(0,nrow=1,ncol=15)
L1[1,15] <- 1
# This give us
L1## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [,15]
## [1,] 1# And the test
contrast(cmodel,L1)## estimate F num df den df guess p-value Wald
## 2.105851e+00 1.825804e+01 1.000000e+00 3.073217e+03 1.987680e-05 1.825804e+01
## df p-value
## 1.000000e+00 1.929096e-05# -- for better looking output
round(contrast(cmodel,L1), digits=2)## estimate F num df den df guess p-value Wald
## 2.11 18.26 1.00 3073.22 0.00 18.26
## df p-value
## 1.00 0.00L <- matrix(0,nrow=5,ncol=15)
L[5,14] <- 1
L[4,13] <- 1
L[3,11] <- 1
L[2,10] <- 1
L[1, 8] <- 1
L## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 0 0 0 0 0 0 0 1 0 0 0 0 0 0
## [2,] 0 0 0 0 0 0 0 0 0 1 0 0 0 0
## [3,] 0 0 0 0 0 0 0 0 0 0 1 0 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 0 1 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0 0 1
## [,15]
## [1,] 0
## [2,] 0
## [3,] 0
## [4,] 0
## [5,] 0round(contrast(cmodel,L),digits=2)## estimate F num df den df guess p-value Wald
## 0.03 -0.03 -0.43 -0.83 0.09 1.73
## df p-value <NA> <NA> <NA> <NA>
## 5.00 156.12 0.13 8.66 5.00 0.12This can take longer depending on method. We did the wald based above and these were pretty fast. The bootstrap CIs an take a very long time
boot.ci <- confint(model1, method="boot", seed=1125, nsim=1000, level=0.99)
round(boot.ci,digits=4)
# for LaTex code
#xtable(boot.ci,digits=4)An alternative are Profile likelihood confidence intervals, which takes less time to compute than bootstrap.
profile.ci<- confint(model1, nsim=1000, level=0.99)## Computing profile confidence intervals ...round(profile.ci,digits=4)## 0.5 % 99.5 %
## .sig01 1.0069 1.5791
## .sigma 5.8628 6.1233
## (Intercept) 11.4653 13.3266
## cses 1.6334 2.1920
## female -1.6490 -0.8309
## minority -3.4051 -2.3696
## meanses 3.1006 4.8300
## zsize 0.0639 0.7697
## sector 1.3430 2.8946This is easy
anova(model1,model2)## Data: hsb
## Models:
## model1: mathach ~ 1 + cses + female + minority + meanses + zsize + sector +
## model1: (1 | id)
## model2: mathach ~ 1 + cses + female + minority + meanses + zsize + sector +
## model2: cses * zsize + (1 | id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model1 9 46307 46369 -23145 46289
## model2 10 46302 46370 -23141 46282 7.7638 1 0.00533 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1For these we’ll use the function “robust” from the source file.
r2 <- robust(model2, hsb$mathach, hsb$id,"between/within")
round(r2, digits=4)## Fixed Est. between/within Model se. Model t p-value Robust se
## (Intercept) 12.3967 156 0.3576 34.6672 0.0000 0.3725
## cses 1.3867 7021 0.2176 6.3727 0.0000 0.2459
## female -1.2422 7021 0.1585 -7.8378 0.0000 0.1726
## minority -2.8687 7021 0.2009 -14.2796 0.0000 0.2378
## meanses 3.9711 156 0.3324 11.9480 0.0000 0.3470
## zsize 0.4161 156 0.1355 3.0702 0.0025 0.1420
## sector 2.1139 156 0.2979 7.0966 0.0000 0.3068
## cses:zsize 0.2909 7021 0.1044 2.7871 0.0053 0.1254
## Robust t p-value
## (Intercept) 33.2755 0.0000
## cses 5.6391 0.0000
## female -7.1993 0.0000
## minority -12.0624 0.0000
## meanses 11.4437 0.0000
## zsize 2.9296 0.0039
## sector 6.8907 0.0000
## cses:zsize 2.3205 0.0203r3 <- robust(model3, hsb$mathach, hsb$id, "between/within")
round(r3, digits=4)## Fixed Est. between/within Model se. Model t p-value Robust se
## (Intercept) 12.0180 155 0.4327 27.7758 0.0000 0.4209
## cses 2.2618 7015 0.3121 7.2461 0.0000 0.3201
## female -0.3582 7015 0.4217 -0.8494 0.3957 0.4061
## minority -4.0871 7015 0.6183 -6.6105 0.0000 0.6494
## meanses 4.2344 155 0.4629 9.1474 0.0000 0.4938
## zsize 0.6932 155 0.1758 3.9429 0.0001 0.1785
## sector 1.7389 155 0.3807 4.5670 0.0000 0.3905
## cses:zsize 0.0321 7015 0.1229 0.2616 0.7937 0.1349
## cses:sector -0.9985 155 0.2538 -3.9349 0.0001 0.2496
## female:meanses -0.0543 7015 0.4433 -0.1226 0.9025 0.4240
## female:zsize -0.4161 7015 0.1677 -2.4815 0.0131 0.1511
## female:sector -0.2368 155 0.3924 -0.6034 0.5472 0.4078
## minority:meanses -0.8046 7015 0.4918 -1.6360 0.1019 0.4895
## minority:zsize 0.0696 7015 0.2213 0.3146 0.7531 0.2348
## minority:sector 2.0697 155 0.4949 4.1824 0.0000 0.5400
## Robust t p-value
## (Intercept) 28.5541 0.0000
## cses 7.0659 0.0000
## female -0.8820 0.3778
## minority -6.2938 0.0000
## meanses 8.5759 0.0000
## zsize 3.8829 0.0002
## sector 4.4526 0.0000
## cses:zsize 0.2383 0.8116
## cses:sector -4.0005 0.0001
## female:meanses -0.1281 0.8981
## female:zsize -2.7540 0.0059
## female:sector -0.5807 0.5623
## minority:meanses -1.6436 0.1003
## minority:zsize 0.2965 0.7668
## minority:sector 3.8329 0.0002# LaTex
# r3.data <- r3[1]$table # turn list into matrix
# xtable(r3.data,digits=2)Compare model and robust using model 3
r4 <- robust(model4, hsb$mathach,hsb$id,"between/within")
round(r4, digits=4)## Fixed Est. between/within Model se. Model t p-value Robust se
## (Intercept) 12.0259 155 0.4692 25.6288 0.0000 0.4269
## female -0.3401 7015 0.4638 -0.7334 0.4634 0.4137
## minority -4.2066 7015 0.6715 -6.2646 0.0000 0.6439
## cses 2.2824 7015 0.3112 7.3346 0.0000 0.3176
## meanses 4.2206 155 0.5004 8.4345 0.0000 0.4961
## zsize 0.6798 155 0.1896 3.5854 0.0005 0.1790
## sector 1.7359 155 0.4174 4.1591 0.0001 0.3978
## cses:zsize 0.0196 7015 0.1226 0.1599 0.8730 0.1342
## cses:sector -1.0033 155 0.2529 -3.9674 0.0001 0.2503
## female:meanses -0.0319 7015 0.4838 -0.0659 0.9474 0.4235
## female:zsize -0.4255 7015 0.1838 -2.3146 0.0207 0.1546
## female:sector -0.3005 155 0.4285 -0.7014 0.4841 0.4150
## minority:meanses -0.7792 7015 0.5392 -1.4451 0.1485 0.4933
## minority:zsize 0.1107 7015 0.2403 0.4608 0.6450 0.2335
## minority:sector 2.1188 155 0.5430 3.9019 0.0001 0.5398
## Robust t p-value
## (Intercept) 28.1682 0.0000
## female -0.8221 0.4111
## minority -6.5327 0.0000
## cses 7.1870 0.0000
## meanses 8.5082 0.0000
## zsize 3.7966 0.0002
## sector 4.3639 0.0000
## cses:zsize 0.1461 0.8838
## cses:sector -4.0083 0.0001
## female:meanses -0.0754 0.9399
## female:zsize -2.7522 0.0059
## female:sector -0.7241 0.4701
## minority:meanses -1.5795 0.1143
## minority:zsize 0.4744 0.6352
## minority:sector 3.9254 0.0001For these we use the bic.html function from the source file
bic.hlm(model1,hsb$id)## n groups total obs haroninc.mean.nj n random n fixed deviance bic.new
## id 160 7185 41.05874 2 7 46289.34 46361.65
## bic.harm bic.ngrps bic.ntot aic
## id 46361.02 46335.02 46369.26 46307.34bic.hlm(model2,hsb$id)## n groups total obs haroninc.mean.nj n random n fixed deviance bic.new
## id 160 7185 41.05874 2 8 46281.58 46362.76
## bic.harm bic.ngrps bic.ntot aic
## id 46362.05 46332.33 46370.37 46301.58bic.hlm(model3,hsb$id)## n groups total obs haroninc.mean.nj n random n fixed deviance bic.new
## id 160 7185 41.05874 4 15 46230.7 46384.2
## bic.harm bic.ngrps bic.ntot aic
## id 46382.85 46327.13 46399.42 46268.7bic.hlm(model4,hsb$id)## n groups total obs haroninc.mean.nj n random n fixed deviance bic.new
## id 160 7185 41.05874 11 15 46220.84 46409.87
## bic.harm bic.ngrps bic.ntot aic
## id 46408.52 46352.8 46451.72 46272.84bic.hlm(model5,hsb$id)## n groups total obs haroninc.mean.nj n random n fixed deviance bic.new
## id 160 7185 41.05874 7 15 46224.49 46393.21
## bic.harm bic.ngrps bic.ntot aic
## id 46391.86 46336.14 46419.84 46268.49# Null or empty HLM
summary(model.a <- lmer(mathach ~ 1 + (1 | id),data=hsb,REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 47121.8 47142.4 -23557.9 47115.8 7182
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.06262 -0.75365 0.02676 0.76070 2.74184
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 8.553 2.925
## Residual 39.148 6.257
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.6371 0.2436 157.6209 51.87 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Random intercept with cses are fixed effect
summary(model.b <- lmer(mathach ~ 1 + cses + (1 | id),data=hsb,REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46728.4 46755.9 -23360.2 46720.4 7181
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.09692 -0.73196 0.01945 0.75739 2.91423
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 8.612 2.935
## Residual 37.005 6.083
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.6494 0.2437 157.7181 51.90 <2e-16 ***
## cses 2.1912 0.1086 7023.0025 20.17 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cses 0.003hlmRsq(hsb, model.b, hsb$id)## harmonic.mean R1sq R2sq
## [1,] 41.05874 0.04491355 0.005671609# Random intercept with more micro level variables
summary(model.c <- lmer(mathach ~ 1 + cses + minority + (1 | id),data=hsb,REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + minority + (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46523.1 46557.5 -23256.5 46513.1 7180
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.13931 -0.72315 0.02938 0.75854 2.94870
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 6.648 2.578
## Residual 36.121 6.010
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 13.4973 0.2243 176.4027 60.18 <2e-16 ***
## cses 1.9372 0.1087 7051.0708 17.82 <2e-16 ***
## minority -3.0806 0.2105 6040.6784 -14.63 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses
## cses -0.038
## minority -0.257 0.160hlmRsq(hsb,model.c, hsb$id)## harmonic.mean R1sq R2sq
## [1,] 41.05874 0.1045467 0.2132329# Random intercept with micro and marco varialbles
summary(model.d <- model.d <- lmer(mathach ~ 1 + cses + minority + female + sector + meanses + zsize + (1 | id),data=hsb,REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + minority + female + sector + meanses + zsize +
## (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46307.3 46369.3 -23144.7 46289.3 7176
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2641 -0.7225 0.0410 0.7599 2.9191
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.61 1.269
## Residual 35.89 5.991
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.4000 0.3576 172.2298 34.672 < 2e-16 ***
## cses 1.9127 0.1084 7076.1669 17.644 < 2e-16 ***
## minority -2.8871 0.2009 3610.0411 -14.373 < 2e-16 ***
## female -1.2402 0.1586 5173.8657 -7.821 6.29e-15 ***
## sector 2.1179 0.2979 153.0838 7.109 4.16e-11 ***
## meanses 3.9634 0.3324 170.3640 11.924 < 2e-16 ***
## zsize 0.4154 0.1355 159.0317 3.065 0.00256 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses minrty female sector meanss
## cses -0.016
## minority -0.041 0.153
## female -0.245 0.050 0.016
## sector -0.653 -0.023 -0.148 -0.006
## meanses 0.122 0.041 0.256 0.048 -0.360
## zsize -0.837 -0.013 -0.089 0.019 0.437 -0.056hlmRsq(hsb,model.d, hsb$id)## harmonic.mean R1sq R2sq
## [1,] 41.05874 0.2148913 0.7403604# Random intercept with some cross-level (fixed) interactions
summary(model.e <- lmer(mathach ~ 1 + cses + minority + female + sector + meanses + zsize
+ cses*sector + female*zsize + minority*sector +(1 | id),data=hsb,REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + minority + female + sector + meanses + zsize +
## cses * sector + female * zsize + minority * sector + (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46258.7 46341.2 -23117.3 46234.7 7173
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1761 -0.7201 0.0375 0.7594 3.0097
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.668 1.291
## Residual 35.594 5.966
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.0251 0.3937 226.0948 30.541 < 2e-16 ***
## cses 2.3261 0.1457 7115.9808 15.966 < 2e-16 ***
## minority -3.7545 0.2894 3552.1880 -12.975 < 2e-16 ***
## female -0.5170 0.3276 4372.6568 -1.578 0.1146
## sector 1.7215 0.3119 168.9717 5.519 1.26e-07 ***
## meanses 3.8984 0.3367 168.2278 11.577 < 2e-16 ***
## zsize 0.7104 0.1620 285.9189 4.385 1.63e-05 ***
## cses:sector -1.0116 0.2174 7078.6332 -4.652 3.34e-06 ***
## female:zsize -0.3783 0.1543 6018.1876 -2.452 0.0142 *
## minority:sector 1.7430 0.3927 2817.0975 4.439 9.39e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses minrty female sector meanss zsize css:sc fml:zs
## cses -0.020
## minority -0.043 0.213
## female -0.456 0.032 0.022
## sector -0.577 0.016 0.090 -0.015
## meanses 0.100 0.051 0.231 0.048 -0.326
## zsize -0.854 -0.025 -0.132 0.463 0.316 -0.042
## cses:sector 0.008 -0.669 -0.141 -0.012 -0.027 -0.025 0.021
## female:zsiz 0.399 -0.007 -0.017 -0.876 0.014 -0.028 -0.520 0.000
## minrty:sctr 0.022 -0.153 -0.720 -0.019 -0.262 -0.077 0.114 0.146 0.020hlmRsq(hsb,model.e, hsb$id)## harmonic.mean R1sq R2sq
## [1,] 41.05874 0.219857 0.7350976# Random intercept with lots of variables and interactions
summary(model.f <- lmer(mathach ~ 1 + cses + minority + female + sector + meanses + zsize
+ cses*sector + female*zsize + minority*sector
+ cses*zsize + sector*female + meanses*female + zsize*minority + meanses*minority + (1 | id),data=hsb,REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + minority + female + sector + meanses + zsize +
## cses * sector + female * zsize + minority * sector + cses *
## zsize + sector * female + meanses * female + zsize * minority +
## meanses * minority + (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46265.4 46382.3 -23115.7 46231.4 7168
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2002 -0.7208 0.0322 0.7612 3.0112
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.631 1.277
## Residual 35.589 5.966
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.01759 0.43260 318.24502 27.780 < 2e-16 ***
## cses 2.25857 0.29858 7085.50794 7.564 4.39e-14 ***
## minority -4.07679 0.61783 3447.62061 -6.599 4.79e-11 ***
## female -0.35494 0.42178 6763.12561 -0.842 0.4001
## sector 1.73533 0.38064 308.27823 4.559 7.42e-06 ***
## meanses 4.25257 0.46236 410.21088 9.198 < 2e-16 ***
## zsize 0.69554 0.17578 359.63362 3.957 9.15e-05 ***
## cses:sector -0.99687 0.24275 7069.53949 -4.106 4.06e-05 ***
## female:zsize -0.41903 0.16771 6588.21823 -2.499 0.0125 *
## minority:sector 2.05946 0.49444 3522.35295 4.165 3.18e-05 ***
## cses:zsize 0.03136 0.11769 7094.85768 0.266 0.7899
## female:sector -0.23692 0.39245 2664.08631 -0.604 0.5461
## female:meanses -0.06222 0.44336 2820.12521 -0.140 0.8884
## minority:zsize 0.06516 0.22111 2517.80727 0.295 0.7682
## minority:meanses -0.79579 0.49139 2206.22376 -1.619 0.1055
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need ithlmRsq(hsb,model.f, hsb$id)## harmonic.mean R1sq R2sq
## [1,] 41.05874 0.2207253 0.7389459summary(model.bplus <- lmer(mathach ~ 1 + cses + meanses + (1 | id),data=hsb,REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + meanses + (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46573.8 46608.2 -23281.9 46563.8 7180
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1675 -0.7257 0.0176 0.7553 2.9445
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 2.647 1.627
## Residual 37.014 6.084
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.6616 0.1484 155.6417 85.31 <2e-16 ***
## cses 2.1912 0.1087 7022.4974 20.16 <2e-16 ***
## meanses 5.8656 0.3594 155.3036 16.32 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses
## cses 0.004
## meanses -0.004 0.000hlmRsq(hsb,model.bplus, hsb$id)## harmonic.mean R1sq R2sq
## [1,] 41.05874 0.1696187 0.6291169summary(model.extra <- lmer(mathach ~ 1 + cses + meanses + female + minority + pracad + sector +(1 | id),data=hsb,REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + meanses + female + minority + pracad + sector +
## (1 | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46300.7 46362.6 -23141.4 46282.7 7176
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2543 -0.7262 0.0348 0.7602 2.9177
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.513 1.230
## Residual 35.888 5.991
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.0039 0.3767 176.4713 31.863 < 2e-16 ***
## cses 1.9129 0.1084 7078.4420 17.648 < 2e-16 ***
## meanses 3.0660 0.4013 165.7240 7.640 1.65e-12 ***
## female -1.2331 0.1582 5042.9657 -7.792 7.94e-15 ***
## minority -2.8873 0.1999 3416.6539 -14.445 < 2e-16 ***
## pracad 3.2688 0.8052 156.3643 4.060 7.75e-05 ***
## sector 0.8826 0.3339 153.4119 2.643 0.00907 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses meanss female minrty pracad
## cses -0.017
## meanses 0.559 0.039
## female -0.243 0.050 0.023
## minority -0.048 0.153 0.249 0.015
## pracad -0.859 -0.009 -0.585 0.030 -0.071
## sector 0.297 -0.010 0.123 -0.031 -0.053 -0.618hlmRsq(hsb, model.extra, hsb$id)## harmonic.mean R1sq R2sq
## [1,] 41.05874 0.2169417 0.7505062# Random intercept and slope for cses
summary(model.bis<- lmer(mathach ~ 1 + cses + meanses + (1 + cses| id),data=hsb,REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + meanses + (1 + cses | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46567.2 46615.3 -23276.6 46553.2 7178
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1778 -0.7260 0.0158 0.7565 2.9403
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 2.6446 1.6262
## cses 0.6701 0.8186 -0.19
## Residual 36.7132 6.0591
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.6588 0.1482 155.6389 85.41 <2e-16 ***
## cses 2.1912 0.1276 156.3840 17.17 <2e-16 ***
## meanses 5.8959 0.3577 154.8624 16.48 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses
## cses -0.083
## meanses -0.004 0.004hlmRsq(hsb, model.bis, hsb$id)## harmonic.mean R1sq R2sq
## [1,] 41.05874 0.1697779 0.6290588summary(model.ris <- lmer(mathach ~ 1 + cses + meanses + female + minority + pracad + sector +(1 +cses| id),data=hsb,REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + meanses + female + minority + pracad + sector +
## (1 + cses | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46299.5 46375.1 -23138.7 46277.5 7174
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2563 -0.7220 0.0341 0.7566 2.9110
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 1.5207 1.2332
## cses 0.4848 0.6963 0.06
## Residual 35.6706 5.9725
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 11.9750 0.3768 176.4770 31.780 < 2e-16 ***
## cses 1.9183 0.1225 156.7347 15.657 < 2e-16 ***
## meanses 3.0487 0.4015 165.7906 7.594 2.15e-12 ***
## female -1.2223 0.1582 5048.8333 -7.728 1.31e-14 ***
## minority -2.8700 0.2004 3219.4296 -14.322 < 2e-16 ***
## pracad 3.2838 0.8054 156.3909 4.077 7.23e-05 ***
## sector 0.9076 0.3340 153.4150 2.717 0.00734 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses meanss female minrty pracad
## cses -0.009
## meanses 0.558 0.033
## female -0.243 0.045 0.023
## minority -0.048 0.134 0.250 0.015
## pracad -0.859 -0.008 -0.585 0.030 -0.072
## sector 0.297 -0.008 0.123 -0.031 -0.053 -0.618hlmRsq(hsb, model.ris, hsb$id)## harmonic.mean R1sq R2sq
## [1,] 41.05874 0.2169097 0.7497793summary(model.ris2 <- lmer(mathach ~ 1 + cses + meanses + minority + female + minority + pracad + sector +(1+cses + minority | id),data=hsb,REML=FALSE))## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mathach ~ 1 + cses + meanses + minority + female + minority +
## pracad + sector + (1 + cses + minority | id)
## Data: hsb
##
## AIC BIC logLik deviance df.resid
## 46292.5 46388.8 -23132.3 46264.5 7171
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2765 -0.7210 0.0317 0.7604 2.9219
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id (Intercept) 1.3833 1.1761
## cses 0.4069 0.6379 0.37
## minority 1.4721 1.2133 -0.04 -0.92
## Residual 35.4840 5.9568
## Number of obs: 7185, groups: id, 160
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 11.9953 0.3782 177.0312 31.718 < 2e-16 ***
## cses 1.8927 0.1203 152.3099 15.737 < 2e-16 ***
## meanses 3.1944 0.4124 166.6987 7.745 8.82e-13 ***
## minority -2.9835 0.2318 130.6876 -12.869 < 2e-16 ***
## female -1.2296 0.1582 5170.4769 -7.775 9.06e-15 ***
## pracad 3.2474 0.8136 158.0909 3.991 0.0001 ***
## sector 0.8111 0.3338 149.7558 2.430 0.0163 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cses meanss minrty female pracad
## cses 0.021
## meanses 0.559 0.024
## minority -0.050 -0.055 0.210
## female -0.248 0.046 0.017 0.013
## pracad -0.865 -0.007 -0.592 -0.058 0.037
## sector 0.309 -0.001 0.133 -0.036 -0.036 -0.623hlmRsq(hsb, model.ris2, hsb$id)## harmonic.mean R1sq R2sq
## [1,] 41.05874 0.2192013 0.7387274Computing p-value for mixture of Chi-squares.
This one is for \(H_o: \tau^2_1 = \tau_{01} = 0\) from page 97 of lecture notes (i.e., testing null model (random intercept) vs random intercept & slope for female):
(p <- .5*(pchisq(5.40,1, lower.tail=FALSE) + pchisq(5.40,2, lower.tail=FALSE))) ## [1] 0.04367113This is the test reported on page 101 of lecture notes.
(p <- .5*(pchisq(1.55, 3, lower.tail=FALSE) + pchisq(1.55, 4, lower.tail=FALSE)) ) ## [1] 0.7442643