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

From: Ben Ward <benjamin.ward_at_bathspa.org>
Date: Mon, 14 Mar 2011 03:24:55 +0000


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.



R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Mon 14 Mar 2011 - 03:31:37 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Mon 14 Mar 2011 - 03:40:21 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive