I see, so I would perhaps use things like plot() and make sure that they 
are ok. In fact I've read a section of the R book, in which he treats 
the dataset as a groupedData object, refits the model, and then uses 
plot() with options to check for various things like heteroscedacity and 
other errors, would this be more valuable to me? I apologise for what 
may seem an elementary question, but this is my first shot at doing this 
kind of stats, it's not covered on my Biology undergraduate course, and 
our statistical help isn't strictly a statistician by profession. And a 
lot of the other lectureres are of a "stick a t-test in it, it'll be 
fine mentality", mostly because they did their work outside of modern 
computers, so ad to settle for less sophisticated analysis, at least 
that what we're told. I'm also going to have a look at the lmer() 
version I mentioned below and see if it comes to radically different 
results, and for completeness and to satisfy myself I will do a nested 
ANOVA, just so I can be happy that, as you say, if my study is balanced 
- which I've done my best to make it so, then they should tell me 
largely the same things.

I'm afraid I do not know what this mcmcsamp() function is, I've had a 
read of 
http://r.789695.n4.nabble.com/lmer-estimation-of-p-values-and-mcmcsamp-td882804.html 
and https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html and had a 
read of ?mcmcsamp, and I get the function is used on the model object, 
and is run as many times as defined by the user as 'n'. I've also read 
of "pvals.fnc" which has an option "withMCMC" and a similar function 
that returns an anova table as well as p-values. Is there a source of 
reading for p-values and such from mcmcsamp approach, that's accesible 
for non-mathematicians like myself?

Thanks,
Ben W.


On 14/03/2011 06:26, John Maindonald wrote:
> One does not expect the dependent variable to be normal, and certainly not
> independently and identically distributed.  Rather, the model is
> constructed using a kitset in which the components are
> 1) 'fixed' effects (in matrix terminology, the X beta part of the model)
> 2)  normal distributions (the random effects, making up the Z b part of the model)
> 3) an exponential family distribution,  in your case normal
> (the e part of the model)
>
> Adding up 1), 2) and 3) to give y = X beta + Z b + e is very unlikely to give
> anything that looks like a normal distribution.  So your test for normality is
> misdirected.  It probably is worth using normal probability plots for a rough
> check that the residuals are are consistent with a normal distribution, and
> that the random effects at the various levels (disks, etc) are pretty much
> consistent with a normal distribution.
>
> You have to extract the relevant quantities from the fitted model object.
>
> For reasons that have been worked over on this list, extensively, it would
> be pointless to use tests to check for normal distributions, even if suitable
> tests were available.  (I do not believe that in general there are good tests.)
>
> If your design is indeed hierarchical and balanced, I think you might find
> it simpler to stick with anova, with an appropriate Error argument.  You
> have as much (or as little!) licence to ignore the degrees of freedom in
> anova as you have not to worry about degrees of freedom for t-statistics
> from your lmer calculations.  Degrees of freedom are part of the deal if
> you want p-values without using the mcmcsamp() approach, or maybe
> profile likelihood.
>
> John Maindonald             email: john.maindonald@anu.edu.au
> phone : +61 2 (6125)3473    fax  : +61 2(6125)5549
> Centre for Mathematics&  Its Applications, Room 1194,
> John Dedman Mathematical Sciences Building (Building 27)
> Australian National University, Canberra ACT 0200.
> http://www.maths.anu.edu.au/~johnm
>
> On 14/03/2011, at 2:24 PM, 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.
>>
>> _______________________________________________
>> R-sig-mixed-models@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>


	[[alternative HTML version deleted]]


