[R-sig-ME] lme4: constrain sigma to 0
Viechtbauer Wolfgang (STAT)
wolfgang.viechtbauer at maastrichtuniversity.nl
Wed Feb 12 17:54:07 CET 2014
If you are really masochistic, you could trick the rma.mv() function from the 'metafor' package to fit such a model (in meta-analyses, you fix the variances of the residuals to known values -- and here, we could set those values to 0).
As an example, let me follow up on what Rolf Turner wanted to do:
> I wanted to be able to fit a fairly simple repeated measures model with
> the covariance over time being an "arbitrary" n-x-n positive definite
> matrix (where "n" is the number of time points). This can't be done
> using lmer() unless one can constrain the "overall" variance to be zero,
> otherwise the diagonal of the covariance matrix is un-identifiable.
##################################################################
library(metafor)
library(nlme)
library(lme4)
data(Orthodont)
V0 <- rep(0, nrow(Orthodont))
res1 <- rma.mv(distance ~ age, V=V0, random = ~ age | Subject, struct="UN", data=Orthodont, control=list(optimizer="nlminb", REMLf=FALSE))
summary(res1)
res2 <- lme(distance ~ age, random = ~ factor(age) - 1 | Subject, data=Orthodont)
summary(res2)
res3 <- gls(distance ~ age, correlation = corSymm(form = ~ 1 | Subject), weights = varIdent(form = ~ 1 | age), data=Orthodont)
summary(res3)
res4 <- lmer(distance ~ age + (factor(age) - 1 | Subject), data=Orthodont)
summary(res4)
# Note that lme() and lmer() also fit this model, but here the residual variance is being estimated.
# As mentioned earlier by Rolf Turner, that model is actually over-parameterized. Marginally, you end
# up with the same fit though:
logLik(res1)
logLik(res2)
logLik(res3)
logLik(res4)
# So, lme() and lmer() are splitting up the variance along the diagonal into variance due to the random
# effects and residual variance. Added up, this yields the same values as with rma.mv() and gls():
round(res1$tau2, 4)
round(as.numeric(VarCorr(res2)[1:4,1]) + res2$sigma^2, 4)
round(unname(c(res3$sigma^2, (res3$sigma * coef(res3$modelStruct, unconstrained=FALSE)[7:9])^2)), 4)
round(unname(diag(VarCorr(res4)[[1]]) + sigma(res4)^2), 4)
# But the split is arbitrary and depends on how the optimization is done:
res2$sigma
sigma(res4)
##################################################################
A few more notes:
1) I had to switch to 'nlminb' for rma.mv() since the default ('optim' with 'Nelder-Mead') did not converge and gave strange results when it did (by increasing the number of iterations).
2) With REMLf=FALSE, a particular constant term is left out from the REML likelihood that is also not included in the likelihood computation for the other functions.
3) rma.mv() is slow as it was never intended for this. But it does work. In other cases, it may not.
Best,
Wolfgang
--
Wolfgang Viechtbauer, Ph.D., Statistician
Department of Psychiatry and Psychology
School for Mental Health and Neuroscience
Faculty of Health, Medicine, and Life Sciences
Maastricht University, P.O. Box 616 (VIJV1)
6200 MD Maastricht, The Netherlands
+31 (43) 388-4170 | http://www.wvbauer.com
> -----Original Message-----
> From: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-
> models-bounces at r-project.org] On Behalf Of Vincent Dorie
> Sent: Wednesday, February 12, 2014 14:35
> To: r-sig-mixed-models
> Subject: Re: [R-sig-ME] lme4: constrain sigma to 0
>
> It was also in an email from Ben a few messages up the chain that blme
> works for this, just not for sigma = 0. e.g.
>
> blmer(formula, data, cov.prior = NULL, resid.prior = point(1.0))
>
>
> On Feb 12, 2014, at 3:37 AM, David Duffy wrote:
>
> > For those who would like to fit these models, you could muck around
> with the OpenMx package, which allows you to fix or constrain any
> parameters you like ;) To fix the residuals to zero, change free to
> FALSE in the
> > "residual variances" bit. Parameters with the same "labels" are
> automatically constrained equal.
> >
> > #
> > # Based on an example in the OpenMX documentation
> > #
> > require(lme4)
> > data(sleepstudy)
> > sleep2 <- reshape(sleepstudy, direction="wide", idvar="Subject",
> timevar="Days")
> > rownames(sleep2) <- sleep2$Subject
> > sleep2 <- sleep2[, -1]
> > names(sleep2) <- gsub("\\.","",names(sleep2)) # Mx no like "."
> >
> > require(OpenMx)
> > growthCurveModel <- mxModel("Sleepstudy as Linear Growth Curve Model",
> > type="RAM",
> > mxData(
> > observed=sleep2,
> > type="raw"
> > ),
> > manifestVars=names(sleep2),
> > latentVars=c("intercept","slope"),
> > # residual variances
> > mxPath(
> > from=names(sleep2),
> > arrows=2,
> > free=TRUE,
> > values = rep(1, ncol(sleep2)),
> > labels=rep("residual", ncol(sleep2))
> > ),
> > # latent variances and covariance
> > mxPath(
> > from=c("intercept","slope"),
> > arrows=2,
> > connect="unique.pairs",
> > free=TRUE,
> > values=c(1, 1, 1),
> > labels=c("vari", "cov", "vars")
> > ),
> > # intercept loadings
> > mxPath(
> > from="intercept",
> > to=names(sleep2),
> > arrows=1,
> > free=FALSE,
> > values = rep(1, ncol(sleep2))
> > ),
> > # slope loadings
> > mxPath(
> > from="slope",
> > to=names(sleep2),
> > arrows=1,
> > free=FALSE,
> > values=seq(0,9)
> > ),
> > # manifest means
> > mxPath(from="one",
> > to=names(sleep2),
> > arrows=1,
> > free=FALSE,
> > values = rep(0, ncol(sleep2))
> > ),
> > # latent means
> > mxPath(from="one",
> > to=c("intercept", "slope"),
> > arrows=1,
> > free=TRUE,
> > values=c(1, 1),
> > labels=c("meani", "means")
> > )
> > ) # close model
> >
> > growthCurveFit <- mxRun(growthCurveModel)
> > print(summary(growthCurveFit))
> > print(growthCurveFit at output$estimate)
> >
> >
> >
> > | David Duffy (MBBS PhD)
> > | email: David.Duffy at qimrberghofer.edu.au ph: INT+61+7+3362-0217 fax:
> -0101
> > | Genetic Epidemiology, QIMR Berghofer Institute of Medical Research
> > | 300 Herston Rd, Brisbane, Queensland 4006, Australia GPG 4D0B994A
> >
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list