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

[Ben Ward]

I would like to add an apology for cross posting on two lists, it was 
only recently pointed out to me that this is frowned upon. Although 
there was some confusion as to which I should direct my question as I 
suppose part of it is not specifically mixed models. In addition I 
should add the packages I used are the ones included with Rcmdr as 
dependencies. In particular, the qqPlot function, as opposed to qqplot, 
is from the car package, and is the command used if one uses the Rcmdr 
gui to construct "Quantile-comparison plot". It's largely the same as 
qqplot, but with red lines and a little bit more information. The 
package e1071 was used for the skewness and kurtosis functions, and I'm 
using the nlme package for the mixed effect model stuff.

Thanks,

Ben W.

On 14/03/2011 03:24, Ben Ward wrote:
> 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.
>