##2020 isSingular fix for exam lmer # to appease lmer these days (these models have run fine for a decade) need to insert a control statement # , control = lmerControl(check.conv.singular = .makeCC(action = "message", tol = 1e-6)) # into the lmer command for both unconditional and conditional models for exam data # binary level 1 predictor, sexM, causes a correlation near -1 for level 1 random effects. # lmer used to not complain about that, now need to placate. # still get the same resultsas in the past R version 3.6.3 (2020-02-29) -- "Holding the Windsock" Copyright (C) 2020 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 (64-bit) ######unconditional model > mixed = read.table("http://web.stanford.edu/~rag/stat209/mExam4347", header = T) > dim(mixed) [1] 2026 6 > library(lme4) Loading required package: Matrix > ggaplmer = lmer(normexam ~ sex + (sex|school), data = mixed) # see models page boundary (singular) fit: see ?isSingular ## this used to work without complaint, we will fix below > summary(ggaplmer) # you can still get the same estimates as in past years Linear mixed model fit by REML ['lmerMod'] Formula: normexam ~ sex + (sex | school) Data: mixed REML criterion at convergence: 5423.9 Scaled residuals: Min 1Q Median 3Q Max -3.8554 -0.6661 0.0006 0.6873 2.7088 Random effects: Groups Name Variance Std.Dev. Corr school (Intercept) 0.196775 0.44359 sexM 0.001263 0.03553 -1.00 Residual 0.814181 0.90232 Number of obs: 2026, groups: school, 33 Fixed effects: Estimate Std. Error t value (Intercept) 0.02467 0.08335 0.296 sexM -0.25886 0.04165 -6.214 Correlation of Fixed Effects: (Intr) sexM -0.393 convergence code: 0 boundary (singular) fit: see ?isSingular > confint(ggaplmer, oldNames = FALSE)#but CI won't work, can bootstrap below ## to see more > ?isSingular #starting httpd help server ... done #> rePCA(ggaplmer) ############2020 version, control statement to appease R ########## unconditional model > ggaplmer = lmer(normexam ~ sex + (sex|school), control = lmerControl(check.conv.singular = .makeCC(action = "message", tol = 1e-6)) , data = mixed) ## we get no complaints > summary(ggaplmer) Linear mixed model fit by REML ['lmerMod'] Formula: normexam ~ sex + (sex | school) Data: mixed Control: lmerControl(check.conv.singular = .makeCC(action = "message", tol = 1e-06)) REML criterion at convergence: 5423.9 Scaled residuals: Min 1Q Median 3Q Max -3.8554 -0.6661 0.0006 0.6873 2.7088 Random effects: Groups Name Variance Std.Dev. Corr school (Intercept) 0.196775 0.44359 sexM 0.001263 0.03553 -1.00 Residual 0.814181 0.90232 Number of obs: 2026, groups: school, 33 Fixed effects: Estimate Std. Error t value (Intercept) 0.02467 0.08335 0.296 sexM -0.25886 0.04165 -6.214 Correlation of Fixed Effects: (Intr) sexM -0.393 #still have to do bootstrap version of CI > confint(ggaplmer, method = "boot", nsim = 4000, boot.type = "perc") Computing bootstrap confidence intervals ... 2170 message(s): boundary (singular) fit: see ?isSingular 2807 warning(s): Model failed to converge with max|grad| = 0.00200353 (tol = 0.002, component 1) (and others) 2.5 % 97.5 % .sig01 0.320914429 0.5662437 .sig02 -1.000000000 0.9999991 .sig03 0.002100234 0.1251865 .sigma 0.873858861 0.9302752 (Intercept) -0.139140852 0.1861349 sexM -0.341577645 -0.1766211 ########## conditional model ## the old way still works, but R complains > ggaplmer2 = lmer(normexam ~ sex*schavg + (sex|school), data = mixed)# do the formal model boundary (singular) fit: see ?isSingular > summary(ggaplmer2) Linear mixed model fit by REML ['lmerMod'] Formula: normexam ~ sex * schavg + (sex | school) Data: mixed REML criterion at convergence: 5406.3 Scaled residuals: Min 1Q Median 3Q Max -3.8320 -0.6769 -0.0049 0.6958 2.7018 Random effects: Groups Name Variance Std.Dev. Corr school (Intercept) 0.098379 0.31365 sexM 0.003121 0.05586 -1.00 Residual 0.815375 0.90298 Number of obs: 2026, groups: school, 33 Fixed effects: Estimate Std. Error t value (Intercept) 0.062091 0.063254 0.982 sexM -0.261676 0.042333 -6.181 schavg 0.958839 0.204804 4.682 sexM:schavg -0.002882 0.141970 -0.020 Correlation of Fixed Effects: (Intr) sexM schavg sexM -0.536 schavg 0.099 -0.027 sexM:schavg -0.029 0.023 -0.563 convergence code: 0 boundary (singular) fit: see ?isSingular ## the 2020 version > ggaplmer2 = lmer(normexam ~ sex*schavg + (sex|school), , control = lmerControl(check.conv.singular = .makeCC(action = "message", tol = 1e-6)), data = mixed) > summary(ggaplmer2) Linear mixed model fit by REML ['lmerMod'] Formula: normexam ~ sex * schavg + (sex | school) Data: mixed Control: lmerControl(check.conv.singular = .makeCC(action = "message", tol = 1e-06)) REML criterion at convergence: 5406.3 Scaled residuals: Min 1Q Median 3Q Max -3.8320 -0.6769 -0.0049 0.6958 2.7018 Random effects: Groups Name Variance Std.Dev. Corr school (Intercept) 0.098379 0.31365 sexM 0.003121 0.05586 -1.00 Residual 0.815375 0.90298 Number of obs: 2026, groups: school, 33 Fixed effects: Estimate Std. Error t value (Intercept) 0.062091 0.063254 0.982 sexM -0.261676 0.042333 -6.181 schavg 0.958839 0.204804 4.682 sexM:schavg -0.002882 0.141970 -0.020 Correlation of Fixed Effects: (Intr) sexM schavg sexM -0.536 schavg 0.099 -0.027 sexM:schavg -0.029 0.023 -0.563 # profile likelihood will give you CI, or bootstrap > confint(ggaplmer2, oldNames = FALSE) Computing profile confidence intervals ... 2.5 % 97.5 % sd_(Intercept)|school 0.21103180 0.4228471 cor_sexM.(Intercept)|school -1.00000000 1.0000000 sd_sexM|school 0.00000000 0.1246263 sigma 0.87526384 0.9313583 (Intercept) -0.05773692 0.1825400 sexM -0.34432594 -0.1787847 schavg 0.57250126 1.3520952 sexM:schavg -0.28174729 0.2736615