[R-sig-ME] mixed-effects model specification question
Kingsford Jones
kingsfordjones at gmail.com
Wed May 7 06:26:10 CEST 2008
Mark,
If you have just 3 observations per rat (and 2 missing values, judging
by the reported sample size), that implies 1 observation within a
tissue sample within a rat. So, unless I'm missing something, the
residual variance and tissue-in-rat variance are confounded and there
are not sufficient degrees of freedom to estimate both. What does the
output of "intervals(myModel)" report? I'm guessing you'll get an
error about a non-positive-definite covariance matrix or a CI for a
variance component that is essentially +/-Inf. Even if the experiment
were conducted perfectly and you could ignore effects such as those
brought up by Antonio, you are trying to fit a model that is too
complex for the available data.
Kingsford Jones
On Tue, May 6, 2008 at 5:56 PM, Mark Kimpel <mwkimpel at gmail.com> wrote:
> The "treatments" are actually two selectively bred lines of rats, with 6
> animals from each line. We need p-values as a filter to determine which of
> ~31,000 genes on the array are significant. The p-values are subjected to
> FDR correction, of course. I've never seen a micro-array paper published
> without p-values.
>
> As for Andrew's comments:
> 1. yes, there are three measurments taken per rat, but they are from
> different brain regions. Since each brain region was processed separately,
> however, the region and batch effects cannot be separated. Having said that,
> I still label this factor as "region" in my model.
> 2. In the past I have done exactly as you suggested and averaged brain
> regions, but I was hoping nlme would allow for more sophistication. Because
> certain genes vary greatly across brain regions, averaging would seem to
> discard some useful information which could be used to improve the inference
> made on treatment effect. If I used averaging, I would just have 12
> measurements made on 12 rats and could do a simple t-test, but I don't think
> that would be preferable, would it?
>
> Here is the output of my model:
>
>
> myModel <- lme(fixed = geneExpression ~ Treatment, random = ~1 |
> Rat/Tissue)
> > myModel
> Linear mixed-effects model fit by REML
> Data: NULL
> Log-restricted-likelihood: 40.60593
> Fixed: geneExpression ~ Treatment
> (Intercept) TreatmentLAD
> 12.628465893 -0.009280945
>
> Random effects:
> Formula: ~1 | Rat
> (Intercept)
> StdDev: 2.752519e-06
>
> Formula: ~1 | Tissue %in% Rat
> (Intercept) Residual
> StdDev: 0.06226097 0.0002623086
>
> Number of Observations: 34
> Number of Groups:
> Rat Tissue %in% Rat
> 12 34
> > summary(myModel)
> Linear mixed-effects model fit by REML
> Data: NULL
> AIC BIC logLik
> -71.21185 -63.88317 40.60593
>
> Random effects:
> Formula: ~1 | Rat
> (Intercept)
> StdDev: 2.752519e-06
>
> Formula: ~1 | Tissue %in% Rat
> (Intercept) Residual
> StdDev: 0.06226097 0.0002623086
>
> Fixed effects: geneExpression ~ Treatment
> Value Std.Error DF t-value p-value
> (Intercept) 12.628466 0.01510064 22 836.2870 0.0000
> TreatmentLAD -0.009281 0.02135553 10 -0.4346 0.6731
> Correlation:
> (Intr)
> TreatmentLAD -0.707
>
> Standardized Within-Group Residuals:
> Min Q1 Med Q3 Max
> -0.0083222376 -0.0027523180 -0.0001242287 0.0034693142 0.0087149017
>
> Number of Observations: 34
> Number of Groups:
> Rat Tissue %in% Rat
> 12 34
>
>
> >
>
> On Tue, May 6, 2008 at 4:42 PM, <Antonio_Paredes at aphis.usda.gov> wrote:
>
> >
> > How were animals allocated to treatments (at random)? What about the
> > housing of the animals. Both of these factors could be influential in
> > selecting the unit of study.
> >
> > You also need to get an ideal at to why a p-value is enough to support the
> > objective of the study. Why not estimation of treatment effects?
> >
> > Tony.
> >
> >
> > *Andrew Robinson <A.Robinson at ms.unimelb.edu.au>*
> > Sent by: r-sig-mixed-models-bounces at r-project.org
> >
> > 05/06/2008 03:19 PM
> > To
> > Mark Kimpel <mwkimpel at gmail.com> cc
> > r-sig-mixed-models at r-project.org
> > Subject
> > Re: [R-sig-ME] mixed-effects model specification question
> >
> >
> >
> >
> > Mark,
> >
> > You should talk to a local statistician about this, but I think that
> > you can probably average across the measurements within each rat, if
> > all you are interested in is a treatment effect. For the analysis of
> > treatment the relevant unit of replication should be the rat, in any
> > case.
> >
> > (Does anyone have any thoughts on why that might be a bad idea?)
> >
> > Also if I understand your design, there are three batches per rat. I
> > suspect that using Rat/Tissue would lead to an over-parametrized
> > model, if my interpretation is correct. My guess is that Rat should
> > be adequate.
> >
> > Final points
> >
> > 1) My speculations would be better-informed if you showed us the
> > output of the model that you proposed - fyi :) .
> >
> > 2) You could try all three configurations and see if it makes any
> > difference to the inference of interest. I suspect that it will
> > not.
> >
> > I hope that this helps,
> >
> > Andrew
> >
> > On Tue, May 06, 2008 at 10:23:08AM -0400, Mark Kimpel wrote:
> > > Perhaps a bit of a newbie question, but I need to get this right. I need
> > to
> > > make sure I am specifying a model correctly. Here is our
> > less-than-perfect
> > > experimental design:
> > >
> > > 36 rats divided into two treatment groups, i.e. 18 per group
> > >
> > > each rat has measurements taken on 3 brain regions, but each brain
> > region
> > > was analyzed in a separate batch (there are strong batch effects) so we
> > > really can't compare the regions per se, but do recognize that the 3
> > regions
> > > from a single rat have variance above and beyond that accounted for by
> > > technical factors. Since, because of the batch effect we are not going
> > to
> > > look at the effect of brain region, I "think" this should be a
> > considered a
> > > random effect.
> > >
> > > So I have rat, treatment, and region(batch) as variables. The only thing
> > I
> > > am interested in getting a p-value for is treatment.
> > >
> > > Is the model below correct and can I squeek by with using nlme to get a
> > > p-value (notwithstanding recent threads on this list)?
> > >
> > > myModel <- lme(fixed = geneExpression ~ Treatment, random = ~1 |
> > > Rat/Tissue)
> > >
> > > Thanks,
> > > Mark
> > > --
> > > Mark W. Kimpel MD ** Neuroinformatics ** Dept. of Psychiatry
> > > Indiana University School of Medicine
> > >
> > > 15032 Hunter Court, Westfield, IN 46074
> > >
> > > (317) 490-5129 Work, & Mobile & VoiceMail
> > > (317) 663-0513 Home (no voice mail please)
> > >
> > > ******************************************************************
> > >
> > > [[alternative HTML version deleted]]
> > >
> > > _______________________________________________
> > > R-sig-mixed-models at r-project.org mailing list
> > > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
> > --
> > Andrew Robinson
> > Department of Mathematics and Statistics Tel: +61-3-8344-6410
> > University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599
> > http://www.ms.unimelb.edu.au/~andrewpr<http://www.ms.unimelb.edu.au/%7Eandrewpr>
>
> > http://blogs.mbs.edu/fishing-in-the-bay/
> >
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
> >
>
>
>
>
> --
> Mark W. Kimpel MD ** Neuroinformatics ** Dept. of Psychiatry
> Indiana University School of Medicine
>
> 15032 Hunter Court, Westfield, IN 46074
>
> (317) 490-5129 Work, & Mobile & VoiceMail
> (317) 663-0513 Home (no voice mail please)
>
> ********************************************************** ********
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list