4
$\begingroup$

Following example and data are completely fabricated:

Suppose I am studying grooming behaviour in apes. I have four cages, 8 apes in each (4 females + males). For 24 hours I did an observations with aim to record the number of events when females were grooming males.


enter image description here


I would like to know, whether I can explain the variability of grooming events [grooming] with these 5 characteristics of males [m_char] (x1, x2,...x5).

enter image description here

My approach was to perform CCA (a.k.a. canonical correspondence analysis) as following:

library(vegan)
my.cca <- cca(grooming ~ x1 + x2 + x3 + x4 + x5,
              data = m_char,
              scale = TRUE)

and test this model by anova to obtain p-values for factors

anova(my.cca, by="terms", permutations=1000)

# Permutation test for cca under reduced model
# Terms added sequentially (first to last)
# Permutation: free
# Number of permutations: 1000
# 
# Model: cca(formula = grooming ~ x1 + x2 + x3 + x4 + x5,
#            data = m_char, scale = TRUE)
#          Df ChiSquare      F   Pr(>F)   
# x1        1    0.5327 1.1054 0.718282   
# x2        1    0.6216 1.2899 0.359640   
# x3        1    0.6341 1.3159 0.314685   
# x4        1    0.6739 1.3984 0.191808   
# x5        1    0.8782 1.8225 0.008991 **
# Residual 10    4.8188                   
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Although, x5 seems to have significant effect on grooming, when I look on plot...

plot(my.cca)

...apes from one cage seems to strongly group together.

I think that the structure of my data causes a problem to cca analysis, because, when you look on table with numbers of interactions...

enter image description here

...only zeros in yellow squares are "true" zeros (i.e. there is no interactions between apes), the rest of zeros outside yellow squares are just a consequence of being in different cage.

Is it possible to somehow tell the CCA to ignore "false zeros"? How to incorporate such data structure appropriately? My main aim is to find whether one or multiple characteristics of males can be used to explain the grooming rates.

Thank you very much for any suggestion!


"grooming" example data:

AF1 <- c(3,7,2,0,rep(0, times = 12))
AF2 <- c(8,0,0,3,rep(0, times = 12))
AF3 <- c(0,2,0,0,rep(0, times = 12))
AF4 <- c(4,0,1,0,rep(0, times = 12))
AF4 <- c(4,0,1,0,rep(0, times = 12))
BF1 <- c(rep(0, times = 4),4,1,5,0,rep(0, times = 8))
BF2 <- c(rep(0, times = 4),3,0,0,1,rep(0, times = 8))
BF3 <- c(rep(0, times = 4),0,2,0,0,rep(0, times = 8))
BF4 <- c(rep(0, times = 4),0,0,7,0,rep(0, times = 8))
CF1 <- c(rep(0, times = 8),4,0,4,0,rep(0, times = 4))
CF2 <- c(rep(0, times = 8),2,0,0,0,rep(0, times = 4))
CF3 <- c(rep(0, times = 8),0,3,0,4,rep(0, times = 4))
CF4 <- c(rep(0, times = 8),0,0,9,0,rep(0, times = 4))
DF1 <- c(rep(0, times = 12),0,2,0,0)
DF2 <- c(rep(0, times = 12),6,0,0,1)
DF3 <- c(rep(0, times = 12),0,0,1,0)
DF4 <- c(rep(0, times = 12),4,0,1,0)

male_id <- c("AM1", "AM2", "AM3", "AM4",
             "BM1", "BM2", "BM3", "BM4",
             "CM1", "CM2", "CM3", "CM4",
             "DM1", "DM2", "DM3", "DM4")

grooming <- data.frame(AF1,AF2,AF3,AF4,
                       BF1,BF2,BF3,BF4,
                       CF1,CF2,CF3,CF4,
                       DF1,DF2,DF3,DF4)
rownames(grooming) <- male_id

"m_char" example data:

x1 <- c(2,5,3,1,8,1,1,6,2,3,5,1,1,6,6,7)
x2 <- c(4,9,1,3,2,1,9,4,3,1,9,2,9,1,4,3)
x3 <- c(1,2,5,4,4,4,5,2,2,1,1,1,2,5,5,9)
x4 <- c(4,1,1,8,8,6,6,6,6,8,1,6,2,2,1,1)
x5 <- c(1,7,3,3,4,1,5,1,3,8,4,8,8,5,9,9)

m_char <- data.frame(x1, x2, x3, x4, x5)
rownames(m_char) <- male_id
$\endgroup$
2
$\begingroup$

I believe the issue is that the permutation test you are using is too liberal; it assumes a Null hypothesis in which all observations are exchangeable. From what you say, observations within a cage are exchangeable, but are not exchangeable between cages.

To use the more restrictive null in the permutation test we can use a restricted permutation design, in which we ask vegan to only permute observations within cages and never between cages. This is most easily done using a blocking factor.

Create the blocking factor:

m_char <- transform(m_char,
                    cage = factor(substring(rownames(m_char), 1, 1)))

Next, define the restricted permutation design (the defaults for how() mean we get randomisation within the levels of the factor passed to blocks)

ctrl <- how(nperm = 1000, blocks = m_char$cage)

Now fit the model and remove the between-cage variation, which is also required for these analyses with blocks

my.cca <- cca(grooming ~ x1 + x2 + x3 + x4 + x5 + Condition(cage),
              data = m_char)

Now do the restricted permutation test

set.seed(10)
anova(my.cca, by = "terms", permutations = ctrl)

This produces:

> anova(my.cca, by = "terms", permutations = ctrl)
Permutation test for cca under reduced model
Terms added sequentially (first to last)
Blocks:  m_char$cage 
Permutation: free
Number of permutations: 1000

Model: cca(formula = grooming ~ x1 + x2 + x3 + x4 + x5 + Condition(cage), data = m_char)
         Df ChiSquare      F  Pr(>F)  
x1        1   0.52470 1.4119 0.15185  
x2        1   0.60901 1.6388 0.03397 *
x3        1   0.21703 0.5840 0.96503  
x4        1   0.36739 0.9886 0.36563  
x5        1   0.83984 2.2600 0.03497 *
Residual  7   2.60129                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

And the triplot shows different behaviour to that you obtained

enter image description here

Does this make sense now?

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.