[R] How to estimate the residual SD for each sample separately in mixed-effects model?

[Ben Bolker]

Michal Figurski <figurski <at> mail.med.upenn.edu> writes:

> I am developing a Mixed-Effects model for a study of immunoassays using 
> 'lme4' library. The study design is as follows: 10 samples were run 
> using 7 different immunoassays, 3 times each, in duplicates. Data 
> attached. I have developed the following model:
> 
> c.lme <- lmer(Result~SPL + (SPL|Assay/Run) -1, data=data)
> 
> This model has excellent predictions - the Rsquared of the predicted vs 
> measured results is almost 1, with very small RMSE. However, I am not 
> interested in the estimates of the mean, but in SDs from the model.
> 
> I access the SDs using b<-VarCorr(c.lme). There:
>   - the 'attr(b$Assay, "stddev")' is the assay-to-assay SD component for 
> each sample (SDaa)
>   - the 'attr(b$Run, "stddev")' is the run-to-run component (SDrr)
>   - the 'attr(b, "sc")' i.e. the residual (SDres), would be the 
> within-run component, but it's a single number for all the samples.
> 
> * The problem:
>   - how to estimate the 'within-run' component (SDres) for each sample 
> separately, as the two other components?
> 
> * Solutions tried:
>   - subtracting SDaa and SDrr from total SD - sometimes produces 
> negative results
>   - adding SDres to SDaa + SDrr is usually greater than total SD

  A couple of thoughts:

* if you're going to add and subtract variance components, you ought
to do it on the variance scale rather than the standard deviation scale.
(i.e. to get the standard deviation of (eps_1 + eps_2) you need
sqrt(Var(eps_1)+Var(eps_2)) rather than SD(eps_1)+SD(eps_2).

* this model assumes that the within-run component is the SAME
for all assays.  If you want to allow the residual variation
to be different for different assays you need to use a model that
allows for heteroscedasticity in the residuals. For the foreseeable
future, this can best be done in the older lme (in the nlme package)
by specifying a weights= argument such as (I think) 
weights=varIdent(form=~1|Assay) [or something like that].

* I would be extremely interested in any pointers from other readers
on ways of specifying that variance in the random-effect variances
(i.e. not just the residual variance) varies among groups, treatments,
etc.. I've poked through various materials (esp. Pinheiro and Bates
2000) and not been able yet to figure out how to do this -- hopefully
there is some simple trick I'm missing.

* Questions on mixed models are best addressed to the mixed models
special-interest mailing list,  r-sig-mixed-models at r-project.org

  Ben Bolker

future)