How to perform CCA with block design in R

Question

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.

---

---

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).

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...

...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

Answer 1

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

Does this make sense now?