How to prove cooperation from behavioural sequences

Question

Situation: Two birds (male and female) protect their eggs in nest against an intruder. Each bird can use either attack or threat for protection, and be either present or absent. There is a pattern emerging from data that behaviour may be complementary - male attacks while female use threat display and vice versa.

My question is: How to statistically prove such cooperation? Or can anybody know some behavioural study which deals with similar analysis? Vast majority of sequential analyses I found are focused on DNA.

Here I provide some dummy data, but my original dataset is composed from dozens of pairs which were recorded exactly 10 minutes while defending their nest. Behavioural sequence of every bird is therefore 600 states long (each second has state). These shorter data should contain pattern similar to whole dataset.

male_seq <- rep(c("absent","present","attack","threat","present","attack",
                  "threat","present","attack","absent"),
                  times = c(3,4,8,2,6,3,2,6,2,1))

female_seq <- rep(c("absent","present","threat","present","threat","present",
                    "threat","attack","present","threat","attack","present",
                    "attack","threat","absent"),
                  times = c(2,6,2,1,2,1,1,3,5,3,1,3,3,2,2))

Answer 1

I post second answer since your last comment

By cooperation I mean "when male is attacking, the female make threats", and I would like to test this hypothesis against an alternative: "when male is attacking, the female do not prefer make threats" (in other words, behaviour of female is independent of male behaviour).

is a game-changer. It seems that the problem can be approach from totally different perspective. First, you are interested only in part of your sample when males are attacking. Second, you are interested if in such cases females make treats more often than we would expect if they made them randomly. To test such hypothesis we can use a permutation test: randomly shuffle either male_seq or female_seq (it doesn't matter) and then count cases where male_seq == "attack" and female_seq == "treat" to obtain null distribution. Next, compare count obtained from your data to counts in the null distribution to obtain $p$-value.

prmfun <- function() {
  sum(female_seq[sample(male_seq) == "attack"] == "threat")
}

mean(replicate(1e5, prmfun()) >= sum(female_seq[male_seq == "attack"] == "threat"))
## [1] 5e-05

You can define your test statistic differently, based on how do you define females' "preference". Permutation test in this case is a direct interpretation of your $H_0$: "behaviour of female is independent of male behaviour", that leads to: "female behaviour is random given male behaviour", so the behaviours are be randomly shuffled under $H_0$.

Moreover, even if you assumed that the behaviours appear in clusters of the same behaviour repeated for some period of time, with permutation test you can shuffle whole clusters:

female_rle <- rle(female_seq)
n_rle <- length(female_rle$values)

prmfun2 <- function() {
  ord <- sample(n_rle)
  sim_female_seq <- rep(female_rle$values[ord], female_rle$lengths[ord])
  sum(sim_female_seq[male_seq == "attack"] == "threat")
}

mean(replicate(1e5, prmfun2()) >= sum(female_seq[male_seq == "attack"] == "threat"))
## [1] 0.00257

In either of the cases, the co-operation patterns in the data you provided seem to be far from random. Notice that in both cases we ignore the autocorrelated nature of this data, we are rather asking: if we picked random point in time when male was attacking, would female be less or more likely to make treats at the same time?

Since you seem to be talking about causality ("when ... then"), while conducting permutation test you may be interested in comparing males behaviour in $t-1$ time to females behaviour at $t$ time (what was females' "reaction" to males behaviour?), but this is something that you have to ask yourself. Permutation tests are flexible and can be easily adapted to the kind of problems you seem to be describing.

Answer 2

You can think of your data in terms of bivariate Markov chain. You have two different variables $X$ for females and $Y$ for males, that describe stochastic process of changes in $X$ and $Y$ at time $t$ to one of four different states. Let's denote by $X_{t-1,i} \rightarrow X_{t,j}$ transition for $X$ from $t-1$ to $t$ time, from $i$-th to $j$-th state. In this case, transition in time to another state is conditional on previous state in $X$ and in $Y$:

$$ \Pr( X_{t-1,i} \rightarrow X_{t,j} ) = \Pr(X_{t,j} | X_{t-1,i},Y_{t-1,k}) \\ \Pr( Y_{t-1,h} \rightarrow Y_{t,k} ) = \Pr(Y_{t,h} | Y_{t-1,k},X_{t-1,i})$$

Transition probabilities can be easily calculated by counting transition histories and normalizing the probabilities afterwards:

states <- c("absent", "present", "attack", "threat")
# data is stored in 3-dimensional array, initialized with
# a very small "default" non-zero count to avoid zeros.
female_counts <- male_counts <- array(1e-16, c(4,4,4), list(states, states, states))
n <- length(male_seq)

for (i in 1:n) {
  male_counts[female_seq[i-1], male_seq[i-1], male_seq[i]] <- male_counts[female_seq[i-1], male_seq[i-1], male_seq[i]] + 1
  female_counts[male_seq[i-1], female_seq[i-1], female_seq[i]] <- female_counts[male_seq[i-1], female_seq[i-1], female_seq[i]] + 1
}

male_counts/sum(male_counts)
female_counts/sum(female_counts)

It can be also easyly simulated using marginal probabilities:

male_sim <- female_sim <- "absent"

for (i in 2:nsim) {
  male_sim[i] <- sample(states, 1, prob = male_counts[female_sim[i-1], male_sim[i-1], ])
  female_sim[i] <- sample(states, 1, prob = female_counts[male_sim[i-1], female_sim[i-1], ])
}

Result of such simulation is plotted below.

Moreover, it can be used to make one-step-ahead predictions:

male_pred <- female_pred <- NULL

for (i in 2:n) {
  curr_m <- male_counts[female_seq[i-1], male_seq[i-1], ]
  curr_f <- female_counts[male_seq[i-1], female_seq[i-1], ]
  male_pred[i] <- sample(names(curr_m)[curr_m == max(curr_m)], 1)
  female_pred[i] <- sample(names(curr_f)[curr_f == max(curr_f)], 1)
}

with 69-86% accuracy on the data you provided:

> mean(male_seq == male_pred, na.rm = TRUE)
[1] 0.8611111
> mean(female_seq == female_pred, na.rm = TRUE)
[1] 0.6944444

If the transitions occurred randomly, the transition probabilities would follow discrete uniform distribution. This is not a proof, but can serve as a way of thinking about your data using a simple model.