Packages

library(glmmTMB)
library(ggplot2)
library(car)
## Loading required package: carData
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.0     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ lubridate 1.9.5     ✔ tibble    3.3.1
## ✔ purrr     1.2.1     ✔ tidyr     1.3.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ dplyr::recode() masks car::recode()
## ✖ purrr::some()   masks car::some()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(emmeans)
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
library(DHARMa)
## Warning: package 'DHARMa' was built under R version 4.5.3
## This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library(knitr)
library(formattable)
## Warning: package 'formattable' was built under R version 4.5.3
library(patchwork)
## 
## Attaching package: 'patchwork'
## 
## The following object is masked from 'package:formattable':
## 
##     area
library(survival)
library(survminer)
## Loading required package: ggpubr
## 
## Attaching package: 'survminer'
## 
## The following object is masked from 'package:survival':
## 
##     myeloma
library(coxme)
## Warning: package 'coxme' was built under R version 4.5.3
## Loading required package: bdsmatrix
## 
## Attaching package: 'bdsmatrix'
## 
## The following object is masked from 'package:base':
## 
##     backsolve

Working Title

Female Avoidance behaviour vs Aggressive males

Authors

Anastasia J. Shavrova, Bruno A. Buzatto, & Michael M. Kasumovic

Description

In the bulb mite, Rhizoglyphus echinopus, two male morphs are present: a fighter and a scrambler male. The females are typically polyamorous, and hte more mates she has the higher her fecundity. To potentially control this, the fighter morph male, if given the chance, will try to eliminate rival males by using his fighter legs to pierce and kill them. There is some evidence that suggests that if a female is mated with a scrambler male this boosts her fecundity, and on the contrary if she is mated with a fighter male. However, recent evidence suggests that in lines selected for fighters, if a female of the fighter lines is mated with a fighter male this boosts her fecundity, and her fecundity is decreased if she is mated with a scrambler male from the scrambler selected lines.

Mating Behaviour

  1. Do females mate faster with scramblers or fighters the first and/or second time around?

Data

Data Description

  • trial = experimental replicate (1-8)
  • mate.order = First mating with first morph, second mating with second morph
  • male.trmnt = male morph treatment as of 4 treatments: “FF”, “FS”, “SS”, “SF”
  • time = time to initiate mating (in seconds)
  • replicate = replicate vial within an experiment (1-8)
  • mate.type = the type of mating pair for 2nd presentation: either the same (FF/SS) or different (SF/FS)
  • male.morph = male morph presented to female: “Fighter” or “Scrambler”
  • female_ID = ID # of female
  • mated.once = whether pair mated at once
  • mated.twice = whether pair mated a second time
  • female.replaced = if the female was replaced for first mating only
  • male.replaced = if male was replaced, generally for second mating
  • male_ID = unique male ID number
mating <- read.csv("mate.csv")
str(mating)
## 'data.frame':    364 obs. of  13 variables:
##  $ trial          : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ mate.order     : int  1 2 1 2 1 2 1 2 1 2 ...
##  $ male.trmnt     : chr  "FF" "FF" "FF" "FF" ...
##  $ time           : int  19 77 89 64 19 19 90 45 39 44 ...
##  $ replicate      : int  1 1 2 2 3 3 4 4 5 5 ...
##  $ mate.type      : chr  "same" "same" "same" "same" ...
##  $ male.morph     : chr  "Fighter" "Fighter" "Fighter" "Fighter" ...
##  $ female_ID      : int  1 1 2 2 3 3 4 4 5 5 ...
##  $ mated.once     : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ mated.twice    : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ female.replaced: int  0 0 0 0 1 1 0 0 0 0 ...
##  $ male.replaced  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ male_ID        : int  1 2 3 4 5 6 7 8 9 10 ...

Summary Stats

female.replaced <- mating %>%
  group_by(male.trmnt) %>%
  count(female.replaced) %>%
  filter(female.replaced == "1")

Mating Latency in first vs second mating with different morphs

First Mating

#first mating
first.mating <- mating %>%
  filter(mate.order == "1")

first.mate.latency <- glmmTMB(log(time)~male.morph + (1|trial), data=first.mating)

plot(simulateResiduals(first.mate.latency)) #residuals all good after log transforming time

Anova(first.mate.latency) #no diff
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: log(time)
##             Chisq Df Pr(>Chisq)
## male.morph 0.0456  1     0.8309
summary(first.mate.latency)
##  Family: gaussian  ( identity )
## Formula:          log(time) ~ male.morph + (1 | trial)
## Data: first.mating
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     331.7     344.4    -161.9     323.7       172 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  trial    (Intercept) 0.05665  0.2380  
##  Residual             0.34370  0.5863  
## Number of obs: 176, groups:  trial, 8
## 
## Dispersion estimate for gaussian family (sigma^2): 0.344 
## 
## Conditional model:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          3.39087    0.10469   32.39   <2e-16 ***
## male.morphScrambler -0.01891    0.08858   -0.21    0.831    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
emmeans(first.mate.latency, ~male.morph, type="response")
##  male.morph response   SE  df asymp.LCL asymp.UCL
##  Fighter        29.7 3.11 Inf      24.2      36.5
##  Scrambler      29.1 3.07 Inf      23.7      35.8
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale
#second mating
second.mating <- mating %>%
  filter(mate.order == "2")

second.mate.latency <- glmmTMB(log(time)~male.trmnt + (1|trial), data=second.mating)
plot(simulateResiduals(second.mate.latency))#residuals all good after log transforming time

Anova(second.mate.latency) #no diff  
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: log(time)
##             Chisq Df Pr(>Chisq)
## male.trmnt 2.2356  3      0.525
second.latency.table <- as.data.frame(
  emmeans(second.mate.latency, ~ male.trmnt, type = "response")
) %>%
  transmute(
    `Mating sequence` = male.trmnt,
    `Mating latency (sec)` = round(response, 2),
    `Standard error` = round(SE, 2)
  )

formattable(
  second.latency.table,
  align = c("c", "c", "c"),
  table.attr = 'style="width: 400px; margin-left: auto; margin-right: auto;"'
)
Mating sequence Mating latency (sec) Standard error
FF 36.56 4.36
FS 34.11 4.04
SF 32.72 3.88
SS 39.67 4.62

Results No differences in mating latency

Plot

p1_cols <- c(
  "Fighter" = "#9E9E9E",   
  "Scrambler" = "#D95F02" 
)

p2_cols <- c(
  "FF" = "#4C83B6",  
  "FS" = "#1B7F3A",  
  "SF" = "#E6D36F",  
  "SS" = "#CC6677"  
)


p1 <- ggplot(first.mating, aes(x = male.morph, y = time, fill = male.morph)) +
  geom_boxplot(outlier.shape = NA, width = 0.45, alpha = 0.9) +
  geom_jitter(
    aes(fill = male.morph),
    shape = 21,
    colour = "black",
    width = 0.12,
    alpha = 0.65,
    size = 2,
    stroke = 0.4
  ) +
  scale_fill_manual(values = p1_cols) +
  labs(
    x = "Pairing",
    y = "Time to Mate (sec)",
    title = "A"
  ) +
  theme_classic(base_size = 14) +
  theme(
    legend.position = "none",
    plot.title = element_text(face = "bold", hjust = -0.08),
    axis.title = element_text(face = "bold")
  )
p2 <- ggplot(second.mating, aes(x = male.trmnt, y = time, fill = male.trmnt)) +
  geom_boxplot(outlier.shape = NA, width = 0.45, alpha = 0.9) +
  geom_jitter(
    aes(fill = male.trmnt),
    shape = 21,
    colour = "black",
    width = 0.12,
    alpha = 0.65,
    size = 2,
    stroke = 0.4
  ) +
  scale_fill_manual(values = p2_cols) +
  scale_x_discrete(
    labels = c(
      "FF" = "Fighter first\nFighter second",
      "FS" = "Fighter first\nScrambler second",
      "SF" = "Scrambler first\nFighter second",
      "SS" = "Scrambler first\nScrambler second"
    )
  ) +
  labs(
    x = "Pairing",
    y = "Time to Mate (sec)",
    title = "B"
  ) +
  theme_classic(base_size = 14) +
  theme(
    legend.position = "none",
    plot.title = element_text(face = "bold", hjust = -0.08),
    axis.title = element_text(face = "bold"),
    axis.text.x = element_text(size = 11)
  )

p1/p2
## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 29 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 29 rows containing missing values or values outside the scale range
## (`geom_point()`).

Monogamous Matings

Data

single.mate.fecund <- read.csv("egg_count_allfem.csv")
str(single.mate.fecund)
## 'data.frame':    673 obs. of  5 variables:
##  $ trial    : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ day      : int  2 2 2 2 2 2 2 2 2 4 ...
##  $ mate.type: chr  "F" "F" "F" "F" ...
##  $ egg.count: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ female_ID: int  1 7 10 18 22 23 25 33 36 1 ...

Any differences in total fecundity?

#create total fecundity data (counting all the eggs per female)
total_fecund_mono <- single.mate.fecund %>%
  group_by(female_ID, mate.type, trial) %>%
  tally(egg.count)

#model 
mod_t_fecund_mono <- glmmTMB(n~mate.type + (1|trial), data=total_fecund_mono)
plot(simulateResiduals(mod_t_fecund_mono)) #all good

Anova(mod_t_fecund_mono) #no diff
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: n
##            Chisq Df Pr(>Chisq)
## mate.type 0.6378  1     0.4245
emmeans(mod_t_fecund_mono, ~mate.type)
##  mate.type emmean   SE  df asymp.LCL asymp.UCL
##  F            251 38.2 Inf       176       325
##  S            220 38.3 Inf       145       295
## 
## Confidence level used: 0.95

Any differences in survival?

# data change to grab last day of females survival
survival.mono <- single.mate.fecund %>%
  group_by(female_ID, mate.type) %>%
  top_n(1, day)

#plot
single.surv <- survfit(Surv(day)~mate.type + (1|female_ID), data=survival.mono)

surv.sing <- ggsurvplot(single.surv, data=survival.mono,
           legend.title="",
           legend.labs=c("Fighter", "Scrambler"),
           conf.int=TRUE,
           conf.int.style="step",
           palette=c("#999999", "#D55E00"),
           legend="bottom",
           xlim=c(0,100),
           ggtheme = theme_pubr(),
           font.x=15,
           font.y=15,
           font.legend=15)

#model
single.surv.mod.allfem <- coxme(Surv(day)~mate.type + egg.count + (1|female_ID) + (1|trial), data=survival.mono)
cox.zph(single.surv.mod.allfem) #model fine
##           chisq df     p
## mate.type 0.634  1 0.426
## egg.count 2.801  1 0.094
## GLOBAL    3.557  2 0.169
summary(single.surv.mod.allfem) #egg count has effect
## Mixed effects coxme model
##  Formula: Surv(day) ~ mate.type + egg.count + (1 | female_ID) + (1 | trial) 
##     Data: survival.mono 
## 
##   events, n = 81, 81
## 
## Random effects:
##       group  variable         sd    variance
## 1 female_ID Intercept 0.41724154 0.174090507
## 2     trial Intercept 0.01477843 0.000218402
##                   Chisq df         p   AIC    BIC
## Integrated loglik 17.51  4 1.541e-03  9.51  -0.07
##  Penalized loglik 39.58 12 8.457e-05 15.57 -13.17
## 
## Fixed effects:
##                coef exp(coef) se(coef)     z        p
## mate.typeS -0.29052   0.74788  0.25299 -1.15    0.251
## egg.count   0.05537   1.05693  0.01270  4.36 1.31e-05
single.surv.mod.grid <- ref_grid(single.surv.mod.allfem)

emmeans(regrid(single.surv.mod.grid), pairwise~mate.type)
## $emmeans
##  mate.type response    SE  df asymp.LCL asymp.UCL
##  F            1.154 0.144 Inf     0.872      1.44
##  S            0.863 0.111 Inf     0.647      1.08
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate    SE  df z.ratio p.value
##  F - S       0.291 0.255 Inf   1.142  0.2533

Any differences in egg laying rate (to peak and from peak)?

#data set up
no_egg <- total_fecund_mono %>%
  filter(n == 0) %>% #these females laid no eggs in their lifetime 
  pull(female_ID)

#removing these females
fecund_females_mono <- single.mate.fecund %>%
  filter(!female_ID %in% no_egg)

#finding each females peak day
peak.day.mono <- fecund_females_mono %>%
  group_by(female_ID, mate.type) %>%
    summarise(
    peak_day = day[which.max(egg.count)][1],
    .groups = "drop"
  )

#finding the "peak" day - day the of highest egg count
fecund_peak <- fecund_females_mono %>%
  group_by(female_ID) %>%
  mutate(
    peak_day  = day[which.max(egg.count)][1],
    peak_eggs = max(egg.count, na.rm = TRUE)
  ) %>%
  ungroup()

#subsetting data to peak
to_peak <- fecund_peak %>%
  group_by(female_ID) %>%
  filter(day >= min(day[egg.count > 0]) & day <= peak_day) %>%
  ungroup()

#setting "peak" value as zero, and days leading up to it in negative - normalizing the data
to_peak <- to_peak %>%
  group_by(female_ID) %>%
  mutate(
    rel_day_to_peak = day - peak_day
  ) %>%
  ungroup()

#subsetting data from peak
from_peak <- fecund_peak %>%
  filter(day >= peak_day) %>%
  ungroup()

from_peak <- from_peak %>%
  group_by(female_ID) %>%
  mutate(
    rel_day_from_peak = day - peak_day
  ) %>%
  ungroup()

Difference to peak

to_peak_mono <- glmmTMB(egg.count~poly(rel_day_to_peak,2) * mate.type + (1|female_ID), ziformula = ~1,
  family = "poisson", data=to_peak) 
plot(simulateResiduals(to_peak_mono))

testDispersion(simulateResiduals(to_peak_mono)) # data not overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.59957, p-value = 0.232
## alternative hypothesis: two.sided

Plot

# to peak
ggplot(to_peak, aes(x = day, y = egg.count, col = mate.type)) + geom_point() + theme_bw()

trying calculating slopes

#to peak
to_peak_summary <- single.mate.fecund %>%
  group_by(female_ID, mate.type, trial) %>%
  arrange(day, .by_group = TRUE) %>%
  summarise(
    first_day = min(day),
    first_eggs = egg.count[which.min(day)],
    peak_day = day[which.max(egg.count)][1],
    peak_eggs = max(egg.count, na.rm = TRUE),
    last_day = max(day),
    days_to_peak = peak_day - first_day,
    slope_to_peak = ifelse(
      days_to_peak == 0,
      NA,
      (peak_eggs - first_eggs) / days_to_peak
    ),
    .groups = "drop"
  )
str(to_peak_summary)
## tibble [81 × 10] (S3: tbl_df/tbl/data.frame)
##  $ female_ID    : int [1:81] 1 2 3 4 5 6 7 7 8 9 ...
##  $ mate.type    : chr [1:81] "F" "F" "F" "F" ...
##  $ trial        : int [1:81] 1 1 1 1 1 1 1 2 1 1 ...
##  $ first_day    : int [1:81] 2 2 2 2 2 2 2 37 2 2 ...
##  $ first_eggs   : int [1:81] 0 30 10 20 10 7 0 0 21 28 ...
##  $ peak_day     : int [1:81] 2 9 9 9 6 9 9 40 14 2 ...
##  $ peak_eggs    : int [1:81] 0 117 90 104 119 78 119 2 108 28 ...
##  $ last_day     : int [1:81] 17 14 20 41 22 29 20 44 25 6 ...
##  $ days_to_peak : int [1:81] 0 7 7 7 4 7 7 3 12 0 ...
##  $ slope_to_peak: num [1:81] NA 12.4 11.4 12 27.2 ...
mod_slope_to_peak <- glmmTMB(
  slope_to_peak ~ mate.type + factor(trial),
  data = to_peak_summary
)
plot(simulateResiduals(mod_slope_to_peak)) #fits all assumptions

summary(mod_slope_to_peak)
##  Family: gaussian  ( identity )
## Formula:          slope_to_peak ~ mate.type + factor(trial)
## Data: to_peak_summary
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     465.5     474.3    -228.7     457.5        63 
## 
## 
## Dispersion estimate for gaussian family (sigma^2): 54.1 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     11.8001     1.4770   7.989 1.36e-15 ***
## mate.typeS       0.1814     1.8216   0.100    0.921    
## factor(trial)2   1.4477     1.8467   0.784    0.433    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#from peak
from_peak_summary <- single.mate.fecund %>%
  group_by(female_ID, mate.type, trial) %>%
  arrange(day, .by_group = TRUE) %>%
  summarise(
    peak_day = day[which.max(egg.count)][1],
    peak_eggs = max(egg.count, na.rm = TRUE),
    last_day = max(day),
    last_eggs = egg.count[which.max(day)],
    days_after_peak = last_day - peak_day,
    decline_rate = ifelse(
      days_after_peak == 0,
      NA,
      (peak_eggs - last_eggs) / days_after_peak
    ),
    .groups = "drop"
  )
from_peak_summary2 <- from_peak_summary %>%
  filter(!is.na(decline_rate), days_after_peak > 0, peak_eggs > 0) #filterring NA's


mod_slope_from_peak <- glmmTMB(
  log(decline_rate+1) ~ mate.type + factor(trial),
  data = from_peak_summary2
)
plot(simulateResiduals(mod_slope_from_peak))

summary(mod_slope_from_peak)
##  Family: gaussian  ( identity )
## Formula:          log(decline_rate + 1) ~ mate.type + factor(trial)
## Data: from_peak_summary2
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     166.4     175.3     -79.2     158.4        64 
## 
## 
## Dispersion estimate for gaussian family (sigma^2): 0.601 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     2.26845    0.15506  14.630   <2e-16 ***
## mate.typeS     -0.14658    0.18986  -0.772    0.440    
## factor(trial)2 -0.05733    0.19214  -0.298    0.765    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot to peak

emm_to_peak <- emmeans(mod_slope_to_peak, ~ mate.type)

emm_to_peak_df <- as.data.frame(emm_to_peak)

ggplot(to_peak_summary, aes(x = mate.type, y = slope_to_peak, color = mate.type)) +
  geom_jitter(width = 0.12, alpha = 0.7, size = 2) +
  geom_pointrange(
    data = emm_to_peak_df,
    aes(
      x = mate.type,
      y = emmean,
      ymin = lower.CL,
      ymax = upper.CL
    ),
    inherit.aes = FALSE,
    color = "black",
    linewidth = 0.8
  ) +
  scale_x_discrete(labels = c("F" = "Fighter", "S" = "Scrambler")) +
  theme_classic() +
  labs(
    x = "Mate type",
    y = "Slope to peak egg laying (egg/per day)",
  ) +
  theme(legend.position = "none")
## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_point()`).

plot from peak

emm_from_peak <- emmeans(mod_slope_from_peak, ~ mate.type)

emm_from_peak_df <- as.data.frame(emm_from_peak) %>%
  mutate(
    decline_est = exp(emmean) - 1,
    decline_lower = exp(lower.CL) - 1,
    decline_upper = exp(upper.CL) - 1
  )

ggplot(from_peak_summary2, aes(x = mate.type, y = decline_rate, color = mate.type)) +
  geom_jitter(width = 0.12, alpha = 0.7, size = 2) +
  geom_pointrange(
    data = emm_from_peak_df,
    aes(
      x = mate.type,
      y = decline_est,
      ymin = decline_lower,
      ymax = decline_upper
    ),
    inherit.aes = FALSE,
    color = "black",
    linewidth = 0.8
  ) +
  scale_x_discrete(labels = c("F" = "Fighter", "S" = "Scrambler")) +
  theme_classic() +
  labs(
    x = "Mate type",
    y = "Post-peak decline rate (eggs/per day)",
  ) +
  theme(legend.position = "none")