[R-sig-ME] Non-constancy of variances in mixed model.

[Ben Ward]

Hi, I've been doing an experiment, measuring the dead-zone-diameters of 
bacteria, when they've been grown with paper diffusion disks of 
antimicrobial. There are two groups, or treatments - one is bacteria 
that have been cultured in said antimicrobial for the past year, the 
other group is of the same species, but lab stock and has not gone had 
any prior contact with the antimicrobial. Because I take four readings 
from each disk, and there are three disks in a petri dish, and three 
petri dishes for each lineage (or replicate), and there are 5 lineages 
for each group/treatment, I want to use lme() so I can add lineage, 
dish, and disk as random effects, and have group as the fixed effect, so 
I'm analysing differences between two groups, a bit like a t.test or 
anova, but without riddiculous degrees of freedom. I've done qqPlots() 
and normality tests on my data, according to shapiro.test() on both my 
subsets (I split the data into the two groups):

 > shapiro.test(DiamEVDettol)

         Shapiro-Wilk normality test

     data:  DiamEVDettol
     W = 0.9823, p-value = 0.0221


 > shapiro.test(DiamNEDettol)

         Shapiro-Wilk normality test

     data:  DiamNEDettol
     W = 0.9788, p-value = 0.007677

 > qqPlot(DiamEVDettol, dist="norm")
 > qqPlot(DiamNEDettol, dist="norm")

Which from the R-Book shows non-normality. qqPlot() I did show a slight 
S shape, and then slightly more profound s-shape respectively 
suggesting, lepikurtosis. Skewness and Kurtosis are:

 > skewness(DiamEVDettol)
     [1] -0.1917142

 > kurtosis(DiamEVDettol)
     [1] -0.7594319

 > skewness(DiamNEDettol)
     [1] 0.1314403

 > kurtosis(DiamNEDettol)
     [1] -0.8365288

(I'm not sure how to get p-values to see if they're significant or not) .

But I'm unsure on how to proceed, I had done:

allrandom <- lme(Diameter~1,random=~1|Group/Lineage/Dish/Disk,data=Dataset)
So as I could do a Variance components analysis:

 > VarCorr(allrandom)
                 Variance     StdDev
     Group =     pdLogChol(1)
     (Intercept) 1.8773750    1.3701734
     Lineage =   pdLogChol(1)
     (Intercept) 0.2648475    0.5146333
     Dish =      pdLogChol(1)
     (Intercept) 0.0601047    0.2451626
     Disk =      pdLogChol(1)
     (Intercept) 0.1456451    0.3816348
     Residual    1.3456346    1.1600149

 > vars<-c(1.8773750,0.2648475,0.0601047,1.3456346)

 > 100*vars/sum(vars)
     [1] 52.914183  7.464779  1.694063 37.926975

Before moving on to:

groupfixed <-lme(Diameter~Group,random=~1|Lineage/Dish/Disk,data=Dataset)

Because Group is what I'm interested in, and is the only one with 
informative factor levels, other than maybe Lineage. I was going to gon 
on and construct similar models with different interactions or terms, 
but this issue of skewness and non normality stopped me, at first 
boxplots show what, in my only 3 years undergraduate experience looked 
like normal data - indeed neater data than I'm used to in biological 
science experiments. However this furthur probing o skewness and 
kurtosis and qqplots and shapiro.tests reveals otherwise. The variance 
is also not the same between the two groups (this was anticipated, 
because of natural selection of the bacteria in one of the groups, that 
variance could be less).

I'm concerned as to what I should do about the data I have, in terms of 
transformations, or whether there's a way for mixed models to take into 
account different distributions of data like a glm so my data does not 
have to be strictly normal or have equal variance.

Another concern of mine is whether I should be using lme as above, or as 
a book I read states, re-coding factor levels, and using lmer, for example:

Treatment<-factor(Treatment)
Liver<-factor(Liver)
Rat<-factor(Rat)
rat<-Treatment:Rat
liver<-Treatment:Rat:Liver
model<-lmer(Glycogen~Treatment+(1|rat)+(1|liver)

so with me it might be:

Group<-factor(Group)
Lineage<-factor(Lineage)
Dish<-factor(Dish)
Disk<-factor(Disk)
lineage<-Group:Lineage
dish<-Group:Lineage:Dish
disk<-Group:Lineage:Dish:Disk
model<-lmer(Diameter~Group+(1|lineage)+(1|dish)+(1|disk)

I'm not clear on what the difference would be here.


Thanks,

Ben W.