Chapter 4 Model prediction

cages=c("C03","C04","C05","C06","C10","C11","C12","C13")

dominance_predictions <- data.frame()
  
for(cage in cages){
    #Filter table by cage
    cage_data <- tube_test_data %>% filter(cage_id == cage)
    #Model with temporal autocorrelation (NULL if the model does not converge)
    M3 <- tryCatch({glmmTMB(win_loss_binomial~1+ar1(time+0|mouse), data = cage_data,family = binomial)}, warning = function(w){NULL})
    #Model without temporal autocorrelation
    M4 <- glmmTMB(win_loss_binomial~1+(1|mouse), data = cage_data,family = binomial)
    
    if(!is.null(M3)){
          cage_model <- data.frame(ggpredict(M3, type = "random", terms = c("mouse","time")), model="M3", cage=cage)
    }else{
          cage_model <- data.frame(ggpredict(M4, type = "random", terms = "mouse"), model="M4", cage=cage)
          cage_model <- bind_rows(cage_model %>% mutate(group = factor(1)),
            cage_model %>% mutate(group = factor(2)),
            cage_model %>% mutate(group = factor(3)),
            cage_model %>% mutate(group = factor(4)),
            cage_model %>% mutate(group = factor(5)),
            cage_model %>% mutate(group = factor(6)))}
    
    # Append predictions to table
    dominance_predictions <- bind_rows(dominance_predictions,cage_model)
}

dominance_predictions <- dominance_predictions %>%
      mutate(group = factor(case_when(
          group == 1 ~ "t1",
          group == 2 ~ "t2",
          group == 3 ~ "t3",
          group == 4 ~ "t4",
          group == 5 ~ "t5",
          group == 6 ~ "t6",
          TRUE ~ NA))) %>%
      rename(mouse=x, dominance=predicted, time=group)

4.1 Dominance visualisation

dominance_predictions %>%
    mutate(time = as.numeric(substr(time, 2, 2))) %>%
    mutate(mouse_number = substr(mouse, 5, 5)) %>%
    ggplot(aes(x=time,y=dominance, group=mouse, color=mouse_number)) +
      geom_line(size=1) +
      scale_color_manual(values=c("#D5CFCF","#86b9d4","#a6897c","#2e2e33","#2e7287"))+
      scale_y_continuous(breaks = c(0,0.25,0.5,0.75,1)) +
      scale_x_continuous(breaks = c(1:6)) +
      facet_wrap(cage ~ ., ncol=4) + 
      theme_minimal()+
      theme(legend.position = "none")

4.2 Hierarchy stability

cages=c("C03","C04","C05","C06","C10","C11","C12","C13")

hierarchy_stability <- c()
  
for(cage_code in cages){
    dominance_predictions_cage <- dominance_predictions %>% 
        filter(cage == cage_code) %>% 
        select(mouse,dominance,time) %>% 
        pivot_wider(names_from="time",values_from="dominance")%>% 
        select(-mouse)
    
    correlation_coefficients <- numeric(length = 5)
    for (i in 1:5) {
      correlation_coefficients[i] <- cor(dominance_predictions_cage[, i], dominance_predictions_cage[, i + 1], method = "spearman")
    }
    hierarchy_stability <- c(hierarchy_stability,mean(correlation_coefficients))
}

names(hierarchy_stability) <- cages

hierarchy_stability %>% 
    as.data.frame() %>% t() %>% as.data.frame() %>%
    tt(digits = 2) |> 
      style_tt(i = 0, color = "white", background = "#a52929")

C03	C04	C05	C06	C10	C11	C12	C13
0.96	0.86	0.94	0.95	1	0.72	1	1

mean(hierarchy_stability)

[1] 0.9292857

sd(hierarchy_stability)

[1] 0.09657396

hierarchy_stability %>%
    as.data.frame() %>% t() %>% t() %>% as.data.frame() %>%
    rownames_to_column(var="cage") %>%
    rename(cage=1,stability=2) %>%
    ggplot(aes(y=fct_rev(cage),x=stability))+
      geom_col(fill="#999999")+
      theme_minimal()

### Between environmental treatments

hierarchy_stability %>% 
    as.data.frame() %>% t() %>% t() %>% as.data.frame() %>%
    rownames_to_column(var="cage") %>%
    rename(cage=1,stability=2) %>%
    mutate(treatment=c("variable","invariable","variable","variable","variable","variable","variable","invariable")) %>%
    group_by(treatment) %>%
    summarise(mean=mean(stability), sd=sd(stability)) %>%
    tt()

treatment	mean	sd
invariable	0.9300000	0.09899495
variable	0.9290476	0.10534169

hierarchy_stability %>% 
    as.data.frame() %>% t() %>% t() %>% as.data.frame() %>%
    rownames_to_column(var="cage") %>%
    rename(cage=1,stability=2) %>%
    mutate(treatment=c("variable","invariable","variable","variable","variable","variable","variable","invariable")) %>%
    ggplot(aes(x=treatment,y=stability,group=treatment))+
      geom_boxplot()

4.2.1 Between sexes

hierarchy_stability %>% 
    as.data.frame() %>% t() %>% t() %>% as.data.frame() %>%
    rownames_to_column(var="cage") %>%
    rename(cage=1,stability=2) %>%
    mutate(sex=c("male","male","male","male","female","female","female","female")) %>%
    group_by(sex) %>%
    summarise(mean=mean(stability), sd=sd(stability)) %>%
    tt()

sex	mean	sd
female	0.9300000	0.14000000
male	0.9285714	0.04648165

hierarchy_stability %>% 
    as.data.frame() %>% t() %>% t() %>% as.data.frame() %>%
    rownames_to_column(var="cage") %>%
    rename(cage=1,stability=2) %>%
    mutate(sex=c("male","male","male","male","female","female","female","female")) %>%
    ggplot(aes(x=sex,y=stability,group=sex))+
      geom_boxplot()