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