Hi-
I want to fit a model with crossed random effects and heteroskedastic
level-1 errors where inferences about fixed effects are of primary
interest. The dimension of the random effects is making the model
computationally prohibitive using lme() where I could model the
heteroskedasticity with the "weights" argument. I am aware that the weights
argument to lmer() cannot be used to estimate a heteroskedastic variance
function with unknown parameters. However my data situation is such that I
am able to assign each unit uniquely to one of four groups and I am willing
to treat the relative error variances of the groups as known. I.e. I would
be able to specify the model I want to fit in lme() using a weights
statement of the form "weights = varIdent(form = ~1 | group, fixed =
knownSDratios)". I am trying to figure out how I can concoct a "weights"
argument for lmer() that will fit this same model. In particular, since the
reported inferences about fixed effects parameters are not invariant to how
the weights passed to lmer() are scaled, I am trying to understand how to
construct a "weights" argument that not only contains the information about
the relative precisions of the errors in the different groups, but also
provides valid inferences about the fixed effects. In a simple example
below with balanced, nested data, scaling the weights to sum to the number
of cases in the data makes lmer() and lme() agree on the inferences about
the fixed effects. However, at the end of the example, it is also clear
that this correspondence does not hold when the data are not balanced. My
question boils down to this: is there a general way to scale lmer()'s
"weights" argument that will cause it to always correspond to what lme()
would report about fixed effects inferences when passed the same fixed
relative precision information?
Thanks for any advice.
J.R.
###############################################################################
## EXAMPLE: generate balanced nested data with heteroskedastic errors of
## variance 0.5, 1, 2, or 4
###############################################################################
ngroup <- 50
npergroup <- 20
n <- ngroup*npergroup
set.seed(5046)
d <- data.frame(unitid = 1:n, groupid = rep(1:ngroup,
each=npergroup), verror = sample(c(0.5, 1, 2, 4), size=n, replace=T),
x = 0.1*rnorm(n))
groupeffx <- data.frame(groupid = 1:ngroup, theta = rnorm(ngroup, sd = 0.25))
d <- merge(d, groupeffx)
d$etrue <- rnorm(n, sd = sqrt(d$verror))
d$y <- 5 + d$x + d$theta + d$etrue
d$verrorf <- factor(paste("v",d$verror,sep=""))
print(tapply(d$etrue, d$verrorf, var))
## function to collect pieces from lme() output
sumLME <- function(o){
tab <- summary(o)$tTable
r <- c(tab[1,1:2], tab[2,1:2], as.numeric(VarCorr(o))[3:4])
names(r) <- c("b0","b0sd","b1","b1sd","sdgroup","sdresid")
return(r)
}
## function to collect pieces from lmer() output
sumLMER <- function(o){
tab <- summary(o)@coefs
r <- c(tab[1,1:2], tab[2,1:2],
as.numeric(attributes(VarCorr(o)$groupid)$stddev),
as.numeric(attr(VarCorr(o), "sc")))
names(r) <- c("b0","b0sd","b1","b1sd","sdgroup","sdresid")
return(r)
}
res <- vector(0, mode="list")
library(nlme)
## lme, ignoring heteroskedasticity
res[["lme.nohet"]] <- sumLME( lme(fixed = y ~ x, data = d, random = ~1
| groupid) )
## lme, heteroskedastic model with variance weights estimated
res[["lme.esthet"]] <- sumLME( lme(fixed = y ~ x, data = d, random =
~1 | groupid, weights = varIdent(form = ~1|verrorf)) )
## lme, heteroskedastic model with true fixed weights
v <- c(v0.5 = 0.5, v1 = 1, v2 = 2, v4 = 4)
sdrats <- sqrt(v / v[1])[-1]
res[["lme.fixedhet"]] <- sumLME( lme(fixed = y ~ x, data = d, random =
~1 | groupid, weights = varIdent(form = ~1|verrorf, fixed = sdrats)) )
detach("package:nlme")
library(lme4)
## lmer, ignoring heteroskedasticity
res[["lmer.nohet"]] <- sumLMER( lmer(y ~ x + (1|groupid), data = d) )
## essentially matches res[["lme.nohet"]], makes sense
## lmer, with fixed weights equal to known precisions
d$w1 <- 1 / d$verror
res[["lmer.w1"]] <- sumLMER( lmer(y ~ x + (1|groupid), data = d, weights = w1) )
## only b0 and b1 matches res[["lme.fixedhet"]]
## lmer, with fixed weights proportional to known precisions but weights sum to
## number of records as they would with the unweighted case.
d$w2 <- nrow(d) * d$w1 / sum(d$w1)
res[["lmer.w2"]] <- sumLMER( lmer(y ~ x + (1|groupid), data = d, weights = w2) )
## this is extremely close to res[["lme.fixedhet"]] except for sdresid
do.call("rbind", res)
## is the case that matched only because of balance?
dsub <- d[sample(1:nrow(d), size= nrow(d)/2, replace=F),]
dsub$w2 <- nrow(dsub) * dsub$w1 / sum(dsub$w1)
sumLMER( lmer(y ~ x + (1|groupid), data = dsub, weights = w2) )
detach("package:lme4")
library(nlme)
sumLME( lme(fixed = y ~ x, data = dsub, random = ~1 | groupid, weights
= varIdent(form = ~1|verrorf, fixed = sdrats)) )
## these no longer match
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
