R Under development (unstable) (2023-11-02 r85465 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> ## use old (<=3.5.2) sample() algorithm if necessary
> if ("sample.kind" %in% names(formals(RNGkind))) {
+     suppressWarnings(RNGkind("Mersenne-Twister", "Inversion", "Rounding"))
+ }
> 
> 
> compFunc <- function(lmeMod, lmerMod, tol = 1e-2){
+     lmeVarCorr <- nlme:::VarCorr(lmeMod)[,"StdDev"]
+     lmeCoef <- summary(lmeMod)$tTable[,-c(3,5)]
+     lmeOut <- c(as.numeric(lmeVarCorr), as.numeric(lmeCoef))
+     keep <- !is.na(lmeOut)
+     lmeOut <- lmeOut[keep]
+     dn <- dimnames(lmeCoef)
+     if(is.null(dn)) dn <- list("", names(lmeCoef))
+     names(lmeOut) <- c(
+         paste(names(lmeVarCorr), "Var"),
+         as.character(do.call(outer, c(dn, list("paste")))))[keep]
+ 
+     ## get nested RE variances in the same order as nlme
+     ## FIXME:  not sure if this works generally
+     vcLmer <- VarCorr(lmerMod)
+     vcLmer <- vcLmer[length(vcLmer):1]
+     ##
+ 
+     lmerVarCorr <- c(sapply(vcLmer, attr, "stddev"),
+                      attr(VarCorr(lmerMod), "sc"))
+     ## differentiate lme4{new} and lme4.0 :
+     lmerCoef <- if(is(lmerMod, "merMod"))
+ 	summary(lmerMod)$coefficients else summary(lmerMod)@coefs
+     lmerOut <- c(lmerVarCorr, as.numeric(lmerCoef))
+     names(lmerOut) <- names(lmeOut)
+ 
+     return(list(target = lmeOut, current = lmerOut, tolerance = tol))
+ }
> 
> if (.Platform$OS.type != "windows") {
+ set.seed(1)
+ nGroups <- 100
+ nObs <- 1000
+ 
+ # explanatory variable with a fixed effect
+ explVar1 <- rnorm(nObs)
+ explVar2 <- rnorm(nObs)
+ 
+ # random intercept among levels of a grouping factor
+ groupFac <- as.factor(rep(1:nGroups,each=nObs/nGroups))
+ randEff0 <- rep(rnorm(nGroups),each=nObs/nGroups)
+ randEff1 <- rep(rnorm(nGroups),each=nObs/nGroups)
+ randEff2 <- rep(rnorm(nGroups),each=nObs/nGroups)
+ 
+ # residuals with heterogeneous variance
+ residSD <- rpois(nObs,1) + 1
+ residError <- rnorm(nObs,sd=residSD)
+ 
+ # response variable
+ respVar <- randEff0 + (1+randEff1)*explVar1 + (1+randEff2)*explVar2 + residError
+ 
+ # rename to fit models on one line
+ y <- respVar
+ x <- explVar1
+ z <- explVar2
+ g <- groupFac
+ v <- residSD^2
+ w <- 1/v
+ 
+ library("nlme")
+ lmeMods <- list(
+     ML1   = lme(y ~ x,   random = ~ 1|g,   weights = varFixed(~v), method = "ML"),
+     REML1 = lme(y ~ x,   random = ~ 1|g,   weights = varFixed(~v), method = "REML"),
+     ML2   = lme(y ~ x,   random = ~ x|g,   weights = varFixed(~v), method = "ML"),
+     REML2 = lme(y ~ x,   random = ~ x|g,   weights = varFixed(~v), method = "REML"),
+     ML1   = lme(y ~ x+z, random = ~ x+z|g, weights = varFixed(~v), method = "ML"),
+     REML2 = lme(y ~ x+z, random = ~ x+z|g, weights = varFixed(~v), method = "REML"))
+ 
+ library("lme4")
+ lmerMods <- list(
+     ML1 =   lmer(y ~ x +  (1|g),      weights = w, REML = FALSE),
+     REML1 = lmer(y ~ x +  (1|g),      weights = w, REML = TRUE),
+     ML2 =   lmer(y ~ x +  (x|g),      weights = w, REML = FALSE),
+     REML2 = lmer(y ~ x +  (x|g),      weights = w, REML = TRUE),
+     ML3 =   lmer(y ~ x + z + (x+z|g), weights = w, REML = FALSE),
+     REML3 = lmer(y ~ x + z + (x+z|g), weights = w, REML = TRUE))
+ 
+ comp <- mapply(compFunc, lmeMods, lmerMods, SIMPLIFY=FALSE)
+ stopifnot(all(sapply(comp, do.call, what = all.equal)))
+ ## Look at the relative differences:
+ sapply(mapply(compFunc, lmeMods, lmerMods, SIMPLIFY=FALSE, tol = 0),
+        do.call, what = all.equal)
+ 
+ ## add simulated weights to the sleepstudy example
+ n <- nrow(sleepstudy)
+ v <- rpois(n,1) + 1
+ w <- 1/v
+ sleepLme <- lme(Reaction ~ Days, random = ~ Days|Subject,
+                 sleepstudy, weights = varFixed(~v),
+                 method = "ML")
+ sleepLmer <- lmer(Reaction ~ Days + (Days|Subject),
+                   sleepstudy, weights = w,
+                   REML = FALSE)
+ sleepComp <- compFunc(sleepLme, sleepLmer)
+ stopifnot(do.call(all.equal, sleepComp))
+ ## look at relative differences:
+ sleepComp$tolerance <- 0
+ do.call(all.equal, sleepComp)
+ 
+ if (require("mlmRev")) {
+     n <- nrow(Chem97)
+     v <- rpois(n,1) + 1
+     w <- 1/v
+     Chem97Lme <- lme(score ~ 1, random = ~ 1|lea/school, Chem97)
+     Chem97Lmer <- lmer(score ~ (1|lea/school), Chem97)
+     Chem97Comp <- compFunc(Chem97Lme, Chem97Lmer)
+     stopifnot(do.call(all.equal, Chem97Comp))
+     ## look at relative differences:
+     Chem97Comp$tolerance <- 0
+     do.call(all.equal, Chem97Comp)
+ }
+ 
+ set.seed(2)
+ n <- 40
+ w <- runif(n)
+ x <- runif(n)
+ g <- factor(sample(1:10,n,replace=TRUE))
+ Z <- model.matrix(~g-1);
+ y <- Z%*%rnorm(ncol(Z)) + x + rnorm(n)/w^.5
+ m <- lmer(y ~ x + (1|g), weights=w, REML = TRUE)
+ 
+ ## CRAN-forbidden:
+ ## has4.0 <- require("lme4.0"))
+ has4.0 <- FALSE
+ if(has4.0) {
+     ## m.0 <- lme4.0::lmer(y ~ x + (1|g), weights=w, REML = TRUE)
+     lmer0 <- get("lmer", envir=asNamespace("lme4.0"))
+     m.0 <- lmer0(y ~ x + (1|g), weights=w, REML = TRUE)
+     dput(fixef(m.0)) # c(-0.73065400610675, 2.02895402562926)
+     dput(sigma(m.0)) # 1.73614301673377
+     dput(VarCorr(m.0)$g[1,1]) # 2.35670451590395
+     dput(unname(coef(summary(m.0))[,"Std. Error"]))
+     ## c(0.95070076853232, 1.37650858268602)
+ }
+ fixef_lme4.0 <- c(-0.7306540061, 2.0289540256)
+ sigma_lme4.0 <- 1.7361430
+ Sigma_lme4.0 <- 2.3567045
+ SE_lme4.0 <- c(0.95070077, 1.37650858)
+ if(has4.0) try(detach("package:lme4.0"))
+ 
+ stopifnot(all.equal(unname(fixef(m)), fixef_lme4.0, tolerance = 1e-3))
+           all.equal(unname(fixef(m)), fixef_lme4.0, tolerance = 0) #-> 1.657e-5
+ 
+ ## but these are not at all equal :
+ (all.equal(sigma(m),                		sigma_lme4.0, tolerance = 10^-3)) # 0.4276
+ (all.equal(as.vector(VarCorr(m)$g), 		Sigma_lme4.0, tolerance = 10^-3)) # 1.038
+ (all.equal(as.vector(summary(m)$coefficients[,2]), SE_lme4.0, tolerance = 10^-3)) # 0.4276
+ ## so, lme4.0 was clearly wrong here
+ 
+ 
+ ##' make sure models that differ only in a constant
+ ##' prior weight have identical deviance:
+ fm <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy,REML=FALSE)
+ fm_wt <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy, weights = rep(5, nrow(sleepstudy)),REML=FALSE)
+ all.equal(deviance(fm), deviance(fm_wt))
+ } ## skip on windows (for speed)
> 
> proc.time()
   user  system elapsed 
   0.31    0.04    0.29