unit tests using nlme instead of lme4. varcomp works with both lme, lmer

kkholst · Dec 20, 2023 · 9b92b51 · 9b92b51
1 parent 8e79ba2
commit 9b92b51
Show file tree

Hide file tree

Showing 2 changed files with 38 additions and 23 deletions.
diff --git a/R/lmers.R b/R/lmers.R
@@ -48,20 +48,32 @@ lmerplot <- function(model,x,id,y,transform,re.form=NULL,varcomp=FALSE,colorbar=
  }
 }
 
-varcomp <- function(x,profile=TRUE,...) {
- cc <- cbind(lme4::fixef(x),diag(as.matrix(vcov(x)))^.5)
- cc <- cbind(cc,cc[,1]-qnorm(0.975)*cc[,2],cc[,1]+qnorm(0.975)*cc[,2],
- 2*(1-pnorm(abs(cc[,1])/cc[,2])))
- pr <- NULL
- if (profile) pr <- confint(x)
- colnames(cc) <- c("Estimate","Std.Err","2.5%","97.5%","p-value")
- vc <- lme4::VarCorr(x)
- res <- structure(list(coef=lme4::fixef(x), vcov=as.matrix(vcov(x)),
- coefmat=cc,
- confint=pr,
- varcomp=vc[[1]][,],
- residual=attributes(vc)$sc^2
- ),
- class="estimate.lmer")
- res
+varcomp <- function(x, profile=TRUE, ...) {
+ if (inherits(x, "lme")) {
+ coefs <- nlme::fixef(x)
+ varcomp <- nlme::getVarCov(x)
+ pr <- nlme::intervals(x)
+ resvar <- summary(x)$sigma
+ } else if (inherits(x, "lmerMod")) {
+ coefs <- lme4::fixef(x)
+ varcomp <- lme4::VarCorr(x)[[1]][, ]
+ resvar <- attributes(VarCorr(x))$sc^2
+ pr <- confint(x) ## Profile likelihood intervals
+ } else {
+ stop("unrecognized model object")
+ }
+ cc <- cbind(coefs, diag(as.matrix(vcov(x)))^.5)
+ cc <- cbind(cc,cc[,1]-qnorm(0.975)*cc[,2],
+ cc[,1]+qnorm(0.975)*cc[,2],
+ 2*(1-pnorm(abs(cc[,1])/cc[,2])))
+ colnames(cc) <- c("Estimate","Std.Err","2.5%","97.5%","p-value")
+ res <- structure(list(coef=coefs,
+ vcov=as.matrix(vcov(x)),
+ coefmat=cc,
+ confint=pr,
+ varcomp=varcomp,
+ varresidual=resvar
+ ),
+ class="estimate.lmer")
+ res
 }
diff --git a/tests/testthat/test-inference.R b/tests/testthat/test-inference.R
@@ -299,22 +299,24 @@ test_that("Prediction with missing data, random intercept", {
 })
 
 
-if (requireNamespace("lme4",quietly = TRUE) && requireNamespace("mets",quietly = TRUE))
+## if (requireNamespace("lme4", quietly = TRUE) && requireNamespace("mets", quietly = TRUE)) {
+if (requireNamespace("mets",quietly = TRUE))
 test_that("Random slope model", {
  set.seed(1)
  m <- lvm()
- regression(m) <- y1 ~ 1*u+1*s
- regression(m) <- y2 ~ 1*u+2*s
- regression(m) <- y3 ~ 1*u+3*s
+ regression(m) <- y1 ~ 1*u + 1*s
+ regression(m) <- y2 ~ 1*u + 2*s
+ regression(m) <- y3 ~ 1*u + 3*s
  latent(m) <- ~u+s
  dw <- sim(m,20)
 
  dd <- mets::fast.reshape(dw)
  dd$num <- dd$num+runif(nrow(dd),0,0.2)
  dd0 <- dd[-c(1:2*3),]
  library(lme4)
- l <- lmer(y~ 1+num +(1+num|id),dd,REML=FALSE)
- sl <- lava:::varcomp(l,profile=FALSE)
+ l <- lme(y ~ 1+num, random=~1+num|id, data=dd, method="ML")
+ ##l0 <- lmer(y~ 1+num +(1+num|id),dd,REML=FALSE)
+ sl0 <- lava:::varcomp(l)
 
  d <- mets::fast.reshape(dd,id="id")
  d0 <- mets::fast.reshape(dd0,id="id")
@@ -337,7 +339,8 @@ test_that("Random slope model", {
  ## missing
  testthat::expect_output(e0 <- estimate(m0,d0,missing=TRUE,param="none",control=list(method="NR",constrain=FALSE,start=coef(e),trace=1)),
  "Iter=")
- l0 <- lmer(y~ 1+num +(1+num|id),dd0,REML=FALSE)
+ ## l0 <- lmer(y ~ 1 + num + (1 + num | id), dd0, REML = FALSE)
+ l0 <- lme(y~ 1+num, random=~1+num|id, data=dd0, method="ML")
  testthat::expect_true((logLik(l0)-logLik(e0))^2<1e-5)
 
  m1 <- lvm(c(y1[0:v],y2[0:v],y3[0:v])~1*u)