[R-sig-eco] adonis and contrasts

[Christine Griffiths]

Dear Jari,

Thank you very much for your help. I have tried both methods, but neither 
of which seem to have worked. This is most likely an oversight on my part 
so I have included an example of the alternative method provided with the 
hope that I am simply missing something.

dataset[1:5,]
  ID  Area Treatment Block BlockID Month X5 X6 X7 X8 X9 X11 X12 X13 X14 X15
1  1 lower 1.Control     a       1    1a 11  2  1  0 58  10  28 137   0   0 

2  2 lower 1.Control     b       2    1a  6 10  2  0 36   1  26 192  65   0 

3  3 lower 1.Control     c       3    1a  2  0  0  0 27  41  33 394  16   0 

4  4 upper 1.Control     d       4    1a 22 54  7  0 26  42  10 275 137   4
5  5 upper 1.Control     e       5    1a  1  5 10  0  0  19  45 599  60   0


X16 X17 X18 X19 X20 X22 X23 X24 X25 X26 X27 X28 X29 X30
87  18  68  10   0   0  43   0 0   0   0   0   0   5
6  35  94  14   0   4  28   2  0   0   0   0   6   8
3  14   4  16   0   2  61   1  0   0   0   0   1  11
22  10  30   1   0  35  19   0  1   3   2   0   0   0
6   3   0   5   0   4  39   0  0   0   0   0   0   6

dat<-dataset[1:16,7:30]
datenv<-dataset[1:16,1:6]

treatment<-as.factor(Treatment)
levels(treatment)
contrasts(treatment)<-cbind(c(-2,1,1),c(0,-1,1))
contrasts(treatment)
treatment<-model.matrix(~treatment)[,-1]

m1<-adonis(dat~treatment,data=datenv,strata=datenv[,"Block"],permutation=999,method="bray")

m1

Call:
adonis(formula = dat ~ treatment, data = datenv, permutations = 999, 
method = "bray", strata = datenv[, "Block"])

                 Df SumsOfSqs   MeanSqs   F.Model     R2 Pr(>F)
treatment  2.000000  0.112451  0.056225  0.680226 0.0947   0.48
Residuals 13.000000  1.074541  0.082657           0.9053
Total     15.000000  1.186991                     1.0000

The above result does not seem to account for the contrast. Am I missing 
something?

Cheers
Christine


--On 26 April 2010 16:41 +0300 Jari Oksanen <jari.oksanen at oulu.fi> wrote:

> On 26/04/10 13:42 PM, "Christine Griffiths"
> <Christine.Griffiths at bristol.ac.uk> wrote:
>
>> Dear All,
>>
>> I have done pairwise comparisons of species compositions between three
>> treatments (1 v 2, 1 v 3, and 2 v 3) which were replicated 6 times (i.e.
>> 6 blocks). To avoid multiple testing I have decided to use contrasts
>> instead (treatments 1 v 2 + 3, and 2 v 3). Oksanen (2009) mentions that
>> this can be done, but when I have tried to incorporate the contrasts,
>> either it doesn't run or the output suggests that the contrast has not
>> been addressed.
>>
>> Below is an example of the model I ran before wanting to add a contrast:
>>
>> m1<-adonis(dat~Treatment,data=datenv,strata=datenv[,"Block"],permutation
>> =1000, method="bray")
>>
>> How do I add the following contrasts into the above model?
>> treatment<-as.factor(Treatment)
>> contrasts(treatment)<-cbind(c(-2,1,1),c(0,-1,1))
>> contrasts(treatment)
>>
>> [,1] [,2]
>> 1.Control   -2    0
>> 2.Radiata    1   -1
>> 3.Aldabra    1    1
>>
>> If this is not possible, is it acceptable to sum the species presence in
>> the grouped treatments (2 + 3) and divide by 2, then combine this data to
>> that of treatment 1 to form a matrix? However, as I would not be able to
>> do more than one contrast, then I don't really overcome the problem of
>> multiple testing.
>>
> Christince,
>
> Everything is possible in R --- though not always easy.
>
> This really is a problem in adonis(). The crux of the problem is adonis()
> handles contrasts in non-standard way: it sets contrasts in function call
> with arguments "contr.ordered" and "contr.unordered". Hank Stevens (the
> adonis author) preferred sum contrasts which are not the standard in R,
> and he wanted to have an easy way to use those within adonis() without
> mucking around with global options. The price is that some of your
> specific options may be difficult to handle. You need a name of a
> contrast setting function as the argument. Writing up your contrast
> functions may not be quite easy. However, the following ad hoc function
> seems to work for this very specific case:
>
> contr.mine <- function(...) cbind(c(-2,1,1), c(0,-1,1))
> adonis(dat~Treatment,data=datenv,strata=datenv[,"Block"],permutation=1000,
> method="bray", contr.unordered = "contr.mine")
>
> In this case you must make a new contrast function for each case (and
> check that they really are meaningful and correct contrasts).
>
> Another alternative is to make explicit model matrices for these
> variables, since adonis() does not re-set contrasts if the argument is a
> matrix which already has a "contrast" attribute (from the help of
> ?contrasts: " The argument ?contrasts¹ is ignored if ?x¹ has a
> matrix?contrasts¹ attribute set.") So the following should work:
>
>## From your example
> treatment<-as.factor(Treatment)
> contrasts(treatment)<-cbind(c(-2,1,1),c(0,-1,1))
>## New line
> treatment <- model.matrix(~ treatment)[,-1]
>## Your code with treatment model.matrix without the (Intercept)
> m1<-adonis(dat~treatment,strata=datenv[,"Block"],permutation=1000,
> method="bray")
>
> BTW, I recommend you use 999 permutations instead of 1000.
>
> It may be that we should have a talk about adonis() contrasts with Hank
> Stevens.
>
> Cheers, Jari
>