Bastimentos Aggression

Preliminaries

Load libraries

library(knitr)       #for R Markdown
library(MASS)        #for negative binomial glm (glm.nb)
library(emmeans)     #for posthoc
library(car)         #for Anova
library(survival)    #for survival analysis
library(MuMIn)       #for model selection (dredge)
library(ggfortify)   #for plotting survival analysis
library(ggsignif)    #for labeling significance in ggplots
library(tidyverse)   #for data processing, put last to avoid function masking

Data processing

Correct variable type and make derived variables

data_raw <- read.csv("Bast_Agg_final.csv",
                 
                 # to avoid reading errors
                 fileEncoding="UTF-8-BOM", na.strings = "")


data_raw <- data_raw %>%
  
  # correct categorical data type
  mutate_at(c("ID","success","pop","experimenter","color","model",
              "perch_type","start_call","conspp","track",
              "approach","challenge","attack","call"), as.factor) %>%
  
  # make a male color variable with just red and yellow
  mutate(color2 = recode(color, Y = "Y", O = "R", R = "R" ))




# make a dataset with only successful trials for analysis 
data <- data_raw %>% filter(success=="S")

Summary Stats

Number of performed assays

data_raw %>% select(color2, model) %>% table()

##       model
## color2  R  Y
##      R 34 27
##      Y 23 25

Number successful assays

data %>% select(color2, model) %>% table()

##       model
## color2  R  Y
##      R 25 22
##      Y 20 21

Percentage that did each behavior

prop <- data.frame(sum(data$approach=="Y")/nrow(data),
                   sum(data$call=="Y")/nrow(data),
                   sum(data$challenge=="Y")/nrow(data),
                   sum(data$attack=="Y")/nrow(data))

colnames(prop) <- c("approach", "call", "challenge", "attack")

kable(prop, digits = 3)

approach	call	challenge	attack
0.92	0.898	0.795	0.58

Main analyses

All model started with main effects and interaction between male color (color2) and model color (model); perch height (perch_height) and conspecific presence (conspp) as covariates. Non-significant covariates and interaction were dropped based on AIC using model selection (MuMIn::dredge). color2 and model were always kept in the final model

Likelihood to approach

Final model: approach ~ color2 * model

However, the SE of estimates is really high, probably because of sample size of frogs that didn’t approach is really small (7/101). Probably don’t want to make anything of this analysis.

# model selection

# mod <- glm(approach ~ color2 * model + perch_height + conspp,
#            data = data, family = binomial(),
#            na.action = "na.fail")
# MuMIn::dredge(mod)


mod <- glm(approach ~ color2 * model ,
           data = data, family = binomial())

summary(mod)

## 
## Call:
## glm(formula = approach ~ color2 * model, family = binomial(), 
##     data = data)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.18993   0.00008   0.43660   0.44740   0.57012  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)
## (Intercept)       19.57    2150.80   0.009    0.993
## color2Y          -17.83    2150.80  -0.008    0.993
## modelY           -17.26    2150.80  -0.008    0.994
## color2Y:modelY    17.78    2150.80   0.008    0.993
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 48.868  on 87  degrees of freedom
## Residual deviance: 43.521  on 84  degrees of freedom
## AIC: 51.521
## 
## Number of Fisher Scoring iterations: 18

Likelihood to call

Final model: call ~ color2 + model no sigificant effects

# model selection

# mod <- glm(call ~ color2 * model + perch_height + conspp,
#            data = data, family = binomial(),
#            na.action = "na.fail")
# MuMIn::dredge(mod)


mod <- glm(call ~ color2 + model ,
           data = data, family = binomial())

summary(mod)

## 
## Call:
## glm(formula = call ~ color2 + model, family = binomial(), data = data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2690   0.3982   0.4417   0.4843   0.5361  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   2.2783     0.6032   3.777 0.000159 ***
## color2Y      -0.4108     0.7090  -0.579 0.562344    
## modelY        0.2167     0.7097   0.305 0.760095    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 58.088  on 87  degrees of freedom
## Residual deviance: 57.672  on 85  degrees of freedom
## AIC: 63.672
## 
## Number of Fisher Scoring iterations: 5

Likelihood to challenge

Final model: challenge ~ color2 * model

# model selection

# mod <- glm(challenge ~ color2 * model + perch_height + conspp,
#            data = data, family = binomial(),
#            na.action = "na.fail")
# MuMIn::dredge(mod)


mod <- glm(challenge ~ color2 * model ,
           data = data, family = binomial())

summary(mod)

## 
## Call:
## glm(formula = challenge ~ color2 * model, family = binomial(), 
##     data = data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5373   0.2857   0.7375   0.7981   0.8446  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)   
## (Intercept)       3.178      1.020   3.115  0.00184 **
## color2Y          -2.331      1.131  -2.061  0.03933 * 
## modelY           -2.197      1.127  -1.950  0.05124 . 
## color2Y:modelY    2.513      1.331   1.888  0.05896 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 89.169  on 87  degrees of freedom
## Residual deviance: 81.666  on 84  degrees of freedom
## AIC: 89.666
## 
## Number of Fisher Scoring iterations: 5

Anova table

Anova(mod, type = 3)

## Analysis of Deviance Table (Type III tests)
## 
## Response: challenge
##              LR Chisq Df Pr(>Chisq)  
## color2         6.0686  1    0.01376 *
## model          5.3817  1    0.02035 *
## color2:model   4.3263  1    0.03753 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Posthoc test: Red males are more likely to challenge red models, and yellow males have no bias (marginal significance)

emmeans(mod, ~model|color2) %>% pairs() %>% test(by = NULL, adjust = "mvt")

##  contrast color2 estimate    SE  df z.ratio p.value
##  R - Y    R         2.197 1.127 Inf   1.950  0.0998
##  R - Y    Y        -0.316 0.708 Inf  -0.446  0.8812
## 
## Results are given on the log odds ratio (not the response) scale. 
## P value adjustment: mvt method for 2 tests

Figure

# make dataset for plotting
data_plot <- data %>% group_by(color2, model) %>%
  summarise(n = n(),
            prop = sum(challenge=="Y")/n(),
            SE = sqrt(prop*(1-prop)/n())
            )

## `summarise()` has grouped output by 'color2'. You can override using the
## `.groups` argument.

# alternative: use model estimates - in this case gives the same result
# emmeans(mod, ~color2|model, type = "response") %>% summary() %>% tibble


plot_challenge <- 
  
  ggplot(data_plot, aes(x = color2, y = prop, fill = model) )+
  
  # bars and error bars
  geom_bar(stat = "identity", width = 0.7, position = position_dodge(0.75),
           color = "black") +
  geom_errorbar(aes(ymin = prop-SE, ymax = prop+SE),
                width = 0.3, position = position_dodge(0.75),
                color = "black") +

  # axes and legends
  labs(x = "Male color", y = "Probability to challenge") +
  
  scale_x_discrete(labels = c("Red", "Yellow")) +
  
  scale_fill_manual(values=c('red3',"goldenrod2"),
                    name = "Decoy color",
                    labels = c("Red", "Yellow")) +
  
  scale_y_continuous(expand = c(0,0), limit = c(0,1),
                     labels = scales::percent) +
  
  # adjust element themes
  theme_classic(base_size = 15)

plot_challenge

Likelihood to attack

Final model: attack ~ color2 * model

# model selection

# mod <- glm(attack ~ color2 * model + perch_height + conspp,
#            data = data, family = binomial(),
#            na.action = "na.fail")
# MuMIn::dredge(mod)


mod <- glm(attack ~ color2 * model ,
           data = data, family = binomial())

summary(mod)

## 
## Call:
## glm(formula = attack ~ color2 * model, family = binomial(), data = data)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6894  -0.9948   0.7409   0.9794   1.5134  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)   
## (Intercept)      1.1527     0.4683   2.461  0.01384 * 
## color2Y         -0.7472     0.6539  -1.143  0.25319   
## modelY          -1.9148     0.6548  -2.924  0.00345 **
## color2Y:modelY   1.9949     0.9160   2.178  0.02942 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 119.76  on 87  degrees of freedom
## Residual deviance: 109.91  on 84  degrees of freedom
## AIC: 117.91
## 
## Number of Fisher Scoring iterations: 4

Anova table

Anova(mod, type = 3)

## Analysis of Deviance Table (Type III tests)
## 
## Response: attack
##              LR Chisq Df Pr(>Chisq)   
## color2         1.3243  1   0.249813   
## model          9.5473  1   0.002002 **
## color2:model   4.8738  1   0.027268 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Posthoc test: Red males are more likely to attack red models, and yellow males have no bias

emmeans(mod, ~model|color2) %>% pairs() %>% test(by = NULL, adjust = "mvt")

##  contrast color2 estimate    SE  df z.ratio p.value
##  R - Y    R          1.91 0.655 Inf   2.924  0.0069
##  R - Y    Y         -0.08 0.641 Inf  -0.125  0.9901
## 
## Results are given on the log odds ratio (not the response) scale. 
## P value adjustment: mvt method for 2 tests

Figure

# make dataset for plotting
data_plot <- data %>% group_by(color2, model) %>%
  summarise(n = n(),
            prop = sum(attack=="Y")/n(),
            SE = sqrt(prop*(1-prop)/n())
            )

## `summarise()` has grouped output by 'color2'. You can override using the
## `.groups` argument.

plot_attack <- 
  
  ggplot(data_plot, aes(x = color2, y = prop, fill = model) )+
  
  # bars and error bars
  geom_bar(stat = "identity", width = 0.7, position = position_dodge(0.75),
           color = "black") +
  geom_errorbar(aes(ymin = prop-SE, ymax = prop+SE),
                width = 0.3, position = position_dodge(0.75),
                color = "black") +
  
    
  #label significance
  geom_signif(y_position = 0.91, xmin = 0.8125, xmax = 1.1875, 
              tip_length = 0.03, textsize = 5,
              annotation = "*")+
  
  # axes and legends
  labs(x = "Male color", y = "Probability to attack") +
  
  scale_x_discrete(labels = c("Red", "Yellow")) +
  
  scale_fill_manual(values=c('red3',"goldenrod2"),
                    name = "Decoy color",
                    labels = c("Red", "Yellow")) +
  
  scale_y_continuous(expand = c(0,0), limit = c(0,1),
                     labels = scales::percent) +
  
  # adjust element themes
  theme_classic(base_size = 15)

plot_attack

Number of calls

final model: call_num ~ color2 * model + conspp

data_call <- data %>% filter(call=="Y")

# model selection
# mod <- glm.nb(call_num ~ color2 * model + perch_height + conspp + offset(log(total_time)),
#            data = data_call,
#            na.action = "na.fail")
# MuMIn::dredge(mod)


mod <- glm.nb(call_num ~ color2 * model + conspp + offset(log(total_time)),
           data = data_call)

summary(mod)

## 
## Call:
## glm.nb(formula = call_num ~ color2 * model + conspp + offset(log(total_time)), 
##     data = data_call, init.theta = 5.262153966, link = log)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.91542  -0.49705  -0.00933   0.59736   1.79144  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -3.3453     0.1028 -32.536  < 2e-16 ***
## color2Y         -0.1932     0.1627  -1.188  0.23503    
## modelY          -0.4560     0.1570  -2.905  0.00367 ** 
## consppY         -0.2833     0.1675  -1.692  0.09068 .  
## color2Y:modelY   0.3974     0.2328   1.707  0.08786 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(5.2622) family taken to be 1)
## 
##     Null deviance: 100.641  on 78  degrees of freedom
## Residual deviance:  87.464  on 74  degrees of freedom
## AIC: 548.5
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  5.26 
##           Std. Err.:  1.22 
## 
##  2 x log-likelihood:  -536.50

Anova(mod, type = 3)

## Analysis of Deviance Table (Type III tests)
## 
## Response: call_num
##              LR Chisq Df Pr(>Chisq)   
## color2         1.3746  1   0.241018   
## model          8.3054  1   0.003953 **
## conspp         2.7545  1   0.096979 . 
## color2:model   2.8700  1   0.090244 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

emmeans(mod, ~model|color2) %>% pairs() %>% test(by = NULL, adjust = "mvt")

##  contrast color2 estimate    SE  df z.ratio p.value
##  R - Y    R        0.4560 0.157 Inf   2.905  0.0073
##  R - Y    Y        0.0586 0.172 Inf   0.340  0.9291
## 
## Results are averaged over the levels of: conspp 
## Results are given on the log (not the response) scale. 
## P value adjustment: mvt method for 2 tests

Intepretation: red males called more at red models, and yellow males showed no bias

Figure

plot_call_num <-
  ggplot(data %>% filter(call == "Y"), 
         aes(x = color2, y = call_num, fill = model )) +
  
  # boxplot and error bars (to add wiskers so they look better)
  stat_boxplot(geom = "errorbar", width = 0.3, 
               position = position_dodge(width = 0.75),
               color = "black") +
  geom_boxplot(color = "black", outlier.shape = NA) +
  geom_point(position = position_jitterdodge(jitter.width = 0.3, dodge.width = 0.75),
             pch=1, size = 1.5, color = "grey40")+
  
  
  #label significance
  geom_signif(y_position = 37, xmin = 0.8125, xmax = 1.1875, 
              tip_length = 0.03, textsize = 5,
              annotation = "**")+
  
  
  # axes and legends
  labs(x = "Male color", y = "Number of call bouts") +
  
  scale_x_discrete(labels = c("Red", "Yellow")) +
  
  scale_fill_manual(values=c('red3',"goldenrod2"),
                    name = "Decoy color",
                    labels = c("Red", "Yellow")) +

  # adjust element themes
  theme_classic(base_size = 15, base_line_size = 1)


plot_call_num

Number of attacks

final model: attack_num ~ color2 + model

no difference in attack number by either model color or male color.

data_attack <- data %>% filter(attack=="Y")

# model selection
# mod <- glm.nb(attack_num ~ color2 * model + perch_height + conspp,
#            data = data.attack ,
#            na.action = "na.fail")
# MuMIn::dredge(mod)

mod <- glm.nb(attack_num ~ color2 + model,
           data = data_attack)

summary(mod)

## 
## Call:
## glm.nb(formula = attack_num ~ color2 + model, data = data_attack, 
##     init.theta = 2.042149627, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.1161  -0.8990   0.1158   0.5021   1.1205  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   3.8539     0.1514  25.460   <2e-16 ***
## color2Y       0.1251     0.2071   0.604    0.546    
## modelY       -0.1145     0.2121  -0.540    0.589    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.0421) family taken to be 1)
## 
##     Null deviance: 56.405  on 50  degrees of freedom
## Residual deviance: 55.828  on 48  degrees of freedom
## AIC: 494.72
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.042 
##           Std. Err.:  0.410 
## 
##  2 x log-likelihood:  -486.716

Anova(mod)

## Analysis of Deviance Table (Type II tests)
## 
## Response: attack_num
##        LR Chisq Df Pr(>Chisq)
## color2  0.37812  1     0.5386
## model   0.29938  1     0.5843

Figure

plot_attack_num <-
  ggplot(data %>% filter(attack == "Y"), 
         aes(x = color2, y = attack_num, fill = model )) +
  
  # boxplot and error bars (to add wiskers so they look better)
  stat_boxplot(geom = "errorbar", width = 0.3, 
               position = position_dodge(width = 0.75),
               color = "black") +
  geom_boxplot(color = "black") +
  geom_point(position = position_jitterdodge(jitter.width = 0.3, dodge.width = 0.75),
             pch=1, size = 1.5, color = "grey40")+
  
  
  # axes and legends
  labs(x = "Male color", y = "Number of attacks") +
  
  scale_x_discrete(labels = c("Red", "Yellow")) +
  
  scale_fill_manual(values=c('red3',"goldenrod2"),
                    name = "Decoy color",
                    labels = c("Red", "Yellow")) +

  # adjust element themes
  theme_classic(base_size = 15, base_line_size = 1)


plot_attack_num

Latency to track

Final model: color2 + model + perch_height + conspp

# recode Y/N behaviors to 1 = Y, 0 = N for survival analysis
data_latency <- data %>% mutate(track = recode(track, "Y" = 1, "N" = 0))

# model selection
# mod <- coxph(Surv(track_lat, track) ~ color2*model + perch_height + conspp,
#                data = data_latency, na.action = "na.fail")
# dredge(mod)


mod <- coxph(Surv(track_lat, track) ~ color2 + model + perch_height + conspp,
               data = data_latency)

summary(mod)

## Call:
## coxph(formula = Surv(track_lat, track) ~ color2 + model + perch_height + 
##     conspp, data = data_latency)
## 
##   n= 88, number of events= 88 
## 
##                   coef exp(coef)  se(coef)      z Pr(>|z|)  
## color2Y       0.270605  1.310758  0.230914  1.172   0.2412  
## modelY       -0.077220  0.925687  0.219640 -0.352   0.7252  
## perch_height -0.002648  0.997355  0.001728 -1.533   0.1254  
## consppY      -0.808958  0.445322  0.341439 -2.369   0.0178 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##              exp(coef) exp(-coef) lower .95 upper .95
## color2Y         1.3108     0.7629    0.8336    2.0610
## modelY          0.9257     1.0803    0.6019    1.4237
## perch_height    0.9974     1.0027    0.9940    1.0007
## consppY         0.4453     2.2456    0.2281    0.8696
## 
## Concordance= 0.576  (se = 0.038 )
## Likelihood ratio test= 9.89  on 4 df,   p=0.04
## Wald test            = 8.56  on 4 df,   p=0.07
## Score (logrank) test = 8.78  on 4 df,   p=0.07

Anova(mod, type = 3)

## Analysis of Deviance Table (Type III tests)
##              LR Chisq Df Pr(>Chisq)  
## color2         1.3640  1    0.24285  
## model          0.1237  1    0.72504  
## perch_height   2.9042  1    0.08835 .
## conspp         6.4416  1    0.01115 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretation: when male’s perch is higher, or when there are conspecifics present, males tracked slower

Latency to approach

Final model: color2 + model

# recode Y/N behaviors to 1 = Y, 0 = N for survival analysis
# assign censored time (i.e. replace NA with trial time)
data_latency <- data %>% mutate(approach = recode(approach, "Y" = 1, "N" = 0),
                                approach_lat = ifelse(approach == 1, approach_lat, total_time))


# model selection
# mod <- coxph(Surv(approach_lat, approach) ~ color2 * model + perch_height + conspp,
#                data = data_latency, na.action = na.fail)
# dredge(mod)

mod <- coxph(Surv(approach_lat, approach) ~ color2 + model ,
               data = data_latency)
summary(mod)

## Call:
## coxph(formula = Surv(approach_lat, approach) ~ color2 + model, 
##     data = data_latency)
## 
##   n= 88, number of events= 81 
## 
##             coef exp(coef) se(coef)      z Pr(>|z|)
## color2Y  0.03890   1.03967  0.22613  0.172    0.863
## modelY  -0.02886   0.97155  0.22413 -0.129    0.898
## 
##         exp(coef) exp(-coef) lower .95 upper .95
## color2Y    1.0397     0.9618    0.6674     1.619
## modelY     0.9716     1.0293    0.6262     1.507
## 
## Concordance= 0.523  (se = 0.037 )
## Likelihood ratio test= 0.05  on 2 df,   p=1
## Wald test            = 0.05  on 2 df,   p=1
## Score (logrank) test = 0.05  on 2 df,   p=1

Latency to challenge

Final model: color2 * model + perch_height

# recode Y/N behaviors to 1 = Y, 0 = N for survival analysis
# assign censored time (i.e. replace NA with trial time)
data_latency <- data %>% mutate(challenge = recode(challenge, "Y" = 1, "N" = 0),
                                challenge_lat = ifelse(challenge == 1, challenge_lat, total_time))


# model selection
# mod <- coxph(Surv(challenge_lat, challenge) ~ color2 * model + perch_height + conspp,
#                data = data_latency, na.action = na.fail)
# dredge(mod)

mod <- coxph(Surv(challenge_lat, challenge) ~ color2 * model + perch_height ,
               data = data_latency)

summary(mod)

## Call:
## coxph(formula = Surv(challenge_lat, challenge) ~ color2 * model + 
##     perch_height, data = data_latency)
## 
##   n= 88, number of events= 70 
## 
##                     coef exp(coef)  se(coef)      z Pr(>|z|)  
## color2Y        -0.636522  0.529129  0.340750 -1.868   0.0618 .
## modelY         -0.580654  0.559532  0.324989 -1.787   0.0740 .
## perch_height   -0.003250  0.996755  0.002076 -1.566   0.1174  
## color2Y:modelY  0.845516  2.329181  0.491369  1.721   0.0853 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## color2Y           0.5291     1.8899    0.2713     1.032
## modelY            0.5595     1.7872    0.2959     1.058
## perch_height      0.9968     1.0033    0.9927     1.001
## color2Y:modelY    2.3292     0.4293    0.8891     6.102
## 
## Concordance= 0.607  (se = 0.039 )
## Likelihood ratio test= 7.2  on 4 df,   p=0.1
## Wald test            = 7.18  on 4 df,   p=0.1
## Score (logrank) test = 7.4  on 4 df,   p=0.1

Anova(mod, type = 3)

## Analysis of Deviance Table (Type III tests)
##              LR Chisq Df Pr(>Chisq)  
## color2         3.6253  1    0.05691 .
## model          3.2711  1    0.07051 .
## perch_height   3.1570  1    0.07560 .
## color2:model   2.9868  1    0.08395 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretation: marginal, but when male perch is higher they challenged slower

Posthoc: marginal, but red males challenged faster when model is red

emmeans(mod, ~model|color2) %>% pairs() %>% test(by = NULL, adjust = "mvt")

##  contrast color2 estimate    SE  df z.ratio p.value
##  R - Y    R         0.581 0.325 Inf   1.787  0.1425
##  R - Y    Y        -0.265 0.368 Inf  -0.719  0.7211
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: mvt method for 2 tests

Figure

data_plot <- 
  survfit(Surv(challenge_lat, challenge) ~ color2 + model,
          data = data_latency) %>%
  fortify(surv.connect = TRUE) %>%
  separate(strata, c("color2","model"),",") %>%
  drop_na()



plot_challenge_lat <- 
  
  ggplot(data_plot, aes(time, 1-surv, color = color2, lty = model)) +
  
  #survival plot
  geom_step(size = 1) +
  
  
  # axes and legends
  labs(x = "Latency to challenge (s)", y = "Proportion challenging") +

  scale_linetype_manual(values = c(1,2),
                        name = "Decoy color",
                        labels = c("Red", "Yellow") )+
  
  scale_color_manual(values=c('red3',"goldenrod2"),
                     name = "Male color",
                     labels = c("Red", "Yellow")) +
  

  
  scale_y_continuous(expand = c(0,0), limit = c(0,1),
                     labels = scales::percent) +
  
  # adjust element themes
  theme_classic(base_size = 15)

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

plot_challenge_lat

Latency to attack

Final model: color2 * model + perch_height

# recode Y/N behaviors to 1 = Y, 0 = N for survival analysis
# assign censored time (i.e. replace NA with trial time)
data_latency <- data %>% mutate(attack = recode(attack , "Y" = 1, "N" = 0),
                                attack_lat = ifelse(attack == 1, attack_lat, total_time))

# model selection
# mod <- coxph(Surv(attack_lat, attack ) ~ color2 * model + perch_height + conspp,
#                data = data_latency, na.action = na.fail)
# dredge(mod)

mod <- coxph(Surv(attack_lat, attack ) ~ color2 * model + perch_height,
               data = data_latency, na.action = na.fail)
summary(mod)

## Call:
## coxph(formula = Surv(attack_lat, attack) ~ color2 * model + perch_height, 
##     data = data_latency, na.action = na.fail)
## 
##   n= 88, number of events= 51 
## 
##                     coef exp(coef)  se(coef)      z Pr(>|z|)  
## color2Y        -0.173405  0.840797  0.373414 -0.464   0.6424  
## modelY         -1.050104  0.349901  0.447808 -2.345   0.0190 *
## perch_height   -0.004471  0.995539  0.003211 -1.392   0.1638  
## color2Y:modelY  1.195178  3.304145  0.600879  1.989   0.0467 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## color2Y           0.8408     1.1893    0.4044    1.7480
## modelY            0.3499     2.8579    0.1455    0.8416
## perch_height      0.9955     1.0045    0.9893    1.0018
## color2Y:modelY    3.3041     0.3027    1.0176   10.7282
## 
## Concordance= 0.62  (se = 0.038 )
## Likelihood ratio test= 10.69  on 4 df,   p=0.03
## Wald test            = 8.56  on 4 df,   p=0.07
## Score (logrank) test = 9.29  on 4 df,   p=0.05

Anova(mod, type = 3)

## Analysis of Deviance Table (Type III tests)
##              LR Chisq Df Pr(>Chisq)  
## color2         0.2179  1    0.64066  
## model          6.1594  1    0.01307 *
## perch_height   2.6516  1    0.10344  
## color2:model   4.1185  1    0.04242 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

posthoc: red males attacked earlier when model is red, and yellow males don’t have a bias

emmeans(mod, ~model|color2) %>% pairs() %>% test(by = NULL, adjust = "mvt")

##  contrast color2 estimate    SE  df z.ratio p.value
##  R - Y    R         1.050 0.448 Inf   2.345  0.0377
##  R - Y    Y        -0.145 0.401 Inf  -0.361  0.9204
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: mvt method for 2 tests

Figure

data_plot <- 
  survfit(Surv(attack_lat, attack) ~ color2 + model,
          data = data_latency) %>%
  fortify(surv.connect = TRUE) %>%
  separate(strata, c("color2","model"),",") %>%
  drop_na()



plot_attack_lat <- 
  
  ggplot(data_plot, aes(time, 1-surv, color = color2, lty = model)) +
  
  #survival plot
  geom_step(size = 1) +
  
  
  # axes and legends
  labs(x = "Time since start of trial (sec)", y = "Proportion attacking") +

  scale_linetype_manual(values = c(1,2),
                        name = "Decoy color",
                        labels = c("Red", "Yellow") )+
  
  scale_color_manual(values=c('red3',"goldenrod2"),
                     name = "Male color",
                     labels = c("Red", "Yellow")) +
  

  
  scale_y_continuous(expand = c(0,0), limit = c(0,1),
                     labels = scales::percent) +
  
  # adjust element themes
  theme_classic(base_size = 15)

plot_attack_lat

Multipanel figure options

attack & challenge

panel <- 
  egg::ggarrange(plot_challenge_lat+ theme(legend.position = "none"),
                 plot_attack_lat + theme(legend.key.width=unit(8, 'mm')),
                 plot_challenge + theme(legend.position = "none"),
                 plot_attack,
               labels = c("A","B","C","D"),
               nrow = 2, ncol = 2)

attack only

panel <- 
  egg::ggarrange(plot_attack_lat + theme(legend.key.width=unit(8, 'mm')),
                 plot_attack,
                 plot_attack_num,
               labels = c("A","B","C"),
               nrow = 3, ncol = 1)

for saving as a high resolution figure

ggsave("Bast_Agg_Fig.jpg", panel, width = 5, height = 10)

R Package References

Bartoń, Kamil. 2022. MuMIn: Multi-Model Inference. https://CRAN.R-project.org/package=MuMIn.

Fox, John, and Sanford Weisberg. 2019. An R Companion to Applied Regression. Third. Thousand Oaks CA: Sage. https://socialsciences.mcmaster.ca/jfox/Books/Companion/.

Fox, John, Sanford Weisberg, and Brad Price. 2022. Car: Companion to Applied Regression. https://CRAN.R-project.org/package=car.

Lenth, Russell V. 2023. Emmeans: Estimated Marginal Means, Aka Least-Squares Means. https://github.com/rvlenth/emmeans.

Ripley, Brian. 2022. MASS: Support Functions and Datasets for Venables and Ripley’s MASS. http://www.stats.ox.ac.uk/pub/MASS4/.

Terry M. Therneau, and Patricia M. Grambsch. 2000. Modeling Survival Data: Extending the Cox Model. New York: Springer.

Therneau, Terry M. 2022. Survival: Survival Analysis. https://github.com/therneau/survival.

Venables, W. N., and B. D. Ripley. 2002. Modern Applied Statistics with s. Fourth. New York: Springer. https://www.stats.ox.ac.uk/pub/MASS4/.

Bastimentos Aggression

Yusan Yang

2024-12-09

Preliminaries

Load libraries

Data processing

Summary Stats

Number of performed assays

Number successful assays

Percentage that did each behavior

Main analyses

Likelihood to approach

Likelihood to call

Likelihood to challenge

Figure

Likelihood to attack

Figure

Number of calls

Figure

Number of attacks

Figure

Latency to track

Latency to approach

Latency to challenge

Figure

Latency to attack

Figure

Multipanel figure options

attack & challenge

attack only

R Package References