tests/priorWeights.R

In lme4: Linear Mixed-Effects Models using 'Eigen' and S4

## 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)