Last updated on 2026-07-06 14:52:53 CEST.
| Flavor | Version | Tinstall | Tcheck | Ttotal | Status | Flags |
|---|---|---|---|---|---|---|
| r-devel-linux-x86_64-debian-clang | 2.0-1 | 84.06 | 782.54 | 866.60 | OK | |
| r-devel-linux-x86_64-debian-gcc | 2.0-1 | 60.88 | 524.57 | 585.45 | OK | |
| r-devel-linux-x86_64-fedora-clang | 2.0-1 | 119.00 | 1195.64 | 1314.64 | ERROR | |
| r-devel-linux-x86_64-fedora-gcc | 2.0-1 | 164.00 | 970.22 | 1134.22 | ERROR | |
| r-devel-windows-x86_64 | 2.0-1 | 89.00 | 700.00 | 789.00 | OK | |
| r-patched-linux-x86_64 | 2.0-1 | 88.53 | 746.02 | 834.55 | OK | |
| r-release-linux-x86_64 | 2.0-1 | 86.57 | 748.47 | 835.04 | OK | |
| r-release-macos-arm64 | 2.0-1 | 18.00 | 163.00 | 181.00 | OK | |
| r-release-macos-x86_64 | 2.0-1 | 57.00 | 801.00 | 858.00 | OK | |
| r-release-windows-x86_64 | 2.0-1 | 92.00 | 687.00 | 779.00 | OK | |
| r-oldrel-macos-arm64 | 2.0-1 | 17.00 | 165.00 | 182.00 | OK | |
| r-oldrel-macos-x86_64 | 2.0-1 | 57.00 | 577.00 | 634.00 | OK | |
| r-oldrel-windows-x86_64 | 2.0-1 | 118.00 | 883.00 | 1001.00 | ERROR |
Version: 2.0-1
Check: examples
Result: ERROR
Running examples in ‘lme4-Ex.R’ failed
The error most likely occurred in:
> ### Name: salamander
> ### Title: Mountain dusky salamander mating
> ### Aliases: salamander
> ### Keywords: datasets
>
> ### ** Examples
>
> ## Making sure Male, Female, and CRoss are treated as factors
> salamander$Male <- factor(salamander$Male)
> salamander$Female <- factor(salamander$Female)
> salamander$Cross <- factor(salamander$Cross)
> ## Fitting the model described in 14.5.3 from McCullagh and Nelder
> sal_mod <- glmer(Mate ~ (1|Female) + (1 | Male) + Cross, data = salamander,
+ family = binomial(link = "logit"))
Error: Downdated VtV is not positive definite
Execution halted
Flavor: r-devel-linux-x86_64-fedora-clang
Version: 2.0-1
Check: tests
Result: ERROR
Running ‘AAAtest-all.R’ [214s/218s]
Running ‘HSAURtrees.R’
Running ‘REMLdev.R’
Running ‘ST.R’
Running ‘agridat_gotway.R’
Running ‘bootMer.R’ [11s/19s]
Running ‘boundary.R’ [29s/30s]
Running ‘confint.R’
Running ‘devCritFun.R’
Running ‘drop.R’
Running ‘drop1contrasts.R’
Running ‘dynload.R’
Running ‘elston.R’
Running ‘evalCall.R’
Running ‘extras.R’
Running ‘falsezero_dorie.R’
Running ‘fewlevels.R’
Running ‘getME.R’
Running ‘glmer-1.R’
Running ‘glmerControlPass.R’ [15s/15s]
Running ‘glmerWarn.R’
Running ‘glmmExt.R’ [12s/13s]
Running ‘glmmWeights.R’ [12s/13s]
Running ‘hatvalues.R’
Running ‘is.R’
Running ‘lmList-tst.R’
Running ‘lme4_nlme.R’
Running ‘lmer-0.R’
Running ‘lmer-1.R’
Running ‘lmer-conv.R’
Running ‘lmer2_ex.R’
Running ‘methods.R’
Running ‘minval.R’
Running ‘modFormula.R’
Running ‘nbinom.R’
Running ‘nlmer-conv.R’
Running ‘nlmer.R’
Running ‘offset.R’
Running ‘optimizer.R’
Running ‘polytomous.R’
Running ‘prLogistic.R’
Running ‘predict_basis.R’
Running ‘predsim.R’
Running ‘priorWeights.R’
Running ‘priorWeightsModComp.R’
Running ‘profile-tst.R’
Running ‘refit.R’
Running ‘resids.R’
Running ‘respiratory.R’
Running ‘simulate.R’
Running ‘test-glmernbref.R’
Running ‘testOptControl.R’
Running ‘testcolonizer.R’
Running ‘testcrab.R’ [14s/15s]
Running ‘throw.R’
Running ‘varcorr.R’
Running ‘vcov-etc.R’
Running the tests in ‘tests/AAAtest-all.R’ failed.
Complete output:
> if (base::require("testthat", quietly = TRUE)) {
+ pkg <- "lme4"
+ require(pkg, character.only=TRUE, quietly=TRUE)
+ if(getRversion() < "3.5.0") { withAutoprint <- identity ; prt <- print } else { prt <- identity }
+ if(Sys.getenv("USER") %in% c("maechler", "bbolker")) withAutoprint({
+ ## for developers' sake:
+ lP <- .libPaths() # ---- .libPaths() : ----
+ prt(lP)
+ ## ---- Entries in .libPaths()[1] : ----
+ prt(list.files(lP[1], include.dirs=TRUE))
+ prt(sessionInfo())
+ prt(packageDescription("Matrix"))
+ ## 'lme4' from packageDescription "file" :
+ prt(attr(packageDescription("lme4"), "file"))
+ })
+ test_check(pkg)
+ ##======== ^^^
+ print(warnings()) # TODO? catch most of these by expect_warning(..)
+ } else {
+ cat( "package 'testthat' not available, cannot run unit tests\n" )
+ }
Saving _problems/test-covariance_structures-398.R
[ FAIL 1 | WARN 0 | SKIP 6 | PASS 1673 ]
══ Skipped tests (6) ═══════════════════════════════════════════════════════════
• On CRAN (1): 'test-eval.R:2:1'
• Skipping (1): 'test-covariance_nlmer.R:22:3'
• empty test (3): , ,
• testLevel < 2 is TRUE (1): 'test-predict.R:689:3'
══ Failed tests ════════════════════════════════════════════════════════════════
── Error ('test-covariance_structures.R:396:1'): (code run outside of `test_that()`) ──
<Rcpp::exception/C++Error/error/condition>
Error: Downdated VtV is not positive definite
Backtrace:
▆
1. └─lme4::glmer(...) at test-covariance_structures.R:396:1
2. └─lme4::optimizeGlmer(...)
3. └─lme4:::optwrap(...)
4. ├─base::withCallingHandlers(...)
5. ├─base::do.call(optfun, arglist)
6. └─lme4 (local) `<fn>`(...)
7. ├─nM$newf(fn(nM$xeval()))
8. │ └─base::stopifnot(length(value <- as.numeric(value)) == 1L)
9. └─lme4 (local) fn(nM$xeval())
10. └─lme4 (local) pwrssUpdate(...)
[ FAIL 1 | WARN 0 | SKIP 6 | PASS 1673 ]
Error:
! Test failures.
Execution halted
Running the tests in ‘tests/glmer-1.R’ failed.
Complete output:
> if (lme4:::testLevel() > 1 || .Platform$OS.type!="windows") withAutoprint({
+
+ ## generalized linear mixed model
+ stopifnot(suppressPackageStartupMessages(require(lme4)))
+ options(show.signif.stars = FALSE)
+
+ source(system.file("test-tools-1.R", package = "Matrix"), keep.source = FALSE)
+ ##
+ ##' Check that coefficient +- "2" * SD contains true value
+ ##'
+ ##' @title Check that confidence interval for coefficients contains true value
+ ##' @param fm fitted model, e.g., from lm(), lmer(), glmer(), ..
+ ##' @param true.coef numeric vector of true (fixed effect) coefficients
+ ##' @param conf.level confidence level for confidence interval
+ ##' @param sd.factor the "2", i.e. default 1.96 factor for the confidence interval
+ ##' @return TRUE or a string of "error"
+ ##' @author Martin Maechler
+ chkFixed <- function(fm, true.coef, conf.level = 0.95,
+ sd.factor = qnorm((1+conf.level)/2))
+ {
+ stopifnot(is.matrix(cf <- coefficients(summary(fm))), ncol(cf) >= 2)
+ cc <- cf[,1]
+ sd <- cf[,2]
+ if(any(out1 <- true.coef < cc - sd.factor*sd))
+ return(sprintf("true coefficient[j], j=%s, is smaller than lower confidence limit",
+ paste(which(out1), collapse=", ")))
+ if(any(out2 <- true.coef > cc + sd.factor*sd))
+ return(sprintf("true coefficient[j], j=%s, is larger than upper confidence limit",
+ paste(which(out2), collapse=", ")))
+ ## else, return
+ TRUE
+ }
+
+
+ ## TODO: (1) move these to ./glmer-ex.R [DONE]
+ ## ---- (2) "rationalize" with ../man/cbpp.Rd
+ #m1e <- glmer1(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ # family = binomial, data = cbpp, doFit = FALSE)
+ ## now
+ #bobyqa(m1e, control = list(iprint = 2L))
+
+ m1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ family = binomial, data = cbpp)
+ m1. <- update(m1, start = getME(m1, c("theta", "fixef")))
+ dm1 <- drop1(m1)
+ stopifnot(all.equal(drop1(m1.), dm1, tol = 1e-10))# Lnx(F28) 64b: 4e-12
+ ## response as a vector of probabilities and usage of argument "weights"
+ m1p <- glmer(incidence / size ~ period + (1 | herd), weights = size,
+ family = binomial, data = cbpp)
+ ## Confirm that these are equivalent:
+ stopifnot(all.equal(fixef(m1), fixef(m1p)),
+ all.equal(ranef(m1), ranef(m1p)),
+ TRUE)
+ ## for(m in c(m1, m1p)) {
+ ## cat("-------\\n\\nCall: ",
+ ## paste(format(getCall(m)), collapse="\\n"), "\\n")
+ ## print(logLik(m)); cat("AIC:", AIC(m), "\\n") ; cat("BIC:", BIC(m),"\\n")
+ ## }
+ stopifnot(all.equal(logLik(m1), logLik(m1p)),
+ all.equal(AIC(m1), AIC(m1p)),
+ all.equal(BIC(m1), BIC(m1p)))
+
+
+ ## changed tolPwrss to 1e-7 to match other default
+ m1b <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ family = binomial, data = cbpp, verbose = 2L,
+ control =
+ glmerControl(optimizer="bobyqa", tolPwrss=1e-7,
+ optCtrl=list(rhobeg=0.2, rhoend=2e-7)))
+
+ ## using nAGQ=9L provides a better evaluation of the deviance
+ m.9 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ family = binomial, data = cbpp, nAGQ = 9)
+
+ ## check with nAGQ = 25
+ m2 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ family = binomial, data = cbpp, nAGQ = 25)
+
+ ## loosened tolerance on parameters
+ stopifnot(is((cm2 <- coef(m2)), "coef.mer"),
+ dim(cm2$herd) == c(15,4),
+ all.equal(fixef(m2),
+ ### lme4a [from an Ubuntu 11.10 amd64 system]
+ c(-1.39922533406847, -0.991407294757321,
+ -1.12782184600404, -1.57946627431248),
+ ##c(-1.3766013, -1.0058773,
+ ## -1.1430128, -1.5922817),
+ tolerance = 5.e-4,
+ check.attributes=FALSE),
+ all.equal(c(-2*logLik(m2)), 100.010030538022, tolerance=1e-9),
+ all.equal(deviance(m2), 73.373, tolerance=1e-5)
+ ## with bobyqa first (AGQ=0), then
+ ##all.equal(deviance(m2), 101.119749563, tolerance=1e-9)
+ )
+
+ ## 32-bit Ubuntu 10.04:
+ coef_m1_lme4.0 <- structure(c(-1.39853505102576,
+ -0.992334712470269, -1.12867541092127,
+ -1.58037389566025),
+ .Names = c("(Intercept)", "period2", "period3",
+ "period4"))
+
+ ## library(glmmADMB)
+ ## mg <- glmmadmb(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ ## family = "binomial", data = cbpp)
+ coef_m1_glmmadmb <- structure(c(-1.39853810064827, -0.99233330126975, -1.12867317840779,
+ -1.58031150854503), .Names = c("(Intercept)", "period2", "period3",
+ "period4"))
+
+ ## library(glmmML)
+ ## mm <- glmmML(cbind(incidence, size - incidence) ~ period,
+ ## cluster=herd,
+ ## family = "binomial", data = cbpp)
+ coef_m1_glmmML <- structure(c(-1.39853234657711, -0.992336901732793, -1.12867036466201,
+ -1.58030977686564), .Names = c("(Intercept)", "period2", "period3",
+ "period4"))
+
+ ## lme4[r 1636], 64-bit ubuntu 11.10:
+ ## c(-1.3788385, -1.0589543,
+ ## -1.1936382, -1.6306271),
+
+ stopifnot(is((cm1 <- coef(m1b)), "coef.mer"),
+ dim(cm1$herd) == c(15,4),
+ all.equal(fixef(m1b),fixef(m1),tolerance=4e-5),
+ is.all.equal4(fixef(m1b),
+ coef_m1_glmmadmb,
+ coef_m1_lme4.0,
+ coef_m1_glmmML,
+ tol = 5e-4)
+ )
+
+
+ ## Deviance for the new algorithm is lower, eventually we should change the previous test
+ ##stopifnot(deviance(m1) <= deviance(m1e))
+
+ showProc.time() #
+
+ if (require('MASS', quietly = TRUE)) {
+ bacteria$wk2 <- bacteria$week > 2
+ contrasts(bacteria$trt) <-
+ structure(contr.sdif(3),
+ dimnames = list(NULL, c("diag", "encourage")))
+ print(fm5 <- glmer(y ~ trt + wk2 + (1|ID),
+ data=bacteria, family=binomial))
+ showProc.time() #
+
+ stopifnot(
+ all.equal(logLik(fm5),
+ ## was -96.127838
+ structure(-96.13069, nobs = 220L, nall = 220L,
+ df = 5L, REML = FALSE,
+ class = "logLik"),
+ tolerance = 5e-4, check.attributes = FALSE)
+ ,
+ all.equal(fixef(fm5),
+ ## was 2.834218798 -1.367099481
+ c("(Intercept)"= 2.831609490, "trtdiag"= -1.366722631,
+ ## now 0.5842291915, -1.599148773
+ "trtencourage"=0.5840147802, "wk2TRUE"=-1.598591346),
+ tolerance = 1e-4 )
+ )
+ }
+
+ ## Failure to specify a random effects term - used to give an obscure message
+ ## Ensure *NON*-translated message; works on Linux,... :
+ if(.Platform$OS.type == "unix") {
+ Sys.setlocale("LC_MESSAGES", "C")
+ tc <- tryCatch(
+ m2 <- glmer(incidence / size ~ period, weights = size,
+ family = binomial, data = cbpp)
+ , error = function(.) .)
+ stopifnot(inherits(tc, "error"),
+ identical(tc$message,
+ "No random effects terms specified in formula"))
+ }
+
+
+ ## glmer - Modeling overdispersion as "mixture" aka
+ ## ----- - *ONE* random effect *PER OBSERVATION" -- example inspired by Ben Bolker:
+
+ ##' <description>
+ ##'
+ ##' <details>
+ ##' @title
+ ##' @param ng number of groups
+ ##' @param nr number of "runs", i.e., observations per groups
+ ##' @param sd standard deviations of group and "Individual" random effects,
+ ##' (\sigma_f, \sigma_I)
+ ##' @param b true beta (fixed effects)
+ ##' @return a data frame (to be used in glmer()) with columns
+ ##' (x, f, obs, eta0, eta, mu, y), where y ~ Pois(lambda(x)),
+ ##' log(lambda(x_i)) = b_1 + b_2 * x + G_{f(i)} + I_i
+ ##' and G_k ~ N(0, \sigma_f); I_i ~ N(0, \sigma_I)
+ ##' @author Ben Bolker and Martin Maechler
+ rPoisGLMMi <- function(ng, nr, sd=c(f = 1, ind = 0.5), b=c(1,2))
+ {
+ stopifnot(nr >= 1, ng >= 1,
+ is.numeric(sd), names(sd) %in% c("f","ind"), sd >= 0)
+ ntot <- nr*ng
+ b.reff <- rnorm(ng, sd= sd[["f"]])
+ b.rind <- rnorm(ntot,sd= sd[["ind"]])
+ x <- runif(ntot)
+ within(data.frame(x,
+ f = factor(rep(LETTERS[1:ng], each=nr)),
+ obs = 1:ntot,
+ eta0 = cbind(1, x) %*% b),
+ {
+ eta <- eta0 + b.reff[f] + b.rind[obs]
+ mu <- exp(eta)
+ y <- rpois(ntot, lambda=mu)
+ })
+ }
+
+ set.seed(1)
+ dd <- rPoisGLMMi(12, 20)
+ m0 <- glmer(y~x + (1|f), family="poisson", data=dd)
+ m1 <- glmer(y~x + (1|f) + (1|obs), family="poisson", data=dd)
+ stopifnot(isTRUE(chkFixed(m0, true.coef = c(1,2))),
+ isTRUE(chkFixed(m1, true.coef = c(1,2))))
+ (a01 <- anova(m0, m1))
+
+ stopifnot(all.equal(a01$Chisq[2], 554.334056, tolerance=1e-5),
+ all.equal(a01$logLik, c(-1073.77193, -796.604902), tolerance=1e-6),
+ a01$ npar == 3:4,
+ na.omit(a01$ Df) == 1)
+
+ if(lme4:::testLevel() > 1) {
+ nsim <- 10
+ set.seed(2)
+ system.time(
+ simR <- lapply(1:nsim, function(i) {
+ cat(i,"", if(i %% 20 == 0)"\n")
+ dd <- rPoisGLMMi(10 + rpois(1, lambda=3),
+ 16 + rpois(1, lambda=5))
+ m0 <- glmer(y~x + (1|f), family="poisson", data=dd)
+ m1 <- glmer(y~x + (1|f) + (1|obs), family="poisson", data=dd)
+ a01 <- anova(m0, m1)
+ stopifnot(a01$ npar == 3:4,
+ na.omit(a01$ Df) == 1)
+ list(chk0 = chkFixed(m0, true.coef = c(1,2)),
+ chk1 = chkFixed(m1, true.coef = c(1,2)),
+ chisq= a01$Chisq[2],
+ lLik = a01$logLik)
+ }))
+
+ ## m0 is the wrong model, so we don't expect much here:
+ table(unlist(lapply(simR, `[[`, "chk0")))
+
+
+ ## If the fixed effect estimates were unbiased and the standard errors correct,
+ ## and N(0,sigma^2) instead of t_{nu} good enough for the fixed effects,
+ ## the confidence interval should contain the true coef in ~95 out of 100:
+ table(unlist(lapply(simR, `[[`, "chk1")))
+
+ ## The tests are all highly significantly in favor of m1 :
+ summary(chi2s <- sapply(simR, `[[`, "chisq"))
+ ## Min. 1st Qu. Median Mean 3rd Qu. Max.
+ ## 158.9 439.0 611.4 698.2 864.3 2268.0
+ stopifnot(chi2s > qchisq(0.9999, df = 1))
+ }
+
+ showProc.time()
+ }) ## skip if windows and testLevel<1
> stopifnot(suppressPackageStartupMessages(require(lme4)))
> options(show.signif.stars = FALSE)
> source(system.file("test-tools-1.R", package = "Matrix"), keep.source = FALSE)
Loading required package: tools
> chkFixed <- function(fm, true.coef, conf.level = 0.95, sd.factor = qnorm((1 +
+ conf.level)/2)) {
+ stopifnot(is.matrix(cf <- coefficients(summary(fm))), ncol(cf) >= 2)
+ cc <- cf[, 1]
+ sd <- cf[, 2]
+ if (any(out1 <- true.coef < cc - sd.factor * sd))
+ return(sprintf("true coefficient[j], j=%s, is smaller than lower confidence limit",
+ paste(which(out1), collapse = ", ")))
+ if (any(out2 <- true.coef > cc + sd.factor * sd))
+ return(sprintf("true coefficient[j], j=%s, is larger than upper confidence limit",
+ paste(which(out2), collapse = ", ")))
+ TRUE
+ }
> m1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ family = binomial, data = cbpp)
> m1. <- update(m1, start = getME(m1, c("theta", "fixef")))
> dm1 <- drop1(m1)
> stopifnot(all.equal(drop1(m1.), dm1, tol = 1e-10))
> m1p <- glmer(incidence/size ~ period + (1 | herd), weights = size, family = binomial,
+ data = cbpp)
> stopifnot(all.equal(fixef(m1), fixef(m1p)), all.equal(ranef(m1), ranef(m1p)),
+ TRUE)
> stopifnot(all.equal(logLik(m1), logLik(m1p)), all.equal(AIC(m1), AIC(m1p)),
+ all.equal(BIC(m1), BIC(m1p)))
> m1b <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ family = binomial, data = cbpp, verbose = 2L, control = glmerControl(optimizer = "bobyqa",
+ tolPwrss = 1e-07, optCtrl = list(rhobeg = 0.2, rhoend = 2e-07)))
npt = 3 , n = 1
rhobeg = 0.2 , rhoend = 2e-07
start par. = 1 fn = 186.7231
rho: 0.020 eval: 4 fn: 184.166 par:0.600000
rho: 0.0020 eval: 7 fn: 184.110 par:0.649419
rho: 0.00020 eval: 10 fn: 184.109 par:0.641956
rho: 2.0e-05 eval: 12 fn: 184.109 par:0.641847
rho: 2.0e-06 eval: 13 fn: 184.109 par:0.641847
rho: 2.0e-07 eval: 15 fn: 184.109 par:0.641839
At return
eval: 18 fn: 184.10869 par: 0.641839
npt = 7 , n = 5
rhobeg = 0.2 , rhoend = 2e-07
start par. = 0.6418386 -1.360476 -0.9761732 -1.111073 -1.559676 fn = 184.1086
rho: 0.020 eval: 8 fn: 184.109 par:0.641839 -1.36048 -0.976173 -1.11107 -1.55968
rho: 0.0020 eval: 15 fn: 184.056 par:0.641943 -1.40262 -0.981786 -1.13822 -1.57895
rho: 0.00020 eval: 30 fn: 184.053 par:0.642120 -1.39845 -0.991559 -1.12809 -1.58002
rho: 2.0e-05 eval: 37 fn: 184.053 par:0.642074 -1.39838 -0.991818 -1.12817 -1.57968
rho: 2.0e-06 eval: 49 fn: 184.053 par:0.642064 -1.39834 -0.991914 -1.12821 -1.57974
rho: 2.0e-07 eval: 57 fn: 184.053 par:0.642064 -1.39833 -0.991924 -1.12821 -1.57975
At return
eval: 69 fn: 184.05313 par: 0.642064 -1.39833 -0.991924 -1.12821 -1.57975
> m.9 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ family = binomial, data = cbpp, nAGQ = 9)
> m2 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ family = binomial, data = cbpp, nAGQ = 25)
> stopifnot(is((cm2 <- coef(m2)), "coef.mer"), dim(cm2$herd) == c(15, 4),
+ all.equal(fixef(m2), c(-1.39922533406847, -0.991407294757321, -1.12782184600404,
+ -1.57946627431248), tolerance = 5e-04, check.attributes = FALSE), all.equal(c(-2 *
+ logLik(m2)), 100.010030538022, tolerance = 1e-09), all.equal(deviance(m2),
+ 73.373, tolerance = 1e-05))
> coef_m1_lme4.0 <- structure(c(-1.39853505102576, -0.992334712470269, -1.12867541092127,
+ -1.58037389566025), .Names = c("(Intercept)", "period2", "period3", "period4"))
> coef_m1_glmmadmb <- structure(c(-1.39853810064827, -0.99233330126975,
+ -1.12867317840779, -1.58031150854503), .Names = c("(Intercept)", "period2", "period3",
+ "period4"))
> coef_m1_glmmML <- structure(c(-1.39853234657711, -0.992336901732793, -1.12867036466201,
+ -1.58030977686564), .Names = c("(Intercept)", "period2", "period3", "period4"))
> stopifnot(is((cm1 <- coef(m1b)), "coef.mer"), dim(cm1$herd) == c(15, 4),
+ all.equal(fixef(m1b), fixef(m1), tolerance = 4e-05), is.all.equal4(fixef(m1b),
+ coef_m1_glmmadmb, coef_m1_lme4.0, coef_m1_glmmML, tol = 5e-04))
> showProc.time()
Time (user system elapsed): 2.176 0.067 2.306
> if (require("MASS", quietly = TRUE)) {
+ bacteria$wk2 <- bacteria$week > 2
+ contrasts(bacteria$trt) <- structure(contr.sdif(3), dimnames = list(NULL, c("diag",
+ "encourage")))
+ print(fm5 <- glmer(y ~ trt + wk2 + (1 | ID), data = bacteria, family = binomial))
+ showProc.time()
+ stopifnot(all.equal(logLik(fm5), structure(-96.13069, nobs = 220L, nall = 220L,
+ df = 5L, REML = FALSE, class = "logLik"), tolerance = 5e-04, check.attributes = FALSE),
+ all.equal(fixef(fm5), c(`(Intercept)` = 2.83160949, trtdiag = -1.366722631,
+ trtencourage = 0.5840147802, wk2TRUE = -1.598591346), tolerance = 1e-04))
+ }
Error in (function (cond) :
error in evaluating the argument 'x' in selecting a method for function 'print': Downdated VtV is not positive definite
Calls: withAutoprint ... <Anonymous> -> <Anonymous> -> stopifnot -> fn -> pwrssUpdate
Execution halted
Running the tests in ‘tests/optimizer.R’ failed.
Complete output:
> library(lme4)
Loading required package: Matrix
> source(system.file("test-tools-1.R", package = "Matrix"), keep.source = FALSE)
Loading required package: tools
> ## N.B. is.all.equal4() and assert.EQ() use 'tol', not 'tolerance'
>
>
> ## should be able to run any example with any bounds-constrained optimizer ...
> ## Nelder_Mead, bobyqa built in; optimx/nlminb, optimx/L-BFGS-B
> ## optimx/Rcgmin will require a bit more wrapping/interface work (requires gradient)
>
> if (.Platform$OS.type != "windows") {
+ fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy) ## Nelder_Mead
+ fm1B <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy,
+ control=lmerControl(optimizer="bobyqa"))
+ stopifnot(all.equal(fixef(fm1),fixef(fm1B)))
+ require(optimx)
+ lmerCtrl.optx <- function(method, ...)
+ lmerControl(optimizer="optimx", ..., optCtrl=list(method=method))
+ glmerCtrl.optx <- function(method, ...)
+ glmerControl(optimizer="optimx", ..., optCtrl=list(method=method))
+
+ (testLevel <- lme4:::testLevel())
+
+ ## FAILS on Windows (on r-forge only, not win-builder)... 'function is infeasible at initial parameters'
+ ## (can we test whether we are on r-forge??)
+ if (.Platform$OS.type != "windows") {
+ fm1C <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy,
+ control=lmerCtrl.optx(method="nlminb"))
+ fm1D <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy,
+ control=lmerCtrl.optx(method="L-BFGS-B"))
+ stopifnot(is.all.equal4(fixef(fm1),fixef(fm1B),fixef(fm1C),fixef(fm1D)))
+
+ fm1E <- update(fm1,control=lmerCtrl.optx(method=c("nlminb","L-BFGS-B")))
+ ## hack equivalence of call and optinfo
+ fm1E@call <- fm1C@call
+ fm1E@optinfo <- fm1C@optinfo
+ assert.EQ(fm1C,fm1E, tol=1e-5, giveRE=TRUE)# prints unless tolerance=0--equality
+ }
+
+ gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ data = cbpp, family = binomial,
+ control=glmerControl(tolPwrss=1e-13))
+ gm1B <- update(gm1, control=glmerControl (tolPwrss=1e-13, optimizer="bobyqa"))
+ gm1C <- update(gm1, control=glmerCtrl.optx(tolPwrss=1e-13, method="nlminb"))
+ gm1D <- update(gm1, control=glmerCtrl.optx(tolPwrss=1e-13, method="L-BFGS-B"))
+ stopifnot(is.all.equal4(fixef(gm1),fixef(gm1B),fixef(gm1C),fixef(gm1D),
+ tol=1e-5))
+
+ if (testLevel > 1) {
+ gm1E <- update(gm1, control=
+ glmerCtrl.optx(tolPwrss=1e-13, method=c("nlminb","L-BFGS-B")))
+ ## hack equivalence of call and optinfo
+ gm1E@call <- gm1C@call
+ gm1E@optinfo <- gm1C@optinfo
+ assert.EQ(gm1E,gm1C, tol=1e-5, giveRE=TRUE)# prints unless tol=0--equality
+ }
+ } ## skip on windows (for speed)
Loading required package: optimx
Error in eval(call, envir = pf) : Downdated VtV is not positive definite
Calls: update ... optimizeGlmer -> optwrap -> deriv12 -> fun -> pwrssUpdate
In addition: Warning messages:
1: In optimx.run(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, :
Gradient not computable after method nlminb
2: In optimx.run(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, :
Gradient not computable after method nlminb
3: In optimx.run(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, :
Gradient not computable after method nlminb
Execution halted
Running the tests in ‘tests/respiratory.R’ failed.
Complete output:
> ## Data originally from Davis 1991 Stat. Med., as packaged in geepack
> ## and transformed (center, id -> factor, idctr created, levels labeled)
> library(lme4)
Loading required package: Matrix
>
> if (.Platform$OS.type != "windows") {
+ load(system.file("testdata","respiratory.RData",package="lme4"))
+ m_glmer_4.L <- glmer(outcome~center+treat+sex+age+baseline+(1|idctr),
+ family=binomial,data=respiratory)
+
+ m_glmer_4.GHQ5 <- glmer(outcome~center+treat+sex+age+baseline+(1|idctr),
+ family=binomial,data=respiratory,nAGQ=5)
+
+ m_glmer_4.GHQ8 <- glmer(outcome~center+treat+sex+age+baseline+(1|idctr),
+ family=binomial,data=respiratory,nAGQ=8)
+
+ m_glmer_4.GHQ16 <- glmer(outcome~center+treat+sex+age+baseline+(1|idctr),
+ family=binomial,data=respiratory,nAGQ=16)
+ } ## skip on windows (for speed)
Error: Downdated VtV is not positive definite
Execution halted
Flavor: r-devel-linux-x86_64-fedora-clang
Version: 2.0-1
Check: tests
Result: ERROR
Running ‘AAAtest-all.R’ [236s/297s]
Running ‘HSAURtrees.R’
Running ‘REMLdev.R’
Running ‘ST.R’
Running ‘agridat_gotway.R’
Running ‘bootMer.R’ [11s/26s]
Running ‘boundary.R’ [29s/42s]
Running ‘confint.R’ [8s/11s]
Running ‘devCritFun.R’
Running ‘drop.R’
Running ‘drop1contrasts.R’
Running ‘dynload.R’ [12s/14s]
Running ‘elston.R’
Running ‘evalCall.R’
Running ‘extras.R’
Running ‘falsezero_dorie.R’
Running ‘fewlevels.R’
Running ‘getME.R’
Running ‘glmer-1.R’
Running ‘glmerControlPass.R’ [14s/14s]
Running ‘glmerWarn.R’
Running ‘glmmExt.R’ [11s/11s]
Running ‘glmmWeights.R’
Running ‘hatvalues.R’
Running ‘is.R’
Running ‘lmList-tst.R’
Running ‘lme4_nlme.R’
Running ‘lmer-0.R’
Running ‘lmer-1.R’
Running ‘lmer-conv.R’
Running ‘lmer2_ex.R’
Running ‘methods.R’
Running ‘minval.R’
Running ‘modFormula.R’
Running ‘nbinom.R’
Running ‘nlmer-conv.R’
Running ‘nlmer.R’
Running ‘offset.R’
Running ‘optimizer.R’
Running ‘polytomous.R’
Running ‘prLogistic.R’
Running ‘predict_basis.R’
Running ‘predsim.R’
Running ‘priorWeights.R’
Running ‘priorWeightsModComp.R’
Running ‘profile-tst.R’
Running ‘refit.R’
Running ‘resids.R’
Running ‘respiratory.R’
Running ‘simulate.R’
Running ‘test-glmernbref.R’
Running ‘testOptControl.R’
Running ‘testcolonizer.R’
Running ‘testcrab.R’
Running ‘throw.R’
Running ‘varcorr.R’
Running ‘vcov-etc.R’
Running the tests in ‘tests/AAAtest-all.R’ failed.
Complete output:
> if (base::require("testthat", quietly = TRUE)) {
+ pkg <- "lme4"
+ require(pkg, character.only=TRUE, quietly=TRUE)
+ if(getRversion() < "3.5.0") { withAutoprint <- identity ; prt <- print } else { prt <- identity }
+ if(Sys.getenv("USER") %in% c("maechler", "bbolker")) withAutoprint({
+ ## for developers' sake:
+ lP <- .libPaths() # ---- .libPaths() : ----
+ prt(lP)
+ ## ---- Entries in .libPaths()[1] : ----
+ prt(list.files(lP[1], include.dirs=TRUE))
+ prt(sessionInfo())
+ prt(packageDescription("Matrix"))
+ ## 'lme4' from packageDescription "file" :
+ prt(attr(packageDescription("lme4"), "file"))
+ })
+ test_check(pkg)
+ ##======== ^^^
+ print(warnings()) # TODO? catch most of these by expect_warning(..)
+ } else {
+ cat( "package 'testthat' not available, cannot run unit tests\n" )
+ }
Saving _problems/test-glmFamily-187.R
Saving _problems/test-isSingular-97.R
[ FAIL 2 | WARN 0 | SKIP 6 | PASS 1766 ]
══ Skipped tests (6) ═══════════════════════════════════════════════════════════
• On CRAN (1): 'test-eval.R:2:1'
• Skipping (1): 'test-covariance_nlmer.R:22:3'
• empty test (3): , ,
• testLevel < 2 is TRUE (1): 'test-predict.R:689:3'
══ Failed tests ════════════════════════════════════════════════════════════════
── Error ('test-glmFamily.R:186:3'): estimated Inverse Gaussian shape is correct ──
<Rcpp::exception/C++Error/error/condition>
Error: Downdated VtV is not positive definite
Backtrace:
▆
1. └─lme4::glmer(...) at test-glmFamily.R:186:3
2. └─lme4::optimizeGlmer(...)
3. └─lme4:::optwrap(...)
4. ├─base::withCallingHandlers(...)
5. ├─base::do.call(optfun, arglist)
6. └─lme4 (local) `<fn>`(...)
7. ├─nM$newf(fn(nM$xeval()))
8. │ └─base::stopifnot(length(value <- as.numeric(value)) == 1L)
9. └─lme4 (local) fn(nM$xeval())
10. └─lme4 (local) pwrssUpdate(...)
── Error ('test-isSingular.R:96:3'): checking singular fit for merMod ──────────
<Rcpp::exception/C++Error/error/condition>
Error: Downdated VtV is not positive definite
Backtrace:
▆
1. ├─base::suppressWarnings(...) at test-isSingular.R:96:3
2. │ └─base::withCallingHandlers(...)
3. └─lme4::glmer(form, family = poisson(link = "log"), data = dat)
4. └─lme4::optimizeGlmer(...)
5. └─lme4:::optwrap(...)
6. ├─base::withCallingHandlers(...)
7. ├─base::do.call(optfun, arglist)
8. └─lme4 (local) `<fn>`(...)
9. ├─nM$newf(fn(nM$xeval()))
10. │ └─base::stopifnot(length(value <- as.numeric(value)) == 1L)
11. └─lme4 (local) fn(nM$xeval())
12. └─lme4 (local) pwrssUpdate(...)
[ FAIL 2 | WARN 0 | SKIP 6 | PASS 1766 ]
Error:
! Test failures.
Execution halted
Running the tests in ‘tests/glmer-1.R’ failed.
Complete output:
> if (lme4:::testLevel() > 1 || .Platform$OS.type!="windows") withAutoprint({
+
+ ## generalized linear mixed model
+ stopifnot(suppressPackageStartupMessages(require(lme4)))
+ options(show.signif.stars = FALSE)
+
+ source(system.file("test-tools-1.R", package = "Matrix"), keep.source = FALSE)
+ ##
+ ##' Check that coefficient +- "2" * SD contains true value
+ ##'
+ ##' @title Check that confidence interval for coefficients contains true value
+ ##' @param fm fitted model, e.g., from lm(), lmer(), glmer(), ..
+ ##' @param true.coef numeric vector of true (fixed effect) coefficients
+ ##' @param conf.level confidence level for confidence interval
+ ##' @param sd.factor the "2", i.e. default 1.96 factor for the confidence interval
+ ##' @return TRUE or a string of "error"
+ ##' @author Martin Maechler
+ chkFixed <- function(fm, true.coef, conf.level = 0.95,
+ sd.factor = qnorm((1+conf.level)/2))
+ {
+ stopifnot(is.matrix(cf <- coefficients(summary(fm))), ncol(cf) >= 2)
+ cc <- cf[,1]
+ sd <- cf[,2]
+ if(any(out1 <- true.coef < cc - sd.factor*sd))
+ return(sprintf("true coefficient[j], j=%s, is smaller than lower confidence limit",
+ paste(which(out1), collapse=", ")))
+ if(any(out2 <- true.coef > cc + sd.factor*sd))
+ return(sprintf("true coefficient[j], j=%s, is larger than upper confidence limit",
+ paste(which(out2), collapse=", ")))
+ ## else, return
+ TRUE
+ }
+
+
+ ## TODO: (1) move these to ./glmer-ex.R [DONE]
+ ## ---- (2) "rationalize" with ../man/cbpp.Rd
+ #m1e <- glmer1(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ # family = binomial, data = cbpp, doFit = FALSE)
+ ## now
+ #bobyqa(m1e, control = list(iprint = 2L))
+
+ m1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ family = binomial, data = cbpp)
+ m1. <- update(m1, start = getME(m1, c("theta", "fixef")))
+ dm1 <- drop1(m1)
+ stopifnot(all.equal(drop1(m1.), dm1, tol = 1e-10))# Lnx(F28) 64b: 4e-12
+ ## response as a vector of probabilities and usage of argument "weights"
+ m1p <- glmer(incidence / size ~ period + (1 | herd), weights = size,
+ family = binomial, data = cbpp)
+ ## Confirm that these are equivalent:
+ stopifnot(all.equal(fixef(m1), fixef(m1p)),
+ all.equal(ranef(m1), ranef(m1p)),
+ TRUE)
+ ## for(m in c(m1, m1p)) {
+ ## cat("-------\\n\\nCall: ",
+ ## paste(format(getCall(m)), collapse="\\n"), "\\n")
+ ## print(logLik(m)); cat("AIC:", AIC(m), "\\n") ; cat("BIC:", BIC(m),"\\n")
+ ## }
+ stopifnot(all.equal(logLik(m1), logLik(m1p)),
+ all.equal(AIC(m1), AIC(m1p)),
+ all.equal(BIC(m1), BIC(m1p)))
+
+
+ ## changed tolPwrss to 1e-7 to match other default
+ m1b <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ family = binomial, data = cbpp, verbose = 2L,
+ control =
+ glmerControl(optimizer="bobyqa", tolPwrss=1e-7,
+ optCtrl=list(rhobeg=0.2, rhoend=2e-7)))
+
+ ## using nAGQ=9L provides a better evaluation of the deviance
+ m.9 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ family = binomial, data = cbpp, nAGQ = 9)
+
+ ## check with nAGQ = 25
+ m2 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ family = binomial, data = cbpp, nAGQ = 25)
+
+ ## loosened tolerance on parameters
+ stopifnot(is((cm2 <- coef(m2)), "coef.mer"),
+ dim(cm2$herd) == c(15,4),
+ all.equal(fixef(m2),
+ ### lme4a [from an Ubuntu 11.10 amd64 system]
+ c(-1.39922533406847, -0.991407294757321,
+ -1.12782184600404, -1.57946627431248),
+ ##c(-1.3766013, -1.0058773,
+ ## -1.1430128, -1.5922817),
+ tolerance = 5.e-4,
+ check.attributes=FALSE),
+ all.equal(c(-2*logLik(m2)), 100.010030538022, tolerance=1e-9),
+ all.equal(deviance(m2), 73.373, tolerance=1e-5)
+ ## with bobyqa first (AGQ=0), then
+ ##all.equal(deviance(m2), 101.119749563, tolerance=1e-9)
+ )
+
+ ## 32-bit Ubuntu 10.04:
+ coef_m1_lme4.0 <- structure(c(-1.39853505102576,
+ -0.992334712470269, -1.12867541092127,
+ -1.58037389566025),
+ .Names = c("(Intercept)", "period2", "period3",
+ "period4"))
+
+ ## library(glmmADMB)
+ ## mg <- glmmadmb(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ ## family = "binomial", data = cbpp)
+ coef_m1_glmmadmb <- structure(c(-1.39853810064827, -0.99233330126975, -1.12867317840779,
+ -1.58031150854503), .Names = c("(Intercept)", "period2", "period3",
+ "period4"))
+
+ ## library(glmmML)
+ ## mm <- glmmML(cbind(incidence, size - incidence) ~ period,
+ ## cluster=herd,
+ ## family = "binomial", data = cbpp)
+ coef_m1_glmmML <- structure(c(-1.39853234657711, -0.992336901732793, -1.12867036466201,
+ -1.58030977686564), .Names = c("(Intercept)", "period2", "period3",
+ "period4"))
+
+ ## lme4[r 1636], 64-bit ubuntu 11.10:
+ ## c(-1.3788385, -1.0589543,
+ ## -1.1936382, -1.6306271),
+
+ stopifnot(is((cm1 <- coef(m1b)), "coef.mer"),
+ dim(cm1$herd) == c(15,4),
+ all.equal(fixef(m1b),fixef(m1),tolerance=4e-5),
+ is.all.equal4(fixef(m1b),
+ coef_m1_glmmadmb,
+ coef_m1_lme4.0,
+ coef_m1_glmmML,
+ tol = 5e-4)
+ )
+
+
+ ## Deviance for the new algorithm is lower, eventually we should change the previous test
+ ##stopifnot(deviance(m1) <= deviance(m1e))
+
+ showProc.time() #
+
+ if (require('MASS', quietly = TRUE)) {
+ bacteria$wk2 <- bacteria$week > 2
+ contrasts(bacteria$trt) <-
+ structure(contr.sdif(3),
+ dimnames = list(NULL, c("diag", "encourage")))
+ print(fm5 <- glmer(y ~ trt + wk2 + (1|ID),
+ data=bacteria, family=binomial))
+ showProc.time() #
+
+ stopifnot(
+ all.equal(logLik(fm5),
+ ## was -96.127838
+ structure(-96.13069, nobs = 220L, nall = 220L,
+ df = 5L, REML = FALSE,
+ class = "logLik"),
+ tolerance = 5e-4, check.attributes = FALSE)
+ ,
+ all.equal(fixef(fm5),
+ ## was 2.834218798 -1.367099481
+ c("(Intercept)"= 2.831609490, "trtdiag"= -1.366722631,
+ ## now 0.5842291915, -1.599148773
+ "trtencourage"=0.5840147802, "wk2TRUE"=-1.598591346),
+ tolerance = 1e-4 )
+ )
+ }
+
+ ## Failure to specify a random effects term - used to give an obscure message
+ ## Ensure *NON*-translated message; works on Linux,... :
+ if(.Platform$OS.type == "unix") {
+ Sys.setlocale("LC_MESSAGES", "C")
+ tc <- tryCatch(
+ m2 <- glmer(incidence / size ~ period, weights = size,
+ family = binomial, data = cbpp)
+ , error = function(.) .)
+ stopifnot(inherits(tc, "error"),
+ identical(tc$message,
+ "No random effects terms specified in formula"))
+ }
+
+
+ ## glmer - Modeling overdispersion as "mixture" aka
+ ## ----- - *ONE* random effect *PER OBSERVATION" -- example inspired by Ben Bolker:
+
+ ##' <description>
+ ##'
+ ##' <details>
+ ##' @title
+ ##' @param ng number of groups
+ ##' @param nr number of "runs", i.e., observations per groups
+ ##' @param sd standard deviations of group and "Individual" random effects,
+ ##' (\sigma_f, \sigma_I)
+ ##' @param b true beta (fixed effects)
+ ##' @return a data frame (to be used in glmer()) with columns
+ ##' (x, f, obs, eta0, eta, mu, y), where y ~ Pois(lambda(x)),
+ ##' log(lambda(x_i)) = b_1 + b_2 * x + G_{f(i)} + I_i
+ ##' and G_k ~ N(0, \sigma_f); I_i ~ N(0, \sigma_I)
+ ##' @author Ben Bolker and Martin Maechler
+ rPoisGLMMi <- function(ng, nr, sd=c(f = 1, ind = 0.5), b=c(1,2))
+ {
+ stopifnot(nr >= 1, ng >= 1,
+ is.numeric(sd), names(sd) %in% c("f","ind"), sd >= 0)
+ ntot <- nr*ng
+ b.reff <- rnorm(ng, sd= sd[["f"]])
+ b.rind <- rnorm(ntot,sd= sd[["ind"]])
+ x <- runif(ntot)
+ within(data.frame(x,
+ f = factor(rep(LETTERS[1:ng], each=nr)),
+ obs = 1:ntot,
+ eta0 = cbind(1, x) %*% b),
+ {
+ eta <- eta0 + b.reff[f] + b.rind[obs]
+ mu <- exp(eta)
+ y <- rpois(ntot, lambda=mu)
+ })
+ }
+
+ set.seed(1)
+ dd <- rPoisGLMMi(12, 20)
+ m0 <- glmer(y~x + (1|f), family="poisson", data=dd)
+ m1 <- glmer(y~x + (1|f) + (1|obs), family="poisson", data=dd)
+ stopifnot(isTRUE(chkFixed(m0, true.coef = c(1,2))),
+ isTRUE(chkFixed(m1, true.coef = c(1,2))))
+ (a01 <- anova(m0, m1))
+
+ stopifnot(all.equal(a01$Chisq[2], 554.334056, tolerance=1e-5),
+ all.equal(a01$logLik, c(-1073.77193, -796.604902), tolerance=1e-6),
+ a01$ npar == 3:4,
+ na.omit(a01$ Df) == 1)
+
+ if(lme4:::testLevel() > 1) {
+ nsim <- 10
+ set.seed(2)
+ system.time(
+ simR <- lapply(1:nsim, function(i) {
+ cat(i,"", if(i %% 20 == 0)"\n")
+ dd <- rPoisGLMMi(10 + rpois(1, lambda=3),
+ 16 + rpois(1, lambda=5))
+ m0 <- glmer(y~x + (1|f), family="poisson", data=dd)
+ m1 <- glmer(y~x + (1|f) + (1|obs), family="poisson", data=dd)
+ a01 <- anova(m0, m1)
+ stopifnot(a01$ npar == 3:4,
+ na.omit(a01$ Df) == 1)
+ list(chk0 = chkFixed(m0, true.coef = c(1,2)),
+ chk1 = chkFixed(m1, true.coef = c(1,2)),
+ chisq= a01$Chisq[2],
+ lLik = a01$logLik)
+ }))
+
+ ## m0 is the wrong model, so we don't expect much here:
+ table(unlist(lapply(simR, `[[`, "chk0")))
+
+
+ ## If the fixed effect estimates were unbiased and the standard errors correct,
+ ## and N(0,sigma^2) instead of t_{nu} good enough for the fixed effects,
+ ## the confidence interval should contain the true coef in ~95 out of 100:
+ table(unlist(lapply(simR, `[[`, "chk1")))
+
+ ## The tests are all highly significantly in favor of m1 :
+ summary(chi2s <- sapply(simR, `[[`, "chisq"))
+ ## Min. 1st Qu. Median Mean 3rd Qu. Max.
+ ## 158.9 439.0 611.4 698.2 864.3 2268.0
+ stopifnot(chi2s > qchisq(0.9999, df = 1))
+ }
+
+ showProc.time()
+ }) ## skip if windows and testLevel<1
> stopifnot(suppressPackageStartupMessages(require(lme4)))
> options(show.signif.stars = FALSE)
> source(system.file("test-tools-1.R", package = "Matrix"), keep.source = FALSE)
Loading required package: tools
> chkFixed <- function(fm, true.coef, conf.level = 0.95, sd.factor = qnorm((1 +
+ conf.level)/2)) {
+ stopifnot(is.matrix(cf <- coefficients(summary(fm))), ncol(cf) >= 2)
+ cc <- cf[, 1]
+ sd <- cf[, 2]
+ if (any(out1 <- true.coef < cc - sd.factor * sd))
+ return(sprintf("true coefficient[j], j=%s, is smaller than lower confidence limit",
+ paste(which(out1), collapse = ", ")))
+ if (any(out2 <- true.coef > cc + sd.factor * sd))
+ return(sprintf("true coefficient[j], j=%s, is larger than upper confidence limit",
+ paste(which(out2), collapse = ", ")))
+ TRUE
+ }
> m1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ family = binomial, data = cbpp)
> m1. <- update(m1, start = getME(m1, c("theta", "fixef")))
> dm1 <- drop1(m1)
> stopifnot(all.equal(drop1(m1.), dm1, tol = 1e-10))
> m1p <- glmer(incidence/size ~ period + (1 | herd), weights = size, family = binomial,
+ data = cbpp)
> stopifnot(all.equal(fixef(m1), fixef(m1p)), all.equal(ranef(m1), ranef(m1p)),
+ TRUE)
> stopifnot(all.equal(logLik(m1), logLik(m1p)), all.equal(AIC(m1), AIC(m1p)),
+ all.equal(BIC(m1), BIC(m1p)))
> m1b <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ family = binomial, data = cbpp, verbose = 2L, control = glmerControl(optimizer = "bobyqa",
+ tolPwrss = 1e-07, optCtrl = list(rhobeg = 0.2, rhoend = 2e-07)))
npt = 3 , n = 1
rhobeg = 0.2 , rhoend = 2e-07
start par. = 1 fn = 186.7231
rho: 0.020 eval: 4 fn: 184.166 par:0.600000
rho: 0.0020 eval: 7 fn: 184.110 par:0.649419
rho: 0.00020 eval: 10 fn: 184.109 par:0.641956
rho: 2.0e-05 eval: 12 fn: 184.109 par:0.641847
rho: 2.0e-06 eval: 13 fn: 184.109 par:0.641847
rho: 2.0e-07 eval: 15 fn: 184.109 par:0.641839
At return
eval: 18 fn: 184.10869 par: 0.641839
npt = 7 , n = 5
rhobeg = 0.2 , rhoend = 2e-07
start par. = 0.6418386 -1.360476 -0.9761732 -1.111073 -1.559676 fn = 184.1086
rho: 0.020 eval: 8 fn: 184.109 par:0.641839 -1.36048 -0.976173 -1.11107 -1.55968
rho: 0.0020 eval: 15 fn: 184.056 par:0.641943 -1.40262 -0.981786 -1.13822 -1.57895
rho: 0.00020 eval: 30 fn: 184.053 par:0.642120 -1.39845 -0.991559 -1.12809 -1.58002
rho: 2.0e-05 eval: 37 fn: 184.053 par:0.642074 -1.39838 -0.991818 -1.12817 -1.57968
rho: 2.0e-06 eval: 49 fn: 184.053 par:0.642064 -1.39834 -0.991914 -1.12821 -1.57974
rho: 2.0e-07 eval: 57 fn: 184.053 par:0.642064 -1.39833 -0.991924 -1.12821 -1.57975
At return
eval: 69 fn: 184.05313 par: 0.642064 -1.39833 -0.991924 -1.12821 -1.57975
> m.9 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ family = binomial, data = cbpp, nAGQ = 9)
> m2 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
+ family = binomial, data = cbpp, nAGQ = 25)
> stopifnot(is((cm2 <- coef(m2)), "coef.mer"), dim(cm2$herd) == c(15, 4),
+ all.equal(fixef(m2), c(-1.39922533406847, -0.991407294757321, -1.12782184600404,
+ -1.57946627431248), tolerance = 5e-04, check.attributes = FALSE), all.equal(c(-2 *
+ logLik(m2)), 100.010030538022, tolerance = 1e-09), all.equal(deviance(m2),
+ 73.373, tolerance = 1e-05))
> coef_m1_lme4.0 <- structure(c(-1.39853505102576, -0.992334712470269, -1.12867541092127,
+ -1.58037389566025), .Names = c("(Intercept)", "period2", "period3", "period4"))
> coef_m1_glmmadmb <- structure(c(-1.39853810064827, -0.99233330126975,
+ -1.12867317840779, -1.58031150854503), .Names = c("(Intercept)", "period2", "period3",
+ "period4"))
> coef_m1_glmmML <- structure(c(-1.39853234657711, -0.992336901732793, -1.12867036466201,
+ -1.58030977686564), .Names = c("(Intercept)", "period2", "period3", "period4"))
> stopifnot(is((cm1 <- coef(m1b)), "coef.mer"), dim(cm1$herd) == c(15, 4),
+ all.equal(fixef(m1b), fixef(m1), tolerance = 4e-05), is.all.equal4(fixef(m1b),
+ coef_m1_glmmadmb, coef_m1_lme4.0, coef_m1_glmmML, tol = 5e-04))
> showProc.time()
Time (user system elapsed): 2.117 0.005 2.143
> if (require("MASS", quietly = TRUE)) {
+ bacteria$wk2 <- bacteria$week > 2
+ contrasts(bacteria$trt) <- structure(contr.sdif(3), dimnames = list(NULL, c("diag",
+ "encourage")))
+ print(fm5 <- glmer(y ~ trt + wk2 + (1 | ID), data = bacteria, family = binomial))
+ showProc.time()
+ stopifnot(all.equal(logLik(fm5), structure(-96.13069, nobs = 220L, nall = 220L,
+ df = 5L, REML = FALSE, class = "logLik"), tolerance = 5e-04, check.attributes = FALSE),
+ all.equal(fixef(fm5), c(`(Intercept)` = 2.83160949, trtdiag = -1.366722631,
+ trtencourage = 0.5840147802, wk2TRUE = -1.598591346), tolerance = 1e-04))
+ }
Error in (function (cond) :
error in evaluating the argument 'x' in selecting a method for function 'print': Downdated VtV is not positive definite
Calls: withAutoprint ... <Anonymous> -> <Anonymous> -> stopifnot -> fn -> pwrssUpdate
Execution halted
Running the tests in ‘tests/glmmWeights.R’ failed.
Complete output:
> if (.Platform$OS.type != "windows") {
+
+ library(lme4)
+ library(testthat)
+
+ source(system.file("testdata/lme-tst-funs.R", package="lme4", mustWork=TRUE))
+ ##-> gSim(), a general simulation function ...
+
+ ## hand-coded Pearson residuals {for sumFun() }
+ mypresid <- function(x) {
+ mu <- fitted(x)
+ (getME(x,"y") - mu) * sqrt(weights(x)) / sqrt(x@resp$family$variance(mu))
+ }
+
+ ## should be equal (up to numerical error) to weights(.,type="working")
+ workingWeights <- function(mod) mod@resp$weights*(mod@resp$muEta()^2)/mod@resp$variance()
+
+ ##' Sum of weighted residuals, 4 ways; the last three are identical
+ sumFun <- function(m) {
+ wrss1 <- m@devcomp$cmp["wrss"]
+ wrss2 <- sum(residuals(m,type="pearson")^2)
+ wrss3 <- sum(m@resp$wtres^2)
+ ## compare to hand-fitted Pearson resids ...
+ wrss4 <- sum(mypresid(m)^2)
+ c(wrss1,wrss2,wrss3,wrss4)
+ }
+ ## The relative "error"/differences of the weights w[] entries
+ rel.diff <- function(w) abs(1 - w[-1]/w[1])
+
+ set.seed(101)
+
+ ## GAMMA
+ g0 <- glmer(y~x+(1|block),data=gSim(),family=Gamma)
+ expect_true(all(rel.diff(sumFun(g0)) < 1e-13))
+ expect_equal(weights(g0, type = "working"), workingWeights(g0),
+ tolerance = 1e-4) ## FIXME: why is such a high tolerance required?
+
+ ## BERNOULLI
+ g1 <- glmer(y~x+(1|block),data=gSim(family=binomial(),nbinom=1),
+ family=binomial)
+ expect_true(all(rel.diff(sumFun(g1)) < 1e-13))
+ expect_equal(weights(g1, type = "working"), workingWeights(g1),
+ tolerance = 1e-5) ## FIXME: why is such a high tolerance required?
+
+
+ ## POISSON
+ (n <- nrow(d.P <- gSim(family=poisson())))
+ g2 <- glmer(y ~ x + (1|block), data = d.P, family=poisson)
+ g2W <- glmer(y ~ x + (1|block), data = d.P, family=poisson, weights = rep(2,n))
+ expect_true(all(rel.diff(sumFun(g2 )) < 1e-13))
+ expect_true(all(rel.diff(sumFun(g2W)) < 1e-13))
+ ## correct
+ expect_equal(weights(g2, type = "working"), workingWeights(g2),
+ tolerance = 1e-5) ## FIXME: why is such a high tolerance required?
+ expect_equal(weights(g2W, type = "working"), workingWeights(g2W),
+ tolerance = 1e-5) ## FIXME: why is such a high tolerance required?
+
+
+ ## non-Bernoulli BINOMIAL
+ g3 <- glmer(y ~ x + (1|block), data= gSim(family=binomial(), nbinom=10),
+ family=binomial)
+ expect_true(all(rel.diff(sumFun(g3)) < 1e-13))
+ expect_equal(weights(g3, type = "working"), workingWeights(g3),
+ tolerance = 1e-4) ## FIXME: why is such a high tolerance required?
+
+
+
+ d.b.2 <- gSim(nperblk = 2, family=binomial())
+ g.b.2 <- glmer(y ~ x + (1|block), data=d.b.2, family=binomial)
+
+ expect_true(all(rel.diff(sumFun(g.b.2 )) < 1e-13))
+
+
+ ## Many blocks of only 2 observations each - (but nicely balanced)
+ ## Want this "as" https://github.com/lme4/lme4/issues/47
+ ## (but it "FAILS" survival already):
+ ##
+ ## n2 = n/2 :
+ n2 <- 2048
+ if(FALSE)
+ n2 <- 100 # for building/testing
+ set.seed(47)
+ dB2 <- gSim(n2, nperblk = 2, x= rep(0:1, each= n2), family=binomial())
+ ## -- -- --- --------
+ gB2 <- glmer(y ~ x + (1|block), data=dB2, family=binomial)
+ expect_true(all(rel.diff(sumFun(gB2)) < 1e-13))
+
+ ## NB: Finite sample bias of \hat\sigma_1 and \hat\beta_1 ("Intercept")
+ ## tend to zero only slowly for n2 -> Inf, e.g., for
+ ## n2 = 2048, b1 ~= 4.3 (instead of 4); s1 ~= 1.3 (instead of 1)
+
+ ## FAILS -----
+ ## library(survival)
+ ## (gSurv.B2 <- clogit(y ~ x + strata(block), data=dB2))
+ ## ## --> Error in Surv(rep(1, 200L), y) : Time and status are different lengths
+ ## summary(gSurv.B2)
+ ## (SE.surf <- sqrt(diag(vcov(gSurv.B2))))
+
+
+
+ g3 <- glmer(y ~ x + (1|block),data=gSim(family=binomial(),nbinom=10),
+ family=binomial)
+ expect_equal(var(sumFun(g3)),0)
+
+ ## check dispersion parameter
+ ## (lowered tolerance to pass checks on my machine -- SCW)
+ expect_equal(sigma(g0)^2, 0.4888248, tolerance=1e-4)
+
+ } ## skip on windows (for speed)
Loading required package: Matrix
Error: Downdated VtV is not positive definite
Execution halted
Running the tests in ‘tests/respiratory.R’ failed.
Complete output:
> ## Data originally from Davis 1991 Stat. Med., as packaged in geepack
> ## and transformed (center, id -> factor, idctr created, levels labeled)
> library(lme4)
Loading required package: Matrix
>
> if (.Platform$OS.type != "windows") {
+ load(system.file("testdata","respiratory.RData",package="lme4"))
+ m_glmer_4.L <- glmer(outcome~center+treat+sex+age+baseline+(1|idctr),
+ family=binomial,data=respiratory)
+
+ m_glmer_4.GHQ5 <- glmer(outcome~center+treat+sex+age+baseline+(1|idctr),
+ family=binomial,data=respiratory,nAGQ=5)
+
+ m_glmer_4.GHQ8 <- glmer(outcome~center+treat+sex+age+baseline+(1|idctr),
+ family=binomial,data=respiratory,nAGQ=8)
+
+ m_glmer_4.GHQ16 <- glmer(outcome~center+treat+sex+age+baseline+(1|idctr),
+ family=binomial,data=respiratory,nAGQ=16)
+ } ## skip on windows (for speed)
Error: Downdated VtV is not positive definite
Execution halted
Flavor: r-devel-linux-x86_64-fedora-gcc
Version: 2.0-1
Check: tests
Result: ERROR
Running 'AAAtest-all.R' [177s]
Running 'HSAURtrees.R' [0s]
Running 'REMLdev.R' [3s]
Running 'ST.R' [3s]
Running 'agridat_gotway.R' [0s]
Running 'bootMer.R' [0s]
Running 'boundary.R' [0s]
Running 'confint.R' [3s]
Running 'devCritFun.R' [0s]
Running 'drop.R' [0s]
Running 'drop1contrasts.R' [0s]
Running 'dynload.R' [1s]
Running 'elston.R' [0s]
Running 'evalCall.R' [0s]
Running 'extras.R' [3s]
Running 'falsezero_dorie.R' [0s]
Running 'fewlevels.R' [0s]
Running 'getME.R' [0s]
Running 'glmer-1.R' [3s]
Running 'glmerControlPass.R' [0s]
Running 'glmerWarn.R' [0s]
Running 'glmmExt.R' [3s]
Running 'glmmWeights.R' [0s]
Running 'hatvalues.R' [0s]
Running 'is.R' [0s]
Running 'lmList-tst.R' [3s]
Running 'lme4_nlme.R' [3s]
Running 'lmer-0.R' [3s]
Running 'lmer-1.R' [3s]
Running 'lmer-conv.R' [3s]
Running 'lmer2_ex.R' [3s]
Running 'methods.R' [0s]
Running 'minval.R' [3s]
Running 'modFormula.R' [0s]
Running 'nbinom.R' [0s]
Running 'nlmer-conv.R' [3s]
Running 'nlmer.R' [3s]
Running 'offset.R' [3s]
Running 'optimizer.R' [3s]
Running 'polytomous.R' [3s]
Running 'prLogistic.R' [3s]
Running 'predict_basis.R' [4s]
Running 'predsim.R' [3s]
Running 'priorWeights.R' [0s]
Running 'priorWeightsModComp.R' [3s]
Running 'profile-tst.R' [3s]
Running 'refit.R' [3s]
Running 'resids.R' [4s]
Running 'respiratory.R' [3s]
Running 'simulate.R' [3s]
Running 'test-glmernbref.R' [5s]
Running 'testOptControl.R' [3s]
Running 'testcolonizer.R' [3s]
Running 'testcrab.R' [3s]
Running 'throw.R' [3s]
Running 'varcorr.R' [3s]
Running 'vcov-etc.R' [3s]
Running the tests in 'tests/AAAtest-all.R' failed.
Complete output:
> if (base::require("testthat", quietly = TRUE)) {
+ pkg <- "lme4"
+ require(pkg, character.only=TRUE, quietly=TRUE)
+ if(getRversion() < "3.5.0") { withAutoprint <- identity ; prt <- print } else { prt <- identity }
+ if(Sys.getenv("USER") %in% c("maechler", "bbolker")) withAutoprint({
+ ## for developers' sake:
+ lP <- .libPaths() # ---- .libPaths() : ----
+ prt(lP)
+ ## ---- Entries in .libPaths()[1] : ----
+ prt(list.files(lP[1], include.dirs=TRUE))
+ prt(sessionInfo())
+ prt(packageDescription("Matrix"))
+ ## 'lme4' from packageDescription "file" :
+ prt(attr(packageDescription("lme4"), "file"))
+ })
+ test_check(pkg)
+ ##======== ^^^
+ print(warnings()) # TODO? catch most of these by expect_warning(..)
+ } else {
+ cat( "package 'testthat' not available, cannot run unit tests\n" )
+ }
Saving _problems/test-resids-82.R
[ FAIL 1 | WARN 0 | SKIP 8 | PASS 1763 ]
══ Skipped tests (8) ═══════════════════════════════════════════════════════════
• On CRAN (1): 'test-eval.R:2:1'
• On Windows (1): 'test-methods.R:345:3'
• Skipping (1): 'test-covariance_nlmer.R:22:3'
• empty test (3): , ,
• getRversion() < "4.6.0" is TRUE (1): 'test-methods.R:899:3'
• testLevel < 2 is TRUE (1): 'test-predict.R:689:3'
══ Failed tests ════════════════════════════════════════════════════════════════
── Failure ('test-resids.R:81:5'): weighted residuals ──────────────────────────
Expected `head(weighted.residuals(gm1), 3)` to equal `structure(...)`.
Differences:
2/3 mismatches (average diff: 0.267)
[2] -5.17 - -4.93 == -0.244
[3] -5.65 - -5.35 == -0.291
[ FAIL 1 | WARN 0 | SKIP 8 | PASS 1763 ]
Error:
! Test failures.
Execution halted
Flavor: r-oldrel-windows-x86_64