1stopifnot(require("testthat"), require("lme4")) 2 3## use old (<=3.5.2) sample() algorithm if necessary 4if ("sample.kind" %in% names(formals(RNGkind))) { 5 suppressWarnings(RNGkind("Mersenne-Twister", "Inversion", "Rounding")) 6} 7 8context("fitting lmer models") 9## is "Nelder_Mead" default optimizer? -- no longer 10(isNM <- formals(lmerControl)$optimizer == "Nelder_Mead") 11 12test_that("lmer", { 13 set.seed(101) 14 d <- data.frame(z=rnorm(200), 15 f=factor(sample(1:10,200, replace=TRUE))) 16 17 ## Using 'method=*' defunct in 2019-05 (after 6 years of deprecation) 18 ## expect_warning(lmer(z~ 1|f, d, method="abc"),"Use the REML argument") 19 ## expect_warning(lmer(z~ 1|f, d, method="Laplace"),"Use the REML argument") 20 ##sp No '...' anymore 21 ##sp expect_warning(lmer(z~ 1|f, d, sparseX=TRUE),"has no effect at present") 22 expect_error(lmer(z~ 1|f, ddd), "bad 'data': object 'ddd' not found") 23 expect_error(lmer(z~ 1|f), "object 'z' not found") 24 expect_error(lmer(z~ 1|f, d[,1:1000]), "bad 'data': undefined columns selected") 25 expect_is(fm1 <- lmer(Yield ~ 1|Batch, Dyestuff), "lmerMod") 26 expect_is(fm1_noCD <- update(fm1,control=lmerControl(calc.derivs=FALSE)), 27 "lmerMod") 28 expect_equal(VarCorr(fm1),VarCorr(fm1_noCD)) 29 ## backward compatibility version {for optimizer="Nelder-Mead" only}: 30 if(isNM) expect_is(fm1.old <- update(fm1,control=lmerControl(use.last.params=TRUE)), 31 "lmerMod") 32 expect_is(fm1@resp, "lmerResp") 33 expect_is(fm1@pp, "merPredD") 34 expect_that(fe1 <- fixef(fm1), is_equivalent_to(1527.5)) 35 expect_that(VarCorr(fm1)[[1]][1,1], ## "bobyqa" : 1764.050060 36 equals(1764.0375195, tolerance = 1e-5)) 37 ## back-compatibility ... 38 if(isNM) expect_that(VarCorr(fm1.old)[[1]][1,1], equals(1764.0726543)) 39 40 expect_that(isREML(fm1), equals(TRUE)) 41 expect_is(REMLfun <- as.function(fm1), "function") 42 expect_that(REMLfun(1), equals(319.792389042002)) 43 expect_that(REMLfun(0), equals(326.023232155879)) 44 expect_that(family(fm1), equals(gaussian())) 45 expect_that(isREML(fm1ML <- refitML(fm1)), equals(FALSE)) 46 expect_that(REMLcrit(fm1), equals(319.654276842342)) 47 expect_that(deviance(fm1ML), equals(327.327059881135)) 48 ## "bobyqa": 49.51009984775 49 expect_that(sigma(fm1), equals(49.5101272946856, tolerance=1e-6)) 50 if(isNM) expect_that(sigma(fm1.old), equals(49.5100503990048)) 51 expect_that(sigma(fm1ML), equals(49.5100999308089)) 52 expect_that(extractAIC(fm1), equals(c(3, 333.327059881135))) 53 expect_that(extractAIC(fm1ML), equals(c(3, 333.327059881135))) 54 ## "bobyqa": 375.71667627943 55 expect_that(vcov(fm1) [1,1], equals(375.714676744, tolerance=1e-5)) 56 if(isNM) expect_that(vcov(fm1.old)[1,1], equals(375.72027872986)) 57 expect_that(vcov(fm1ML) [1,1], equals(313.09721874266, tolerance=1e-7)) 58 # was 313.0972246957 59 expect_is(fm2 <- refit(fm1, Dyestuff2$Yield), "lmerMod") 60 expect_that(fixef(fm2), is_equivalent_to(5.6656)) 61 expect_that(VarCorr(fm2)[[1]][1,1], is_equivalent_to(0)) 62 expect_that(getME(fm2, "theta"), is_equivalent_to(0)) 63 expect_that(X <- getME(fm1, "X"), is_equivalent_to(array(1, c(1, 30)))) 64 expect_is(Zt <- getME(fm1, "Zt"), "dgCMatrix") 65 expect_that(dim(Zt), equals(c(6L, 30L))) 66 expect_that(Zt@x, equals(rep.int(1, 30L))) 67 expect_equal(dimnames(Zt), 68 list(levels(Dyestuff$Batch), 69 rownames(Dyestuff))) 70 ## "bobyqa": 0.8483237982 71 expect_that(theta <- getME(fm1, "theta"), equals(0.84832031, tolerance=6e-6, check.attributes=FALSE)) 72 if(isNM) expect_that(getME(fm1.old, "theta"), is_equivalent_to(0.848330078)) 73 expect_is(Lambdat <- getME(fm1, "Lambdat"), "dgCMatrix") 74 expect_that(as(Lambdat, "matrix"), is_equivalent_to(diag(theta, 6L, 6L))) 75 expect_is(fm3 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject), sleepstudy), 76 "lmerMod") 77 expect_that(getME(fm3,"n_rtrms"), equals(2L)) 78 expect_that(getME(fm3,"n_rfacs"), equals(1L)) 79 80 expect_equal(getME(fm3, "lower"), c(`Subject.(Intercept)` = 0, Subject.Days = 0)) 81 82 expect_error(fm4 <- lmer(Reaction ~ Days + (1|Subject), 83 subset(sleepstudy,Subject==levels(Subject)[1])), "must have > 1") 84 expect_warning(fm4 <- lFormula(Reaction ~ Days + (1|Subject), 85 subset(sleepstudy,Subject==levels(Subject)[1]), 86 control=lmerControl(check.nlev.gtr.1="warning")), "must have > 1") 87 expect_warning(fm4 <- lmer(Reaction ~ Days + (1|Subject), 88 subset(sleepstudy,Subject %in% levels(Subject)[1:4]), 89 control=lmerControl(check.nlev.gtreq.5="warning")), 90 "< 5 sampled levels") 91 sstudy9 <- subset(sleepstudy, Days == 1 | Days == 9) 92 expect_error(lmer(Reaction ~ 1 + Days + (1 + Days | Subject), 93 data = sleepstudy, subset = (Days == 1 | Days == 9)), 94 "number of observations \\(=36\\) <= number of random effects \\(=36\\)") 95 expect_error(lFormula(Reaction ~ 1 + Days + (1 + Days | Subject), 96 data = sleepstudy, subset = (Days == 1 | Days == 9)), 97 "number of observations \\(=36\\) <= number of random effects \\(=36\\)") 98 ## with most recent Matrix (1.1-1), should *not* flag this 99 ## for insufficient rank 100 dat <- readRDS(system.file("testdata", "rankMatrix.rds", package="lme4")) 101 expect_is(lFormula(y ~ (1|sample) + (1|day) + (1|day:sample) + 102 (1|operator) + (1|day:operator) + (1|sample:operator) + 103 (1|day:sample:operator), 104 data = dat, 105 control = lmerControl(check.nobs.vs.rankZ = "stop")), 106 "list") 107 ## check scale 108 ss <- within(sleepstudy, Days <- Days*1e6) 109 expect_warning(lmer(Reaction ~ Days + (1|Subject), data=ss), 110 "predictor variables are on very different scales") 111 112 ## Promote warning to error so that warnings or errors will stop the test: 113 options(warn=2) 114 expect_is(lmer(Yield ~ 1|Batch, Dyestuff, REML=TRUE), "lmerMod") 115 expect_is(lmer(Yield ~ 1|Batch, Dyestuff, start=NULL), "lmerMod") 116 expect_is(lmer(Yield ~ 1|Batch, Dyestuff, verbose=0L), "lmerMod") 117 expect_is(lmer(Yield ~ 1|Batch, Dyestuff, subset=TRUE), "lmerMod") 118 expect_is(lmer(Yield ~ 1|Batch, Dyestuff, weights=rep(1,nrow(Dyestuff))), "lmerMod") 119 expect_is(lmer(Yield ~ 1|Batch, Dyestuff, na.action="na.exclude"), "lmerMod") 120 expect_is(lmer(Yield ~ 1|Batch, Dyestuff, offset=rep(0,nrow(Dyestuff))), "lmerMod") 121 expect_is(lmer(Yield ~ 1|Batch, Dyestuff, contrasts=NULL), "lmerMod") 122 expect_is(lmer(Yield ~ 1|Batch, Dyestuff, devFunOnly=FALSE), "lmerMod") 123 expect_is(lmer(Yield ~ 1|Batch, Dyestuff, control=lmerControl(optimizer="Nelder_Mead")), "lmerMod") 124 expect_is(lmer(Yield ~ 1|Batch, Dyestuff, control=lmerControl()), "lmerMod") 125 ## avoid _R_CHECK_LENGTH_1_LOGIC2_ errors ... 126 if (getRversion() < "3.6.0" || packageVersion("optimx")>"2018.7.10") { 127 expect_error(lmer(Yield ~ 1|Batch, Dyestuff, control=lmerControl(optimizer="optimx")),"must specify") 128 expect_is(lmer(Yield ~ 1|Batch, Dyestuff, 129 control=lmerControl(optimizer="optimx", 130 optCtrl=list(method="L-BFGS-B"))), 131 "lmerMod") 132 } 133 expect_error(lmer(Yield ~ 1|Batch, Dyestuff, control=lmerControl(optimizer="junk")), "couldn't find optimizer function") 134 ## disable test ... should be no warning 135 expect_is(lmer(Reaction ~ 1 + Days + (1 + Days | Subject), 136 data = sleepstudy, subset = (Days == 1 | Days == 9), 137 control=lmerControl(check.nobs.vs.rankZ="ignore", 138 check.nobs.vs.nRE="ignore", 139 check.conv.hess="ignore", 140 ## need to ignore relative gradient check too; 141 ## surface is flat so *relative* gradient gets large 142 check.conv.grad="ignore")), 143 "merMod") 144 expect_is(lmer(Reaction ~ 1 + Days + (1|obs), 145 data = transform(sleepstudy,obs=seq(nrow(sleepstudy))), 146 control=lmerControl(check.nobs.vs.nlev="ignore", 147 check.nobs.vs.nRE="ignore", 148 check.nobs.vs.rankZ="ignore")), 149 "merMod") 150 expect_error(lmer(Reaction ~ 1 + Days + (1|obs), 151 data = transform(sleepstudy,obs=seq(nrow(sleepstudy))), 152 "number of levels of each grouping factor")) 153 154 ## check for errors with illegal input checking options 155 flags <- lme4:::.get.checkingOpts(names(formals(lmerControl))) 156 .t <- lapply(flags, function(OPT) { 157 ## set each to invalid string: 158 ## cat(OPT,"\n") 159 expect_error(lFormula(Reaction~1+Days+(1|Subject), data = sleepstudy, 160 control = do.call(lmerControl, 161 ## Deliberate: fake typo 162 ## vvv 163 setNames(list("warnign"), OPT))), 164 "invalid control level") 165 }) 166 ## disable warning via options 167 options(lmerControl=list(check.nobs.vs.rankZ="ignore",check.nobs.vs.nRE="ignore")) 168 expect_is(fm4 <- lmer(Reaction ~ Days + (1|Subject), 169 subset(sleepstudy,Subject %in% levels(Subject)[1:4])), "merMod") 170 expect_is(lmer(Reaction ~ 1 + Days + (1 + Days | Subject), 171 data = sleepstudy, subset = (Days == 1 | Days == 9), 172 control=lmerControl(check.conv.hess="ignore", 173 check.conv.grad="ignore")), 174 "merMod") 175 options(lmerControl=NULL) 176 ## check for when ignored options are set 177 options(lmerControl=list(junk=1,check.conv.grad="ignore")) 178 expect_warning(lmer(Reaction ~ Days + (1|Subject),sleepstudy), 179 "some options") 180 options(lmerControl=NULL) 181 options(warn=0) 182 expect_error(lmer(Yield ~ 1|Batch, Dyestuff, junkArg=TRUE), "unused argument") 183 expect_warning(lmer(Yield ~ 1|Batch, Dyestuff, control=list()), 184 "passing control as list is deprecated") 185 if(FALSE) ## Hadley broke this 186 expect_warning(lmer(Yield ~ 1|Batch, Dyestuff, control=glmerControl()), 187 "passing control as list is deprecated") 188 189 ss <- transform(sleepstudy,obs=factor(seq(nrow(sleepstudy)))) 190 expect_warning(lmer(Reaction ~ 1 + (1|obs), data=ss, 191 control=lmerControl(check.nobs.vs.nlev="warning", 192 check.nobs.vs.nRE="ignore")), 193 "number of levels of each grouping factor") 194 195 ## test deparsing of very long terms inside mkReTrms 196 set.seed(101) 197 longNames <- sapply(letters[1:25], 198 function(x) paste(rep(x,8),collapse="")) 199 tstdat <- data.frame(Y=rnorm(10), 200 F=factor(1:10), 201 matrix(runif(250),ncol=25, 202 dimnames=list(NULL, 203 longNames))) 204 expect_is(lFormula(Y~1+(aaaaaaaa+bbbbbbbb+cccccccc+dddddddd+ 205 eeeeeeee+ffffffff+gggggggg+hhhhhhhh+ 206 iiiiiiii+jjjjjjjj+kkkkkkkk+llllllll|F), 207 data=tstdat, 208 control=lmerControl(check.nobs.vs.nlev="ignore", 209 check.nobs.vs.nRE="ignore", 210 check.nobs.vs.rankZ="ignore")),"list") 211 212 ## do.call(new,...) bug 213 new <- "foo" 214 expect_is(refit(fm1),"merMod") 215 rm("new") 216 217 ## test subset-with-( printing from summary 218 fm1 <- lmer(z~1|f,d,subset=(z<1e9)) 219 expect_equal(sum(grepl("Subset: \\(",capture.output(summary(fm1)))),1) 220 221 ## test messed-up Hessian 222 fm1 <- lmer(z~ as.numeric(f) + 1|f, d) 223 fm1@optinfo$derivs$Hessian[2,2] <- NA 224 expect_warning(lme4:::checkConv(fm1@optinfo$derivs, 225 coefs=c(1,1), 226 ctrl=lmerControl()$checkConv,lbound=0), 227 "Problem with Hessian check") 228 229 ## test ordering of Ztlist names 230 ## this is a silly model, just using it for a case 231 ## where nlevs(RE term 1) < nlevs(RE term 2)x 232 data(cbpp) 233 cbpp <- transform(cbpp,obs=factor(1:nrow(cbpp))) 234 fm0 <- lmer(incidence~1+(1|herd)+(1|obs),cbpp, 235 control=lmerControl(check.nobs.vs.nlev="ignore", 236 check.nobs.vs.rankZ="ignore", 237 check.nobs.vs.nRE="ignore", 238 check.conv.grad="ignore", 239 check.conv.singular="ignore", 240 check.conv.hess="ignore")) 241 fm0B <- update(fm0, .~1+(1|obs)+(1|herd)) 242 expect_equal(names(getME(fm0,"Ztlist")), 243 c("obs.(Intercept)", "herd.(Intercept)")) 244 ## stable regardless of order in formula 245 expect_equal(getME(fm0,"Ztlist"),getME(fm0B,"Ztlist")) 246 ## no optimization (GH #408) 247 fm_noopt <- lmer(z~1|f,d, 248 control=lmerControl(optimizer=NULL)) 249 expect_equal(unname(unlist(getME(fm_noopt,c("theta","beta")))), 250 c(0.244179074357121, -0.0336616441209862)) 251 expect_error(lmer(z~1|f,d, 252 control=lmerControl(optimizer="none")), 253 "deprecated use") 254 my_opt <- function(fn,par,lower,upper,control) { 255 opt <- optim(fn=fn,par=par,lower=lower, 256 upper=upper,control=control,,method="L-BFGS-B") 257 return(list(par=opt$par,fval=opt$value,conv=opt$convergence)) 258 } 259 expect_is(fm_noopt <- lmer(z~1|f,d, 260 control=lmerControl(optimizer=my_opt)),"merMod") 261 262 ## test verbose option for nloptwrap 263 cc <- capture.output(lmer(Reaction~1+(1|Subject), 264 data=sleepstudy, 265 control=lmerControl(optimizer="nloptwrap", 266 optCtrl=list(xtol_abs=1e-6, ftol_abs=1e-6)), 267 verbose=5)) 268 expect_equal(sum(grepl("^iteration:",cc)),14) 269 270}) ## test_that(..) 271 272test_that("coef_lmer", { 273 ## test coefficient extraction in the case where RE contain 274 ## terms that are missing from the FE ... 275 set.seed(101) 276 d <- data.frame(resp=runif(100), 277 var1=factor(sample(1:5,size=100,replace=TRUE)), 278 var2=runif(100), 279 var3=factor(sample(1:5,size=100,replace=TRUE))) 280 library(lme4) 281 mix1 <- lmer(resp ~ 0 + var1 + var1:var2 + (1|var3), data=d) 282 c1 <- coef(mix1) 283 expect_is(c1, "coef.mer") 284 cd1 <- c1$var3 285 expect_is (cd1, "data.frame") 286 n1 <- paste0("var1", 1:5) 287 nn <- c(n1, paste(n1, "var2", sep=":")) 288 expect_identical(names(cd1), c("(Intercept)", nn)) 289 expect_equal(fixef(mix1), 290 setNames(c(0.2703951, 0.3832911, 0.451279, 0.6528842, 0.6109819, 291 0.4949802, 0.1222705, 0.08702069, -0.2856431, -0.01596725), 292 nn), tolerance= 6e-6)# 64-bit: 6.73e-9 293}) 294 295test_that("getCall", { 296 ## GH #535 297 getClass <- function() "foo" 298 expect_is(glmer(round(Reaction) ~ 1 + (1|Subject), sleepstudy, 299 family=poisson), "glmerMod") 300 rm(getClass) 301}) 302 303test_that("better info about optimizer convergence", 304{ 305 set.seed(14) 306 cbpp$var <- rnorm(nrow(cbpp), 10, 10) 307 308 suppressWarnings(gm2 <- 309 glmer(cbind(incidence, size - incidence) ~ period * var + (1 | herd), 310 data = cbpp, family = binomial, 311 control=glmerControl(optimizer=c("bobyqa","Nelder_Mead"))) 312 ) 313 314 ## FIXME: with new update, suppressWarnings(update(gm2)) will give 315 ## Error in as.list.environment(X[[i]], ...) : 316 ## promise already under evaluation: recursive default argument reference or earlier problems? 317 op <- options(warn=-1) 318 gm3 <- update(gm2, 319 control=glmerControl(optimizer="bobyqa", 320 optCtrl=list(maxfun=2))) 321 options(op) 322 323 cc <-capture.output(print(summary(gm2))) 324 expect_equal(tail(cc,3)[1], 325 "optimizer (Nelder_Mead) convergence code: 0 (OK)") 326}) 327 328context("convergence warnings etc.") 329 330fm1 <- lmer(Reaction~ Days + (Days|Subject), sleepstudy) 331suppressMessages(fm0 <- lmer(Reaction~ Days + (Days|Subject), sleepstudy[1:20,])) 332 333msg_in_output <- function(x, str) { 334 cc <- capture.output(.prt.warn(x)) 335 any(grepl(str , cc)) 336} 337 338test_that("convergence warnings from limited evals", { 339 expect_warning(fm1B <- update(fm1, control=lmerControl(optCtrl=list(maxeval=3))), 340 "convergence code 5") 341 expect_true(msg_in_output(fm1B@optinfo, "convergence code: 5")) 342 expect_warning(fm1C <- update(fm1, control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=3))), 343 "maximum number of function evaluations exceeded") 344 expect_true(msg_in_output(fm1C@optinfo, 345 "maximum number of function evaluations exceeded")) 346 ## one extra (spurious) warning here ... 347 expect_warning(fm1D <- update(fm1, control=lmerControl(optimizer="Nelder_Mead",optCtrl=list(maxfun=3))), 348 "failure to converge in 3 evaluations") 349 expect_true(msg_in_output(fm1D@optinfo, 350 "failure to converge in 3 evaluations")) 351 expect_message(fm0D <- update(fm0, control=lmerControl(optimizer="Nelder_Mead")), 352 "boundary") 353 expect_true(msg_in_output(fm0D@optinfo, 354 "(OK)")) 355}) 356 357## GH 533 358test_that("test for zero non-NA cases", { 359 data_bad <- sleepstudy 360 data_bad$Days <- NA_real_ 361 expect_error(lmer(Reaction ~ Days + (1| Subject), data_bad), 362 "0 \\(non-NA\\) cases") 363}) 364