1library(lme4)
2library(testthat)
3
4## use old (<=3.5.2) sample() algorithm if necessary
5if ("sample.kind" %in% names(formals(RNGkind))) {
6    suppressWarnings(RNGkind("Mersenne-Twister", "Inversion", "Rounding"))
7}
8
9context("lmList")
10test_that("basic lmList", {
11    set.seed(17)
12    fm1. <- lmList(Reaction ~ Days | Subject, sleepstudy, pool=FALSE)
13    fm1  <- lmList(Reaction ~ Days | Subject, sleepstudy)
14    cf.fm1 <- data.frame(
15        `(Intercept)` =
16        c(244.19267, 205.05495, 203.48423, 289.68509, 285.73897, 264.25161,
17          275.01911, 240.16291, 263.03469, 290.10413, 215.11177, 225.8346,
18          261.14701, 276.37207, 254.96815, 210.44909, 253.63604, 267.0448),
19    Days =
20        c(21.764702, 2.2617855, 6.1148988, 3.0080727, 5.2660188, 9.5667679,
21          9.1420455, 12.253141, -2.8810339, 19.025974, 13.493933, 19.504017,
22          6.4334976, 13.566549, 11.348109, 18.056151, 9.1884448, 11.298073))
23    expect_equal(signif(coef(fm1), 8), cf.fm1,
24		    tolerance = 1e-7, check.attributes=FALSE)
25    expect_equal(coef(fm1.), coef(fm1))
26    expect_true(inherits(formula(fm1), "formula")) ## <- had been wrong till 2015-04-09
27
28    sm1. <- summary(fm1.)
29    sm1 <- summary(fm1)
30    expect_equal(sm1$RSE, 25.5918156267, tolerance = 1e-10)
31    cf1 <- confint(fm1)
32
33
34    ## Calling the plot.lmList4.confint() method :
35    expect_true(inherits(pcf1 <- plot(cf1), "trellis"))
36
37})
38
39test_that("orthodont", {
40    data(Orthodont, package="nlme")
41    fm2 <- lmList(distance ~ age | Subject, Orthodont)
42    fe2 <- fixef(fm2)
43    expect_equal(fe2, c("(Intercept)" = 16.7611111111111,
44                               age = 0.660185185185185))
45    expect_true(inherits(pairs(fm2), "trellis"))
46})
47
48test_that("simulated", {
49    set.seed(12)
50    d <- data.frame(
51        g = sample(c("A","B","C","D","E"), 250, replace=TRUE),
52        y1 = runif(250, max=100),
53        y2 = sample(c(0,1), 250, replace=TRUE)
54    )
55
56    fm3.1 <- lmList(y1 ~ 1 | g, data=d)
57    expect_equal(coef(fm3.1),
58    structure(list(`(Intercept)` = c(45.8945525606396, 50.1127995110841,
59                 49.5320538515225, 52.4286874305165, 48.7716343882989)),
60              .Names = "(Intercept)",
61              row.names = c("A",
62                            "B", "C", "D", "E"), class = "data.frame",
63              label = "Coefficients", effectNames = "(Intercept)",
64              standardized = FALSE))
65
66    cf31 <- confint(fm3.1)
67    expect_true(inherits(plot(cf31), "trellis"))
68
69    fm3.2 <- lmList(y2 ~ 1 | g, data=d, family=binomial)
70    ##                                ^^^^^^^^ "glmList"
71    cf32 <- suppressMessages(confint(fm3.2,quiet=TRUE))
72    expect_identical(dim(cf32), c(5L,2:1))
73    expect_true(inherits(plot(cf32), "trellis"))
74    expect_equal(unname(getDataPart(signif(drop(cf32), 6))),
75	    cbind(c(-0.400041, -0.311489, -1.07774, -0.841075, -0.273828),
76	  c( 0.743188,  0.768538, 0.0723138, 0.274392,  0.890795)))
77
78
79})
80
81test_that("cbpp", {
82    ## "glmList" (2) -- here,  herd == 8 has only one observation => not estimable
83    expect_warning(fm4 <- lmList(cbind(incidence, size - incidence) ~ period | herd,
84                                 family=binomial, data=cbpp),
85                   "Fitting failed for ")
86
87    cf4 <- coef(fm4) # with some 5 NA's
88    ## match NA locations
89    expect_equal(dim(cf4),c(15,4))
90    expect_identical(which(is.na(cf4)),
91                     sort(as.integer(c(8+15*(0:3), 47))))
92
93    expect_warning(fm4B <- lme4::lmList(incidence ~ period | herd,
94                                        data=cbpp),
95                   "Fitting failed")
96
97
98
99if(FALSE) {
100        ## FIXME: this is actually an nlme bug ...
101        ## https://bugs.r-project.org/bugzilla/show_bug.cgi?id=16542
102        try(summary(fm4))
103        ## Error in `[<-`(`*tmp*`, use, use, ii, value = lst[[ii]]) :
104        ##   subscript out of bounds
105        library(nlme)
106        data("cbpp",package="lme4")
107        fm6 <- nlme::lmList(incidence ~ period | herd, data=cbpp)
108        try(coef(fm6))  ## coef does *not* work here
109        try(summary(fm6))
110
111        ## this is a slightly odd example because the residual df from
112        ##  these fits are in fact zero ...  so pooled.SD fails, as it should
113
114    }
115})
116
117test_that("NA,weights,offsets", {
118
119    ## from GH #320
120    set.seed(101)
121    x <- 1:8
122    y <- c(2,2,5,4,3,1,2,1)
123    g <- c(1,1,1,2,2,3,3,3)
124    dat <- data.frame(x=x, y=y, g=g)
125    m1 <- lmList(y ~ x | g, data=dat)
126    expect_false(any(is.na(coef(m1))))
127    w <- runif(nrow(sleepstudy))
128    m2 <- lmList(Reaction ~ Days | Subject,
129                 weights=w, sleepstudy)
130    ss <- subset(sleepstudy,Subject==levels(Subject)[1])
131    m2X <- lm(Reaction ~ Days, ss, weights=w[1:nrow(ss)])
132    expect_equal(coef(m2X),as.matrix(coef(m2))[1,])
133    m3 <- lmList(Reaction ~ Days | Subject, sleepstudy)
134    m4 <- lmList(Reaction ~ Days | Subject,
135                 offset=w, sleepstudy)
136    m4X <- lm(Reaction ~ Days, ss, offset=w[1:nrow(ss)])
137    expect_equal(coef(m4X),as.matrix(coef(m4))[1,])
138    expect_false(identical(m2,m3))
139    expect_false(identical(m4,m3))
140    m5 <- lmList(Reaction ~ Days + offset(w) | Subject, sleepstudy)
141    expect_equal(coef(m5),coef(m4))
142
143    ## more from GH 320
144
145    dat2 <- data.frame(dat,xx=c(NA,NA,NA,1:4,NA))
146    m5 <- lmList(y ~ x | g, data=dat2)
147    expect_equal(unlist(coef(m5)[1,]),
148                 coef(lm(y~x,subset=(g==1))))
149    expect_equal(unlist(coef(m5)[3,]),
150                 coef(lm(y~x,subset=(g==3))))
151})
152
153test_that("pooled",
154{
155    ## GH #26
156    fm_lme4 <- lme4:::lmList(Reaction ~ Days | Subject, sleepstudy)
157    fm_nlme <- nlme:::lmList(Reaction ~ Days | Subject, sleepstudy)
158    fm_nlme_nopool <- nlme:::lmList(Reaction ~ Days | Subject, sleepstudy, pool=FALSE)
159    ci_lme4_pooled <- confint(fm_lme4,pool=TRUE) #get low and high CI estimates and pooled sd
160    ci_nlme_pooled <- nlme:::intervals(fm_nlme,pool=TRUE)
161    expect_equal(unname(ci_lme4_pooled[,,1]),unname(ci_nlme_pooled[,c(1,3),1]))
162    ci_lme4_nopool1 <- confint(fm_lme4,pool=FALSE)
163    ci_lme4_nopool2 <- confint(fm_lme4)
164    expect_identical(ci_lme4_nopool1,ci_lme4_nopool2)
165    ## BUG in nlme::intervals ... ? can't get CIs on unpooled fits
166    ## nlme::intervals(fm_nlme,pool=FALSE)
167    ## nlme::intervals(fm_nlme_nopool)
168    expect_equal(ci_lme4_nopool1[1:3,,1],
169                 structure(c(179.433862895996, 193.026448122379, 186.785722998616,
170                             308.951475285822, 217.083442786712, 220.182727910474),
171                           .Dim = c(3L, 2L), .Dimnames = list(c("308", "309", "310"),
172                                                              c("2.5 %", "97.5 %"))))
173})
174
175test_that("derived variables",
176          {
177              fm_lme4 <- lme4:::lmList(log(Reaction) ~ Days | Subject, sleepstudy)
178              fm_nlme <- nlme:::lmList(log(Reaction) ~ Days | Subject, sleepstudy)
179              expect_equal(c(coef(fm_lme4)),c(coef(fm_nlme)),tolerance=1e-5)
180          })
181
182test_that("subset", {
183    data(MathAchieve, package="nlme")
184    data(MathAchSchool, package="nlme")
185    RB <- merge(MathAchieve, MathAchSchool[, c("School", "Sector")],
186                by="School")
187    names(RB) <- tolower(names(RB))
188    RB$cses <- with(RB, ses - meanses)
189    cat.list.nlme <- nlme::lmList(mathach ~ cses | school,
190                  subset = sector=="Catholic",
191                  data=RB)
192    cat.list.lme4 <- lme4::lmList(mathach ~ cses | school,
193                                  subset = sector=="Catholic", data=RB)
194    expect_equal(c(coef(cat.list.lme4)),
195                 c(coef(cat.list.nlme)),tolerance=1e-5)
196})
197
198