1if (.Platform$OS.type != "windows") { 2 3 require(lme4) 4 source(system.file("test-tools-1.R", package = "Matrix"))# identical3() etc 5 6 ## use old (<=3.5.2) sample() algorithm if necessary 7 if ("sample.kind" %in% names(formals(RNGkind))) { 8 suppressWarnings(RNGkind("Mersenne-Twister", "Inversion", "Rounding")) 9 } 10 11 ## Check that quasi families throw an error 12 assertError(lmer(cbind(incidence, size - incidence) ~ period + (1|herd), 13 data = cbpp, family = quasibinomial)) 14 assertError(lmer(incidence ~ period + (1|herd), 15 data = cbpp, family = quasipoisson)) 16 assertError(lmer(incidence ~ period + (1|herd), 17 data = cbpp, family = quasi)) 18 19 ## check bug found by Kevin Buhr 20 set.seed(7) 21 n <- 10 22 X <- data.frame(y=runif(n), x=rnorm(n), z=sample(c("A","B"), n, TRUE)) 23 fm <- lmer(log(y) ~ x | z, data=X) ## ignore grouping factors with 24 ## gave error inside model.frame() 25 stopifnot(all.equal(unname(fixef(fm)), -0.8345, tolerance=.01)) 26 27 ## is "Nelder_Mead" default optimizer? 28 isNM <- formals(lmerControl)$optimizer == "Nelder_Mead" 29 isOldB <- formals(lmerControl)$optimizer == "bobyqa" 30 isOldTol <- environment(nloptwrap)$defaultControl$xtol_abs == 1e-6 31 ## check working of Matrix methods on vcov(.) etc ---------------------- 32 fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy) 33 V <- vcov(fm) 34 V1 <- vcov(fm1) 35 TOL <- 0 # to show the differences below 36 TOL <- 1e-5 # for the check 37 stopifnot( 38 all.equal(diag(V), if(isNM) 0.176076 else if(isOldB) 0.176068575 else if (isOldTol) 0.1761714 else 0.1760782, 39 tolerance = TOL) 40 , 41 all.equal(as.numeric(chol(V)), if(isNM) 0.4196165 else if(isOldB) 0.41960526 else if(isOldTol) 0.4197278 else 0.4196167, 42 tolerance=TOL) 43 , 44 all.equal(diag(V1), c(46.5639, 2.39), tolerance = 40*TOL)# (for "all" algos) 45 , 46 dim(C1 <- chol(V1)) == c(2,2) 47 , 48 all.equal(as.numeric(C1), 49 c(6.82377, 0, -0.212575, 1.53127), tolerance=20*TOL)# ("all" algos) 50 , 51 dim(chol(crossprod(getME(fm1, "Z")))) == 36 52 , TRUE) 53 ## printing 54 signif(chol(crossprod(getME(fm,"Z"))), 4)# -> simple 4 x 4 sparse 55 56 showProc.time() # 57 58 ## From: Stephane Laurent 59 ## To: r-sig-mixed-models@.. 60 ## "crash with the latest update of lme4" 61 ## 62 ## .. example for which lmer() crashes with the last update of lme4 ...{R-forge}, 63 ## .. but not with version CRAN version (0.999999-0) 64 lsDat <- data.frame( 65 Operator = as.factor(rep(1:5, c(3,4,8,8,8))), 66 Part = as.factor( 67 c(2L, 3L, 5L, 68 1L, 1L, 2L, 3L, 69 1L, 1L, 2L, 2L, 3L, 3L, 4L, 5L, 70 1L, 2L, 3L, 3L, 4L, 4L, 5L, 5L, 71 1L, 2L, 2L, 3L, 3L, 4L, 5L, 5L)), 72 y = 73 c(0.34, -1.23, -2.46, 74 -0.84, -1.57,-0.31, -0.18, 75 -0.94, -0.81, 0.77, 0.4, -2.37, -2.78, 1.29, -0.95, 76 -1.58, -2.06, -3.11,-3.2, -0.1, -0.49,-2.02, -0.75, 77 1.71, -0.85, -1.19, 0.13, 1.35, 1.92, 1.04, 1.08)) 78 79 xtabs( ~ Operator + Part, data=lsDat) # --> 4 empty cells, quite a few with only one obs.: 80 ## Part 81 ## Operator 1 2 3 4 5 82 ## 1 0 1 1 0 1 83 ## 2 2 1 1 0 0 84 ## 3 2 2 2 1 1 85 ## 4 1 1 2 2 2 86 ## 5 1 2 2 1 2 87 lsD29 <- lsDat[1:29, ] 88 89 ## FIXME: rank-Z test should probably not happen in this case: 90 (sm3 <- summary(m3 <- lm(y ~ Part*Operator, data=lsDat)))# ok: some interactions not estimable 91 stopifnot(21 == nrow(coef(sm3)))# 21 *are* estimable 92 sm4 <- summary(m4 <- lm(y ~ Part*Operator, data=lsD29)) 93 stopifnot(20 == nrow(coef(sm4)))# 20 *are* estimable 94 lf <- lFormula(y ~ (1|Part) + (1|Operator) + (1|Part:Operator), data = lsDat) 95 dim(Zt <- lf$reTrms$Zt)## 31 x 31 96 c(rankMatrix(Zt)) ## 21 97 c(rankMatrix(Zt,method="qr")) ## 31 || 29 (64 bit Lnx), then 21 (!) 98 c(rankMatrix(t(Zt),method="qr")) ## 30, then 21 ! 99 nrow(lsDat) 100 fm3 <- lmer(y ~ (1|Part) + (1|Operator) + (1|Part:Operator), data = lsDat, 101 control=lmerControl(check.nobs.vs.rankZ="warningSmall")) 102 103 lf29 <- lFormula(y ~ (1|Part) + (1|Operator) + (1|Part:Operator), data = lsD29) 104 (fm4 <- update(fm3, data=lsD29)) 105 fm4. <- update(fm4, REML=FALSE, 106 control=lmerControl(optimizer="nloptwrap", 107 optCtrl=list(ftol_abs=1e-6, 108 xtol_abs=1e-6))) 109 ## summary(fm4.) 110 stopifnot( 111 all.equal(as.numeric(formatVC(VarCorr(fm4.))[,"Std.Dev."]), 112 c(1.04066, 0.63592, 0.52914, 0.48248), tol = 1e-4) 113 ) 114 115 116 showProc.time() 117 cat('Time elapsed: ', proc.time(),'\n') # for ``statistical reasons'' 118 119} ## skip on windows (for speed) 120