[R-sig-ME] Specifying fixed residual variances for multilevel models

[Viechtbauer Wolfgang (STAT)]

In meta-analytic models, we commonly fix the residual variances to (at least approximately) known values. Following your example, you could fit the desired models with the 'metafor' package with:

library(metafor)

rma.mv(y, V, random = ~ 1 | group)
rma.mv(y, V2, random = ~ 1 | group)

This will really fix those residual variances to V and V2 (not just up to that proportionality constant).

At the moment, rma.mv() isn't optimized for speed, so it takes a bit to get the model to converge. The next version of metafor (to be released in the near future) will be faster.

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 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
> observations
> > 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.021891/.5899889)^2
> [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
> > (.7569/.437)^2
> [1] 2.999951
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models