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)


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@aphis.usda.gov

	[[alternative HTML version deleted]]

