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-mod...@r-project.org Ben Bolker future) ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
