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