lme4: constrain sigma to 0
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