-----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 Zaslavsky, Alan M.
Sent: Sunday, February 23, 2014 00:02
To: r-sig-mixed-models at r-project.org
Subject: [R-sig-ME] Specifying fixed residual variances for multilevel
models
We have a bunch of analyses in which we would like to specify residual
variances that are calculated outside the multilevel modeling function.
We were proceeding on the assumption that using varFixed() for the
weights argument to lme() would do this, but after doing a few
simulations it became evident that this was only fixing the residual
variances up to a proportionality constant and in fact rescaling the
weights had no effect on estimates of random effects variances.
I've spent a few hours scanning the documentation, list comments and FAQs
without tracking down a solution. This seems to be more a feature than a
bug, although comments I see on this list and Github, among other places
suggest that the handling of residual variances or weights is the subject
of current development efforts (lme4a?). Not surprising --- in fact the
terms 'weights' has at least as many meanings as 'fixed effects', with
potentially different analytic implications (especially for mixed
models). I would urge talking explicitly about residual variance models,
which have a clear interpretation.
Any thoughts?
Alan Zaslavsky
zaslavsk at hcp.med.harvard.edu
NLME EXAMPLE
group=rep(1:40,10) ## 40 groups of 10
V=401:800/400 ## residual variances
y=rnorm(40)[group]+rnorm(400)*sqrt(V)
lme(fixed=y~1,random=~1|group) ## unweighted
Linear mixed-effects model fit by REML
Data: NULL
Log-restricted-likelihood: -696.8343
Fixed: y ~ 1
(Intercept)
0.3656883
Random effects:
Formula: ~1 | group
(Intercept) Residual
StdDev: 0.9231727 1.257827
Number of Observations: 400
Number of Groups: 40
lme(fixed=y~1,random=~1|group,weight=varFixed(~V))
Linear mixed-effects model fit by REML
Data: NULL
Log-restricted-likelihood: -692.0942
Fixed: y ~ 1
(Intercept)
0.3513349
Random effects:
Formula: ~1 | group
(Intercept) Residual
StdDev: 0.9261319 1.021891
Variance function:
Structure: fixed weights
Formula: ~V
Number of Observations: 400
Number of Groups: 40
## the one above specified the residual variances that were used in
generating the
## the data and there estimated scale factor close to 1, as it should
V2=3*V
lme(fixed=y~1,random=~1|group,weight=varFixed(~V2))
Linear mixed-effects model fit by REML
Data: NULL
Log-restricted-likelihood: -692.0942
Fixed: y ~ 1
(Intercept)
0.3513349
Random effects:
Formula: ~1 | group
(Intercept) Residual
StdDev: 0.9261319 0.5899889
Variance function:
Structure: fixed weights
Formula: ~V2
Number of Observations: 400
Number of Groups: 40
### here we increased the specified residual variance by a factor of 3,
but it was ignored except that the residual scale factor went down by a
factor of sqrt(3) as shown below. All very good if you only know
*relative* residual variances but not helpful if you know them
*absolutely*
[1] 3.000001
LMER (same data)
lmer(formula=y~1+(1|group),weights=V)
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ 1 + (1 | group)
REML criterion at convergence: 1570.624
Random effects:
Groups Name Std.Dev.
group (Intercept) 0.7569
Residual 1.2824
Number of obs: 400, groups: group, 40
Fixed Effects:
(Intercept)
0.3788
lmer(formula=y~1+(1|group),weights=V2)
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ 1 + (1 | group)
REML criterion at convergence: 2010.069
Random effects:
Groups Name Std.Dev.
group (Intercept) 0.437
Residual 1.282
Number of obs: 400, groups: group, 40
Fixed Effects:
(Intercept)
0.3788
############ even stranger -- rescaling residual variance leaves
residual var estimate alone, but changes level-2 variance by factor of 3