[BioC] technical replicates (again!): a summary
Naomi Altman
naomi at stat.psu.edu
Wed Mar 31 23:12:40 CEST 2004
Dear Ramon,
I did not look through all of your R-code, but here is what statistical
theory says:
If you treat the biological subject as a random effect, then the mixed
model ANOVA tests of all the treatment and other effects are identical to
what you get if you average the technical reps (within treatment). (This
assumes a balanced dye-swap design, with no missing observations.)
If you run lme on the experiment, and set everything up correctly, this is
what you will get and the df will be the same.
If you use SAS PROC MIXED with a Type III analysis, you will also get this.
If you use SAS PROC MIXED with the REML (default) option, you may get
something different because of the way the variance of the random effects
are estimated. When the variances match those obtained from the Type III
analysis, then you will get the same ANOVA as the other methods. However,
when the variance are estimated to be 0, then the biological subject
effect is (in effect) dropped from the model, so that the technical and
biological replicates are treated identically. (This does not happen in
lme, which also uses REML, because lme estimates the log(variance) which
cannot ever be 0.)
After struggling with this (and huge amounts of imbalance due to bad spots,
unequal spot replication, etc) I am pretty much settling for doing the
averaging (over all technical reps) before sending my data to the ANOVA
routine. (I should do a weighted analysis is account for this, but so far
I have not.) To do this, I have to do a lot of slow
preprocessing. However, because I am working in a non-medical environment
in which the number of replicates is generally small, the number of
treatments is generally large, and investigator patience is reasonably
good, this seems to be the best solution for our group. The biggest
draw-back is that to date I have not fully automated the procedure, so I
cannot off-load the data-processing to a research assistant (or publish
it to this list).
Regarding randomized blocks:
We have a RCB design if each biological sample gets all of the treatments
(e.g. cancer cells and normal cells from the same subject; different
organs from the same mouse). In this case, the "error" for testing
treatment effects is biological sample * treatment, and the technical
replicates appear as "pure error". The d.f. should be correct. If I have
time, I will look carefully at your design and model over the weekend to
see if I can determine why you are not getting the right d.f.
I hope that this is helpful, and apologize for not giving this all of the
attention needed.
--Naomi
At 10:49 AM 3/31/2004, Ramon Diaz-Uriarte wrote:
>Dear Gordon, Naomi, and BioC list,
>
>The issue of how to deal with technical replicates (such as those obtained
>when we do dye-swaps of the same biological samples in cDNA arrays) has come
>up in the BioC list several times. What follows is a short summary, with
>links to mails in BioC plus some questions/comments.
>
>
>There seem to be three major ways of approaching the issue:
>
>
>1. Treat the technical replicates as ordinary replicates
>*************************************************************
>E.g., Gordon Smyth in sept. 2003
>(https://www.stat.math.ethz.ch/pipermail/bioconductor/2003-September/002405.html)
>
>However, this makes me (and Naomi Altman ---e.g.,
>https://www.stat.math.ethz.ch/pipermail/bioconductor/2003-December/003340.html)
>
>uneasy (tech. reps. are not independent biological reps. which leads to the
>usual inflation of dfs and deflation of se).
>
>I guess part of the key to Gordon's suggestion is his comment that even if
>the
>s.e. are slightly underestimated, the ordering is close to being the optimal
>one. But I don't see why the ordering out to be much worse if we use the
>means of technical replicates as in 3. below. (Haven't done the math, but it
>seems that, specially in the pressence of strong tech. reps. covariance and
>small number of independent samples we ought to be better of using the means
>of the tech. reps).
>
>
>2. Mixed effects models with subject as random effect (e.g., via lme).
>******************************************************************************
>
>In late August of 2003 I asked about these issues, and Gordon seemed to agree
>that trying the lme approach could be a way to go.
>(https://www.stat.math.ethz.ch/pipermail/bioconductor/2003-August/002224.html).
>
>However, in September, I posted an aswer and included code that, at least for
>some cases, shows potential problems with using lme when the number of
>technical replicates is small.
>(https://www.stat.math.ethz.ch/pipermail/bioconductor/2003-September/002430.html)
>
>Nevertheless, Naomi Altman reports using lme/mixed models in a couple of
>emails
>(https://www.stat.math.ethz.ch/pipermail/bioconductor/2003-December/003191.html;
>
>(https://www.stat.math.ethz.ch/pipermail/bioconductor/2004-January/003481.html).
>
>After reading about randomizedBlock (package statmod) in a message in BioC (I
>think from Gordon), I have tried aggain the mixed models approach, since with
>tech. reps and no other random effects, we should be able to use
>randomizedBlock. Details in 5. below:
>
>
>3. Take the average of the technical replicates
>****************************************************
>Except for being possibly conservative (and not estimating tech. reps.
>variance component), I think this is a "safe" procedure.
>This is what I have ended up doing routinely after my disappointing tries
>with
>lme
>(https://www.stat.math.ethz.ch/pipermail/bioconductor/2003-September/002430.html)
>and what Bill Kenworthy seemed to end up doing
>(https://www.stat.math.ethz.ch/pipermail/bioconductor/2004-January/003493.html).
>
>I think this is also what is done at least some times in literature (e.g.,
>Huber et al., 2002, Bioinformatics, 18, S96--S104 [the vsn paper]).
>
>*********
>
>4. Dealing with replicates in future versions of limma
>***********************************************************
>
>Now, in Sept. 2004 Gordon mentioned that an explicit treatment of tech. reps.
>would be in a future version of limma
>(
>https://www.stat.math.ethz.ch/pipermail/bioconductor/2003-September/002411.html)
>and I wonder if Gordon meant via mixed-effects models, or some other way, or
>if there has been some progress in this area.
>
>
>
>5. Using randomizedBlock
>*****************************
>In a simple set up of control and treatment with dye-swaps, I have done some
>numerical comparisons of the outcome of a t-test on the mean of the technical
>replicates with lme and with randomizedBlock. [The function is attached]. The
>outcome of the t-test of the means of replicates and randomizedBlock, in
>terms of the t-statistic, tends to be the same (if we "positivize" the dye
>swaps). The only difference, then, lies in the df we then use to put a
>p-value on the statistic. But I don't see how we can use the dfs from
>randomizedBlock: they seem way too large. Where am I getting lost?
>
>
>Best,
>
>
>R.
>
>
>
>--
>Ramón Díaz-Uriarte
>Bioinformatics Unit
>Centro Nacional de Investigaciones Oncológicas (CNIO)
>(Spanish National Cancer Center)
>Melchor Fernández Almagro, 3
>28029 Madrid (Spain)
>Fax: +-34-91-224-6972
>Phone: +-34-91-224-6900
>
>http://bioinfo.cnio.es/~rdiaz
>PGP KeyID: 0xE89B3462
>(http://bioinfo.cnio.es/~rdiaz/0xE89B3462.asc)
>
>
>
>_______________________________________________
>Bioconductor mailing list
>Bioconductor at stat.math.ethz.ch
>https://www.stat.math.ethz.ch/mailman/listinfo/bioconductor
Naomi S. Altman 814-865-3791 (voice)
Associate Professor
Bioinformatics Consulting Center
Dept. of Statistics 814-863-7114 (fax)
Penn State University 814-865-1348 (Statistics)
University Park, PA 16802-2111
More information about the Bioconductor
mailing list