R help - Apr 2010 - How to estimate the residual SD for each sample separately in mixed-effects model?

Dear R-helpers,

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
  - ...

I have no idea how to do this in a formally acceptable way. Any help 
would be appreciated.

Kind regards,

-- 
Michal J. Figurski, PhD
HUP, Pathology & Laboratory Medicine
Biomarker Research Laboratory
3400 Spruce St. 7 Maloney
Philadelphia, PA 19104
tel. (215) 662-3413
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: Data.csv
URL:
<https://stat.ethz.ch/pipermail/r-help/attachments/20100429/d6b88cf9/attachment.pl>

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)