library("testthat") library("lme4") (testLevel <- if (nzchar(s <- Sys.getenv("LME4_TEST_LEVEL"))) as.numeric(s) else 1) ## use old (<=3.5.2) sample() algorithm if necessary if ("sample.kind" %in% names(formals(RNGkind))) suppressWarnings(RNGkind("Mersenne-Twister", "Inversion", "Rounding")) L <- load(system.file("testdata", "lme-tst-fits.rda", package="lme4", mustWork=TRUE)) ## FIXME: should test for old R versions, skip reloading test data in that ## case? fm0 <- fit_sleepstudy_0 fm1 <- fit_sleepstudy_1 fm2 <- fit_sleepstudy_2 gm1 <- fit_cbpp_1 gm2 <- fit_cbpp_2 gm3 <- fit_cbpp_3 ## More objects to use in all contexts : set.seed(101) dNA <- data.frame( xfac= sample(letters[1:10], 100, replace=TRUE), xcov= runif(100), resp= rnorm(100)) dNA[sample(1:100, 10), "xcov"] <- NA CI.boot <- function(fm, nsim=10, seed=101, ...) suppressWarnings(confint(fm, method="boot", nsim=nsim, quiet=TRUE, seed=seed, ...)) ## rSimple <- function(rep = 2, m.u = 2, kinds = c('fun', 'boring', 'meh')) { stopifnot(is.numeric(rep), rep >= 1, is.numeric(m.u), m.u >= 1, is.character(kinds), (nk <- length(kinds)) >= 1) nobs <- rep * m.u * nk data.frame(kind= rep(kinds, each=rep*m.u), unit = gl(m.u, 1, nobs), y = round(50*rnorm(nobs))) } d12 <- rSimple() data("Pixel", package="nlme") nPix <- nrow(Pixel) fmPix <- lmer(pixel ~ day + I(day^2) + (day | Dog) + (1 | Side/Dog), data = Pixel) test_that("summary", { ## test for multiple-correlation-warning bug and other 'correlation = *' variants ## Have 2 summary() versions, each with 3 print(.) ==> 6 x capture.output(.) sf.aa <- summary(fit_agridat_archbold) msg1 <- "Correlation.* not shown by default" ## message => *not* part of capture.*(.) expect_message(c1 <- capture.output(sf.aa), msg1) # correlation = NULL - default cF <- capture.output(print(sf.aa, correlation=FALSE)) ## TODO? ensure the above gives *no* message/warning/error expect_identical(c1, cF) expect_message( cT <- capture.output(print(sf.aa, correlation=TRUE)) , "Correlation.* could have been required in summary()") expect_identical(cF, cT[seq_along(cF)]) sfT.aa <- summary(fit_agridat_archbold, correlation=TRUE) ## no message any more ## expect_message(cT2 <- capture.output(sfT.aa), msg1) ## expect_identical(cF, cT2) cT3 <- capture.output(print(sfT.aa, correlation=TRUE)) expect_identical(cT, cT3) cF2 <- capture.output(print(sfT.aa, correlation=FALSE)) expect_identical(cF, cF2) }) test_that("lmer anova", { aa <- suppressMessages(anova(fm0,fm1)) expect_that(aa, is_a("anova")) expect_equal(names(aa), c("npar", "AIC", "BIC", "logLik", "deviance", "Chisq", "Df", "Pr(>Chisq)")) expect_warning(suppressMessages(do.call(anova,list(fm0,fm1))), "assigning generic names") ## dat <- data.frame(y = 1:5, u = c(rep("A",2), rep("B",3)), t = c(rep("A",3), rep("B",2))) datfun <- function(x) dat aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa <- dat expect_is(stats::anova(lmer(y ~ u + (1 | t), dat = aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa, REML=FALSE), lmer(y ~ 1 + (1 | t), dat = aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa, REML=FALSE)), "anova") expect_equal(rownames(stats::anova(lmer(y ~ u + (1 | t), dat = aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa, REML=FALSE), lmer(y ~ 1 + (1 | t), dat = aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa, REML=FALSE), model.names=c("a","b"))), c("b","a")) ff <- function(form) { lmer(form, dat=dat, REML=FALSE, control=lmerControl(check.conv.singular="ignore")) } expect_error(rownames(stats::anova(ff(y ~ u + (1 | t)), ff(y ~ 1 + (1 | t)), model.names=c("a","b","c"))), "different lengths") z <- 1 ## output not tested (but shouldn't fail) ss <- stats::anova(lmer(y ~ u + (1 | t), data = datfun(z), REML=FALSE), lmer(y ~ 1 + (1 | t), data = datfun(z), REML=FALSE)) ## ## from Roger Mundry via Roman Lustrik full <- lmer(resp ~ xcov + (1|xfac), data=dNA) null <- lmer(resp ~ 1 + (1|xfac), data=dNA) expect_error(anova(null,full), "models were not all fitted to the same size of dataset") }) if (requireNamespace("merDeriv")) { test_that("summary with merDeriv", { library(merDeriv) cc <- capture.output(print(summary(fm1))) expect_true(any(grepl("Correlation of Fixed Effects", cc))) ## WARNING, this will detach package but may not undo ## method loading ... detach("package:merDeriv") }) } ## Github issue #256 from Jonas Lindeløv -- issue is *not* specific for this dataset test_that("Two models with subset() within lmer()", { full3 <- lmer(y ~ kind + (1|unit), subset(d12, kind != 'boring'), REML=FALSE) null3 <- update(full3, .~. - kind) op <- options(warn = 2) # no warnings! ano3 <- anova(full3, null3)## issue #256: had warning in data != data[[1]] : ... o3 <- capture.output(ano3) # now prints with only one 'Data:' expect_equal(1, grep("^Data:", o3)) d12sub <- subset(d12, kind != 'boring') expect_is(full3s <- lmer(y ~ kind + (1|unit), d12sub, REML=FALSE), "lmerMod") expect_is(null3s <- update(full3s, .~. - kind), "lmerMod") expect_is(ano3s <- anova(full3s, null3s), "anova") expect_equal(ano3, ano3s, check.attributes=FALSE) options(op) }) test_that("anova() of glmer+glm models", { dat <<- data.frame(y = 1:5, u = c(rep("A",2), rep("B",3)), t = c(rep("A",3), rep("B",2))) cs <- glmerControl(check.conv.singular = "ignore") ## ignore singular fits gm1 <- glmer(y~(1|u), data=dat[1:4,], family=poisson, control = cs) gm0 <- glm(y~1, data=dat[1:4,], family=poisson) gm2 <- glmer(y~(1|u), data=dat[1:4,], family=poisson,nAGQ=2, control = cs) aa <- anova(gm1,gm0) expect_equal(aa[2,"Chisq"],0) expect_error(anova(gm2,gm0),"incommensurate") }) test_that("anova() of lmer+glm models", { dat2 <- dat set.seed(101) dat2$y <- rnorm(5) fm1 <- lmer(y~(1|u),data=dat2,REML=FALSE) fm0 <- lm(y~1,data=dat2) aa2 <- anova(fm1,fm0) expect_equal(aa2[2,"Chisq"],0) expect_warning(anova(fm1,type="III"),"additional arguments ignored") }) test_that("set p-values to NA for equivalent models: #583", { fm0B <- fm0 aa <- suppressMessages(anova(fm0B,fm0)) expect_true(all(is.na(aa[["Pr(>Chisq)"]]))) }) test_that("long names", { ## GH names(sleepstudy) <- c("Reaction", "Days", "Subject_xxxxxxxxxxxxxxxxxxxxxxxxxxx") fm1 <- lmer(Reaction ~ Days + (Days | Subject_xxxxxxxxxxxxxxxxxxxxxxxxxxx), sleepstudy) fm2 <- lmer(Reaction ~ Days + (Days || Subject_xxxxxxxxxxxxxxxxxxxxxxxxxxx), sleepstudy) expect_equal(length(attributes(suppressMessages(anova(fm1,fm2)))$heading),4) }) if (testLevel>1) { context("bootMer confint()") set.seed(47) test_that("bootMer", { ## testing bug-fix for ordering of sd/cor components in sd/cor matrix with >2 rows ## FIXME: This model makes no sense [and CI.boot() fails for "nloptwrap"!] dd <- expand.grid(A=factor(1:3),B=factor(1:10),rep=1:10) dd$y <- suppressMessages(simulate(~1 + (A|B), newdata=dd, newparams=list(beta=1,theta=rep(1,6), sigma=1), family=gaussian, seed=101))[[1]] m1 <- lmer(y ~ 1 + (A|B), data=dd, control=lmerControl(calc.deriv=FALSE)) ci <- CI.boot(m1,seed=101) ci2 <- CI.boot(m1,seed=101) expect_equal(ci,ci2) ci_50 <- CI.boot(m1,level=0.5,seed=101) expect_true(all(ci_50[,"25 %"]>ci[,"2.5 %"])) expect_true(all(ci_50[,"75 %"]1 context("confint_other") test_that("confint", { load(system.file("testdata", "gotway_hessianfly.rda", package = "lme4")) ## generated via: ## gotway_hessianfly_fit <- glmer(cbind(y, n-y) ~ gen + (1|block), ## data=gotway.hessianfly, family=binomial, ## control=glmerControl(check.nlev.gtreq.5="ignore")) ## gotway_hessianfly_prof <- profile(gotway_hessianfly_fit,which=1) ## save(list=ls(pattern="gotway"),file="gotway_hessianfly.rda") expect_equal(confint(gotway_hessianfly_prof)[1,1],0) ## FIXME: should add tests for {-1,1} bounds on correlations as well expect_equal(c(confint(fm1,method="Wald",parm="beta_")), c(232.301892,8.891041,270.508318,12.043531), tolerance=1e-5) ## Wald gives NA for theta values expect_true(all(is.na(confint(fm1,method="Wald",parm="theta_")))) ## check names ci1.p <- suppressWarnings(confint(fm1,quiet=TRUE)) ci1.w <- confint(fm1,method="Wald") ci1.b <- CI.boot(fm1, nsim=2) expect_equal(dimnames(ci1.p), list(c(".sig01", ".sigma", "(Intercept)", "Days"), c("2.5 %", "97.5 %"))) expect_equal(dimnames(ci1.p),dimnames(ci1.w)) expect_equal(dimnames(ci1.p),dimnames(ci1.b)) ci1.p.n <- suppressWarnings(confint(fm1,quiet=TRUE,oldNames=FALSE)) ci1.w.n <- confint(fm1,method="Wald", oldNames=FALSE) ci1.b.n <- CI.boot(fm1, nsim=2, oldNames=FALSE) expect_equal(dimnames(ci1.p.n), list(c("sd_(Intercept)|Subject", "sigma", "(Intercept)", "Days"), c("2.5 %", "97.5 %"))) expect_equal(dimnames(ci1.p.n),dimnames(ci1.w.n)) expect_equal(dimnames(ci1.p.n),dimnames(ci1.b.n)) ## test case of slightly wonky (spline fit fails) but monotonic profiles: ## simfun <- function(J,n_j,g00,g10,g01,g11,sig2_0,sig01,sig2_1){ N <- sum(rep(n_j,J)) x <- rnorm(N) z <- rnorm(J) mu <- c(0,0) sig <- matrix(c(sig2_0,sig01,sig01,sig2_1),ncol=2) u <- MASS::mvrnorm(J,mu=mu,Sigma=sig) b_0j <- g00 + g01*z + u[,1] b_1j <- g10 + g11*z + u[,2] y <- rep(b_0j,each=n_j)+rep(b_1j,each=n_j)*x + rnorm(N,0,sqrt(0.5)) sim_data <- data.frame(Y=y,X=x,Z=rep(z,each=n_j), group=rep(1:J,each=n_j)) } set.seed(102) dat <- simfun(10,5,1,.3,.3,.3,(1/18),0,(1/18)) fit <- lmer(Y~X+Z+X:Z+(X||group),data=dat) if (Sys.info()["sysname"] != "SunOS" && .Platform$OS.type != "windows") { ## doesn't produce warnings on Solaris, or win-builder ... expect_warning(pp <- profile(fit,"theta_"), "non-monotonic profile") expect_warning(cc <- confint(pp),"falling back to linear interpolation") ## very small/unstable problem, needs large tolerance expect_equal(unname(cc[2,]), c(0, 0.509), tolerance=0.09) # "bobyqa" had 0.54276 } badprof <- readRDS(system.file("testdata","badprof.rds", package="lme4")) expect_warning(cc <- confint(badprof), "falling back to linear") expect_equal(cc, array(c(0, -1, 2.50856219044636, 48.8305727797906, NA, NA, 33.1204478717389, 1, 7.33374326592662, 68.7254711217912, -6.90462047196017, NA), dim = c(6L, 2L), dimnames = list(c(".sig01", ".sig02", ".sig03", ".sigma", "(Intercept)", "cYear"), c("2.5 %", "97.5 %"))), tolerance=1e-3) }) test_that("refit", { s1 <- simulate(fm1) expect_is(refit(fm1,s1), "merMod") s2 <- simulate(fm1,2) expect_error(refit(fm1,s2), "refit not implemented .* lists") data(Orthodont,package = "nlme") fmOrth <- fm <- lmer(distance ~ I(age - 11) + (I(age - 11) | Subject), data = Orthodont) expect_equal(s1 <- simulate(fm,newdata = Orthodont,seed = 101), s2 <- simulate(fm,seed = 101)) ## works *without* offset ... m5 <- glmer(round(Reaction) ~ Days + (1|Subject), data = sleepstudy, family=poisson, offset=rep(0,nrow(sleepstudy))) m5R <- refit(m5) ## lots of fussy details make expect_equal() on the whole object difficult expect_equal(coef(m5),coef(m5R),tolerance=3e-6) expect_equal(VarCorr(m5),VarCorr(m5R),tolerance=1e-6) expect_equal(logLik(m5),logLik(m5R)) }) if (testLevel>1) { context("predict method") test_that("predict", { d1 <- with(cbpp, expand.grid(period = unique(period), herd = unique(herd))) d2 <- data.frame(period = "1", herd = unique(cbpp$herd)) d3 <- expand.grid(period = as.character(1:3), herd = unique(cbpp$herd)) p0 <- predict(gm1) p1 <- predict(gm1,d1) p2 <- predict(gm1,d2) p3 <- predict(gm1,d3) expect_equal(p0[1], p1[1]) expect_equal(p0[1], p2[1]) expect_equal(p0[1], p3[1]) expect_warning(predict(gm1, ReForm=NA), "is deprecated") ## matrix-valued predictors: Github #201 from Fabian S. sleepstudy$X <- cbind(1, sleepstudy$Days) m <- lmer(Reaction ~ -1 + X + (Days | Subject), sleepstudy) pm <- predict(m, newdata=sleepstudy) expect_is(pm, "numeric") expect_equal(quantile(pm, names = FALSE), c(211.0108, 260.9496, 296.873, 328.6378, 458.1584), tol=1e-5) op <- options(warn = 2) # there should be no warnings! if (require("MEMSS",quietly=TRUE)) { ## test spurious warning with factor as response variable data("Orthodont", package = "MEMSS") # (differently "coded" from the 'default' "nlme" one) silly <- glmer(Sex ~ distance + (1|Subject), data = Orthodont, family = binomial) sillypred <- data.frame(distance = c(20, 25)) ps <- predict(silly, sillypred, re.form=NA, type = "response") expect_is(ps, "numeric") expect_equal(unname(ps), c(0.999989632, 0.999997201), tolerance=1e-6) detach("package:MEMSS") } ## a case with interactions (failed in one temporary version): expect_warning(fmPixS <<- update(fmPix, .~. + Side), "nearly unidentifiable|unable to evaluate scaled gradient|failed to converge") ## (1|2|3); 2 and 3 seen (as Error??) on CRAN's Windows 32bit options(op) set.seed(1); ii <- sample(nrow(Pixel), 16) expect_equal(predict(fmPix, newdata = Pixel[ii,]), fitted(fmPix )[ii]) expect_equal(predict(fmPixS, newdata = Pixel[ii,]), fitted(fmPixS)[ii]) set.seed(7); n <- 100; y <- rnorm(n) dd <- data.frame(id = factor(sample(10, n, replace = TRUE)), x1 = 1, y = y, x2 = rnorm(n, mean = sign(y))) expect_message(m <- lmer(y ~ x1 + x2 + (1 | id), data = dd), "fixed-effect model matrix is rank deficient") expect_is(summary(m),"summary.merMod") ii <- sample(n, 16) expect_equal(predict(m, newdata = dd[ii,]), fitted(m)[ii]) ## predict(*, new..) gave Error in X %*% fixef(object) - now also drops col. ## predict(*, new..) with NA in data {and non-simple model}, issue #246: m1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy) sleepst.NA <- sleepstudy ; sleepst.NA$Days[2] <- NA m2 <- update(fm1, data = sleepst.NA) ## maybe tricky for evaluation; fm1 was defined elsewhere, so data expect_equal(length(predict(m2, sleepst.NA[1:4,])),4) ## Wrong 'b' constructed in mkNewReTrms() -- issue #257 data(Orthodont,package="nlme") Orthodont <- within(Orthodont, nsex <- as.numeric(Sex == "Male")) m3 <- lmer(distance ~ age + (age|Subject) + (0 + Sex |Subject), data=Orthodont, control=lmerControl(check.conv.hess="ignore", check.conv.grad="ignore")) m4 <- lmer(distance ~ age + (age|Subject) + (0 + nsex|Subject), data=Orthodont) expect_equal(p3 <- predict(m3, Orthodont), fitted(m3), tolerance=1e-14) expect_equal(p4 <- predict(m4, Orthodont), fitted(m4), tolerance=1e-14) ## related to GH #275 (*passes*), ss <- sleepstudy set.seed(1) ss$noiseChar <- ifelse(runif(nrow(sleepstudy)) > 0.8, "Yes", "No") ss$noiseFactor <- factor(ss$noiseChar) fm4 <- lmer(Reaction ~ Days + noiseChar + (Days | Subject), ss) expect_equal(predict(fm4, newdata = model.frame(fm4)[2:3, ])[2], predict(fm4, newdata = model.frame(fm4)[3, ])) fm3 <- lmer(Reaction ~ Days + noiseFactor + (Days | Subject), ss) expect_equal(predict(fm3, newdata = model.frame(fm3)[2:3, ])[2], predict(fm3, newdata = model.frame(fm3)[3, ])) ## complex-basis functions in RANDOM effect fm5 <- lmer(Reaction~Days+(poly(Days,2)|Subject),sleepstudy) expect_equal(predict(fm5,sleepstudy[1,]),fitted(fm5)[1]) ## complex-basis functions in FIXED effect fm6 <- lmer(Reaction~poly(Days,2)+(1|Subject),sleepstudy) expect_equal(predict(fm6,sleepstudy[1,]),fitted(fm6)[1]) ## GH #414: no warning about dropping contrasts on random effects op <- options(warn = 2) # there should be no warnings! set.seed(1) dat <- data.frame( fac = factor(rep(c("a", "b"), 100)), grp = rep(1:25, each = 4)) dat$y <- 0 contr <- 0.5 * contr.sum(2) rownames(contr) <- c("a", "b") colnames(contr) <- "a" contrasts(dat$fac) <- contr m1_contr <- lmer(y~fac+(fac|grp),dat) pp <- predict(m1_contr,newdata=dat) options(op) }) ## testLevel>1 test_that("simulate", { ## simulate() will look for data in environment of formula, find ## unmodified version of cbpp -- need to re-add observation-level factor ee <- environment(formula(gm2)) ee$cbpp$obs <- factor(seq(nrow(ee$cbpp))) expect_is(simulate(gm2), "data.frame") expect_warning(simulate(gm2, ReForm = NA), "is deprecated") expect_warning(simulate(gm2, REForm = NA), "is deprecated") p1 <- simulate(gm2, re.form = NULL, seed = 101) p2 <- simulate(gm2, re.form = ~0, seed = 101) p3 <- simulate(gm2, re.form = NA, seed = 101) p4 <- simulate(gm2, re.form = NULL, seed = 101) expect_warning(p5 <- simulate(gm2, ReForm = ~0, seed = 101), "is deprecated") p6 <- simulate(gm2, re.form = NA, seed = 101) expect_warning(p7 <- simulate(gm2, REForm = NULL, seed = 101), "is deprecated") p8 <- simulate(gm2, re.form = ~0, seed = 101) p9 <- simulate(gm2, re.form = NA, seed = 101) p10 <- simulate(gm2,use.u = FALSE, seed = 101) p11 <- simulate(gm2,use.u = TRUE, seed = 101) ## minimal check of content: expect_identical(colSums(p1[,1]), c(incidence = 95, 747)) expect_identical(colSums(p2[,1]), c(incidence = 109, 733)) ## equivalences: ## group ~0: expect_equal(p2,p3) expect_equal(p2,p5) expect_equal(p2,p6) expect_equal(p2,p8) expect_equal(p2,p9) expect_equal(p2,p10) ## group 1: expect_equal(p1,p4) expect_equal(p1,p7) expect_equal(p1,p11) expect_error(simulate(gm2,use.u = TRUE, re.form = NA), "should specify only one") ## ## hack: test with three REs p1 <- lmer(diameter ~ (1|plate) + (1|plate) + (1|sample), Penicillin, control = lmerControl(check.conv.hess = "ignore", check.conv.grad = "ignore")) expect_is(sp1 <- simulate(p1, seed=123), "data.frame") expect_identical(dim(sp1), c(nrow(Penicillin), 1L)) expect_equal(fivenum(sp1[,1]), c(20.864, 22.587, 23.616, 24.756, 28.599), tolerance=0.01) ## Pixel example expect_identical(dim(simulate(fmPixS)), c(nPix, 1L)) expect_identical(dim(simulate(fmPix )), c(nPix, 1L)) ## simulation with newdata smaller/larger different from original fm <- lmer(diameter ~ 1 + (1|plate) + (1|sample), Penicillin) expect_is(simulate(fm,newdata=Penicillin[1:10,],allow.new.levels=TRUE),"data.frame") expect_is(simulate(fm,newdata=do.call(rbind,replicate(4,Penicillin,simplify=FALSE))),"data.frame") ## negative binomial sims set.seed(101) dd <- data.frame(f=factor(rep(1:10,each=20)), x=runif(200), y=rnbinom(200,size=2,mu=2)) g1 <- glmer.nb(y ~ x + (1|f), data=dd) th.g1 <- getME(g1, "glmer.nb.theta") ## changed to setting seed internally ts1 <- table(s1 <- simulate(g1,seed=101)[,1]) ## ts1B <- table(s1 <- simulate(g1,seed=101)[,1]) expect_equal(fixef(g1), c("(Intercept)" = 0.630067, x = -0.0167248), tolerance = 1e-4) ## ?? Travis is getting hung up here/ignoring tolerance spec?? expect_equal(th.g1, 2.013, tolerance = 1e-4) expect_equal(th.g1, g1@call$family[["theta"]])# <- important for pkg{effects} eval() expect_identical(sum(s1), 413) expect_identical(as.vector(ts1[as.character(0:5)]), ## c(51L, 54L, 36L, 21L, 14L, 9L)) c(49L,56L,32L,25L,11L,9L)) ## de novo NB simulation ... s2 <- simulate(~x + (1|f),seed=101, family=MASS::negative.binomial(theta=th.g1), newparams=getME(g1,c("theta","beta")), newdata=dd)[,1] expect_equal(s1,s2) ## Simulate with newdata with *new* RE levels: d <- sleepstudy[-1] # droping the response ("Reaction") ## d$Subject <- factor(rep(1:18, each=10)) ## Add 18 new subjects: d <- rbind(d, d) d$Subject <- factor(rep(1:36, each=10)) d$simulated <- simulate(fm1, seed=1, newdata = d, re.form=NULL, allow.new.levels = TRUE)[,1] expect_equal(mean(d$simulated), 299.9384608) ## Simulate with weights: newdata <- with(cbpp, expand.grid(period=unique(period), herd=unique(herd))) ss <- simulate(gm1, newdata=newdata[1:3,], weights=20, seed=101)[[1]] expect_equal(ss, matrix(c(4,2,0,16,18,20),nrow=3, dimnames=list(NULL,c("incidence","")))) ss <- simulate(gm3, newdata=newdata[1:3,], weights=20, seed=101)[[1]] expect_equal(ss,c(0.2,0.1,0.0)) ss <- simulate(gm1, newdata=newdata[1,], weights=20, seed=101)[[1]] expect_equal(unname(ss),matrix(c(4,16),nrow=1)) ## simulate Gamma, from function and de novo set.seed(102) dd <- data.frame(x=rep(seq(-2,2,length=15),10), f=factor(rep(1:10,each=15))) u <- rnorm(10) dd$y <- with(dd, rgamma(nrow(dd),shape=2, scale=exp(2+1*x+u[as.numeric(f)])/2)) g1 <- glmer(y~x+(1|f),family=Gamma(link="log"),dd) s1 <- simulate(g1,seed=101) s2 <- suppressMessages(simulate(~x+(1|f), family=Gamma(link="log"), seed=101, newdata=dd, newparams=getME(g1,c("theta","beta","sigma")))) expect_equal(s1,s2) dd$y2 <- s2[[1]] g2 <- glmer(y2~x+(1|f),family=Gamma(link="log"),dd) expect_equal(fixef(g2), tolerance = 4e-7, # 32-bit windows showed 1.34e-7 c("(Intercept)" = 2.81887136759369, x= 1.06543222163626)) ## simulate with re.form = NULL and derived/offset components in formula fm7 <- lmer(Reaction ~ Days + offset(Days) + (1|Subject), sleepstudy) s7 <- simulate(fm7, seed = 101, re.form = NULL) ## thought this would break but it doesn't ??? f_wrap <- function() { Reaction ~ Days + offset(Days) + (1|Subject) } fm8 <- lmer(f_wrap(), sleepstudy) s8 <- simulate(fm8, seed = 101, re.form = NULL) expect_identical(s7, s8) ## harder: insert NA values in the offset and see if it handles this OK?? }) context("misc") test_that("misc", { expect_equal(df.residual(fm1),176) if (suppressWarnings(require(ggplot2))) { ## ggplot calls sample() [for silly start-up messages ## throws warning because we're using backward-compatible RNGkind expect_is(fortify.merMod(fm1), "data.frame") expect_is(fortify.merMod(gm1), "data.frame") } expect_is(as.data.frame(VarCorr(fm1)), "data.frame") }) } ## testLevel>1 context("plot") test_that("plot", { ## test getData() within plot function: reported by Dieter Menne doFit <- function(){ data(Orthodont,package = "nlme") data1 <- Orthodont lmer(distance ~ age + (age|Subject), data = data1) } data(Orthodont, package = "nlme") fm0 <- lmer(distance ~ age + (age|Subject), data = Orthodont) expect_is(plot(fm0), "trellis") suppressWarnings(rm("Orthodont")) fm <- doFit() pp <- plot(fm, resid(., scaled = TRUE) ~ fitted(.) | Sex, abline = 0) expect_is(pp, "trellis") ## test qqmath/getIDLabels() expect_is(q1 <- lattice::qqmath(fm,id=0.05),"trellis") cake2 <- transform(cake,replicate=as.numeric(replicate), recipe=as.numeric(recipe)) fm2 <- lmer(angle ~ recipe + temp + (1|recipe:replicate), cake2, REML= FALSE) expect_is(lattice::qqmath(fm2, id=0.05), "trellis") expect_is(lattice::qqmath(fm2, id=0.05, idLabels=~recipe), "trellis") expect_is(lattice::qqmath(fm2, 0.05, ~recipe), "trellis") # was expect_warning() expect_is(lattice::qqmath(fm2, 0.05), "trellis") # was expect_warning() }) context("misc") test_that("summary", { ## test that family() works when $family element is weird ## FIXME: is convergence warning here a false positive? gnb <- suppressWarnings(glmer(TICKS~1+(1|BROOD), family=MASS::negative.binomial(theta=2), data=grouseticks)) expect_is(family(gnb),"family") }) if (testLevel>1) { context("profile") test_that("profile", { ## FIXME: can we deal with convergence warning messages here ... ? ## fit profile on default sd/cor scale ... p1 <- suppressWarnings(profile(fm1,which="theta_")) ## and now on var/cov scale ... p2 <- suppressWarnings(profile(fm1,which="theta_", prof.scale="varcov")) ## because there are no correlations, squaring the sd results ## gives the same result as profiling on the variance scale ## in the first place expect_equal(confint(p1)^2,confint(p2), tolerance=1e-5) ## or via built-in varianceProf() function expect_equal(unname(confint(varianceProf(p1))), unname(confint(p2)), tolerance=1e-5) p3 <- profile(fm2,which=c(1,3,4)) p4 <- suppressWarnings(profile(fm2,which="theta_",prof.scale="varcov", signames=FALSE)) ## compare only for sd/var components, not corr component ## FAILS on r-patched-solaris-x86 2018-03-30 ??? ## 2/6 mismatches (average diff: 4.62) ## [1] 207 - 216 == -9.23697 ## [4] 1422 - 1422 == -0.00301 if (Sys.info()["sysname"] != "SunOS") { expect_equal(unname(confint(p3)^2), unname(confint(p4)[c(1,3,4),]), tolerance=1e-3) } ## check naming convention properly adjusted expect_equal(as.character(unique(p4$.par)), c("var_(Intercept)|Subject", "cov_Days.(Intercept)|Subject", "var_Days|Subject", "sigma")) }) test_that("densityplot is robust", { p <- readRDS(system.file("testdata","harmel_profile.rds", package="lme4")) expect_warning(lattice::densityplot(p), "unreliable profiles for some variables") }) } ## testLevel>1 context("model.frame") test_that("model.frame", { ## non-syntactic names d <- sleepstudy names(d)[1] <- "Reaction Time" ee <- function(m,nm) { expect_equal(names(model.frame(m, fixed.only=TRUE)),nm) } m <- lmer(Reaction ~ 1 + (1 | Subject), sleepstudy) ee(m,"Reaction") m2 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy) ee(m2,c("Reaction","Days")) m3 <- lmer(`Reaction Time` ~ Days + (1 | Subject), d) ee(m3, c("Reaction Time","Days")) m4 <- lmer(Reaction ~ log(1+Days) + (1 | Subject), sleepstudy) ee(m4, c("Reaction","log(1 + Days)")) }) context("influence measures") d <- as.data.frame(ChickWeight) colnames(d) <- c("y", "x", "subj", "tx") dNAs <- d dNAs$y[c(1, 3, 5)] <- NA fitNAs <- lmer(y ~ tx*x + (x | subj), data = dNAs, na.action=na.exclude) test_that("influence/hatvalues works", { ifm1 <- influence(fm1, do.coef=FALSE) expect_equal(unname(head(ifm1$hat)), c(0.107483311203734, 0.102096105816528, 0.0980557017761242, 0.0953620990825215, 0.0940152977357202, 0.0940152977357202), tolerance=1e-6) expect_equal(nrow(dNAs),length(hatvalues(fitNAs))) }) test_that("influence OK with tibbles", { if (requireNamespace("tibble")) { ## make small data set/example so influence() isn't too slow ... ss <- tibble::as_tibble(sleepstudy[1:60,]) smallfit <- lmer(Reaction ~ 1 + (1 | Subject), data = ss) i1 <- influence(smallfit, ncores = 1) expect_equal(head(i1[["fixed.effects[-case]"]]), structure(c(286.35044481665, 286.179896199062, 286.327301507498, 285.014692121823, 284.36060419176, 283.297551183126), dim = c(6L, 1L), dimnames = list(c("1", "2", "3", "4", "5", "6"), "(Intercept)")), tolerance = 1e-6) } }) test_that("rstudent", { rfm1 <- rstudent(fm1) expect_equal(unname(head(rfm1)), c(-1.45598270922089, -1.49664543508657, -2.11747425025103, -0.0729690066951975, 0.772716397142335, 2.37859408861768), tolerance=1e-6) expect_equal(nrow(dNAs),length(rstudent(fitNAs))) }) test_that("cooks distance", { expect_equal( unname(head(cooks.distance(fm1))), c(0.127645976734753, 0.127346548123793, 0.243724627125036, 0.000280638917214881, 0.0309804642689636, 0.293554225380831), tolerance=1e-6) expect_equal(nrow(dNAs),length(cooks.distance(fitNAs))) }) test_that("cooks distance on subject-level influence", { ifm1S <- influence(fm1, "Subject", ncores=1) expect_equal( unname(head(cooks.distance(ifm1S),2)), c(0.33921460279262, 0.290309061006305), tolerance = 1e-6) }) test_that("cooks distance on glmer models", { inf <- influence(gm1) inf.h <- influence(gm1, "herd", ncores=1) cook <- cooks.distance(inf) expect_equal(unname(head(cook, 3)), c(0.0532998800033037, 0.0405931172763581, 0.252608337928438), tolerance = 1e-6) cook.h <- cooks.distance(inf.h) expect_equal(unname(head(cook.h, 3)), c(0.256630560723611, 0.00525856231971531, 0.103355658099396), tolerance = 1e-6) }) ## tweaked example so estimated var = 0 zerodat <- data.frame(x=seq(0,1,length.out=120), f=rep(1:3,each=40)) zerodat$y1 <- simulate(~x+(1|f), family=gaussian, seed=102, newparams=list(beta=c(1,1), theta=c(0.001), sigma=1), newdata=zerodat)[[1]] zerodat$y2 <- simulate(~x+(1|f), family=poisson, seed=102, newparams=list(beta=c(1,1), theta=c(0.001)), newdata=zerodat)[[1]] test_that("rstudent matches for zero-var cases", { lmer_zero <- lmer(y1~x+(1|f), data=zerodat) glmer_zero <- glmer(y2~x+(1|f),family=poisson, data=zerodat) lm_zero <- lm(y1~x, data=zerodat) glm_zero <- glm(y2~x,family=poisson, data=zerodat) expect_equal(suppressWarnings(rstudent(glmer_zero)), rstudent(glm_zero), tolerance=0.01) expect_equal(suppressWarnings(rstudent(lmer_zero)), rstudent(lm_zero),tolerance=0.01) }) if (testLevel>1) { ## n.b. influence() doesn't work under system.time(); ## weird evaluation stuff ? ## FIXME: work on timing some more i1 <- influence(fm1, ncores=1) test_that("full version of influence", { expect_equal(c(head(i1[["fixed.effects[-case]"]],1)), c(252.323536264131, 10.3222704729148)) }) cd <- cooks.distance(i1) expect_equal(unname(head(cd,2)), c(0.016503344184025, 0.0106634053477361)) if (parallel::detectCores() > 1) { test_that("parallel influence", { i2 <- suppressMessages(influence(fm1, ncores=2)) ## if (packageVersion("Matrix") != "1.4.2") ## fow now,as they differ str(i1) str(i2) print(all.equal(i2, i2)) # to see diff # expect_equal(i1, i2) ## <<<<-------------- FAILS (4 MM) }) } } ## car method testing: influence timing with ncores > 1 ... ## car version 3.0.10. ## L <- load(system.file("testdata", "lme-tst-fits.rda", ## package="lme4", mustWork=TRUE)) ## data("sleepstudy", package="lme4") ## library(lme4) ## library(car) ## WANT warning about S3 method overwrite ... ## fm1 <- fit_sleepstudy_1 ## library(pracma) ## because system.time() is weird ## tic(); i1 <- influence(fm1); toc() ## 2+ seconds ## tic(); i2 <- influence(fm1, ncores=8); toc() ## 3.4 seconds test_that("influence with nAGQ=0", { gm1Q0 <- update(gm1, nAGQ=0) expect_is(influence(gm1Q0), "influence.merMod") }) if (testLevel > 1) withAutoprint({ test_that("cook's distance comparison", { ## generate data with zero variance set.seed(101) n <- 50 dd <- data.frame(x = rnorm(n), f = factor(rep(1:2, each = n/2))) suppressMessages(dd$y <- simulate(~ x + (1|f), newdata = dd, newparams = list(beta = c(2,2), theta = 0, sigma = 5), family = gaussian)[[1]]) (fm2 <- lmer(y~x + (1|f), dd, REML = FALSE)) (fm2L <- lm(y~x , dd)) i2 <- influence(fm2) ## hatvalues version does **not** match exactly ... expect_equal(cooks.distance(i2), cooks.distance(fm2L)) expect_equal(cooks.distance(fm2), cooks.distance(fm2L), tolerance = 1e-2) }) }) ## testLevel > 1