[R-sig-eco] Adonis and Random Effects
Erin Nuccio
enuccio at gmail.com
Mon Mar 11 00:09:11 CET 2013
Hi again,
OK, figuring out if it's possible to use Unifrac with nested.npmanova may be necessary....
I just realized my test comparing nested.npmanova and adonis on the same model had no strata for adonis. When I add the strata GrasslandPlot to adonis, my p values are equal to 1. So adonis with no strata gives me similar values to nested.npmanova for the following model: community_data ~ Grassland + GrasslandPlot.
So, (community_data ~ Grassland + GrasslandPlot) approximates the correct statistics, but since this ignores all strata, I'm not sure if it's justified. Thoughts?
Thanks,
Erin
On Mar 10, 2013, at 3:42 PM, Erin Nuccio wrote:
> Thanks Steve, that is helpful.
>
> However, I've run into a small problem with nested.npmanova. It appears that I cannot supply my own distance matrix, and need to supply the raw species data. I am using Unifrac distances, which is not an option for vegdist. Anyone know if there is a workaround here?
>
> I did compare nested.npmanova to adonis with bray distance using the same model (community_data ~ Grassland + GrasslandPlot), and it looks like the F values are similar for Grassland (F values: 3.6 vs. 3.3), and the same for GrasslandPlot. The R2 values seem to stay the same no matter what I do in adonis, and the p values are all ~ 0.001.
>
> So, in case there is no way to use Unifrac distances with nested.npmanova, my backup plan would be to perform two adonis functions, and use the second function to get the approximate F value for Grassland and correct F value for GrasslandPlot:
> adonis(community_data ~ Treatment*Grassland + GrasslandPlot, strata=GrasslandPlot)
> adonis(community_data ~ Grassland + GrasslandPlot, strata=GrasslandPlot)
>
> Does this seem reasonable? Of course, the best thing would be to use the Unifrac distances with nested.npmanova if it's possible.
>
> Thank you,
> Erin
>
>
>
>
>
> On Mar 10, 2013, at 8:17 AM, JOHN S BREWER wrote:
>
>> Erin,
>>
>> Please check the February 25 post I made called "Permanova with nested data." It explains how to test whole plot and split plot effects correctly in adonis. But to answer your question, even if you treat Grassland as a fixed-plot effect (which seems perfectly reasonable), Grassland is a whole-plot effect. Using the model formula given and strata, adonis uses the split-plot error term (i.e., the residual error term) to test all effects. That's wrong because Grassland needs to be tested with the whole-plot error term. In the post I referred to, I describe how you can do a separate test for the whole plot using the BiodiversityR package and the nested.npmanova function. In this case, you would only include Grassland and GrasslandPlot as terms in the model. It's just doing a two-way nested manova. The whole-plot effect of Grassland will be tested correctly using the GrasslandPlot term. GrasslandPlot will be tested with the residual error term, which will be wrong, but you can ignore that. I've tried it with my own data and it works. One cautionary note. See the posts by Jari Oksanen and others about the versions of BiodiversityR and R used.
>>
>> Hope this helps
>>
>> Steve
>> ________________________________________
>> From: Erin Nuccio [enuccio at gmail.com]
>> Sent: Saturday, March 09, 2013 9:09 PM
>> To: JOHN S BREWER
>> Cc: r-sig-ecology at r-project.org
>> Subject: Re: [R-sig-eco] Adonis and Random Effects
>>
>> Hi Steve and R list,
>>
>> I was hoping you could clarify something you mentioned in previous post.
>>
>> A quick recap... I have a split-plot design where I determined the microbial communities at 3 grasslands (see post script for design). I am trying quantify the how much of my community can be explained by Treatment or Grassland effect. After talking with a statistician, it seems like treating Grassland as a Fixed effect would be reasonable (because I have such a small number of grasslands).
>>
>> You mentioned that if I treat Grassland as a Fixed effect, and use the following formula, the Grassland effect would not be tested correctly:
>>
>> adonis(formula = community_distance_matrix ~ Treatment*Grassland + GrasslandPlot, strata = GrasslandPlot)
>>
>> Why is this? Is there any way to remedy this?
>>
>> Thanks for your feedback,
>> Erin
>>
>>
>> Experimental design:
>> 4 split plots * 2 Treatments * 3 Grasslands = 24 observations
>> Treatment: 2 levels (each within 1 split plot)
>> Grassland: 3 levels
>> GrasslandPlot: 12 levels (4 split plots nested in 3 Grasslands)
>>
>>
>>
>>
>>
>> On Feb 4, 2013, at 6:22 AM, Steve Brewer wrote:
>>
>>> Erin,
>>>
>>> There have been a lot of similar queries (e.g., repeated measures, nested
>>> permanova). Jari can correct me if I am wrong, but as far as I know, no
>>> one has developed a way to define multiple error terms in adonis.
>>>
>>>
>>> You can use adonis, however, to get the split-plot effects. If you want to
>>> make a grassland a random effect, use the following statement
>>>
>>> adonis(formula = community_distance_matrix ~ Treatment + Grassland +
>>> GrasslandPlot, strata = GrasslandPlot)
>>>
>>>
>>> The treatment effect will be correct because the residual error term
>>> (which is equivalent to treatment x GrasslandPlot interaction nested
>>> within Grassland) is the correct error term. The Grassland effect,
>>> however, will not be tested correctly because it is using the residual
>>> error term when it should be using GrasslandPLot as the error term. You
>>> can determine what the F stat for Grassland should be, however, using the
>>> Ms Grassland and MS GrasslandPlot from the anova table to construct the F
>>> test. You just won't get a p-value for the test.
>>>
>>> If you want to treat Grassland as a fixed effect, the model is similar but
>>> defines the interaction
>>>
>>> adonis(formula = community_distance_matrix ~ Treatment*Grassland +
>>> GrasslandPlot, strata = GrasslandPlot)
>>>
>>>
>>> In this case, the treatment x grassland interaction will be tested
>>> correctly, as will the treatment effect, but not the Grassland effect.
>>>
>>> Unfortunately, you cannot just take averages of abundances across the
>>> treatment and control in each plot and then do a separate analysis of
>>> Grassland and GrasslandPLot (unless you're using Euclidean distances). I
>>> suspect you're not using Euclidean distances.
>>>
>>> Hope this helps some.
>>>
>>> Good luck,
>>> Steve
>>>
>>>
>>> J. Stephen Brewer
>>> Professor
>>> Department of Biology
>>> PO Box 1848
>>> University of Mississippi
>>> University, Mississippi 38677-1848
>>> Brewer web page - http://home.olemiss.edu/~jbrewer/
>>> FAX - 662-915-5144
>>> Phone - 662-915-1077
>>>
>>>
>>>
>>>
>>> On 2/4/13 1:14 AM, "Erin Nuccio" <enuccio at gmail.com> wrote:
>>>
>>>> Hello List,
>>>>
>>>> Is adonis capable of modeling random effects? I'm analyzing the impact
>>>> of a treatment on the microbial community in a split-plot design (2
>>>> treatments per plot, 4 plots per grassland, 3 grasslands total). I would
>>>> like to quantify how much of the variance is due to the Treatment versus
>>>> the Grassland. It seems like Grassland should be a random effect, since
>>>> there are thousands of grasslands, and I'm only looking at 3.
>>>>
>>>> Thanks for your help,
>>>> Erin
>>>>
>>>>
>>>> Here are my factors:
>>>> 'data.frame': 24 obs. of 4 variables:
>>>> $ Treatment : Factor w/ 2 levels "T1","T2": 1 1 1 1 1 2 2 2 1 1 ...
>>>> $ Grassland : Factor w/ 3 levels "G1","G2","G3": 3 3 1 1 1 2 2 1 2 2
>>>> ...
>>>> $ Plot : Factor w/ 4 levels "P1","P2","P3","P4": 1 2 2 3 4 1 3 2
>>>> 1 2 ...
>>>> $ GrasslandPlot: Factor w/ 12 levels "G1:P1","G1:P2","G1:P3"..: 9 10 2 3
>>>> 4 5 7 2 5 6 ...
>>>>
>>>> And here's the error message:
>>>> Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
>>>> contrasts can be applied only to factors with 2 or more levels
>>>> In addition: Warning messages:
>>>> 1: In Ops.factor(1, Grassland) : | not meaningful for factors
>>>> 2: In Ops.factor(1, GrasslandPlot) : | not meaningful for factors
>>>>
>>>> _______________________________________________
>>>> R-sig-ecology mailing list
>>>> R-sig-ecology at r-project.org
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>>>
>>>
>>
>
More information about the R-sig-ecology
mailing list