[R-sig-ME] Specifying the random effects term for class vs continuous variables

Douglas Bates bates at stat.wisc.edu
Wed Apr 6 20:05:58 CEST 2011 

On Tue, Mar 22, 2011 at 6:32 PM,  <Susan.A.Shriner at aphis.usda.gov> wrote:
> Hello list,
>
> I am returning to R and statistical modeling after a several year hiatus
> so I'm a bit rusty.  I am trying to determine the best way to specify my
> random effect.
>
> Briefly, my data are from an experimental infection study of rats with
> influenza viruses. We inoculated rats with one of five different subtypes
> (strains) of the virus, and then each day for seven days we sampled 3 rats
> per subtype by  collecting tissues from nasal turbinates, trachea, cranial
> lung, and caudal lung. We then tested the tissues via PCR to determine the
> concentration of virus in the tissue (our response variable).
>
> My set up is balanced except for having unequal numbers of males and
> females (we had 48 males and 57 females).
>
> After doing some reading and perusing the list serve I came up with three
> specifications to compare. But subsequently I came across a post that
> discriminated how you would specify a continuous variable vs. a class
> variable (I have included the relevant part of the post at the end of my
> message).  Out of the set ups that I compared, the specification that came
> up as best by AIC was one that the poster indicated was appropriate for
> continuous data.
>
> **So, my question is, is the specification invalid for a class variable or
> do I need to discard it? And, does anyone have a better suggestion?**
>
> As background on the effect of individual rats on virus replication, there
> is individual heterogeneity in rat immune responses so some individuals
> have high viral replication rates and some much lower. But, if one tissue
> type has comparatively high replication rates, then the corresponding
> tissues associated with that animal are also likely to show relatively
> high replication rates (but there are large differences between tissues).
>
>
> Here are the original models I evaluated:
> lmer.Rat1 <- lmer(logQuantity ~ Subtype + Tissue + Sex + DPI + Subtype*DPI
> + Sex*DPI + (1|Rat), data = Rat)
> lmer.Rat2 <- lmer(logQuantity ~ Subtype + Tissue + Sex + DPI + Subtype*DPI
> + Sex*DPI + (1|Rat) + (1|Tissue), data = Rat)
> lmer.Rat3 <- lmer(logQuantity ~ Subtype + Tissue + Sex + DPI + Subtype*DPI
> + Sex*DPI + (Tissue|Rat), data = Rat)
> anova(lmer.Rat1, lmer.Rat2, lmer.Rat3)

I think you want the random effects in the second model to be (1|Rat)
+ (1|Tissue:Rat) .  An equivalent specification is (1|Rat/Tissue) but
I prefer the former.

This is the model that is more directly comparable to model 3 with
(Tissue|Rat) for the random effects.  Typically I would write the
random effects for that model as (0+Tissue|Rat) because then the
variance components are easier to interpret.  The two formulations
produce equivalent models.

>
> with the following result:
> Data: Rat
> Models:
> lmer.Rat1: logQuantity ~ Subtype + Tissue + Sex + DPI + Subtype * DPI +
> lmer.Rat1:     Sex * DPI + (1 | Rat)
> lmer.Rat2: logQuantity ~ Subtype + Tissue + Sex + DPI + Subtype * DPI +
> lmer.Rat2:     Sex * DPI + (1 | Rat) + (1 | Tissue)
> lmer.Rat3: logQuantity ~ Subtype + Tissue + Sex + DPI + Subtype * DPI +
> lmer.Rat3:     Sex * DPI + (Tissue | Rat)
>          Df     AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)
> lmer.Rat1 17 1085.25 1151.3 -525.63
> lmer.Rat2 18 1092.61 1162.6 -528.30   0.00      1          1
> lmer.Rat3 26  978.04 1079.1 -463.02 130.57      8     <2e-16 ***
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
>
> The ANOVA shows model 3 to be the best by AIC; here's its summary:
> Linear mixed model fit by REML
> Formula: logQuantity ~ Subtype + Tissue + Sex + DPI + Subtype * DPI + Sex
> * DPI + (Tissue | Rat)
>   Data: Rat
>  AIC  BIC logLik deviance REMLdev
>  1025 1126 -486.3      926   972.5
> Random effects:
>  Groups   Name                   Variance Std.Dev. Corr
>  Rat      (Intercept)            0.69534  0.83387
>          TissueCranial lung     0.36080  0.60067   0.365
>          TissueNasal turbinates 2.55347  1.59796  -0.589  0.081
>          TissueTrachea          0.68626  0.82841  -0.909 -0.051  0.419
>  Residual                        0.23086  0.48048
> Number of obs: 360, groups: Rat, 90
>
> Fixed effects:
>                       Estimate Std. Error t value
> (Intercept)             1.52152    0.28794   5.284
> SubtypeH3N8             1.12680    0.32798   3.436
> SubtypeH4N6            -0.63473    0.33786  -1.879
> SubtypeH4N8            -1.18605    0.35007  -3.388
> SubtypeH6N2            -1.29543    0.33388  -3.880
> TissueCranial lung      0.39796    0.09560   4.163
> TissueNasal turbinates  1.12349    0.18304   6.138
> TissueTrachea          -0.16199    0.11294  -1.434
> SexM                   -0.56213    0.23896  -2.352
> DPI                    -0.25194    0.06559  -3.841
> SubtypeH3N8:DPI        -0.18220    0.08426  -2.162
> SubtypeH4N6:DPI         0.17228    0.08713   1.977
> SubtypeH4N8:DPI         0.19240    0.08801   2.186
> SubtypeH6N2:DPI         0.24590    0.08673   2.835
> SexM:DPI                0.09556    0.06633   1.441
>
> Correlation of Fixed Effects:
>            (Intr) SbH3N8 SbH4N6 SbH4N8 SbH6N2 TssCrl TssNst TssTrc SexM
> DPI    SH3N8: SH4N6: SH4N8: SH6N2:
> SubtypeH3N8 -0.573
> SubtypeH4N6 -0.639  0.490
> SubtypeH4N8 -0.708  0.470  0.507
> SubtypeH6N2 -0.624  0.495  0.522  0.498
> TissCrnllng -0.019  0.000  0.000  0.000  0.000
> TssNsltrbnt -0.214  0.000  0.000  0.000  0.000  0.196
> TissueTrach -0.293  0.000  0.000  0.000  0.000  0.212  0.422
> SexM        -0.471  0.000  0.103  0.342  0.072  0.000  0.000  0.000
>
> DPI         -0.862  0.578  0.614  0.687  0.605  0.000  0.000  0.000  0.417
>
> SbtH3N8:DPI  0.504 -0.897 -0.427 -0.417 -0.434  0.000  0.000  0.000  0.000
> -0.639
> SbtH4N6:DPI  0.537 -0.440 -0.898 -0.425 -0.471  0.000  0.000  0.000  0.005
> -0.631  0.473
> SbtH4N8:DPI  0.628 -0.430 -0.449 -0.905 -0.445  0.000  0.000  0.000 -0.291
> -0.736  0.478  0.463
> SbtH6N2:DPI  0.530 -0.442 -0.474 -0.421 -0.897  0.000  0.000  0.000  0.019
> -0.627  0.476  0.531  0.461
> SexM:DPI     0.347  0.010  0.003 -0.273  0.011  0.000  0.000  0.000 -0.899
> -0.352  0.017 -0.118  0.257 -0.122
>
>
>
> Here is an excerpt from the post: ([R-sig-ME] Simple lme/lmer random
> effects questions, Aug 11, 2008); Note: temp is a fixed effect and sleeve
> is the random effect
>
> ### Assuming temp is continuous...
>>
> m1 <- lmer(response ~ sex*temp + (1|sleeve.in.temp) + (temp-1|
> sleeve.in.temp), data) # enforcing independence of slope and intercept
> m2 <- lmer(response ~ sex*temp + (temp|sleeve.in.temp), data) #
> enforcing dependence of slope and intercept
> m3 <- lmer(response ~ sex*temp + (1|sleeve.in.temp), data) # No
> difference in slopes; merely difference in intercept
> anova(m1,m2,m3)
> ###
> ### Assuming temp is categorical...
> m1 <- lmer(response ~ sex*temp + (1|sleeve.in.temp) + (1|temp), data)
> m2 <- lmer(response ~ sex*temp + (1|sleeve.in.temp), data) #
> enforcing dependence of slope and intercept
> m3 <- lmer(response ~ sex*temp + (1|temp), data) # No difference in
> slopes; merely difference in intercept
> anova(m1,m2,m3)
>
> (Response from Dr. Hank Stevens, Associate Professor
> 338 Pearson Hall
> Botany Department
> Miami University
> Oxford, OH 45056)
>
>
>
> Susan Shriner
> susan.a.shriner at aphis.usda.gov
>
>        [[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