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

[Zaslavsky, Alan M.]

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