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

[Vincent Dorie]

I'm not familiar with how lme works, but if the lmer formulation is sufficient then you can fit models of that sort using blmer in the blme package. Simply fix the common residual variance to one and supply the weights, at which point they'll behave as desired.

For example,

library(blme)
blmer(formula = y ~ 1 + (1 | group), weights = V,
      resid.prior = point(1.0), cov.prior = NULL)

Vince

On Feb 22, 2014, at 6:02 PM, Zaslavsky, Alan M. wrote:

> 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