Non-constancy of variances in mixed model.

Ben Ward

Sun, Mar 13, 2011 10:43 PM

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.