Figures and Statistic
Johannes Krietsch
2023-09-13
- Data available relative to clutch initiation
- Mate guarding intensity in relation to daily patterns
- Mate guarding intensity in relation year
- Mate guarding intensity in relation to breeding state
- Mate guarding intensity in relation to breeding state breeders vs. random pairs
- Female moves away
- Number of separating flights
- Distance moved away by sex
- Nest attendance by sex
- Mate guarding intensity and extra-pair paternity
- Split events and extra-pair paternity
- Mate guarding intensity of polyandrous females
- Behavioural observations
#==============================================================================================================Data and code from “Mutual mate guarding and limited sexual conflict
in a sex-role reversed shorebird” Contributor: Johannes Krietsch
📍 This script runs relative to the project’s root directory and
describes how I created all figures (stored in ./OUTPUTS/FIGURES) and
performed all statistical analysis (summary tables stored
in./OUTPUTS/ESM).
#==============================================================================================================
# Packages
sapply( c('data.table', 'magrittr', 'ggplot2', 'viridis', 'auksRuak', 'foreach', 'sf', 'knitr',
'stringr', 'ggnewscale', 'doFuture', 'patchwork', 'activity', 'glmmTMB', 'effects', 'broomExtra',
'flextable', 'officer', 'dplyr', 'performance', 'ggh4x', 'DHARMa'),
require, character.only = TRUE)## data.table magrittr ggplot2 viridis auksRuak foreach sf knitr stringr ggnewscale doFuture patchwork activity glmmTMB
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## effects broomExtra flextable officer dplyr performance ggh4x DHARMa
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE# Lines to run to create html output
opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
# rmarkdown::render('./R/7_figures_and_statistic.R', output_dir = './OUTPUTS/R_COMPILED')
# Data
dID = fread('./DATA/NANO_TAGS_UNIQUE_BY_DAY.txt', sep = '\t', header = TRUE) %>% data.table
dp = fread('./DATA/PAIR_WISE_INTERACTIONS_BREEDING_PAIRS.txt', sep = '\t', header = TRUE) %>% data.table
dr = fread('./DATA/PAIR_WISE_INTERACTIONS_BREEDING_PAIRS_RANDOM.txt', sep = '\t', header = TRUE) %>% data.table
dn = fread('./DATA/NESTS.txt', sep = '\t', header = TRUE, nThread = 20) %>% data.table
st_transform_DT(dn) # change projection
dn[, initiation_y := as.POSIXct(format(initiation, format = '%m-%d %H:%M:%S'),
format = '%m-%d %H:%M:%S', tz = 'UTC')]
dn = dn[year_ > 2017]
# plot settings
margin_ = unit(c(2, 2, 2, 2), 'pt')
margin_top = unit(c(2, 2, 6, 2), 'pt')
sample_size_label = 2.5
egg_laying_color = 'grey85'
male_col = 'steelblue4'
female_col = 'indianred3'# 'firebrick3'
col1 = '#7aa048' # '#7aa048'
col2 = '#E69F00'
# assign periods
dp[, period := fcase(
date_rel_pair >= -5 & date_rel_pair <= -1, "[-5,-1]",
date_rel_pair >= 0 & date_rel_pair <= 3, "[0,3]",
date_rel_pair >= 4 & date_rel_pair <= 10, "[4,10]"
)]
# set order
setorder(dp, pairID, nestID, datetime_1)
# subset data 10 days around clutch initiation
dm = dp[date_rel_pair >= -10 & date_rel_pair <= 10]
dmr = dr[date_rel_pair >= -10 & date_rel_pair <= 10]
# sample size
dm[, .N, pairID] |> nrow()## [1] 64## [1] 68# mean for each day
dpmi = dp[interaction == TRUE, .(N_int = .N), by = .(pairID, nestID, date_rel_pair)]
dpm = unique(dp, by = c('pairID', 'nestID', 'date_rel_pair', 'initiation_rel', 'N', 'year_'))
dpm = merge(dpm, dpmi, by = c('pairID', 'nestID', 'date_rel_pair'), all.x = TRUE)
dpm[is.na(N_int), N_int := 0]
dpm[, interaction_per_day := N_int / N]
dpm[, period := fcase(
date_rel_pair >= -5 & date_rel_pair <= -1, "[-5,-1]",
date_rel_pair >= 0 & date_rel_pair <= 3, "[0,3]",
date_rel_pair >= 4 & date_rel_pair <= 10, "[4,10]"
)]
setorder(dpm, pairID, nestID, date_rel_pair)
# random pairs
dr[, N := .N, by = .(pairID, nestID, date_rel_pair)]
drmi = dr[interaction == TRUE, .(N_int = .N), by = .(pairID, nestID, date_rel_pair)]
drm = unique(dr, by = c('pairID', 'nestID', 'date_rel_pair', 'initiation_rel', 'N', 'year_'))
drm = merge(drm, drmi, by = c('pairID', 'nestID', 'date_rel_pair'), all.x = TRUE)
drm[is.na(N_int), N_int := 0]
drm[, interaction_per_day := N_int / N]
drm[, period := fcase(
date_rel_pair >= -5 & date_rel_pair <= -1, "[-5,-1]",
date_rel_pair >= 0 & date_rel_pair <= 3, "[0,3]",
date_rel_pair >= 4 & date_rel_pair <= 10, "[4,10]"
)]
setorder(drm, pairID, nestID, date_rel_pair)
# nest method
dmn = unique(dm, by = 'nestID')
xn = dmn$nestID
dns = dn[nestID %in% xn]
dns = dns[, .N, by = initiation_method]
dns[, N_percentage_total := round(N / nrow(dmn) * 100, 0)]
# start word file for ESM
ESM = read_docx()
# parameter names
pn = fread("parname; parameter
(Intercept); Intercept
sin(hh2rad(HH)); Sinus time
cos(hh2rad(HH)); Cosinus time
anyEPYTRUE; EPP (yes)
sexM; Sex (male)
initiation_rel; Clutch initiation date
poly(initiation_rel, 2)1; Clutch initiation date (linear)
poly(initiation_rel, 2)2; Clutch initiation date (quadratic)
date_rel_pair; Day relative to clutch initiation
poly(datetime_rel_pair, 2)1; Day relative to clutch initiation (linear)
poly(datetime_rel_pair, 2)2; Day relative to clutch initiation (quadratic)
poly(date_rel_pair, 2)1; Day relative to clutch initiation (linear)
poly(date_rel_pair, 2)2; Day relative to clutch initiation (quadratic)
year_2019; Year (2019)
f_polyandrous_firstTRUE; First clutch of polyandrous female (yes)
typerandomization; Data type (random pairs)
date_rel_pair:typerandomization; Day relative to clutch initiation:data type (random pairs)
poly(datetime_rel_pair0, 2)1:typerandomization; Day relative to clutch initiation (linear):data type (random pairs)
poly(datetime_rel_pair0, 2)2:typerandomization; Day relative to clutch initiation (quadratic):data type (random pairs)
typerandomization:initiation_rel; Clutch initiation date:data type (random pairs)
typerandomization:poly(initiation_rel, 2)1; Clutch initiation date (linear):data type (random pairs)
typerandomization:poly(initiation_rel, 2)2; Clutch initiation date (quadratic):data type (random pairs)
sd__(Intercept); Random intercept
r2marg; R² marginal
r2cond; R² conditional
", sep = ';')Data available relative to clutch initiation
dIDs = unique(dID[date_rel_pair >= -10 & date_rel_pair <= 10],
by = c('sex', 'nestID', 'date_rel_pair'))
dIDs = dIDs[, .N, by = .(sex, date_rel_pair, both_tagged_overlapping)]
dIDs## sex date_rel_pair both_tagged_overlapping N
## 1: F -10 TRUE 7
## 2: F -9 TRUE 11
## 3: F -8 TRUE 15
## 4: F -7 TRUE 22
## 5: F -6 TRUE 27
## 6: M -10 TRUE 5
## 7: M -9 TRUE 8
## 8: M -8 TRUE 13
## 9: M -7 TRUE 21
## 10: M -6 TRUE 28
## 11: M -5 TRUE 36
## 12: M -4 TRUE 38
## 13: M -3 TRUE 46
## 14: M -8 FALSE 10
## 15: M -7 FALSE 11
## 16: M -6 FALSE 14
## 17: M -5 FALSE 17
## 18: M -4 FALSE 18
## 19: F -5 TRUE 34
## 20: F -4 TRUE 34
## 21: F -3 TRUE 42
## 22: F -2 TRUE 43
## 23: F -1 TRUE 48
## 24: F 0 TRUE 52
## 25: F 1 TRUE 59
## 26: F 2 TRUE 58
## 27: M -2 TRUE 46
## 28: M -1 TRUE 51
## 29: M 0 TRUE 55
## 30: M 1 TRUE 60
## 31: M 2 TRUE 61
## 32: M 3 TRUE 59
## 33: F -9 FALSE 1
## 34: F -8 FALSE 1
## 35: F -7 FALSE 1
## 36: F -6 FALSE 2
## 37: F -5 FALSE 3
## 38: F -4 FALSE 4
## 39: F -3 FALSE 5
## 40: F -2 FALSE 5
## 41: F -1 FALSE 6
## 42: F 0 FALSE 8
## 43: F 1 FALSE 11
## 44: F 2 FALSE 11
## 45: F 3 FALSE 11
## 46: F 3 TRUE 52
## 47: F 4 TRUE 48
## 48: M 4 TRUE 55
## 49: M 5 TRUE 57
## 50: M 6 TRUE 55
## 51: M 7 TRUE 51
## 52: M 8 TRUE 50
## 53: M 9 TRUE 47
## 54: M 10 TRUE 45
## 55: F 5 TRUE 43
## 56: M -3 FALSE 18
## 57: F 4 FALSE 9
## 58: F 5 FALSE 7
## 59: F 6 FALSE 6
## 60: F 7 FALSE 5
## 61: F 8 FALSE 5
## 62: F 9 FALSE 5
## 63: F 10 FALSE 5
## 64: F 6 TRUE 34
## 65: F 7 TRUE 31
## 66: F 8 TRUE 30
## 67: F 9 TRUE 27
## 68: F 10 TRUE 23
## 69: M -2 FALSE 18
## 70: M -1 FALSE 17
## 71: M 0 FALSE 19
## 72: M 1 FALSE 19
## 73: M 2 FALSE 19
## 74: M 3 FALSE 19
## 75: M 4 FALSE 18
## 76: M 5 FALSE 17
## 77: M 6 FALSE 17
## 78: M 7 FALSE 16
## 79: M 8 FALSE 16
## 80: M 9 FALSE 15
## 81: M 10 FALSE 15
## 82: M -10 FALSE 6
## 83: M -9 FALSE 6
## sex date_rel_pair both_tagged_overlapping N# pairwise sample size
du = unique(dp, by = c('pairID', 'nestID', 'date_rel_pair'))
dss = unique(du[date_rel_pair >= -10 & date_rel_pair <= 10],
by = c('nestID', 'date_rel_pair'))
dss = dss[, .N, by = date_rel_pair]
dss## date_rel_pair N
## 1: -10 3
## 2: -9 6
## 3: -8 11
## 4: -7 17
## 5: -6 23
## 6: -5 30
## 7: -4 33
## 8: -3 42
## 9: -2 40
## 10: -1 44
## 11: 0 48
## 12: 1 53
## 13: 2 52
## 14: 3 45
## 15: 4 36
## 16: 5 36
## 17: 6 27
## 18: 7 23
## 19: 8 22
## 20: 9 19
## 21: 10 17
## date_rel_pair N## [1] 68## [1] 64## [1] 63## [1] 60# plot data available
dIDs[, sex := factor(sex, levels = c('F', 'M', 'pair'))]
# plot by nest
dp[, datetime_rel_pair_min := min(datetime_rel_pair), by = nestID]
du = unique(dp, by = c('nestID'))
setorder(du, by = datetime_rel_pair_min)
dp[, nestID := factor(nestID, levels = c(du$nestID))]
# last interaction
dp[interaction == TRUE, max(datetime_rel_pair)]## [1] 18.25826ggplot(data = dp) +
geom_rect(aes(xmin = 0, xmax = 3, ymin = -Inf, ymax = Inf), fill = egg_laying_color) +
geom_tile(aes(datetime_rel_pair, nestID, fill = interaction), width = 0.5, show.legend = TRUE) +
scale_fill_manual(values = c('TRUE' = col1, 'FALSE' = col2), name = '',
labels = c('Not together', 'Together'), drop = FALSE) +
geom_vline(aes(xintercept = 0), color = 'black', size = 0.2) +
geom_vline(aes(xintercept = 3), color = 'black', size = 0.2) +
xlab('Day relative to clutch initiation (= 0)') + ylab('Nest') +
scale_x_continuous(limits = c(-12, 19), breaks = seq(-12, 19, 1), expand = expansion(add = c(0.2, 0.2))) +
theme_classic(base_size = 10) +
theme(legend.position = c(0.07, 0.95), legend.background = element_blank(), plot.margin = margin_,
legend.spacing.y = unit(-0.2, "cm"), legend.title = element_blank())## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.## Warning: Removed 2399 rows containing missing values (`geom_tile()`).
Mate guarding intensity in relation to daily patterns
# Function to deal with time
dt2hh <- function(x) {
h <- as.POSIXlt(x)
h$hour + h$min / 60 + h$sec / 3600
}
hh2rad <- function(x) {
x * pi / 12
}
dm[, HH := dt2hh(datetime_1)]
### before clutch initiation
dx = dm[period == "[-5,-1]"]
# sample size
unique(dx$nestID) |> length() # N nests## [1] 51## [1] 49## [1] 48## [1] 18289# model
m <- glmmTMB(interaction ~ sin(hh2rad(HH)) + cos(hh2rad(HH)) +
poly(initiation_rel, 2) + poly(datetime_rel_pair, 2) +
(datetime_rel_pair | nestID),
family = binomial(link = "logit"),
data = dx,
control = glmmTMBControl(parallel = 15)
)
plot(allEffects(m))
## Family: binomial ( logit )
## Formula: interaction ~ sin(hh2rad(HH)) + cos(hh2rad(HH)) + poly(initiation_rel, 2) + poly(datetime_rel_pair, 2) + (datetime_rel_pair | nestID)
## Data: dx
##
## AIC BIC logLik deviance df.resid
## 8668.3 8746.4 -4324.1 8648.3 18279
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## nestID (Intercept) 18.299 4.278
## datetime_rel_pair 2.506 1.583 0.93
## Number of obs: 18289, groups: nestID, 51
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.82472 0.25516 11.071 < 2e-16 ***
## sin(hh2rad(HH)) 0.31843 0.04469 7.126 1.03e-12 ***
## cos(hh2rad(HH)) 0.18539 0.04194 4.421 9.84e-06 ***
## poly(initiation_rel, 2)1 55.91407 30.77934 1.817 0.069277 .
## poly(initiation_rel, 2)2 -100.62627 28.96458 -3.474 0.000513 ***
## poly(datetime_rel_pair, 2)1 130.97631 46.02796 2.846 0.004433 **
## poly(datetime_rel_pair, 2)2 -83.67335 5.59856 -14.946 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.65799, p-value < 2.2e-16
## alternative hypothesis: two.sided
# check effect of time of the day
e = effect("sin(hh2rad(HH))", m, xlevels = 100) |>
data.frame() |>
setDT()
ggplot(e, aes(y = fit, x = HH)) +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.5) +
geom_line(linewidth = 1) +
theme_classic(base_size = 10) +
xlab("Hour") +
ylab("Proportion of time together") 
## [1] 0.9751856## [1] 0.9879668## [1] 1.278118## HH
## 1: 16## HH
## 1: 3.88# create clean summary table
y = tidy(m) |> data.table()
x = r2(m) |> data.table()
setnames(x, c('estimate'))
x[, estimate := as.numeric(estimate)]
x[, term := c('r2cond', 'r2marg')]
y = rbindlist(list(y, x), use.names = TRUE, fill = TRUE)
y[, row_order := rownames(y) |> as.numeric()]
y = merge(y, pn, by.x = 'term', by.y = 'parname')
setorder(y, row_order)
y = y[, .(Parameter = parameter, Estimate = estimate, SE = std.error, Statistic = statistic, p = p.value)]
y = y %>% mutate_if(is.numeric, ~round(., 3)) # round all numeric columns
# save table in word
ft = flextable(y) |> autofit()
ft = bold(ft, bold = TRUE, part = "header")
ESM = ESM |> body_add_par(paste0('Table S1. GLMM together and time -5 to -1')) |> body_add_par('') |>
body_add_flextable(ft)
ESM = ESM |> body_add_break(pos = 'after')
### during egg-laying
dx = dm[period == "[0,3]"]
# sample size
unique(dx$nestID) |> length() # N nests## [1] 56## [1] 53## [1] 50## [1] 19311# model
m <- glmmTMB(interaction ~ sin(hh2rad(HH)) + cos(hh2rad(HH)) +
poly(initiation_rel, 2) + poly(datetime_rel_pair, 2) +
(datetime_rel_pair | nestID),
family = binomial(link = "logit"),
data = dx,
control = glmmTMBControl(parallel = 15)
)
plot(allEffects(m))
## Family: binomial ( logit )
## Formula: interaction ~ sin(hh2rad(HH)) + cos(hh2rad(HH)) + poly(initiation_rel, 2) + poly(datetime_rel_pair, 2) + (datetime_rel_pair | nestID)
## Data: dx
##
## AIC BIC logLik deviance df.resid
## 19936.2 20014.8 -9958.1 19916.2 19301
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## nestID (Intercept) 1.1178 1.057
## datetime_rel_pair 0.5624 0.750 -0.54
## Number of obs: 19311, groups: nestID, 56
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.27736 0.14572 1.903 0.0570 .
## sin(hh2rad(HH)) -0.02875 0.02582 -1.114 0.2654
## cos(hh2rad(HH)) -0.21515 0.02504 -8.591 < 2e-16 ***
## poly(initiation_rel, 2)1 17.23941 18.94586 0.910 0.3629
## poly(initiation_rel, 2)2 -74.80918 17.16283 -4.359 1.31e-05 ***
## poly(datetime_rel_pair, 2)1 -175.55415 16.44745 -10.674 < 2e-16 ***
## poly(datetime_rel_pair, 2)2 -5.89177 3.04347 -1.936 0.0529 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 1.0444, p-value = 0.248
## alternative hypothesis: two.sided
# check effect of time of the day
e = effect("sin(hh2rad(HH))", m, xlevels = 100) |>
data.frame() |>
setDT()
ggplot(e, aes(y = fit, x = HH)) +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.5) +
geom_line(linewidth = 1) +
theme_classic(base_size = 10) +
xlab("Hour") +
ylab("Proportion of time together") 
## [1] 61.0468## [1] 70.7514## [1] 9.7046## HH se
## 1: 0.486 0.04124426## HH se
## 1: 12.6 0.03582247# create clean summary table
y = tidy(m) |> data.table()
x = r2(m) |> data.table()
setnames(x, c('estimate'))
x[, estimate := as.numeric(estimate)]
x[, term := c('r2cond', 'r2marg')]
y = rbindlist(list(y, x), use.names = TRUE, fill = TRUE)
y[, row_order := rownames(y) |> as.numeric()]
y = merge(y, pn, by.x = 'term', by.y = 'parname')
setorder(y, row_order)
y = y[, .(Parameter = parameter, Estimate = estimate, SE = std.error, Statistic = statistic, p = p.value)]
y = y %>% mutate_if(is.numeric, ~round(., 3)) # round all numeric columns
# save table in word
ft = flextable(y) |> autofit()
ft = bold(ft, bold = TRUE, part = "header")
ESM = ESM |> body_add_par(paste0('Table S2. GLMM together and time 0 to 3')) |> body_add_par('') |>
body_add_flextable(ft)
ESM = ESM |> body_add_break(pos = 'after')Mate guarding intensity in relation year
### before clutch initiation
dx = dpm[period == "[-5,-1]"]
# sample size
unique(dx$nestID) |> length() # N nests## [1] 51## [1] 49## [1] 48## [1] 189## year_ N
## 1: 2018 7
## 2: 2019 44## year_ N
## 1: 2018 28
## 2: 2019 161# beta models only accept proportion in the (0,1) interval
dx[interaction_per_day == 1, interaction_per_day := 0.9999]
dx[interaction_per_day == 0, interaction_per_day := 0.0001]
dx[, year_ := as.character(year_)]
# model
m <- glmmTMB(interaction_per_day ~ poly(date_rel_pair, 2) + poly(initiation_rel, 2) + year_ +
(date_rel_pair | nestID),
family = beta_family(link = "logit"), data = dx, REML = TRUE,
control = glmmTMBControl(parallel = 15)
)
plot(allEffects(m))
## Family: beta ( logit )
## Formula: interaction_per_day ~ poly(date_rel_pair, 2) + poly(initiation_rel, 2) + year_ + (date_rel_pair | nestID)
## Data: dx
##
## AIC BIC logLik deviance df.resid
## -768.2 -735.8 394.1 -788.2 185
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## nestID (Intercept) 2.8940 1.7012
## date_rel_pair 0.3978 0.6307 0.94
## Number of obs: 189, groups: nestID, 51
##
## Dispersion parameter for beta family (): 4.61
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.3902 0.4556 5.247 1.55e-07 ***
## poly(date_rel_pair, 2)1 5.1467 2.5120 2.049 0.0405 *
## poly(date_rel_pair, 2)2 -4.9467 1.2217 -4.049 5.14e-05 ***
## poly(initiation_rel, 2)1 1.4986 1.5900 0.943 0.3459
## poly(initiation_rel, 2)2 -6.7147 1.7232 -3.897 9.75e-05 ***
## year_2019 -0.4179 0.4022 -1.039 0.2988
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.88419, p-value = 0.448
## alternative hypothesis: two.sided
## Warning: Model's distribution-specific variance is negative. Results are not reliable.setnames(x, c('estimate'))
x[, estimate := as.numeric(estimate)]
x[, term := c('r2cond', 'r2marg')]
y = rbindlist(list(y, x), use.names = TRUE, fill = TRUE)
y[, row_order := rownames(y) |> as.numeric()]
y = merge(y, pn, by.x = 'term', by.y = 'parname')
setorder(y, row_order)
y = y[, .(Parameter = parameter, Estimate = estimate, SE = std.error, Statistic = statistic, p = p.value)]
y = y %>% mutate_if(is.numeric, ~round(., 3)) # round all numeric columns
# save table in word
ft = flextable(y) |> autofit()
ft = bold(ft, bold = TRUE, part = "header")
ESM = ESM |> body_add_par(paste0('Table S3. GLMM together and year -5 to -1')) |> body_add_par('') |>
body_add_flextable(ft)
ESM = ESM |> body_add_break(pos = 'after')
# check effect of time of the day
e = effect("year_", m, xlevels = 2) |>
data.frame() |>
setDT()
### during egg-laying
dx = dpm[period == "[0,3]"]
# sample size
unique(dx$nestID) |> length() # N nests## [1] 56## [1] 53## [1] 50## [1] 198## year_ N
## 1: 2019 51
## 2: 2018 5## year_ N
## 1: 2019 181
## 2: 2018 17# beta models only accept proportion in the (0,1) interval
dx[interaction_per_day == 1, interaction_per_day := 0.9999]
dx[interaction_per_day == 0, interaction_per_day := 0.0001]
dx[, year_ := as.character(year_)]
# model
m <- glmmTMB(interaction_per_day ~ date_rel_pair + poly(initiation_rel, 2) + year_ +
(date_rel_pair | nestID),
family = beta_family(link = "logit"), data = dx, REML = TRUE,
control = glmmTMBControl(parallel = 15)
)
plot(allEffects(m))
## Family: beta ( logit )
## Formula: interaction_per_day ~ date_rel_pair + poly(initiation_rel, 2) + year_ + (date_rel_pair | nestID)
## Data: dx
##
## AIC BIC logLik deviance df.resid
## -185.1 -155.5 101.6 -203.1 194
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## nestID (Intercept) 0.3000 0.5477
## date_rel_pair 0.1192 0.3453 0.69
## Number of obs: 198, groups: nestID, 56
##
## Dispersion parameter for beta family (): 6.2
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.31478 0.43245 3.040 0.002363 **
## date_rel_pair -0.87742 0.07931 -11.063 < 2e-16 ***
## poly(initiation_rel, 2)1 5.37171 1.91786 2.801 0.005096 **
## poly(initiation_rel, 2)2 -5.71412 1.65056 -3.462 0.000536 ***
## year_2019 0.09730 0.44618 0.218 0.827376
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.83629, p-value = 0.12
## alternative hypothesis: two.sided
## Warning: mu of 2.3 is too close to zero, estimate of random effect variances may be unreliable.
## Warning: Model's distribution-specific variance is negative. Results are not reliable.setnames(x, c('estimate'))
x[, estimate := as.numeric(estimate)]
x[, term := c('r2cond', 'r2marg')]
y = rbindlist(list(y, x), use.names = TRUE, fill = TRUE)
y[, row_order := rownames(y) |> as.numeric()]
y = merge(y, pn, by.x = 'term', by.y = 'parname')
setorder(y, row_order)
y = y[, .(Parameter = parameter, Estimate = estimate, SE = std.error, Statistic = statistic, p = p.value)]
y = y %>% mutate_if(is.numeric, ~round(., 3)) # round all numeric columns
# save table in word
ft = flextable(y) |> autofit()
ft = bold(ft, bold = TRUE, part = "header")
ESM = ESM |> body_add_par(paste0('Table S4. GLMM together and year 0 to 3')) |> body_add_par('') |>
body_add_flextable(ft)
ESM = ESM |> body_add_break(pos = 'after')
# Plot time together breeding pairs split by year
# pairwise sample size 2018
dus = unique(dp, by = c('pairID', 'nestID', 'date_rel_pair'))
dss_2018 = unique(dus[year_ == 2018 & date_rel_pair >= -10 & date_rel_pair <= 10],
by = c('nestID', 'date_rel_pair'))
dss_2018 = dss_2018[, .N, by = date_rel_pair]
# pairwise sample size 2019
dus = unique(dp, by = c('pairID', 'nestID', 'date_rel_pair'))
dss_2019 = unique(dus[year_ == 2019 & date_rel_pair >= -10 & date_rel_pair <= 10],
by = c('nestID', 'date_rel_pair'))
dss_2019 = dss_2019[, .N, by = date_rel_pair]
du[, year_ := as.character(year_)]
# merge
dss = merge(dss_2018[, .(N_2018 = N, date_rel_pair)],
dss_2019[, .(N_2019 = N, date_rel_pair)], by = 'date_rel_pair', all = TRUE)
dss[is.na(N_2018), N_2018 := 0]
dss[, N_year_label := paste0(N_2018, '/', N_2019)]
dpm[, year_cha := as.character(year_)]
p1 =
ggplot() +
geom_text(data = dss, aes(date_rel_pair, Inf, label = N_year_label), vjust = 1, size = sample_size_label) +
geom_rect(aes(xmin = -0.5, xmax = 3.5, ymin = -0.01, ymax = 1), fill = egg_laying_color) +
geom_boxplot(data = dpm,
aes(date_rel_pair, interaction_per_day, group = interaction(date_rel_pair, year_cha),
color = year_cha),
lwd = 0.3, outlier.size = 0.7, outlier.alpha = 0,
position = position_dodge2(preserve = "single")) +
geom_point(data = dpm,
aes(date_rel_pair, interaction_per_day, group = interaction(date_rel_pair, year_cha),
color = year_cha), # shape = data_quality
position=position_jitterdodge(), size = 0.2) +
scale_color_manual(values = c(col1, col2), name = '',
labels = c('2018', '2019'), drop = FALSE) +
scale_x_continuous(limits = c(-10.4, 10.4), breaks = seq(-10, 10, 1),
labels = c('-10', '', '-8', '', '-6', '', '-4', '', '-2', '', '0',
'', '2', '', '4', '', '6', '', '8', '', '10'),
expand = expansion(add = c(0.2, 0.2))) +
scale_y_continuous(limits = c(-0.01, 1.01), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0, 0.05))) +
theme_classic(base_size = 10) +
theme(legend.position = c(0.93, 0.85), legend.background = element_blank(), plot.margin = margin_,
legend.spacing.y = unit(-0.2, "cm"), legend.title = element_blank()) +
ylab('Proportion of time together') +
xlab('Day relative to clutch initiation (= 0)')
p1## Warning: Removed 105 rows containing missing values (`stat_boxplot()`).## Warning: Removed 105 rows containing missing values (`geom_point()`).
# clutch initiation distribution
di = dn[!is.na(year_) & data_type == 'study_site', .(initiation_mean = mean(initiation, na.rm = TRUE)),
by = year_]
dn = merge(dn, di, by = 'year_', all.x = TRUE)
dn[, initiation_rel_ := difftime(initiation, initiation_mean, units = 'days') %>% as.numeric %>% round(., 0)]
dns = dn[data_type == 'study_site' & !is.na(initiation_rel_)]
dns[, year_ := as.character(year_)]
# Sample size
dss = dns[, .(median = median(initiation_rel_), q25 = quantile(initiation_rel_, probs = c(0.25)),
q75 = quantile(initiation_rel_, probs = c(0.75)), .N, max = max(initiation_rel_)), by = year_]
dss[, year_ := as.character(year_)]
# x axis ticks
x = data.table(seq1 = seq(-16, 16, 1))
y = data.table(seq1 = seq(-16, 16, 2), seq2 = as.character(seq(-16, 16, 2)))
xy = merge(x, y, by = 'seq1', all.x = TRUE)
xy[is.na(seq2), seq2 := '']
p2 =
ggplot(data = dns) +
geom_violin(aes(initiation_rel_, year_, color = year_, fill = year_), alpha = 0.2, show.legend = FALSE) +
geom_point(data = dss, aes(median, year_, color = year_), size = 1.5) +
geom_linerange(data = dss, aes(y = year_, xmin = q75, xmax = q25, color = year_), size = 0.3) +
geom_text(data = dss, aes(Inf, year_, label = N), hjust = 1, size = sample_size_label) +
scale_color_manual(values = c(col1, col2), name = '',
labels = c('2018', '2019'), drop = FALSE) +
scale_fill_manual(values = c(col1, col2), name = '',
labels = c('2018', '2019'), drop = FALSE) +
scale_x_continuous(limits = c(-16, 16), breaks = seq(-16, 16, 1),
labels = xy$seq2,
expand = expansion(add = c(0.5, 0.5))) +
ylab('Year') + xlab('Clutch initiation date (standardized)') +
theme_classic(base_size = 10) +
theme(legend.position = "none")
p2
# merge plots
p1 + p2 +
plot_layout(nrow = 2, heights = c(3, 1)) +
plot_annotation(tag_levels = 'a')## Warning: Removed 105 rows containing missing values (`stat_boxplot()`).
## Removed 105 rows containing missing values (`geom_point()`).
Mate guarding intensity in relation to breeding state
# plot time together
# merge data
du = rbindlist(list(dpm[, .(pairID, nestID, year_, date_rel_pair, prop = interaction_per_day,
type = 'm_f_together')],
drm[, .(pairID, nestID, year_, date_rel_pair, prop = interaction_per_day,
type = 'm_f_together_randomized')]
))
# pairwise sample size
ds = unique(dpm, by = c('pairID', 'nestID', 'date_rel_pair'))
dss = unique(ds[date_rel_pair >= -10 & date_rel_pair <= 10],
by = c('nestID', 'date_rel_pair'))
dss = dss[, .N, by = date_rel_pair]
dss## date_rel_pair N
## 1: -10 3
## 2: -9 6
## 3: -8 11
## 4: -7 17
## 5: -6 23
## 6: -5 30
## 7: -4 33
## 8: -3 42
## 9: -2 40
## 10: -1 44
## 11: 0 48
## 12: 1 53
## 13: 2 52
## 14: 3 45
## 15: 4 36
## 16: 5 36
## 17: 6 27
## 18: 7 23
## 19: 8 22
## 20: 9 19
## 21: 10 17
## date_rel_pair N# Plot time together breeding pairs vs. randomized pairs
pa =
ggplot() +
geom_text(data = dss, aes(date_rel_pair, Inf, label = N), vjust = 1, size = sample_size_label) +
geom_rect(aes(xmin = -0.5, xmax = 3.5, ymin = -0.01, ymax = 1), fill = egg_laying_color) +
geom_boxplot(data = du,
aes(date_rel_pair, prop, group = interaction(date_rel_pair, type), color = type),
lwd = 0.3, outlier.size = 0.7, outlier.alpha = 0) +
geom_point(data = du,
aes(date_rel_pair, prop, group = interaction(date_rel_pair, type), color = type),
position=position_jitterdodge(), size = 0.2) +
# scale_shape_manual(values=c(20, 19)) +
scale_color_manual(values = c(col1, col2), name = '',
labels = c('Breeding pair', 'Random pair'), drop = FALSE) +
scale_x_continuous(limits = c(-10.4, 10.4), breaks = seq(-10, 10, 1),
labels = c('-10', '', '-8', '', '-6', '', '-4', '', '-2', '', '0',
'', '2', '', '4', '', '6', '', '8', '', '10'),
expand = expansion(add = c(0.2, 0.2))) +
scale_y_continuous(limits = c(-0.01, 1.01), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0, 0.05))) +
theme_classic(base_size = 10) +
theme(legend.position = c(0.9, 0.85), legend.background = element_blank(), plot.margin = margin_top,
legend.spacing.y = unit(-0.2, "cm"), legend.title = element_blank()) +
ylab('Proportion of time together') +
xlab('Day relative to clutch initiation (= 0)')
pa## Warning: Removed 449 rows containing missing values (`stat_boxplot()`).## Warning: Removed 449 rows containing missing values (`geom_point()`).
### statistic
### before clutch initiation
dx = dpm[period == "[-5,-1]"]
# sample size
unique(dx$nestID) |> length() # N nests## [1] 51## [1] 49## [1] 48## [1] 189# beta models only accept proportion in the (0,1) interval
dx[interaction_per_day == 1, interaction_per_day := 0.9999]
dx[interaction_per_day == 0, interaction_per_day := 0.0001]
# model
m <- glmmTMB(interaction_per_day ~ poly(date_rel_pair, 2) + poly(initiation_rel, 2) +
(date_rel_pair | nestID),
family = beta_family(link = "logit"), data = dx, REML = TRUE,
control = glmmTMBControl(parallel = 15)
)
plot(allEffects(m))
## Family: beta ( logit )
## Formula: interaction_per_day ~ poly(date_rel_pair, 2) + poly(initiation_rel, 2) + (date_rel_pair | nestID)
## Data: dx
##
## AIC BIC logLik deviance df.resid
## -769.2 -740.0 393.6 -787.2 185
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## nestID (Intercept) 2.5565 1.599
## date_rel_pair 0.3625 0.602 0.95
## Number of obs: 189, groups: nestID, 51
##
## Dispersion parameter for beta family (): 4.11
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.9965 0.1637 12.197 < 2e-16 ***
## poly(date_rel_pair, 2)1 4.7641 2.1455 2.221 0.0264 *
## poly(date_rel_pair, 2)2 -4.9250 1.1172 -4.408 1.04e-05 ***
## poly(initiation_rel, 2)1 1.5356 1.4499 1.059 0.2896
## poly(initiation_rel, 2)2 -6.9725 1.5702 -4.441 8.97e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.90071, p-value = 0.512
## alternative hypothesis: two.sided
## Warning: Model's distribution-specific variance is negative. Results are not reliable.setnames(x, c('estimate'))
x[, estimate := as.numeric(estimate)]
x[, term := c('r2cond', 'r2marg')]
y = rbindlist(list(y, x), use.names = TRUE, fill = TRUE)
y[, row_order := rownames(y) |> as.numeric()]
y = merge(y, pn, by.x = 'term', by.y = 'parname')
setorder(y, row_order)
y = y[, .(Parameter = parameter, Estimate = estimate, SE = std.error, Statistic = statistic, p = p.value)]
y = y %>% mutate_if(is.numeric, ~round(., 3)) # round all numeric columns
# save table in word
ft = flextable(y) |> autofit()
ft = bold(ft, bold = TRUE, part = "header")
ESM = ESM |> body_add_par(paste0('Table S5. GLMM together -5 to -1')) |> body_add_par('') |>
body_add_flextable(ft)
ESM = ESM |> body_add_break(pos = 'after')
# descriptive part
x = effect("poly(initiation_rel,2)", m, xlevels = 19) |>
data.frame() |>
setDT() |>
print()## initiation_rel fit se lower upper
## 1: -6 0.5455966 0.12734293 0.3035408 0.7678628
## 2: -5 0.6765325 0.08842456 0.4851747 0.8227487
## 3: -4 0.7733287 0.05536427 0.6465862 0.8641669
## 4: -3 0.8390069 0.03412081 0.7599632 0.8955975
## 5: -2 0.8817228 0.02267445 0.8291883 0.9196652
## 6: -1 0.9089738 0.01721518 0.8688323 0.9377114
## 7: 0 0.9260696 0.01473584 0.8912122 0.9503799
## 8: 1 0.9363634 0.01360764 0.9036332 0.9584877
## 9: 2 0.9418126 0.01318128 0.9096683 0.9629839
## 10: 3 0.9434168 0.01327367 0.9107812 0.9645788
## 11: 4 0.9414636 0.01391899 0.9071637 0.9635993
## 12: 5 0.9356022 0.01533224 0.8978948 0.9600046
## 13: 6 0.9247544 0.01801982 0.8805707 0.9534560
## 14: 7 0.9068497 0.02303440 0.8504044 0.9434143
## 15: 8 0.8783705 0.03233304 0.7989899 0.9291816
## 16: 9 0.8337927 0.04900487 0.7140400 0.9097356
## 17: 10 0.7654406 0.07665636 0.5842835 0.8834077
## 18: 11 0.6653149 0.11618049 0.4152441 0.8476735
## 19: 12 0.5313916 0.15848562 0.2441706 0.7992177## initiation_rel fit se lower upper
## 1: -1 0.9089738 0.01721518 0.8688323 0.9377114
## 2: 0 0.9260696 0.01473584 0.8912122 0.9503799
## 3: 1 0.9363634 0.01360764 0.9036332 0.9584877
## 4: 2 0.9418126 0.01318128 0.9096683 0.9629839
## 5: 3 0.9434168 0.01327367 0.9107812 0.9645788
## 6: 4 0.9414636 0.01391899 0.9071637 0.9635993
## 7: 5 0.9356022 0.01533224 0.8978948 0.9600046
## 8: 6 0.9247544 0.01801982 0.8805707 0.9534560
## 9: 7 0.9068497 0.02303440 0.8504044 0.9434143## fit se
## 1: 74.32375 6.55854## fit se
## 1: 73.48621 8.653208## date_rel_pair fit se lower upper
## 1: -5 0.7854713 0.05273188 0.6638354 0.8716075
## 2: -4 0.9020369 0.02040895 0.8537595 0.9355793
## 3: -3 0.9378609 0.01364559 0.9048460 0.9599282
## 4: -2 0.9416035 0.01372853 0.9078566 0.9634879
## 5: -1 0.9182172 0.02282730 0.8604016 0.9533854## fit se
## 1: 89.7038 2.466845## date_rel_pair fit se lower upper
## 1: -4 0.9020369 0.02040895 0.8537595 0.9355793
## 2: -3 0.9378609 0.01364559 0.9048460 0.9599282
## 3: -2 0.9416035 0.01372853 0.9078566 0.9634879
## 4: -1 0.9182172 0.02282730 0.8604016 0.9533854## [1] 34## [1] 40# extract effect from model for plot
e = effect("poly(initiation_rel,2)", m, xlevels = 100) |>
data.frame() |>
setDT()
# data for points
dms = dm[period == "[-5,-1]"]
dms = dms[, N_ini := .N, by = .(pairID, nestID)]
du = unique(dms, by = c('pairID', 'nestID', 'initiation_rel'))
du = du[!is.na(N_ini)]
du[, .(min(N_ini), max(N_ini))] # check min and max## V1 V2
## 1: 27 686## V1 V2
## 1: -6 12dms = dms[interaction == TRUE & period == "[-5,-1]", .(N_int = .N), by = .(pairID, nestID, initiation_rel)]
du = merge(du, dms, by = c('pairID', 'nestID', 'initiation_rel'), all.x = TRUE)
du[is.na(N_int), N_int := 0]
du[, int_prop := N_int / N_ini]
d0 = copy(du)
# point sizes range
du[, .(min(N_ini), max(N_ini))]## V1 V2
## 1: 27 686pb =
ggplot() +
geom_text(aes(-7.8, Inf, label = 'Day -5 to -1'), vjust = 1, hjust = 0, size = 3.3) +
geom_point(data = du, aes(initiation_rel, int_prop, size = N_ini), shape = 1, color = col1) +
geom_line(data = e, aes(y = fit, x = initiation_rel), size = 0.8, color = col1) +
geom_ribbon(data = e, aes(y = fit, x = initiation_rel, ymin = lower, ymax = upper), alpha = 0.2,
fill = col1) +
scale_x_continuous(limits = c(-8, 12), breaks = seq(-8, 12, 1),
labels = c('-8', '', '-6', '', '-4', '', '-2', '', '0',
'', '2', '', '4', '', '6', '', '8', '', '10', '', '12'),
expand = expansion(add = c(0.2, 0.2))) +
scale_y_continuous(limits = c(-0.01, 1.01), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0.05, 0.05))) +
scale_size_area(max_size = 4, breaks=c(100, 300, 500)) +
theme_classic(base_size = 10) +
theme(legend.position = "none", legend.background = element_blank(), plot.margin = margin_,
legend.spacing.y = unit(-0.2, "cm"), legend.title = element_blank()) +
ylab('Proportion of time together') +
xlab('Clutch initiation date (standardized)')
pb
### during egg-laying
dx = dpm[period == "[0,3]"]
# sample size
unique(dx$nestID) |> length() # N nests## [1] 56## [1] 53## [1] 50## [1] 198# beta models only accept proportion in the (0,1) interval
dx[interaction_per_day == 1, interaction_per_day := 0.9999]
dx[interaction_per_day == 0, interaction_per_day := 0.0001]
dx[, year_ := as.character(year_)]
# model
m <- glmmTMB(interaction_per_day ~ date_rel_pair + poly(initiation_rel, 2) +
(date_rel_pair | nestID),
family = beta_family(link = "logit"), data = dx, REML = TRUE,
control = glmmTMBControl(parallel = 15)
)
plot(allEffects(m))
## Family: beta ( logit )
## Formula: interaction_per_day ~ date_rel_pair + poly(initiation_rel, 2) + (date_rel_pair | nestID)
## Data: dx
##
## AIC BIC logLik deviance df.resid
## -186.9 -160.5 101.4 -202.9 194
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## nestID (Intercept) 0.2872 0.5359
## date_rel_pair 0.1185 0.3442 0.71
## Number of obs: 198, groups: nestID, 56
##
## Dispersion parameter for beta family (): 6.18
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.40471 0.13026 10.784 < 2e-16 ***
## date_rel_pair -0.87830 0.07912 -11.100 < 2e-16 ***
## poly(initiation_rel, 2)1 5.34767 1.90590 2.806 0.005018 **
## poly(initiation_rel, 2)2 -5.63191 1.61186 -3.494 0.000476 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.85907, p-value = 0.208
## alternative hypothesis: two.sided
## Warning: mu of 2.3 is too close to zero, estimate of random effect variances may be unreliable.## Warning: Model's distribution-specific variance is negative. Results are not reliable.setnames(x, c('estimate'))
x[, estimate := as.numeric(estimate)]
x[, term := c('r2cond', 'r2marg')]
y = rbindlist(list(y, x), use.names = TRUE, fill = TRUE)
y[, row_order := rownames(y) |> as.numeric()]
y = merge(y, pn, by.x = 'term', by.y = 'parname')
setorder(y, row_order)
y = y[, .(Parameter = parameter, Estimate = estimate, SE = std.error, Statistic = statistic, p = p.value)]
y = y %>% mutate_if(is.numeric, ~round(., 3)) # round all numeric columns
# save table in word
ft = flextable(y) |> autofit()
ft = bold(ft, bold = TRUE, part = "header")
ESM = ESM |> body_add_par(paste0('Table S6. GLMM together and year 0 to 3')) |> body_add_par('') |>
body_add_flextable(ft)
ESM = ESM |> body_add_break(pos = 'after')
# descriptive part
effect("date_rel_pair", m, xlevels = 4) |>
data.frame() |>
setDT() |>
print() * 100## date_rel_pair fit se lower upper
## 1: 0 0.8432725 0.02078156 0.7978097 0.8800499
## 2: 1 0.6909330 0.03293572 0.6225287 0.7518819
## 3: 2 0.4815553 0.04695502 0.3906101 0.5737393
## 4: 3 0.2784609 0.04900092 0.1926146 0.3843540## date_rel_pair fit se lower upper
## 1: 0 84.32725 2.078156 79.78097 88.00499
## 2: 100 69.09330 3.293572 62.25287 75.18819
## 3: 200 48.15553 4.695502 39.06101 57.37393
## 4: 300 27.84609 4.900092 19.26146 38.43540## initiation_rel fit se lower upper
## 1: -8.00 0.1188473 0.05521910 0.04550561 0.2761908
## 2: -7.00 0.1710422 0.06114130 0.08100722 0.3256826
## 3: -6.10 0.2265871 0.06230326 0.12687634 0.3713329
## 4: -5.10 0.2944725 0.05907007 0.19239586 0.4223819
## 5: -4.20 0.3572608 0.05351710 0.25981421 0.4681404
## 6: -3.20 0.4244436 0.04689673 0.33555293 0.5185063
## 7: -2.30 0.4794718 0.04243316 0.39711938 0.5629552
## 8: -1.30 0.5322917 0.04001219 0.45322861 0.6097659
## 9: -0.38 0.5721997 0.03969050 0.49276112 0.6480824
## 10: 0.57 0.6045098 0.04037406 0.52279087 0.6807796
## 11: 1.50 0.6276571 0.04145708 0.54297246 0.7051705
## 12: 2.50 0.6436558 0.04291416 0.55533561 0.7231781
## 13: 3.40 0.6505374 0.04468024 0.55818037 0.7328292
## 14: 4.40 0.6500552 0.04764191 0.55137538 0.7373682
## 15: 5.30 0.6422928 0.05179469 0.53511907 0.7369063
## 16: 6.30 0.6252752 0.05880573 0.50425083 0.7324319
## 17: 7.20 0.6020369 0.06782391 0.46397254 0.7255735
## 18: 8.20 0.5669508 0.08113706 0.40555758 0.7152872
## 19: 9.10 0.5267975 0.09578817 0.34286251 0.7037354
## 20: 10.00 0.4787146 0.11195805 0.27487221 0.6899008
## 21: 11.00 0.4171046 0.12944062 0.20025546 0.6715838
## 22: 12.00 0.3494116 0.14290573 0.13453013 0.6498161
## initiation_rel fit se lower upper# extract effect from model for plot
e = effect("poly(initiation_rel,2)", m, xlevels = 100) |>
data.frame() |>
setDT()
# data for points
dms = dm[period == "[0,3]"]
dms = dms[, N_ini := .N, by = .(pairID, nestID)]
du = unique(dms, by = c('pairID', 'nestID', 'initiation_rel'))
du = du[!is.na(N_ini)]
du[, .(min(N_ini), max(N_ini))] # check min and max## V1 V2
## 1: 14 559## V1 V2
## 1: -8 12dms = dms[interaction == TRUE & period == "[0,3]", .(N_int = .N), by = .(pairID, nestID, initiation_rel)]
du = merge(du, dms, by = c('pairID', 'nestID', 'initiation_rel'), all.x = TRUE)
du[is.na(N_int), N_int := 0]
du[, int_prop := N_int / N_ini]
d0 = copy(du)
# point sizes range
du[, .(min(N_ini), max(N_ini))]## V1 V2
## 1: 14 559pc =
ggplot() +
geom_text(aes(-7.8, Inf, label = 'Day 0 to 3'), vjust = 1, hjust = 0, size = 3.3) +
geom_point(data = du, aes(initiation_rel, int_prop, size = N_ini), shape = 1, color = col1) +
geom_line(data = e, aes(y = fit, x = initiation_rel), size = 0.8, color = col1) +
geom_ribbon(data = e, aes(y = fit, x = initiation_rel, ymin = lower, ymax = upper), alpha = 0.2,
fill = col1) +
scale_x_continuous(limits = c(-8, 12), breaks = seq(-8, 12, 1),
labels = c('-8', '', '-6', '', '-4', '', '-2', '', '0',
'', '2', '', '4', '', '6', '', '8', '', '10', '', '12'),
expand = expansion(add = c(0.4, 0.4))) +
scale_y_continuous(limits = c(-0.01, 1.01), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0.05, 0.05))) +
scale_size_area(max_size = 4, breaks=c(100, 300, 500)) +
theme_classic(base_size = 10) +
theme(legend.position = "none", legend.background = element_blank(), plot.margin = margin_,
legend.spacing.y = unit(-0.2, "cm"), legend.title = element_blank()) +
ylab('') +
xlab('Clutch initiation date (standardized)')
pc
## Warning: Removed 449 rows containing missing values (`stat_boxplot()`).## Warning: Removed 449 rows containing missing values (`geom_point()`).
Mate guarding intensity in relation to breeding state breeders vs. random pairs
# Statistic together breeding pairs vs. random pairs
# assign random pairs nestID as pairID
dmr[, nestID := pairID]
# merge data
dprm = rbindlist(list(dpm[, .(pairID, nestID, ID1, ID2, year_, date_rel_pair, initiation_rel, interaction_per_day,
period, type = 'breeder')],
drm[, .(pairID, nestID, ID1, ID2, year_, date_rel_pair, initiation_rel, interaction_per_day,
period, type = 'randomization')]
))
### before clutch initiation
dx = dprm[period == "[-5,-1]"]
# sample size
unique(dx$pairID) |> length() # N nests## [1] 228## [1] 57## [1] 54## [1] 428## type N
## 1: breeder 51
## 2: randomization 44## type N
## 1: breeder 189
## 2: randomization 239## type N
## 1: breeder 50
## 2: randomization 178## type N
## 1: breeder 189
## 2: randomization 239# beta models only accept proportion in the (0,1) interval
dx[interaction_per_day == 1, interaction_per_day := 0.9999]
dx[interaction_per_day == 0, interaction_per_day := 0.0001]
# model
m <- glmmTMB(interaction_per_day ~ poly(date_rel_pair, 2) * type + poly(initiation_rel, 2) * type +
(date_rel_pair | nestID),
family = beta_family(link = "logit"), data = dx, REML = TRUE,
control = glmmTMBControl(parallel = 15)
)
plot(allEffects(m))
## Family: beta ( logit )
## Formula: interaction_per_day ~ poly(date_rel_pair, 2) * type + poly(initiation_rel, 2) * type + (date_rel_pair | nestID)
## Data: dx
##
## AIC BIC logLik deviance df.resid
## -1345.0 -1288.2 686.5 -1373.0 424
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## nestID (Intercept) 1.0597 1.0294
## date_rel_pair 0.2139 0.4624 0.98
## Number of obs: 428, groups: nestID, 59
##
## Dispersion parameter for beta family (): 5.63
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.0286 0.1148 17.664 < 2e-16 ***
## poly(date_rel_pair, 2)1 7.6921 2.4854 3.095 0.001969 **
## poly(date_rel_pair, 2)2 -7.8189 1.4913 -5.243 1.58e-07 ***
## typerandomization -4.1261 0.1497 -27.570 < 2e-16 ***
## poly(initiation_rel, 2)1 9.2074 2.1132 4.357 1.32e-05 ***
## poly(initiation_rel, 2)2 -14.2893 2.0705 -6.901 5.15e-12 ***
## poly(date_rel_pair, 2)1:typerandomization -4.1259 2.3597 -1.749 0.080376 .
## poly(date_rel_pair, 2)2:typerandomization 6.8911 2.0248 3.403 0.000666 ***
## typerandomization:poly(initiation_rel, 2)1 -7.9776 2.5538 -3.124 0.001785 **
## typerandomization:poly(initiation_rel, 2)2 13.2206 2.4283 5.444 5.20e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.9235, p-value = 0.656
## alternative hypothesis: two.sided
## Warning in fitTMB(TMBStruc): Model convergence problem; singular convergence (7). See vignette('troubleshooting')## Warning: mu of 1.1 is too close to zero, estimate of random effect variances may be unreliable.## Warning: Model's distribution-specific variance is negative. Results are not reliable.setnames(x, c('estimate'))
x[, estimate := as.numeric(estimate)]
x[, term := c('r2cond', 'r2marg')]
y = rbindlist(list(y, x), use.names = TRUE, fill = TRUE)
y[, row_order := rownames(y) |> as.numeric()]
y = merge(y, pn, by.x = 'term', by.y = 'parname')
setorder(y, row_order)
y = y[, .(Parameter = parameter, Estimate = estimate, SE = std.error, Statistic = statistic, p = p.value)]
y = y %>% mutate_if(is.numeric, ~round(., 3)) # round all numeric columns
# save table in word
ft = flextable(y) |> autofit()
ft = bold(ft, bold = TRUE, part = "header")
ESM = ESM |> body_add_par(paste0('Table S7. GLMM together vs. randomized -5 to -1')) |> body_add_par('') |>
body_add_flextable(ft)
ESM = ESM |> body_add_break(pos = 'after')
# descriptive part
effect("type", m, xlevels = 2) |>
data.frame() |>
setDT() |>
print() * 100## NOTE: type is not a high-order term in the model## type fit se lower upper
## 1: breeder 0.9499649 0.007163943 0.93385284 0.9623107
## 2: randomization 0.1182364 0.014231329 0.09299892 0.1491961## Warning in Ops.factor(left, right): '*' ist nicht sinnvoll für Faktoren## type fit se lower upper
## 1: NA 94.99649 0.7163943 93.385284 96.23107
## 2: NA 11.82364 1.4231329 9.299892 14.91961### during egg-laying
dx = dprm[period == "[0,3]"]
# sample size
unique(dx$pairID) |> length() # N nests## [1] 198## [1] 58## [1] 54## [1] 384## type N
## 1: breeder 56
## 2: randomization 45## type N
## 1: breeder 198
## 2: randomization 186## type N
## 1: breeder 53
## 2: randomization 145## type N
## 1: breeder 198
## 2: randomization 186# beta models only accept proportion in the (0,1) interval
dx[interaction_per_day == 1, interaction_per_day := 0.9999]
dx[interaction_per_day == 0, interaction_per_day := 0.0001]
# model
m <- glmmTMB(interaction_per_day ~ date_rel_pair * type + poly(initiation_rel, 2) * type +
(date_rel_pair | nestID),
family = beta_family(link = "logit"), data = dx, REML = TRUE,
control = glmmTMBControl(parallel = 15)
)## Warning in fitTMB(TMBStruc): Model convergence problem; singular convergence (7). See vignette('troubleshooting')
## Family: beta ( logit )
## Formula: interaction_per_day ~ date_rel_pair * type + poly(initiation_rel, 2) * type + (date_rel_pair | nestID)
## Data: dx
##
## AIC BIC logLik deviance df.resid
## -586.6 -539.2 305.3 -610.6 380
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## nestID (Intercept) 0.23271 0.4824
## date_rel_pair 0.03029 0.1740 1.00
## Number of obs: 384, groups: nestID, 61
##
## Dispersion parameter for beta family (): 5.99
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.35189 0.12258 11.029 < 2e-16 ***
## date_rel_pair -0.85909 0.06176 -13.910 < 2e-16 ***
## typerandomization -3.39124 0.17227 -19.686 < 2e-16 ***
## poly(initiation_rel, 2)1 5.18865 2.11600 2.452 0.014202 *
## poly(initiation_rel, 2)2 -7.04362 1.94987 -3.612 0.000303 ***
## date_rel_pair:typerandomization 0.81175 0.08775 9.250 < 2e-16 ***
## typerandomization:poly(initiation_rel, 2)1 0.18971 2.60351 0.073 0.941911
## typerandomization:poly(initiation_rel, 2)2 6.30402 2.64381 2.384 0.017105 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.9083, p-value = 0.448
## alternative hypothesis: two.sided
## Warning: mu of 0.5 is too close to zero, estimate of random effect variances may be unreliable.setnames(x, c('estimate'))
x[, estimate := as.numeric(estimate)]
x[, term := c('r2cond', 'r2marg')]
y = rbindlist(list(y, x), use.names = TRUE, fill = TRUE)
y[, row_order := rownames(y) |> as.numeric()]
y = merge(y, pn, by.x = 'term', by.y = 'parname')
setorder(y, row_order)
y = y[, .(Parameter = parameter, Estimate = estimate, SE = std.error, Statistic = statistic, p = p.value)]
y = y %>% mutate_if(is.numeric, ~round(., 3)) # round all numeric columns
# save table in word
ft = flextable(y) |> autofit()
ft = bold(ft, bold = TRUE, part = "header")
ESM = ESM |> body_add_par(paste0('Table S8. GLMM together vs. randomized 0 to 3')) |> body_add_par('') |>
body_add_flextable(ft)
ESM = ESM |> body_add_break(pos = 'after')
# descriptive part
effect("type", m, xlevels = 2) |>
data.frame() |>
setDT() |>
print() * 100## NOTE: type is not a high-order term in the model## type fit se lower upper
## 1: breeder 0.5791374 0.03224322 0.51477540 0.6409174
## 2: randomization 0.1106463 0.01525525 0.08401653 0.1443862## Warning in Ops.factor(left, right): '*' ist nicht sinnvoll für Faktoren## type fit se lower upper
## 1: NA 57.91374 3.224322 51.477540 64.09174
## 2: NA 11.06463 1.525525 8.401653 14.43862### after egg-laying
dx = dprm[period == "[4,10]"]
# sample size
unique(dx$nestID) |> length() # N nests## [1] 49## [1] 48## [1] 46## [1] 489## type N
## 1: breeder 38
## 2: randomization 42## type N
## 1: breeder 180
## 2: randomization 309## type N
## 1: breeder 38
## 2: randomization 174## type N
## 1: breeder 180
## 2: randomization 309# beta models only accept proportion in the (0,1) interval
dx[interaction_per_day == 1, interaction_per_day := 0.9999]
dx[interaction_per_day == 0, interaction_per_day := 0.0001]
# model
m <- glmmTMB(interaction_per_day ~ date_rel_pair * type + initiation_rel * type +
(date_rel_pair | nestID),
family = beta_family(link = "logit"), data = dx, REML = TRUE,
control = glmmTMBControl(parallel = 15)
)
plot(allEffects(m))
## Family: beta ( logit )
## Formula: interaction_per_day ~ date_rel_pair * type + initiation_rel * type + (date_rel_pair | nestID)
## Data: dx
##
## AIC BIC logLik deviance df.resid
## -1346.4 -1304.5 683.2 -1366.4 485
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## nestID (Intercept) 0.553349 0.7439
## date_rel_pair 0.005791 0.0761 -0.71
## Number of obs: 489, groups: nestID, 49
##
## Dispersion parameter for beta family (): 4.69
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.29720 0.29399 -4.412 1.02e-05 ***
## date_rel_pair -0.14772 0.04431 -3.334 0.000857 ***
## typerandomization -0.56922 0.37637 -1.512 0.130434
## initiation_rel -0.15517 0.03120 -4.973 6.59e-07 ***
## date_rel_pair:typerandomization 0.12122 0.05494 2.206 0.027365 *
## typerandomization:initiation_rel 0.15830 0.03040 5.207 1.91e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 1.0094, p-value = 0.952
## alternative hypothesis: two.sided
## Warning: mu of 0.1 is too close to zero, estimate of random effect variances may be unreliable.setnames(x, c('estimate'))
x[, estimate := as.numeric(estimate)]
x[, term := c('r2cond', 'r2marg')]
y = rbindlist(list(y, x), use.names = TRUE, fill = TRUE)
y[, row_order := rownames(y) |> as.numeric()]
y = merge(y, pn, by.x = 'term', by.y = 'parname')
setorder(y, row_order)
y = y[, .(Parameter = parameter, Estimate = estimate, SE = std.error, Statistic = statistic, p = p.value)]
y = y %>% mutate_if(is.numeric, ~round(., 3)) # round all numeric columns
# save table in word
ft = flextable(y) |> autofit()
ft = bold(ft, bold = TRUE, part = "header")
ESM = ESM |> body_add_par(paste0('Table S9. GLMM together vs. randomized 4 to 10')) |> body_add_par('') |>
body_add_flextable(ft)
ESM = ESM |> body_add_break(pos = 'after')
# descriptive part
effect("type", m, xlevels = 2) |>
data.frame() |>
setDT() |>
print() * 100## NOTE: type is not a high-order term in the model## type fit se lower upper
## 1: breeder 0.1031824 0.01109636 0.08332741 0.1271122
## 2: randomization 0.1142535 0.01188839 0.09289099 0.1397718## Warning in Ops.factor(left, right): '*' ist nicht sinnvoll für Faktoren## type fit se lower upper
## 1: NA 10.31824 1.109636 8.332741 12.71122
## 2: NA 11.42535 1.188839 9.289099 13.97718Female moves away
# Proportion of split events
dms = dm[split == TRUE, .(N_split = .N), by = .(pairID, nestID, date_rel_pair)]
du = unique(dm[split == TRUE], by = c('pairID', 'nestID', 'date_rel_pair'))
du = merge(du, dms, by = c('pairID', 'nestID', 'date_rel_pair'), all.x = TRUE)
du[is.na(N_split), N_split := 0]
du[, split_prop := N_split / N]
d1 = copy(du)
dm = merge(dm, du[, .(pairID, nestID, date_rel_pair, N_split)], by = c('pairID', 'nestID', 'date_rel_pair'),
all.x = TRUE)
# Times females split
dms = dm[split == TRUE & IDsplitting == 'ID2', .(N_f_split = .N), by = .(pairID, nestID, date_rel_pair)]
du = unique(dm[split == TRUE], by = c('pairID', 'nestID', 'date_rel_pair'))
du = merge(du, dms, by = c('pairID', 'nestID', 'date_rel_pair'), all.x = TRUE)
du[is.na(N_f_split), N_f_split := 0]
du[, f_split_prop := N_f_split /N_split]
# for plot later with EPY
dusm = copy(du)
# pairwise sample size
dus = unique(dp[split == TRUE], by = c('pairID', 'nestID', 'date_rel_pair'))
dss = unique(dus[date_rel_pair >= -10 & date_rel_pair <= 10],
by = c('nestID', 'date_rel_pair'))
dss = dss[, .N, by = date_rel_pair]
dss## date_rel_pair N
## 1: -10 3
## 2: -7 13
## 3: -5 23
## 4: -4 23
## 5: 0 47
## 6: 1 51
## 7: 2 48
## 8: 3 39
## 9: -3 24
## 10: -2 23
## 11: -1 33
## 12: 4 27
## 13: 5 23
## 14: 6 18
## 15: 7 16
## 16: 8 13
## 17: 9 13
## 18: 10 10
## 19: -8 6
## 20: -6 20
## 21: -9 3
## date_rel_pair N# plot splits females
pa =
ggplot() +
geom_text(data = dss, aes(date_rel_pair, Inf, label = N), vjust = 1, size = sample_size_label) +
geom_rect(aes(xmin = -0.5, xmax = 3.5, ymin = -0.01, ymax = 1), fill = egg_laying_color) +
geom_boxplot(data = du,
aes(date_rel_pair, f_split_prop, group = interaction(date_rel_pair)), colour = col1,
lwd = 0.4, outlier.size = 0.7, outlier.alpha = 0) +
geom_jitter(data = du, aes(date_rel_pair, f_split_prop), colour = col1, size = 0.5) +
scale_x_continuous(limits = c(-10.4, 10.4), breaks = seq(-10, 10, 1),
labels = c('-10', '', '-8', '', '-6', '', '-4', '', '-2', '', '0',
'', '2', '', '4', '', '6', '', '8', '', '10'),
expand = expansion(add = c(0.2, 0.2))) +
scale_y_continuous(limits = c(-0.01, 1.01), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0, 0.05))) +
theme_classic(base_size = 11) +
theme(legend.position = c(0.9, 0.15), legend.background = element_blank(), plot.margin = margin_top) +
ylab('Proportion female moves away') +
xlab('Day relative to clutch initiation (= 0)')
pa
### Models
# Statistic females moving away
### before clutch initiation
dx = dp[split == TRUE & period == "[-5,-1]"]
dx[, IDsplitting := ifelse(IDsplitting == 'ID1', 0, 1)] # males = 0
# sample size
unique(dx$nestID) |> length() # N nests## [1] 46## [1] 46## [1] 44## [1] 355# model
m <- glmmTMB(IDsplitting ~ date_rel_pair + initiation_rel +
(date_rel_pair | nestID),
family = binomial(link = "logit"), data = dx, REML = TRUE,
control = glmmTMBControl(parallel = 15)
)## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')## Warning in fitTMB(TMBStruc): Model convergence problem; singular convergence (7). See vignette('troubleshooting')
## Family: binomial ( logit )
## Formula: IDsplitting ~ date_rel_pair + initiation_rel + (date_rel_pair | nestID)
## Data: dx
##
## AIC BIC logLik deviance df.resid
## NA NA NA NA 352
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## nestID (Intercept) 0.568258 0.75383
## date_rel_pair 0.004046 0.06361 1.00
## Number of obs: 355, groups: nestID, 46
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.30403 0.28709 -1.059 0.28959
## date_rel_pair -0.08690 0.08475 -1.025 0.30518
## initiation_rel -0.12401 0.04799 -2.584 0.00976 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.97876, p-value = 0.624
## alternative hypothesis: two.sided
# create clean summary table
y = tidy(m) |> data.table()
x = r2(m) |> data.table()
setnames(x, c('estimate'))
x[, estimate := as.numeric(estimate)]
x[, term := c('r2cond', 'r2marg')]
y = rbindlist(list(y, x), use.names = TRUE, fill = TRUE)
y[, row_order := rownames(y) |> as.numeric()]
y = merge(y, pn, by.x = 'term', by.y = 'parname')
setorder(y, row_order)
y = y[, .(Parameter = parameter, Estimate = estimate, SE = std.error, Statistic = statistic, p = p.value)]
y = y %>% mutate_if(is.numeric, ~round(., 3)) # round all numeric columns
# save table in word
ft = flextable(y) |> autofit()
ft = bold(ft, bold = TRUE, part = "header")
ESM = ESM |> body_add_par(paste0('Table S10. GLMM female moving away -5 to -1')) |> body_add_par('') |>
body_add_flextable(ft)
ESM = ESM |> body_add_break(pos = 'after')
# descriptive part
x = effect("date_rel_pair", m, xlevels = 4) |>
data.frame() |>
setDT() |>
print() ## date_rel_pair fit se lower upper
## 1: -5 0.4979504 0.05819991 0.3859224 0.6101846
## 2: -4 0.4762425 0.04429049 0.3909982 0.5628939
## 3: -2 0.4331753 0.04201998 0.3533538 0.5166224
## 4: -1 0.4119740 0.05401498 0.3115623 0.5202898## fit se
## 1: 45.48355 4.963134## initiation_rel fit se lower upper
## 1: -6 0.6639446 0.08036431 0.4937456 0.8000911
## 2: -5 0.6357374 0.07348464 0.4837470 0.7647421
## 3: -4 0.6065656 0.06590450 0.4729339 0.7259494
## 4: -3 0.5766135 0.05799943 0.4608918 0.6844992
## 5: -2 0.5460874 0.05032813 0.4469327 0.6417162
## 6: -1 0.5152109 0.04370025 0.4299492 0.5995964
## 7: 0 0.4842179 0.03917629 0.4083979 0.5607714
## 8: 1 0.4533458 0.03775115 0.3809013 0.5278216
## 9: 2 0.4228283 0.03965974 0.3475814 0.5018352
## 10: 3 0.3928888 0.04412442 0.3105063 0.4818536
## 11: 4 0.3637332 0.04997871 0.2724037 0.4660688
## 12: 5 0.3355454 0.05625287 0.2354623 0.4529649
## 13: 6 0.3084829 0.06229789 0.2010724 0.4415589
## 14: 7 0.2826746 0.06771400 0.1699821 0.4312600
## 15: 8 0.2582191 0.07226966 0.1424915 0.4217138
## 16: 9 0.2351857 0.07584677 0.1186047 0.4127026## fit se
## 1: 59.06932 6.196354## fit se
## 1: 35.37122 5.450715# extract effect from model for plot
e = effect("initiation_rel", m, xlevels = 100) |>
data.frame() |>
setDT()
# data for points
dms = dm[split == TRUE & period == "[-5,-1]"]
dms[, N_splits_season := .N, by = .(pairID, nestID, initiation_rel)]
du = unique(dms[!is.na(N_splits_season)], by = c('pairID', 'nestID', 'initiation_rel'))
du[, .(min(N_splits_season), max(N_splits_season))] # check min and max## V1 V2
## 1: 1 30## V1 V2
## 1: -6 9# Proportion of split events
dms = dm[split == TRUE & IDsplitting == 'ID2' & period == "[-5,-1]",
.(N_split_female = .N), by = .(pairID, nestID, initiation_rel)]
du = merge(du, dms, by = c('pairID', 'nestID', 'initiation_rel'), all.x = TRUE)
du[is.na(N_split_female), N_split_female := 0]
du[, split_prop := N_split_female / N_splits_season]
d1 = copy(du)
# point sizes range
du[, .(min(N_splits_season), max(N_splits_season))]## V1 V2
## 1: 1 30pb =
ggplot() +
geom_text(aes(-7.8, Inf, label = 'Day -5 to -1'), vjust = 1, hjust = 0, size = 3.3) +
geom_hline(yintercept = 0.5, color = 'black', linetype = 'dashed') +
geom_point(data = du, aes(initiation_rel, split_prop, size = N_splits_season), shape = 1,
color = col1) +
geom_line(data = e, aes(y = fit, x = initiation_rel), size = 0.8, color = col1) +
geom_ribbon(data = e, aes(y = fit, x = initiation_rel, ymin = lower, ymax = upper), alpha = 0.2,
fill = col1) +
scale_x_continuous(limits = c(-8, 12), breaks = seq(-8, 12, 1),
labels = c('-8', '', '-6', '', '-4', '', '-2', '', '0',
'', '2', '', '4', '', '6', '', '8', '', '10', '', '12'),
expand = expansion(add = c(0.4, 0.4))) +
scale_y_continuous(limits = c(-0.01, 1.01), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0.05, 0.05))) +
scale_size_area(max_size = 4, breaks=c(10, 20, 30)) +
theme_classic(base_size = 10) +
theme(legend.position = 'none', legend.background = element_blank(), plot.margin = margin_,
legend.spacing.y = unit(-0.2, "cm"), legend.title = element_blank()) +
ylab('Proportion female moves away') +
xlab('Clutch initiation date (standardized)')
pb
### during egg-laying
dx = dp[split == TRUE & period == "[0,3]"]
dx[, IDsplitting := ifelse(IDsplitting == 'ID1', 0, 1)] # males = 0
# sample size
unique(dx$nestID) |> length() # N nests## [1] 56## [1] 53## [1] 50## [1] 906# model
m <- glmmTMB(IDsplitting ~ date_rel_pair + initiation_rel +
(date_rel_pair | nestID),
family = binomial(link = "logit"), data = dx, REML = TRUE,
control = glmmTMBControl(parallel = 15)
)
plot(allEffects(m))
## Family: binomial ( logit )
## Formula: IDsplitting ~ date_rel_pair + initiation_rel + (date_rel_pair | nestID)
## Data: dx
##
## AIC BIC logLik deviance df.resid
## 1252.7 1281.6 -620.4 1240.7 903
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## nestID (Intercept) 0.50412 0.7100
## date_rel_pair 0.02779 0.1667 -0.94
## Number of obs: 906, groups: nestID, 56
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.05260 0.16403 -0.321 0.748
## date_rel_pair 0.08633 0.07868 1.097 0.273
## initiation_rel 0.04534 0.02759 1.643 0.100
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 1.0018, p-value = 0.856
## alternative hypothesis: two.sided
# create clean summary table
y = tidy(m) |> data.table()
x = r2(m) |> data.table()
setnames(x, c('estimate'))
x[, estimate := as.numeric(estimate)]
x[, term := c('r2cond', 'r2marg')]
y = rbindlist(list(y, x), use.names = TRUE, fill = TRUE)
y[, row_order := rownames(y) |> as.numeric()]
y = merge(y, pn, by.x = 'term', by.y = 'parname')
setorder(y, row_order)
y = y[, .(Parameter = parameter, Estimate = estimate, SE = std.error, Statistic = statistic, p = p.value)]
y = y %>% mutate_if(is.numeric, ~round(., 3)) # round all numeric columns
# save table in word
ft = flextable(y) |> autofit()
ft = bold(ft, bold = TRUE, part = "header")
ESM = ESM |> body_add_par(paste0('Table S11. GLMM female moving away 0 to 3')) |> body_add_par('') |>
body_add_flextable(ft)
ESM = ESM |> body_add_break(pos = 'after')
# descriptive part
x = effect("initiation_rel", m, xlevels = 16) |>
data.frame() |>
setDT() |>
print() ## initiation_rel fit se lower upper
## 1: -8.0 0.4277411 0.06061992 0.3150828 0.5484266
## 2: -6.7 0.4422266 0.05311757 0.3420438 0.5473424
## 3: -5.3 0.4579364 0.04519115 0.3715739 0.5469028
## 4: -4.0 0.4726000 0.03823406 0.3988252 0.5475911
## 5: -2.7 0.4873110 0.03208413 0.4249552 0.5500642
## 6: -1.3 0.5031775 0.02714329 0.4501418 0.5561418
## 7: 0.0 0.5179057 0.02508641 0.4687226 0.5667442
## 8: 1.3 0.5326028 0.02607878 0.4813261 0.5831999
## 9: 2.7 0.5483667 0.03012947 0.4889098 0.6064742
## 10: 4.0 0.5629189 0.03567601 0.4922061 0.6311636
## 11: 5.3 0.5773636 0.04214007 0.4933716 0.6571085
## 12: 6.7 0.5927721 0.04957493 0.4932303 0.6852392
## 13: 8.0 0.6069192 0.05661682 0.4923115 0.7108503
## 14: 9.3 0.6208890 0.06361410 0.4909096 0.7355590
## 15: 11.0 0.6388561 0.07251393 0.4886007 0.7660985
## 16: 12.0 0.6492499 0.07755384 0.4870675 0.7829989## fit se
## 1: 46.51654 4.273169## fit se
## 1: 58.47844 4.789844## date_rel_pair fit se lower upper
## 1: 0 0.4890901 0.04086899 0.4099418 0.5687892
## 2: 1 0.5106699 0.02786685 0.4561416 0.5649454
## 3: 2 0.5322100 0.02551729 0.4820410 0.5817361
## 4: 3 0.5536307 0.03577435 0.4829120 0.6222434## fit se
## 1: 52.14002 3.250687# extract effect from model for plot
e = effect("initiation_rel", m, xlevels = 100) |>
data.frame() |>
setDT()
# data for points
dms = dm[split == TRUE & period == "[0,3]"]
dms[, N_splits_season := .N, by = .(pairID, nestID, initiation_rel)]
du = unique(dms[!is.na(N_splits_season)], by = c('pairID', 'nestID', 'initiation_rel'))
du[, .(min(N_splits_season), max(N_splits_season))] # check min and max## V1 V2
## 1: 1 37## V1 V2
## 1: -8 12# Proportion of split events
dms = dm[split == TRUE & IDsplitting == 'ID2' & period == "[0,3]",
.(N_split_female = .N), by = .(pairID, nestID, initiation_rel)]
du = merge(du, dms, by = c('pairID', 'nestID', 'initiation_rel'), all.x = TRUE)
du[is.na(N_split_female), N_split_female := 0]
du[, split_prop := N_split_female / N_splits_season]
d1 = copy(du)
# point sizes range
du[, .(min(N_splits_season), max(N_splits_season))]## V1 V2
## 1: 1 37pc =
ggplot() +
geom_text(aes(-7.8, Inf, label = 'Day 0 to 3'), vjust = 1, hjust = 0, size = 3.3) +
geom_hline(yintercept = 0.5, color = 'black', linetype = 'dashed') +
geom_point(data = du, aes(initiation_rel, split_prop, size = N_splits_season), shape = 1,
color = col1) +
geom_line(data = e, aes(y = fit, x = initiation_rel), size = 0.8, color = col1) +
geom_ribbon(data = e, aes(y = fit, x = initiation_rel, ymin = lower, ymax = upper), alpha = 0.2,
fill = col1) +
scale_x_continuous(limits = c(-8, 12), breaks = seq(-8, 12, 1),
labels = c('-8', '', '-6', '', '-4', '', '-2', '', '0',
'', '2', '', '4', '', '6', '', '8', '', '10', '', '12'),
expand = expansion(add = c(0.4, 0.4))) +
scale_y_continuous(limits = c(-0.01, 1.01), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0.05, 0.05))) +
scale_size_area(max_size = 4, breaks=c(10, 20, 30)) +
theme_classic(base_size = 10) +
theme(legend.position = 'none', legend.background = element_blank(), plot.margin = margin_,
legend.spacing.y = unit(-0.2, "cm"), legend.title = element_blank()) +
ylab('') +
xlab('Clutch initiation date (standardized)')
pc
Number of separating flights
# subset data
dms = unique(dm, by = c('pairID', 'nestID', 'date_rel_pair'))
dms[is.na(N_split), N_split := 0]
# descriptive statistic
dms[period == "[-5,-1]", .(mean = mean(N_split), min = min(N_split), max = max(N_split))]## mean min max
## 1: 1.878307 0 9## mean min max
## 1: 4.575758 0 17# pairwise sample size
dus = unique(dp, by = c('pairID', 'nestID', 'date_rel_pair'))
dss = unique(dus[date_rel_pair >= -10 & date_rel_pair <= 10],
by = c('nestID', 'date_rel_pair'))
dss = dss[, .N, by = date_rel_pair]
dss## date_rel_pair N
## 1: -10 3
## 2: -9 6
## 3: -8 11
## 4: -7 17
## 5: -6 23
## 6: -5 30
## 7: -4 33
## 8: -3 42
## 9: -2 40
## 10: -1 44
## 11: 0 48
## 12: 1 53
## 13: 2 52
## 14: 3 45
## 15: 4 36
## 16: 5 36
## 17: 6 27
## 18: 7 23
## 19: 8 22
## 20: 9 19
## 21: 10 17
## date_rel_pair Nggplot() +
geom_rect(aes(xmin = -0.5, xmax = 3.5, ymin = -0.3, ymax = 16.1), fill = egg_laying_color) +
geom_text(data = dss, aes(date_rel_pair, Inf, label = N), vjust = 1, size = sample_size_label) +
geom_boxplot(data = dms,
aes(date_rel_pair, N_split, group = interaction(date_rel_pair)),
lwd = 0.4, outlier.size = 0.7, outlier.alpha = 0, color = col1) +
geom_point(data = dms,
aes(date_rel_pair, N_split, group = interaction(date_rel_pair)), color = col1,
position=position_jitter(height = 0), size = 0.2) +
scale_x_continuous(limits = c(-10.4, 10.4), breaks = seq(-10, 10, 1),
labels = c('-10', '', '-8', '', '-6', '', '-4', '', '-2', '', '0',
'', '2', '', '4', '', '6', '', '8', '', '10'),
expand = expansion(add = c(0.2, 0.2))) +
scale_y_continuous(limits = c(-0.3, 16.7), breaks = seq(0, 16, 1),
labels = c('0', '', '2', '', '4', '', '6', '', '8', '', '10', '', '12', '', '14', '', '16'),
expand = expansion(add = c(0, 0))) +
theme_classic(base_size = 10) +
theme(legend.position = c(0.87, 0.94), legend.background = element_blank(), plot.margin = margin_) +
ylab('Number of separating flights') +
xlab('Day relative to clutch initiation (= 0)')## Warning: Removed 1 rows containing non-finite values (`stat_boxplot()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
Distance moved away by sex
# moved away
dms = dm[split == TRUE]
# merge male and female data for plot
dms_m = dms[IDsplitting == 'ID1', .(pairID, nestID, ID1, ID2, date_rel_pair, period, initiation_rel, sex = sex1,
split_distance = distance1_before, stay_distance = distance2_before)]
dms_f = dms[IDsplitting == 'ID2', .(pairID, nestID, ID1, ID2, date_rel_pair, period, initiation_rel, sex = sex2,
split_distance = distance2_before, stay_distance = distance1_before)]
dms = rbindlist(list(dms_m, dms_f))
# how often did the other ID stay within 30 m?
ggplot(data = dms) +
geom_histogram(aes(x = stay_distance)) +
theme_classic(base_size = 10)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## [1] 14.83902## [1] 71.94929# how much was the differences when the other ID moved too?
dms[, delta_split := split_distance - stay_distance]
dms[stay_distance > 30, median(delta_split)]## [1] 64.09144## [1] 0.2025201## [1] 4034.475ggplot(data = dms[stay_distance > 30]) +
geom_histogram(aes(x = delta_split)) +
theme_classic(base_size = 10)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# pairwise sample size
du = unique(dm[split == TRUE], by = c('pairID', 'nestID', 'date_rel_pair'))
dss_m = unique(du[IDsplitting == 'ID1'], by = c('nestID', 'date_rel_pair'))
dss_m = dss_m[, .N, by = date_rel_pair]
# pairwise sample size
dss_f = unique(du[IDsplitting == 'ID2'], by = c('nestID', 'date_rel_pair'))
dss_f = dss_f[, .N, by = date_rel_pair]
# merge
dss = merge(dss_m[, .(N_m = N, date_rel_pair)], dss_f[, .(N_f = N, date_rel_pair)], by = 'date_rel_pair',
all.x = TRUE)
dss[, N_label := paste0(N_f, '/', N_m)]
# adjust distance above 1000 m
dms[, split_distance1000 := split_distance]
dms[split_distance > 1000, split_distance1000 := 1000]
ggplot() +
geom_rect(aes(xmin = -0.5, xmax = 3.5, ymin = -0.01, ymax = 1005), fill = egg_laying_color) +
geom_text(data = dss, aes(date_rel_pair, Inf, label = N_label), vjust = 1, size = sample_size_label) +
geom_boxplot(data = dms,
aes(date_rel_pair, split_distance1000, group = interaction(date_rel_pair, sex),
color = sex),
lwd = 0.4, outlier.size = 0.7, outlier.alpha = 0) +
geom_point(data = dms,
aes(date_rel_pair, split_distance1000, group = interaction(date_rel_pair, sex),
color = sex), position=position_jitterdodge(), size = 0.2) +
scale_color_manual(values = c(female_col, male_col), name = '',
labels = c('Female moves away', 'Male moves away'), drop = FALSE) +
scale_x_continuous(limits = c(-10.4, 10.4), breaks = seq(-10, 10, 1),
labels = c('-10', '', '-8', '', '-6', '', '-4', '', '-2', '', '0',
'', '2', '', '4', '', '6', '', '8', '', '10'),
expand = expansion(add = c(0.2, 0.2))) +
scale_y_continuous(expand = expansion(add = c(0, 50))) +
theme_classic(base_size = 10) +
theme(legend.position = c(0.87, 0.92), legend.background = element_blank(), plot.margin = margin_) +
ylab('Distance moved when separated (m)') +
xlab('Day relative to clutch initiation (= 0)')
# ggsave('./OUTPUTS/FIGURES/male_female_split_events_distance_moved.tiff', plot = last_plot(), width = 177,
# height = 89, units = c('mm'), dpi = 'print')
# statistic
### before clutch initiation
dx = dms[period == "[-5,-1]"]
# sample size
unique(dx$nestID) |> length() # N nests## [1] 46## [1] 46## [1] 44## [1] 355## sex N
## 1: M 40
## 2: F 43## sex N
## 1: M 196
## 2: F 159m <- glmmTMB(split_distance ~ sex + date_rel_pair + initiation_rel + (date_rel_pair | nestID),
family = gaussian, data = dx, REML = TRUE,
control = glmmTMBControl(parallel = 15)
)## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
## Family: gaussian ( identity )
## Formula: split_distance ~ sex + date_rel_pair + initiation_rel + (date_rel_pair | nestID)
## Data: dx
##
## AIC BIC logLik deviance df.resid
## NA NA NA NA 351
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## nestID (Intercept) 201.8 14.20
## date_rel_pair 304.9 17.46 -1.00
## Residual 20245.5 142.29
## Number of obs: 355, groups: nestID, 46
##
## Dispersion estimate for gaussian family (sigma^2): 2.02e+04
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 176.404 20.272 8.702 <2e-16 ***
## sexM -29.511 16.525 -1.786 0.0741 .
## date_rel_pair 4.284 6.687 0.641 0.5218
## initiation_rel 4.556 3.440 1.325 0.1853
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.9405, p-value = 0.504
## alternative hypothesis: two.sided
# create clean summary table
y = tidy(m) |> data.table()
x = r2(m) |> data.table()
setnames(x, c('estimate'))
x[, estimate := as.numeric(estimate)]
x[, term := c('r2cond', 'r2marg')]
y = rbindlist(list(y, x), use.names = TRUE, fill = TRUE)
y[, row_order := rownames(y) |> as.numeric()]
y = merge(y, pn, by.x = 'term', by.y = 'parname')
setorder(y, row_order)
y = y[, .(Parameter = parameter, Estimate = estimate, SE = std.error, Statistic = statistic, p = p.value)]
y = y %>% mutate_if(is.numeric, ~round(., 3)) # round all numeric columns
# save table in word
ft = flextable(y) |> autofit()
ft = bold(ft, bold = TRUE, part = "header")
ESM = ESM |> body_add_par(paste0('Table S12. GLMM split distance -5 to -1')) |> body_add_par('') |>
body_add_flextable(ft)
ESM = ESM |> body_add_break(pos = 'after')
# descriptive part
x = effect("sex", m, xlevels = 2) |>
data.frame() |>
setDT() |>
print()## sex fit se lower upper
## 1: F 169.6253 15.77300 138.6038 200.6468
## 2: M 140.1139 15.21246 110.1948 170.0329## date_rel_pair fit se lower upper
## 1: -5.0 143.7870 23.59875 97.37422 190.1997
## 2: -4.5 145.9288 20.78034 105.05916 186.7984
## 3: -4.0 148.0706 18.14108 112.39173 183.7495
## 4: -3.5 150.2124 15.77119 119.19453 181.2304
## 5: -3.0 152.3543 13.81005 125.19343 179.5151
## 6: -2.5 154.4961 12.45229 130.00561 178.9866
## 7: -2.0 156.6379 11.90617 133.22151 180.0543
## 8: -1.5 158.7798 12.28046 134.62722 182.9323
## 9: -1.0 160.9216 13.49880 134.37287 187.4703### during egg-laying
dx = dms[period == "[0,3]"]
# sample size
unique(dx$nestID) |> length() # N nests## [1] 56## [1] 53## [1] 50## [1] 906## sex N
## 1: M 55
## 2: F 52## sex N
## 1: M 438
## 2: F 468m <- glmmTMB(split_distance ~ sex + date_rel_pair + initiation_rel + (date_rel_pair | nestID),
family = gaussian, data = dx, REML = TRUE,
control = glmmTMBControl(parallel = 15)
)## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Family: gaussian ( identity )
## Formula: split_distance ~ sex + date_rel_pair + initiation_rel + (date_rel_pair | nestID)
## Data: dx
##
## AIC BIC logLik deviance df.resid
## NA NA NA NA 902
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## nestID (Intercept) 2.877e+04 1.696e+02
## date_rel_pair 9.320e-07 9.654e-04 0.30
## Residual 4.814e+04 2.194e+02
## Number of obs: 906, groups: nestID, 56
##
## Dispersion estimate for gaussian family (sigma^2): 4.81e+04
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 222.160 28.174 7.885 3.13e-15 ***
## sexM -21.889 15.438 -1.418 0.156
## date_rel_pair -4.178 7.743 -0.540 0.590
## initiation_rel -2.538 6.249 -0.406 0.685
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.79839, p-value = 0.024
## alternative hypothesis: two.sided
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced## Warning: Can't compute random effect variances. Some variance components equal zero. Your model may suffer from singularity (see `?lme4::isSingular` and
## `?performance::check_singularity`).
## Solution: Respecify random structure! You may also decrease the `tolerance` level to enforce the calculation of random effect variances.## Random effect variances not available. Returned R2 does not account for random effects.setnames(x, c('estimate'))
x[, estimate := as.numeric(estimate)]
x[, term := c('r2cond', 'r2marg')]
y = rbindlist(list(y, x), use.names = TRUE, fill = TRUE)
y[, row_order := rownames(y) |> as.numeric()]
y = merge(y, pn, by.x = 'term', by.y = 'parname')
setorder(y, row_order)
y = y[, .(Parameter = parameter, Estimate = estimate, SE = std.error, Statistic = statistic, p = p.value)]
y = y %>% mutate_if(is.numeric, ~round(., 3)) # round all numeric columns
# save table in word
ft = flextable(y) |> autofit()
ft = bold(ft, bold = TRUE, part = "header")
ESM = ESM |> body_add_par(paste0('Table S13. GLMM split distance 0 to 3')) |> body_add_par('') |>
body_add_flextable(ft)
ESM = ESM |> body_add_break(pos = 'after')
# descriptive part
x = effect("sex", m, xlevels = 2) |>
data.frame() |>
setDT() |>
print()## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced## sex fit se lower upper
## 1: F 215.6459 25.75691 165.0955 266.1964
## 2: M 193.7566 25.78813 143.1449 244.3683## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced
## Warning in sqrt(as.numeric(object$diag.cov.random)): NaNs produced## date_rel_pair fit se lower upper
## 1: 0.00 211.0768 26.87443 158.3331 263.8204
## 2: 0.38 209.4892 25.82757 158.8001 260.1783
## 3: 0.75 207.9434 25.09943 158.6833 257.2035
## 4: 1.10 206.4812 24.69842 158.0081 254.9542
## 5: 1.50 204.8100 24.60034 156.5295 253.0906
## 6: 1.90 203.1389 24.89029 154.2893 251.9885
## 7: 2.20 201.8856 25.35508 152.1237 251.6474
## 8: 2.60 200.2144 26.28317 148.6312 251.7977
## 9: 3.00 198.5433 27.53026 144.5125 252.5741Nest attendance by sex
# assign parameters
dm[, m_at_nest := at_nest1 == TRUE | at_nest2 == TRUE & interaction == TRUE]
dm[, f_at_nest := at_nest2 == TRUE | at_nest1 == TRUE & interaction == TRUE]
dm[, both_at_nest := at_nest1 == TRUE & interaction == TRUE | at_nest2 == TRUE & interaction == TRUE]
dm[, m_alone_at_nest := at_nest1 == TRUE & interaction == FALSE]
dm[, f_alone_at_nest := at_nest2 == TRUE & interaction == FALSE]
# First nest location visit
ds1 = dm[, .(first_pairwise_location = min(datetime_rel_pair)), by = nestID]
ds2 = dm[m_at_nest == TRUE | f_at_nest == TRUE, .(first_nest_visit = min(datetime_rel_pair)), by = nestID]
ds = merge(ds1, ds2, by = 'nestID')
# how many days with data before?
ds[, days_data_before_nest_visit := first_pairwise_location - first_nest_visit]
# subset pairs with at least one day of data before
dss = ds[days_data_before_nest_visit < -1]
# summary (relative to clutch initiation)
dss |> nrow()## [1] 28## mean min max
## 1: -2.581659 -6.768299 0.127662# Male and female together
dms = dm[interaction == TRUE, .(N_int = .N), by = .(pairID, nestID, date_rel_pair)]
du = unique(dm, by = c('pairID', 'nestID', 'date_rel_pair'))
du = merge(du, dms, by = c('pairID', 'nestID', 'date_rel_pair'), all.x = TRUE)
du[is.na(N_int), N_int := 0]
du[, int_prop := N_int / N]
d0 = copy(du)
# Proportion of time males at nest
dms = dm[m_at_nest == TRUE, .(N_m_at_nest = .N), by = .(pairID, nestID, date_rel_pair)]
du = unique(dm, by = c('pairID', 'nestID', 'date_rel_pair'))
du = merge(du, dms, by = c('pairID', 'nestID', 'date_rel_pair'), all.x = TRUE)
du[is.na(N_m_at_nest), N_m_at_nest := 0]
du[, m_at_nest_prop := N_m_at_nest / N]
d1 = copy(du)
# Proportion of time females at nest
dms = dm[f_at_nest == TRUE, .(N_f_at_nest = .N), by = .(pairID, nestID, date_rel_pair)]
du = unique(dm, by = c('pairID', 'nestID', 'date_rel_pair'))
du = merge(du, dms, by = c('pairID', 'nestID', 'date_rel_pair'), all.x = TRUE)
du[is.na(N_f_at_nest), N_f_at_nest := 0]
du[, f_at_nest_prop := N_f_at_nest / N]
d2 = copy(du)
# Proportion of time both at nest
dms = dm[both_at_nest == TRUE, .(N_both_at_nest = .N), by = .(pairID, nestID, date_rel_pair)]
du = unique(dm, by = c('pairID', 'nestID', 'date_rel_pair'))
du = merge(du, dms, by = c('pairID', 'nestID', 'date_rel_pair'), all.x = TRUE)
du[is.na(N_both_at_nest), N_both_at_nest := 0]
du[, both_at_nest_prop := N_both_at_nest / N]
d3 = copy(du)
# Proportion of time male alone at nest
dms = dm[m_alone_at_nest == TRUE, .(N_m_alone_at_nest = .N), by = .(pairID, nestID, date_rel_pair)]
du = unique(dm, by = c('pairID', 'nestID', 'date_rel_pair'))
du = merge(du, dms, by = c('pairID', 'nestID', 'date_rel_pair'), all.x = TRUE)
du[is.na(N_m_alone_at_nest), N_m_alone_at_nest := 0]
du[, m_alone_at_nest_prop := N_m_alone_at_nest / N]
d4 = copy(du)
# Proportion of time female alone at nest
dms = dm[f_alone_at_nest == TRUE, .(N_f_alone_at_nest = .N), by = .(pairID, nestID, date_rel_pair)]
du = unique(dm, by = c('pairID', 'nestID', 'date_rel_pair'))
du = merge(du, dms, by = c('pairID', 'nestID', 'date_rel_pair'), all.x = TRUE)
du[is.na(N_f_alone_at_nest), N_f_alone_at_nest := 0]
du[, f_alone_at_nest_prop := N_f_alone_at_nest / N]
d5 = copy(du)
# male alone and not at the nest
d6 = merge(d0[, .(pairID, nestID, ID1, ID2, date_rel_pair, int_prop)],
d4[, .(pairID, nestID, ID1, ID2, date_rel_pair, m_alone_at_nest_prop)],
by = c('pairID', 'nestID', 'ID1', 'ID2','date_rel_pair'), all.x = TRUE)
d6[, m_alone_prop := 1 - c(int_prop + m_alone_at_nest_prop)]
d6[, total := m_alone_prop + int_prop + m_alone_at_nest_prop]
# female alone and not at the nest
d7 = merge(d0[, .(pairID, nestID, ID1, ID2, date_rel_pair, int_prop)],
d5[, .(pairID, nestID, ID1, ID2, date_rel_pair, f_alone_at_nest_prop)],
by = c('pairID', 'nestID', 'ID1', 'ID2', 'date_rel_pair'), all.x = TRUE)
d7[, f_alone_prop := 1 - c(int_prop + f_alone_at_nest_prop)]
d7[, total := f_alone_prop + int_prop + f_alone_at_nest_prop]
# merge data
du = rbindlist(list(d0[, .(pairID, nestID, ID1, ID2, date_rel_pair, prop = int_prop, type = 'm_f_together')],
d1[, .(pairID, nestID, ID1, ID2, date_rel_pair, prop = m_at_nest_prop, type = 'm_at_nest_prop')],
d2[, .(pairID, nestID, ID1, ID2, date_rel_pair, prop = f_at_nest_prop, type = 'f_at_nest_prop')],
d3[, .(pairID, nestID, ID1, ID2, date_rel_pair, prop = both_at_nest_prop, type = 'both_at_nest_prop')],
d4[, .(pairID, nestID, ID1, ID2, date_rel_pair, prop = m_alone_at_nest_prop, type = 'm_alone_at_nest_prop')],
d5[, .(pairID, nestID, ID1, ID2, date_rel_pair, prop = f_alone_at_nest_prop, type = 'f_alone_at_nest_prop')],
d6[, .(pairID, nestID, ID1, ID2, date_rel_pair, prop = m_alone_prop, type = 'm_alone_prop')],
d7[, .(pairID, nestID, ID1, ID2, date_rel_pair, prop = f_alone_prop, type = 'f_alone_prop')]
))
# descriptive statistic
# day before clutch initiation
du[type == 'f_at_nest_prop' & date_rel_pair == -1, quantile(prop, c(0.5, 0.25, 0.75), na.rm = TRUE)] * 100## 50% 25% 75%
## 10.072464 6.700938 18.914625du[type == 'm_at_nest_prop' & date_rel_pair == -1, quantile(prop, c(0.5, 0.25, 0.75), na.rm = TRUE)] * 100## 50% 25% 75%
## 11.084746 6.719141 18.658088# at clutch initiation
du[type == 'f_at_nest_prop' & date_rel_pair == 0, quantile(prop, c(0.5, 0.25, 0.75), na.rm = TRUE)] * 100## 50% 25% 75%
## 15.53886 10.00000 23.25013du[type == 'm_at_nest_prop' & date_rel_pair == 0, quantile(prop, c(0.5, 0.25, 0.75), na.rm = TRUE)] * 100## 50% 25% 75%
## 20.39513 14.38519 26.82687# at clutch last day of egg laying
du[type == 'f_at_nest_prop' & date_rel_pair == 3, quantile(prop, c(0.5, 0.25, 0.75), na.rm = TRUE)] * 100## 50% 25% 75%
## 6.9767442 0.7194245 25.3333333du[type == 'm_at_nest_prop' & date_rel_pair == 3, quantile(prop, c(0.5, 0.25, 0.75), na.rm = TRUE)] * 100## 50% 25% 75%
## 65.74074 47.36842 75.75758## without partner at nest
# clutch initiation
du[type == 'f_alone_at_nest_prop' & date_rel_pair == 0, quantile(prop, c(0.5, 0.25, 0.75), na.rm = TRUE)] * 100## 50% 25% 75%
## 0.000000 0.000000 1.471429du[type == 'm_alone_at_nest_prop' & date_rel_pair == 0, quantile(prop, c(0.5, 0.25, 0.75), na.rm = TRUE)] * 100## 50% 25% 75%
## 4.375000 1.969462 8.352804# at clutch last day of egg laying
du[type == 'f_alone_at_nest_prop' & date_rel_pair == 3, quantile(prop, c(0.5, 0.25, 0.75), na.rm = TRUE)] * 100## 50% 25% 75%
## 0 0 0du[type == 'm_alone_at_nest_prop' & date_rel_pair == 3, quantile(prop, c(0.5, 0.25, 0.75), na.rm = TRUE)] * 100## 50% 25% 75%
## 46.42857 35.86957 69.87952### without partner away from nest
# clutch initiation
du[type == 'f_alone_prop' & date_rel_pair == 0, quantile(prop, c(0.5, 0.25, 0.75), na.rm = TRUE)] * 100## 50% 25% 75%
## 12.038867 8.034188 19.743590du[type == 'm_alone_prop' & date_rel_pair == 0, quantile(prop, c(0.5, 0.25, 0.75), na.rm = TRUE)] * 100## 50% 25% 75%
## 8.874074 3.502747 13.965116# at clutch last day of egg laying
du[type == 'f_alone_prop' & date_rel_pair == 3, quantile(prop, c(0.5, 0.25, 0.75), na.rm = TRUE)] * 100## 50% 25% 75%
## 70.83333 54.34783 93.52518du[type == 'm_alone_prop' & date_rel_pair == 3, quantile(prop, c(0.5, 0.25, 0.75), na.rm = TRUE)] * 100## 50% 25% 75%
## 18.54839 10.94891 27.90698# pairwise sample size
dus = unique(dp, by = c('pairID', 'nestID', 'date_rel_pair'))
dss = unique(dus[datetime_rel_pair >= -10 & datetime_rel_pair <= 10],
by = c('nestID', 'date_rel_pair'))
dss = dss[, .N, by = date_rel_pair]
dss## date_rel_pair N
## 1: -10 2
## 2: -9 6
## 3: -8 11
## 4: -7 17
## 5: -6 23
## 6: -5 30
## 7: -4 33
## 8: -3 42
## 9: -2 40
## 10: -1 44
## 11: 0 48
## 12: 1 53
## 13: 2 52
## 14: 3 45
## 15: 4 36
## 16: 5 36
## 17: 6 27
## 18: 7 23
## 19: 8 22
## 20: 9 19
## 21: 10 14
## date_rel_pair N# order
dus = du[type == 'm_at_nest_prop' | type == 'f_at_nest_prop']
dus[, type := factor(type, levels = c('m_at_nest_prop', 'f_at_nest_prop'))]
# Males and females at the nest
pa =
ggplot() +
geom_text(data = dss, aes(date_rel_pair, Inf, label = N), vjust = 1, size = sample_size_label) +
geom_rect(aes(xmin = -0.5, xmax = 3.5, ymin = -0.01, ymax = 1), fill = egg_laying_color) +
geom_boxplot(data = dus,
aes(date_rel_pair, prop, group = interaction(date_rel_pair, type), color = type),
lwd = 0.3, outlier.size = 0.7, outlier.alpha = 0) +
geom_point(data = dus,
aes(date_rel_pair, prop, group = interaction(date_rel_pair, type), color = type),
position=position_jitterdodge(), size = 0.2) +
scale_color_manual(values = c(male_col, female_col), name = '',
labels = c('Male', 'Female'), drop = FALSE) +
scale_x_continuous(limits = c(-10.4, 10.4), breaks = seq(-10, 10, 1),
labels = c('-10', '', '-8', '', '-6', '', '-4', '', '-2', '', '0',
'', '2', '', '4', '', '6', '', '8', '', '10'),
expand = expansion(add = c(0.2, 0.2))) +
scale_y_continuous(limits = c(-0.01, 1.01), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0, 0.05))) +
theme_classic(base_size = 10) +
theme(legend.position = c(0.07, 0.9), legend.background = element_blank(), plot.margin = margin_top) +
ylab('Proportion of time at nest') +
xlab('Day relative to clutch initiation (= 0)')
pa
# Males alone and alone at nest
pd =
ggplot() +
geom_text(aes(-5.2, Inf, label = 'Male'), vjust = 1, hjust = 0, size = 3.3) +
geom_rect(aes(xmin = -0.5, xmax = 3.5, ymin = -0.01, ymax = 1), fill = egg_laying_color) +
geom_boxplot(data = du[type == 'm_alone_prop' | type == 'm_alone_at_nest_prop'],
aes(date_rel_pair, prop, group = interaction(date_rel_pair, type), color = type),
lwd = 0.3, outlier.size = 0.7, outlier.alpha = 0) +
geom_point(data = du[type == 'm_alone_prop' | type == 'm_alone_at_nest_prop'],
aes(date_rel_pair, prop, group = interaction(date_rel_pair, type), color = type),
position=position_jitterdodge(), size = 0.2) +
scale_color_manual(values = c(male_col, col2), name = '',
labels = c('at nest', 'not at nest'), drop = FALSE) +
scale_x_continuous(limits = c(-5.4, 5.4), breaks = seq(-5, 5, 1),
labels = c('', '-4', '', '-2', '', '0',
'', '2', '', '4', ''),
expand = expansion(add = c(0.2, 0.2))) +
scale_y_continuous(limits = c(-0.01, 1.01), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0, 0.05))) +
theme_classic(base_size = 10) +
theme(legend.position = c(0.18, 0.9), legend.background = element_blank(), plot.margin = margin_top) +
ylab('Proportion of time alone') +
xlab('Day relative to clutch initiation (= 0)')
pd## Warning: Removed 336 rows containing missing values (`stat_boxplot()`).## Warning: Removed 336 rows containing missing values (`geom_point()`).
# Female alone and alone at nest
pe =
ggplot() +
geom_text(aes(-5.2, Inf, label = 'Female'), vjust = 1, hjust = 0, size = 3.3) +
geom_rect(aes(xmin = -0.5, xmax = 3.5, ymin = -0.01, ymax = 1), fill = egg_laying_color) +
geom_boxplot(data = du[type == 'f_alone_prop' | type == 'f_alone_at_nest_prop'],
aes(date_rel_pair, prop, group = interaction(date_rel_pair, type), color = type),
lwd = 0.3, outlier.size = 0.7, outlier.alpha = 0) +
geom_point(data = du[type == 'f_alone_prop' | type == 'f_alone_at_nest_prop'],
aes(date_rel_pair, prop, group = interaction(date_rel_pair, type), color = type),
position=position_jitterdodge(), size = 0.2) +
scale_color_manual(values = c(female_col, col2), name = '',
labels = c('at nest', 'not at nest'), drop = FALSE) +
scale_x_continuous(limits = c(-5.4, 5.4), breaks = seq(-5, 5, 1),
labels = c('', '-4', '', '-2', '', '0',
'', '2', '', '4', ''),
expand = expansion(add = c(0.2, 0.2))) +
scale_y_continuous(limits = c(-0.01, 1.01), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0, 0.05))) +
theme_classic(base_size = 10) +
theme(legend.position = c(0.18, 0.9), legend.background = element_blank(), plot.margin = margin_top) +
ylab('') +
xlab('Day relative to clutch initiation (= 0)')
pe## Warning: Removed 336 rows containing missing values (`stat_boxplot()`).
## Removed 336 rows containing missing values (`geom_point()`).
# statistic
# merge du with clutch initiation date
dni = unique(dp[, .(nestID, initiation_rel)], by = 'nestID')
du = merge(du, dni, by = 'nestID', all.x = TRUE)
### during egg-laying male
dx = du[type == 'm_at_nest_prop' & date_rel_pair >= 0 & date_rel_pair <= 3]
# sample size
unique(dx$nestID) |> length() # N nests## [1] 56## [1] 53## [1] 50## [1] 198# beta models only accept proportion in the (0,1) interval
dx[prop == 1, prop := 0.9999]
dx[prop == 0, prop := 0.0001]
# model
m <- glmmTMB(prop ~ date_rel_pair + initiation_rel +
(date_rel_pair | nestID),
family = beta_family(link = "logit"), data = dx, REML = TRUE,
control = glmmTMBControl(parallel = 15)
)
plot(allEffects(m))
## Family: beta ( logit )
## Formula: prop ~ date_rel_pair + initiation_rel + (date_rel_pair | nestID)
## Data: dx
##
## AIC BIC logLik deviance df.resid
## -187.3 -164.3 100.7 -201.3 194
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## nestID (Intercept) 0.2125 0.4610
## date_rel_pair 0.1805 0.4248 -0.11
## Number of obs: 198, groups: nestID, 56
##
## Dispersion parameter for beta family (): 14.9
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.31409 0.09689 -13.562 < 2e-16 ***
## date_rel_pair 0.55305 0.07133 7.753 8.97e-15 ***
## initiation_rel 0.03929 0.02280 1.723 0.0849 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.64663, p-value < 2.2e-16
## alternative hypothesis: two.sided
## Warning: mu of 0.4 is too close to zero, estimate of random effect variances may be unreliable.setnames(x, c('estimate'))
x[, estimate := as.numeric(estimate)]
x[, term := c('r2cond', 'r2marg')]
y = rbindlist(list(y, x), use.names = TRUE, fill = TRUE)
y[, row_order := rownames(y) |> as.numeric()]
y = merge(y, pn, by.x = 'term', by.y = 'parname')
setorder(y, row_order)
y = y[, .(Parameter = parameter, Estimate = estimate, SE = std.error, Statistic = statistic, p = p.value)]
y = y %>% mutate_if(is.numeric, ~round(., 3)) # round all numeric columns
# save table in word
ft = flextable(y) |> autofit()
ft = bold(ft, bold = TRUE, part = "header")
ESM = ESM |> body_add_par(paste0('Table S14. GLMM male at nest 0 to 3')) |> body_add_par('') |>
body_add_flextable(ft)
ESM = ESM |> body_add_break(pos = 'after')
# descriptive part
x = effect("initiation_rel", m, xlevels = 22) |>
data.frame() |>
setDT() |>
print() ## initiation_rel fit se lower upper
## 1: -8.00 0.3072924 0.04692959 0.2231132 0.4066091
## 2: -7.00 0.3157188 0.04340265 0.2369031 0.4067783
## 3: -6.10 0.3234078 0.04024640 0.2495714 0.4072345
## 4: -5.10 0.3320643 0.03682160 0.2638064 0.4081910
## 5: -4.20 0.3399536 0.03388728 0.2766082 0.4095918
## 6: -3.20 0.3488247 0.03090763 0.2905805 0.4119639
## 7: -2.30 0.3568994 0.02861248 0.3026513 0.4150833
## 8: -1.30 0.3659675 0.02668240 0.3151081 0.4200017
## 9: -0.38 0.3743950 0.02566880 0.3252949 0.4262246
## 10: 0.57 0.3831780 0.02554497 0.3342029 0.4346448
## 11: 1.50 0.3918506 0.02639521 0.3411920 0.4449508
## 12: 2.50 0.4012526 0.02832834 0.3468768 0.4581727
## 13: 3.40 0.4097773 0.03084445 0.3506002 0.4716877
## 14: 4.40 0.4193128 0.03432674 0.3535130 0.4881130
## 15: 5.30 0.4279469 0.03793651 0.3552929 0.5038470
## 16: 6.30 0.4375918 0.04234209 0.3565760 0.5220791
## 17: 7.20 0.4463130 0.04657245 0.3572625 0.5389479
## 18: 8.20 0.4560419 0.05148884 0.3576364 0.5580033
## 19: 9.10 0.4648269 0.05605548 0.3577057 0.5752943
## 20: 10.00 0.4736337 0.06071733 0.3575792 0.5926088
## 21: 11.00 0.4834381 0.06597192 0.3572584 0.6117657
## 22: 12.00 0.4932553 0.07127192 0.3567878 0.6307354
## initiation_rel fit se lower upper## fit se
## 1: 34.05026 3.479543## fit se
## 1: 43.75707 4.444587## date_rel_pair fit se lower upper
## 1: 0 0.2136641 0.01614057 0.1835551 0.2472166
## 2: 1 0.3208376 0.02006158 0.2826203 0.3616178
## 3: 2 0.4509405 0.03312643 0.3868102 0.5167430
## 4: 3 0.5881159 0.04691009 0.4935642 0.6765822## fit se
## 1: 39.33895 2.905967# extract effect from model for plot
e = effect("initiation_rel", m, xlevels = 100) |>
data.frame() |>
setDT()
# data for points
dms = dm[period == "[0,3]"]
dms[, N_nest_season := .N, by = .(pairID, nestID, initiation_rel)]
dus = unique(dms[!is.na(N_nest_season)], by = c('pairID', 'nestID', 'initiation_rel'))
dus[, .(min(N_nest_season), max(N_nest_season))] # check min and max## V1 V2
## 1: 14 559## V1 V2
## 1: -8 12dms = dm[m_at_nest == TRUE & period == "[0,3]",
.(N_m_at_nest = .N), by = .(pairID, nestID, initiation_rel)]
dus = merge(dus, dms, by = c('pairID', 'nestID', 'initiation_rel'), all.x = TRUE)
dus[is.na(N_m_at_nest), N_m_at_nest := 0]
dus[, m_at_nest_prop := N_m_at_nest / N_nest_season]
d1 = copy(dus)
# point sizes range
dus[, .(min(N_nest_season), max(N_nest_season))]## V1 V2
## 1: 14 559pb =
ggplot() +
geom_text(aes(-7.8, Inf, label = 'Day 0 to 3 (Male)'), vjust = 1, hjust = 0, size = 3.3) +
geom_point(data = dus, aes(initiation_rel, m_at_nest_prop, size = N_nest_season), shape = 1,
color = male_col) +
geom_line(data = e, aes(y = fit, x = initiation_rel), size = 0.8, color = male_col) +
geom_ribbon(data = e, aes(y = fit, x = initiation_rel, ymin = lower, ymax = upper), alpha = 0.2,
fill = male_col) +
scale_x_continuous(limits = c(-8, 12), breaks = seq(-8, 12, 1),
labels = c('-8', '', '-6', '', '-4', '', '-2', '', '0',
'', '2', '', '4', '', '6', '', '8', '', '10', '', '12'),
expand = expansion(add = c(0.4, 0.4))) +
scale_y_continuous(limits = c(-0.01, 1.01), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0.05, 0.05))) +
scale_size_area(max_size = 4, breaks=c(100, 300, 500)) +
theme_classic(base_size = 10) +
theme(legend.position = 'none', legend.background = element_blank(), plot.margin = margin_,
legend.spacing.y = unit(-0.2, "cm"), legend.title = element_blank()) +
ylab('Proportion of time at nest') +
xlab('Clutch initiation date (standardized)')
pb
### during egg-laying female
dx = du[type == 'f_at_nest_prop' & date_rel_pair >= 0 & date_rel_pair <= 3]
# sample size
unique(dx$nestID) |> length() # N nests## [1] 56## [1] 53## [1] 50## [1] 198# beta models only accept proportion in the (0,1) interval
dx[prop == 1, prop := 0.9999]
dx[prop == 0, prop := 0.0001]
# model
m <- glmmTMB(prop ~ date_rel_pair + poly(initiation_rel, 2) +
(date_rel_pair | nestID),
family = beta_family(link = "logit"), data = dx, REML = TRUE,
control = glmmTMBControl(parallel = 15)
)
plot(allEffects(m))
## Family: beta ( logit )
## Formula: prop ~ date_rel_pair + poly(initiation_rel, 2) + (date_rel_pair | nestID)
## Data: dx
##
## AIC BIC logLik deviance df.resid
## -500.1 -473.8 258.0 -516.1 194
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## nestID (Intercept) 0.1898 0.4357
## date_rel_pair 0.1220 0.3493 0.25
## Number of obs: 198, groups: nestID, 56
##
## Dispersion parameter for beta family (): 14.3
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.59367 0.10112 -15.761 < 2e-16 ***
## date_rel_pair -0.25078 0.06831 -3.671 0.000241 ***
## poly(initiation_rel, 2)1 4.68538 1.43945 3.255 0.001134 **
## poly(initiation_rel, 2)2 -3.94660 1.31607 -2.999 0.002711 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.69202, p-value = 0.12
## alternative hypothesis: two.sided
## Warning: mu of 0.2 is too close to zero, estimate of random effect variances may be unreliable.setnames(x, c('estimate'))
x[, estimate := as.numeric(estimate)]
x[, term := c('r2cond', 'r2marg')]
y = rbindlist(list(y, x), use.names = TRUE, fill = TRUE)
y[, row_order := rownames(y) |> as.numeric()]
y = merge(y, pn, by.x = 'term', by.y = 'parname')
setorder(y, row_order)
y = y[, .(Parameter = parameter, Estimate = estimate, SE = std.error, Statistic = statistic, p = p.value)]
y = y %>% mutate_if(is.numeric, ~round(., 3)) # round all numeric columns
# save table in word
ft = flextable(y) |> autofit()
ft = bold(ft, bold = TRUE, part = "header")
ESM = ESM |> body_add_par(paste0('Table S15. GLMM female at nest 0 to 3')) |> body_add_par('') |>
body_add_flextable(ft)
ESM = ESM |> body_add_break(pos = 'after')
# descriptive part
x = effect("poly(initiation_rel,2)", m, xlevels = 22) |>
data.frame() |>
setDT() |>
print() ## initiation_rel fit se lower upper
## 1: -8.00 0.02692465 0.01112447 0.01183428 0.06008728
## 2: -7.00 0.03653727 0.01224275 0.01874102 0.07002663
## 3: -6.10 0.04693322 0.01286012 0.02717278 0.07988372
## 4: -5.10 0.06033331 0.01310202 0.03911158 0.09196774
## 5: -4.20 0.07382885 0.01304209 0.05188287 0.10403915
## 6: -3.20 0.08997902 0.01297415 0.06745903 0.11905710
## 7: -2.30 0.10503700 0.01324833 0.08162817 0.13417807
## 8: -1.30 0.12164419 0.01417781 0.09633238 0.15248455
## 9: -0.38 0.13612313 0.01556571 0.10823123 0.16983404
## 10: 0.57 0.14957880 0.01727286 0.11860389 0.18692757
## 11: 1.50 0.16067413 0.01894130 0.12672479 0.20161864
## 12: 2.50 0.16973554 0.02055354 0.13294672 0.21418969
## 13: 3.40 0.17496910 0.02185119 0.13594851 0.22230799
## 14: 4.40 0.17726425 0.02331197 0.13584360 0.22798264
## 15: 5.30 0.17606877 0.02490214 0.13218599 0.23064773
## 16: 6.30 0.17117898 0.02724155 0.12390201 0.23172397
## 17: 7.20 0.16377235 0.02993627 0.11289114 0.23159917
## 18: 8.20 0.15259314 0.03341431 0.09761105 0.23063014
## 19: 9.10 0.14034291 0.03664432 0.08230101 0.22909940
## 20: 10.00 0.12655623 0.03956163 0.06675604 0.22690120
## 21: 11.00 0.11014350 0.04195900 0.05051708 0.22357616
## 22: 12.00 0.09337895 0.04304786 0.03640644 0.21922415
## initiation_rel fit se lower upper## initiation_rel fit se lower upper
## 1: 1.5 0.1606741 0.01894130 0.12672479 0.2016186
## 2: 2.5 0.1697355 0.02055354 0.13294672 0.2141897
## 3: 3.4 0.1749691 0.02185119 0.13594851 0.2223080
## 4: 4.4 0.1772642 0.02331197 0.13584360 0.2279826
## 5: 5.3 0.1760688 0.02490214 0.13218599 0.2306477
## 6: 6.3 0.1711790 0.02724155 0.12390201 0.2317240
## 7: 7.2 0.1637724 0.02993627 0.11289114 0.2315992
## 8: 8.2 0.1525931 0.03341431 0.09761105 0.2306301## fit se
## 1: 8.469194 1.356103## fit se
## 1: 17.39896 2.190557## fit se
## 1: 13.68523 3.597213## date_rel_pair fit se lower upper
## 1: 0 0.1980059 0.01817284 0.16458406 0.2362951
## 2: 1 0.1611645 0.01615073 0.13179168 0.1956087
## 3: 2 0.1300661 0.01781148 0.09878116 0.1693969
## 4: 3 0.1042229 0.01972076 0.07124198 0.1500062## fit se
## 1: 14.83649 1.796395# extract effect from model for plot
e = effect("poly(initiation_rel,2)", m, xlevels = 100) |>
data.frame() |>
setDT()
# data for points
dms = dm[date_rel_pair >= 0 & date_rel_pair <= 3]
dms[, N_nest_season := .N, by = .(pairID, nestID, initiation_rel)]
dus = unique(dms[!is.na(N_nest_season)], by = c('pairID', 'nestID', 'initiation_rel'))
dus[, .(min(N_nest_season), max(N_nest_season))] # check min and max## V1 V2
## 1: 14 559## V1 V2
## 1: -8 12dms = dm[f_at_nest == TRUE & date_rel_pair >= 0 & date_rel_pair <= 3,
.(N_f_at_nest = .N), by = .(pairID, nestID, initiation_rel)]
dus = merge(dus, dms, by = c('pairID', 'nestID', 'initiation_rel'), all.x = TRUE)
dus[is.na(N_f_at_nest), N_f_at_nest := 0]
dus[, f_at_nest_prop := N_f_at_nest / N_nest_season]
d1 = copy(du)
# point sizes range
dus[, .(min(N_nest_season), max(N_nest_season))]## V1 V2
## 1: 14 559pc =
ggplot() +
geom_text(aes(-7.8, Inf, label = 'Day 0 to 3 (Female)'), vjust = 1, hjust = 0, size = 3.3) +
geom_point(data = dus, aes(initiation_rel, f_at_nest_prop, size = N_nest_season), shape = 1,
color = female_col) +
geom_line(data = e, aes(y = fit, x = initiation_rel), size = 0.8, color = female_col) +
geom_ribbon(data = e, aes(y = fit, x = initiation_rel, ymin = lower, ymax = upper), alpha = 0.2,
fill = female_col) +
scale_x_continuous(limits = c(-8, 12), breaks = seq(-8, 12, 1),
labels = c('-8', '', '-6', '', '-4', '', '-2', '', '0',
'', '2', '', '4', '', '6', '', '8', '', '10', '', '12'),
expand = expansion(add = c(0.4, 0.4))) +
scale_y_continuous(limits = c(-0.01, 1.01), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0.05, 0.05))) +
scale_size_area(max_size = 4, breaks=c(100, 300, 500)) +
theme_classic(base_size = 10) +
theme(legend.position = 'none', legend.background = element_blank(), plot.margin = margin_,
legend.spacing.y = unit(-0.2, "cm"), legend.title = element_blank()) +
ylab('') +
xlab('Clutch initiation date (standardized)')
pc
# merge plots
pa + pb + pc + pd + pe +
plot_layout(design = "
11
23
45
") +
# plot_layout(heights = c(1, 4, 4)) +
plot_annotation(tag_levels = 'a')## Warning: Removed 336 rows containing missing values (`stat_boxplot()`).
## Removed 336 rows containing missing values (`geom_point()`).## Warning: Removed 336 rows containing missing values (`stat_boxplot()`).## Warning: Removed 336 rows containing missing values (`geom_point()`).
Mate guarding intensity and extra-pair paternity
# pairwise sample size
du = unique(dp, by = c('pairID', 'nestID', 'date_rel_pair'))
dss = unique(du[anyEPY == FALSE & date_rel_pair >= -10 & date_rel_pair <= 10],
by = c('nestID', 'date_rel_pair'))
dss = dss[, .N, by = date_rel_pair]
# pairwise sample size
dss_epy = unique(du[anyEPY == TRUE & date_rel_pair >= -10 & date_rel_pair <= 10],
by = c('nestID', 'date_rel_pair'))
dss_epy = dss_epy[, .N, by = date_rel_pair]
dss_epy## date_rel_pair N
## 1: -7 1
## 2: -6 1
## 3: -5 2
## 4: -4 3
## 5: -3 4
## 6: -2 4
## 7: -1 5
## 8: 0 5
## 9: 1 6
## 10: 2 6
## 11: 3 6
## 12: 4 4
## 13: 5 4
## 14: 6 2
## 15: 7 2
## 16: 8 2
## 17: 9 2
## 18: 10 2# merge
dss = merge(dss, dss_epy[, .(N_epy = N, date_rel_pair)], by = 'date_rel_pair', all.x = TRUE)
dss[, N_epy_label := paste0(N, '/', N_epy)]
# Proportion of time together breeders
dps = dp[interaction == TRUE, .(N_int = .N), by = .(pairID, nestID, date_rel_pair)]
du = unique(dp, by = c('pairID', 'nestID', 'date_rel_pair'))
du = merge(du, dps, by = c('pairID', 'nestID', 'date_rel_pair'), all.x = TRUE)
du[is.na(N_int), N_int := 0]
du[, int_prop := N_int / N]
# N
du[, .N, by = pairID] |> nrow()## [1] 64## [1] 68du = du[date_rel_pair >= -10 & date_rel_pair <= 10]
# median per day
dmd = du[, .(int_prop_median = median(int_prop)), by = date_rel_pair]
# order
du[, any_EPY_plot := ifelse(anyEPY == TRUE, 'EPY', 'No EPY')]
du[, any_EPY_plot := factor(any_EPY_plot, levels = c('No EPY', 'EPY'))]
### plot proportion of time together
pa =
ggplot() +
geom_text(data = dss, aes(date_rel_pair, Inf, label = N_epy_label), vjust = 1, size = sample_size_label) +
geom_rect(aes(xmin = -0.5, xmax = 3.5, ymin = -0.01, ymax = 1), fill = egg_laying_color) +
geom_boxplot(data = du[!is.na(anyEPY)],
aes(date_rel_pair, int_prop, group = interaction(date_rel_pair, any_EPY_plot),
color = any_EPY_plot),
lwd = 0.3, outlier.size = 0.7, outlier.alpha = 0) +
geom_point(data = du[!is.na(anyEPY)],
aes(date_rel_pair, int_prop, group = interaction(date_rel_pair, any_EPY_plot),
color = any_EPY_plot),
position=position_jitterdodge(), size = 0.2) +
scale_color_manual(values = c(col1, col2), name = '',
labels = c('No EPP', 'EPP'), drop = FALSE) +
scale_x_continuous(limits = c(-5.4, 5.4), breaks = seq(-5, 5, 1),
labels = c('', '-4', '', '-2', '', '0',
'', '2', '', '4', ''),
expand = expansion(add = c(0.2, 0.2))) +
scale_y_continuous(limits = c(-0.01, 1.01), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0, 0.05))) +
theme_classic(base_size = 10) +
theme(legend.position = c(0.18, 0.14), legend.background = element_blank(), plot.margin = margin_) +
ylab('Proportion of time together') +
xlab('Day relative to clutch initiation (= 0)')
pa## Warning: Removed 168 rows containing missing values (`stat_boxplot()`).## Warning: Removed 10 rows containing missing values (`geom_text()`).## Warning: Removed 168 rows containing missing values (`geom_point()`).
# descriptive statistic
# fertile period
du[anyEPY == TRUE & date_rel_pair >= -5 & date_rel_pair <= 2, quantile(int_prop, c(0.5, 0.25, 0.75), na.rm = TRUE)]## 50% 25% 75%
## 0.7826087 0.6575682 0.9580420du[anyEPY == FALSE & date_rel_pair >= -5 & date_rel_pair <= 2, quantile(int_prop, c(0.5, 0.25, 0.75), na.rm = TRUE)]## 50% 25% 75%
## 0.8787879 0.6331382 0.9744107# before clutch initiation period
du[anyEPY == TRUE & date_rel_pair >= -5 & date_rel_pair <= -1, quantile(int_prop, c(0.5, 0.25, 0.75), na.rm = TRUE)]## 50% 25% 75%
## 0.9027559 0.6818182 0.9950000du[anyEPY == FALSE & date_rel_pair >= -5 & date_rel_pair <= -1, quantile(int_prop, c(0.5, 0.25, 0.75), na.rm = TRUE)]## 50% 25% 75%
## 0.9583333 0.8897059 1.0000000# after clutch initiation period
du[anyEPY == TRUE & date_rel_pair >= 0 & date_rel_pair <= 3, quantile(int_prop, c(0.5, 0.25, 0.75), na.rm = TRUE)]## 50% 25% 75%
## 0.6612903 0.4802294 0.8231225du[anyEPY == FALSE & date_rel_pair >= 0 & date_rel_pair <= 3, quantile(int_prop, c(0.5, 0.25, 0.75), na.rm = TRUE)]## 50% 25% 75%
## 0.5977011 0.2888849 0.8176772### before clutch initiation
dx = du[!is.na(anyEPY) & period == "[-5,-1]"]
dx[, prop := int_prop]
# sample size
unique(dx$nestID) |> length() # N nests## [1] 49## [1] 47## [1] 46## [1] 183## anyEPY N
## 1: FALSE 44
## 2: TRUE 5## anyEPY N
## 1: FALSE 165
## 2: TRUE 18# beta models only accept proportion in the (0,1) interval
dx[prop == 1, prop := 0.9999]
dx[prop == 0, prop := 0.0001]
# model
m <- glmmTMB(prop ~ poly(date_rel_pair, 2) + poly(initiation_rel, 2) + anyEPY +
(date_rel_pair | nestID),
family = beta_family(link = "logit"), data = dx, REML = TRUE,
control = glmmTMBControl(parallel = 15)
)
plot(allEffects(m))
## Family: beta ( logit )
## Formula: prop ~ poly(date_rel_pair, 2) + poly(initiation_rel, 2) + anyEPY + (date_rel_pair | nestID)
## Data: dx
##
## AIC BIC logLik deviance df.resid
## -746.9 -714.8 383.5 -766.9 179
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## nestID (Intercept) 2.9049 1.704
## date_rel_pair 0.3981 0.631 0.94
## Number of obs: 183, groups: nestID, 49
##
## Dispersion parameter for beta family (): 4.26
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.0252 0.1793 11.295 < 2e-16 ***
## poly(date_rel_pair, 2)1 4.7759 2.2607 2.113 0.0346 *
## poly(date_rel_pair, 2)2 -4.7494 1.1244 -4.224 2.40e-05 ***
## poly(initiation_rel, 2)1 1.0120 1.5770 0.642 0.5211
## poly(initiation_rel, 2)2 -6.7902 1.6612 -4.088 4.36e-05 ***
## anyEPYTRUE -0.2078 0.4092 -0.508 0.6116
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.88534, p-value = 0.488
## alternative hypothesis: two.sided
## Warning: Model's distribution-specific variance is negative. Results are not reliable.setnames(x, c('estimate'))
x[, estimate := as.numeric(estimate)]
x[, term := c('r2cond', 'r2marg')]
y = rbindlist(list(y, x), use.names = TRUE, fill = TRUE)
y[, row_order := rownames(y) |> as.numeric()]
y = merge(y, pn, by.x = 'term', by.y = 'parname')
setorder(y, row_order)
y = y[, .(Parameter = parameter, Estimate = estimate, SE = std.error, Statistic = statistic, p = p.value)]
y = y %>% mutate_if(is.numeric, ~round(., 3)) # round all numeric columns
# save table in word
ft = flextable(y) |> autofit()
ft = bold(ft, bold = TRUE, part = "header")
ESM = ESM |> body_add_par(paste0('Table S16. GLMM together and EPY -5 to -1')) |> body_add_par('') |>
body_add_flextable(ft)
ESM = ESM |> body_add_break(pos = 'after')
# extract effect from model for plot
e1 = effect("anyEPY", m, xlevels = 2) |>
data.frame() |>
setDT()
### during egg-laying fertile period
dx = du[!is.na(anyEPY) & date_rel_pair >= 0 & date_rel_pair <= 2]
dx[, prop := int_prop]
# sample size
unique(dx$nestID) |> length() # N nests## [1] 54## [1] 51## [1] 48## [1] 147## anyEPY N
## 1: FALSE 48
## 2: TRUE 6## anyEPY N
## 1: FALSE 130
## 2: TRUE 17# beta models only accept proportion in the (0,1) interval
dx[prop == 1, prop := 0.9999]
dx[prop == 0, prop := 0.0001]
# model
m <- glmmTMB(prop ~ date_rel_pair + poly(initiation_rel, 2) + anyEPY +
(date_rel_pair | nestID),
family = beta_family(link = "logit"), data = dx, REML = TRUE,
control = glmmTMBControl(parallel = 15)
)## Warning in fitTMB(TMBStruc): Model convergence problem; singular convergence (7). See vignette('troubleshooting')
## Family: beta ( logit )
## Formula: prop ~ date_rel_pair + poly(initiation_rel, 2) + anyEPY + (date_rel_pair | nestID)
## Data: dx
##
## AIC BIC logLik deviance df.resid
## -140.6 -113.7 79.3 -158.6 143
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## nestID (Intercept) 0.38363 0.6194
## date_rel_pair 0.09925 0.3150 1.00
## Number of obs: 147, groups: nestID, 54
##
## Dispersion parameter for beta family (): 10.1
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.44422 0.14059 10.272 < 2e-16 ***
## date_rel_pair -0.93535 0.08982 -10.414 < 2e-16 ***
## poly(initiation_rel, 2)1 4.09314 1.56979 2.607 0.009122 **
## poly(initiation_rel, 2)2 -5.31044 1.47013 -3.612 0.000304 ***
## anyEPYTRUE 0.50542 0.41272 1.225 0.220720
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.92423, p-value = 0.552
## alternative hypothesis: two.sided
## Warning: mu of 3.7 is too close to zero, estimate of random effect variances may be unreliable.## Warning: Model's distribution-specific variance is negative. Results are not reliable.setnames(x, c('estimate'))
x[, estimate := as.numeric(estimate)]
x[, term := c('r2cond', 'r2marg')]
y = rbindlist(list(y, x), use.names = TRUE, fill = TRUE)
y[, row_order := rownames(y) |> as.numeric()]
y = merge(y, pn, by.x = 'term', by.y = 'parname')
setorder(y, row_order)
y = y[, .(Parameter = parameter, Estimate = estimate, SE = std.error, Statistic = statistic, p = p.value)]
y = y %>% mutate_if(is.numeric, ~round(., 3)) # round all numeric columns
# save table in word
ft = flextable(y) |> autofit()
ft = bold(ft, bold = TRUE, part = "header")
ESM = ESM |> body_add_par(paste0('Table S17. GLMM together and EPY 0 to 2')) |> body_add_par('') |>
body_add_flextable(ft)
ESM = ESM |> body_add_break(pos = 'after')
# extract effect from model for plot
e2 = effect("anyEPY", m, xlevels = 2) |>
data.frame() |>
setDT()
# plot EPY effect
# merge data
e1[, type := '-5 to -1']
e2[, type := '0 to 2']
e = rbind(e1, e2)
# pairwise sample size
du = unique(dp, by = c('pairID', 'nestID', 'date_rel_pair'))
du = du[date_rel_pair >= -5 & date_rel_pair <= 2]
du[, type := ifelse(date_rel_pair >= 0, '0 to 2', '-5 to -1')]
dsss = unique(du[anyEPY == FALSE], by = c('nestID', 'type'))
dsss = dsss[, .N, by = type]
dsss_epy = unique(du[anyEPY == TRUE], by = c('nestID', 'type'))
dsss_epy = dsss_epy[, .N, by = type]
# merge
dsss = merge(dsss, dsss_epy[, .(N_epy = N, type)], by = 'type', all.x = TRUE)
dsss[, N_epy_label := paste0(N, '/', N_epy)]
# p values
dsp = data.table(type = c( '-5 to -1', '0 to 2'),
p_value = c('p = 0.61', 'p = 0.22'))
p1 =
ggplot() +
geom_point(data = e, aes(type, fit, group = interaction(anyEPY, type), color = anyEPY),
position = position_dodge(width = 0.5)) +
geom_linerange(data = e, aes(x = type, ymin = upper, ymax = lower, color = anyEPY), size = 0.3,
position = position_dodge(width = 0.5)) +
scale_color_manual(values = c(col1, col2), name = '',
labels = c('No EPP', 'EPP'), drop = FALSE) +
geom_text(data = dsss, aes(type, Inf, label = N_epy_label), vjust = 1, size = sample_size_label) +
geom_text(data = dsp, aes(type, c(0.8, 0.55), label = p_value), vjust = 1, size = sample_size_label) +
scale_y_continuous(limits = c(-0.01, 1.01), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0, 0.05))) +
theme_classic(base_size = 10) +
theme(legend.position = c(0.3, 0.15), legend.background = element_blank(),
plot.margin = unit(c(2, 2, 0, 2), 'pt'), axis.title.x = element_blank()) +
ylab('Proportion of time together') +
xlab('')
p1
Split events and extra-pair paternity
# pairwise sample size
du = unique(dp[split == TRUE], by = c('pairID', 'nestID', 'date_rel_pair'))
dss = unique(du[anyEPY == FALSE & date_rel_pair >= -10 & date_rel_pair <= 10],
by = c('nestID', 'date_rel_pair'))
dss = dss[, .N, by = date_rel_pair]
# pairwise sample size
dss_epy = unique(du[anyEPY == TRUE & date_rel_pair >= -10 & date_rel_pair <= 10],
by = c('nestID', 'date_rel_pair'))
dss_epy = dss_epy[, .N, by = date_rel_pair]
dss_epy## date_rel_pair N
## 1: -6 1
## 2: -5 1
## 3: -4 3
## 4: -3 4
## 5: 0 5
## 6: 1 6
## 7: 2 6
## 8: 3 4
## 9: 4 2
## 10: 5 3
## 11: 6 2
## 12: 7 2
## 13: 8 2
## 14: 9 2
## 15: 10 2
## 16: -2 2
## 17: -1 3# merge
dss = merge(dss, dss_epy[, .(N_epy = N, date_rel_pair)], by = 'date_rel_pair', all.x = TRUE)
dss[, N_epy_label := paste0(N, '/', N_epy)]
# plot splits and merges females
dus = dusm[!is.na(anyEPY)]
# order
dus[, any_EPY_plot := ifelse(anyEPY == TRUE, 'EPY', 'No EPY')]
dus[, any_EPY_plot := factor(any_EPY_plot, levels = c('No EPY', 'EPY'))]
pb =
ggplot() +
geom_text(data = dss, aes(date_rel_pair, Inf, label = N_epy_label), vjust = 1, size = sample_size_label) +
geom_rect(aes(xmin = -0.5, xmax = 3.5, ymin = -0.01, ymax = 1), fill = egg_laying_color) +
geom_boxplot(data = dus,
aes(date_rel_pair, f_split_prop, group = interaction(date_rel_pair, any_EPY_plot),
color = any_EPY_plot),
lwd = 0.3, outlier.size = 0.7, outlier.alpha = 0) +
geom_point(data = dus,
aes(date_rel_pair, f_split_prop, group = interaction(date_rel_pair, any_EPY_plot),
color = any_EPY_plot),
position=position_jitterdodge(), size = 0.2) +
scale_color_manual(values = c(col1, col2), name = '',
labels = c('No EPP', 'EPP'), drop = FALSE) +
scale_x_continuous(limits = c(-5.4, 5.4), breaks = seq(-10, 10, 1),
labels = c('-10', '', '-8', '', '-6', '', '-4', '', '-2', '', '0',
'', '2', '', '4', '', '6', '', '8', '', '10'),
expand = expansion(add = c(0.2, 0.2))) +
scale_y_continuous(limits = c(-0.01, 1.01), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0, 0.05))) +
theme_classic(base_size = 10) +
theme(legend.position = c(10.86, 10.12), legend.background = element_blank(), plot.margin = margin_) +
ylab('Proportion female moves away') +
xlab('Day relative to clutch initiation (= 0)')
# merge plots
pa + pb +
plot_layout(ncol = 2) +
plot_annotation(tag_levels = 'a')## Warning: Removed 168 rows containing missing values (`stat_boxplot()`).## Warning: Removed 10 rows containing missing values (`geom_text()`).## Warning: Removed 168 rows containing missing values (`geom_point()`).## Warning: Removed 115 rows containing missing values (`stat_boxplot()`).## Warning: Removed 10 rows containing missing values (`geom_text()`).## Warning: Removed 115 rows containing missing values (`geom_point()`).
# ggsave('./OUTPUTS/FIGURES/male_female_together_epy.tiff', plot = last_plot(), width = 177, height = 89,
# units = c('mm'), dpi = 'print')
# descriptive statistic
# movements away fertile period
dus[anyEPY == TRUE & date_rel_pair >= -5 & date_rel_pair <= 2, quantile(f_split_prop, c(0.5, 0.25, 0.75), na.rm = TRUE)]## 50% 25% 75%
## 0.500 0.400 0.875## [1] 0.6111111dus[anyEPY == FALSE & date_rel_pair >= -5 & date_rel_pair <= 2, quantile(f_split_prop, c(0.5, 0.25, 0.75), na.rm = TRUE)]## 50% 25% 75%
## 0.5000000 0.1666667 0.7329545## [1] 0.4764009# movements before clutch initiation
dus[anyEPY == TRUE & date_rel_pair >= -5 & date_rel_pair <= -1, quantile(f_split_prop, c(0.5, 0.25, 0.75), na.rm = TRUE)]## 50% 25% 75%
## 0.6666667 0.5000000 0.8750000## [1] 0.6576923dus[anyEPY == FALSE & date_rel_pair >= -5 & date_rel_pair <= -1, quantile(f_split_prop, c(0.5, 0.25, 0.75), na.rm = TRUE)]## 50% 25% 75%
## 0.5000000 0.0000000 0.6666667## [1] 0.4302763# movements away after clutch initiation
dus[anyEPY == TRUE & date_rel_pair >= 0 & date_rel_pair <= 3, quantile(f_split_prop, c(0.5, 0.25, 0.75), na.rm = TRUE)]## 50% 25% 75%
## 0.5000000 0.3333333 0.7500000## [1] 0.5484127dus[anyEPY == FALSE & date_rel_pair >= -5 & date_rel_pair <= 3, quantile(f_split_prop, c(0.5, 0.25, 0.75), na.rm = TRUE)]## 50% 25% 75%
## 0.5000000 0.1666667 0.7500000## [1] 0.4835055### Models
# Statistic females moving away
### before clutch initiation
dx = dp[!is.na(anyEPY) & split == TRUE & period == "[-5,-1]"]
dx[, IDsplitting := ifelse(IDsplitting == 'ID1', 0, 1)] # males = 0
# sample size
unique(dx$nestID) |> length() # N nests## [1] 44## [1] 44## [1] 42## [1] 346## anyEPY N
## 1: FALSE 39
## 2: TRUE 5## anyEPY N
## 1: FALSE 290
## 2: TRUE 56# model
m <- glmmTMB(IDsplitting ~ date_rel_pair + initiation_rel + anyEPY +
(date_rel_pair | nestID),
family = binomial(link = "logit"), data = dx, REML = TRUE,
control = glmmTMBControl(parallel = 15)
)## Warning in fitTMB(TMBStruc): Model convergence problem; singular convergence (7). See vignette('troubleshooting')
## Family: binomial ( logit )
## Formula: IDsplitting ~ date_rel_pair + initiation_rel + anyEPY + (date_rel_pair | nestID)
## Data: dx
##
## AIC BIC logLik deviance df.resid
## 466.3 493.2 -226.1 452.3 343
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## nestID (Intercept) 0.646337 0.80395
## date_rel_pair 0.008897 0.09432 1.00
## Number of obs: 346, groups: nestID, 44
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.30856 0.30216 -1.021 0.3072
## date_rel_pair -0.05577 0.08811 -0.633 0.5268
## initiation_rel -0.12042 0.05087 -2.367 0.0179 *
## anyEPYTRUE 0.69207 0.45497 1.521 0.1282
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 0.99005, p-value = 0.808
## alternative hypothesis: two.sided
# create clean summary table
y = tidy(m) |> data.table()
x = r2(m) |> data.table()
setnames(x, c('estimate'))
x[, estimate := as.numeric(estimate)]
x[, term := c('r2cond', 'r2marg')]
y = rbindlist(list(y, x), use.names = TRUE, fill = TRUE)
y[, row_order := rownames(y) |> as.numeric()]
y = merge(y, pn, by.x = 'term', by.y = 'parname')
setorder(y, row_order)
y = y[, .(Parameter = parameter, Estimate = estimate, SE = std.error, Statistic = statistic, p = p.value)]
y = y %>% mutate_if(is.numeric, ~round(., 3)) # round all numeric columns
# save table in word
ft = flextable(y) |> autofit()
ft = bold(ft, bold = TRUE, part = "header")
ESM = ESM |> body_add_par(paste0('Table S18. GLMM female moves away and EPY -5 to -1')) |> body_add_par('') |>
body_add_flextable(ft)
ESM = ESM |> body_add_break(pos = 'after')
# extract effect from model for plot
e1 = effect("anyEPY", m, xlevels = 2) |>
data.frame() |>
setDT()
### during egg-laying
dx = dp[!is.na(anyEPY) & split == TRUE & date_rel_pair >= 0 & date_rel_pair <= 2]
dx[, IDsplitting := ifelse(IDsplitting == 'ID1', 0, 1)] # males = 0
# sample size
unique(dx$nestID) |> length() # N nests## [1] 53## [1] 50## [1] 48## [1] 715## anyEPY N
## 1: FALSE 47
## 2: TRUE 6## anyEPY N
## 1: FALSE 624
## 2: TRUE 91# model
m <- glmmTMB(IDsplitting ~ date_rel_pair + initiation_rel + anyEPY +
(date_rel_pair | nestID),
family = binomial(link = "logit"), data = dx, REML = TRUE,
control = glmmTMBControl(parallel = 15)
)
plot(allEffects(m))
## Family: binomial ( logit )
## Formula: IDsplitting ~ date_rel_pair + initiation_rel + anyEPY + (date_rel_pair | nestID)
## Data: dx
##
## AIC BIC logLik deviance df.resid
## 1003.7 1035.7 -494.9 989.7 712
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## nestID (Intercept) 0.56470 0.7515
## date_rel_pair 0.09633 0.3104 -0.89
## Number of obs: 715, groups: nestID, 53
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.070258 0.184321 0.381 0.703
## date_rel_pair 0.002277 0.111585 0.020 0.984
## initiation_rel 0.023183 0.030247 0.766 0.443
## anyEPYTRUE 0.112732 0.313026 0.360 0.719
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated
##
## data: simulationOutput
## dispersion = 1.0036, p-value = 0.664
## alternative hypothesis: two.sided
# create clean summary table
y = tidy(m) |> data.table()
x = r2(m) |> data.table()
setnames(x, c('estimate'))
x[, estimate := as.numeric(estimate)]
x[, term := c('r2cond', 'r2marg')]
y = rbindlist(list(y, x), use.names = TRUE, fill = TRUE)
y[, row_order := rownames(y) |> as.numeric()]
y = merge(y, pn, by.x = 'term', by.y = 'parname')
setorder(y, row_order)
y = y[, .(Parameter = parameter, Estimate = estimate, SE = std.error, Statistic = statistic, p = p.value)]
y = y %>% mutate_if(is.numeric, ~round(., 3)) # round all numeric columns
# save table in word
ft = flextable(y) |> autofit()
ft = bold(ft, bold = TRUE, part = "header")
ESM = ESM |> body_add_par(paste0('Table S19. GLMM female moves away and EPY 0 to 2')) |> body_add_par('') |>
body_add_flextable(ft)
ESM = ESM |> body_add_break(pos = 'after')
# extract effect from model for plot
e2 = effect("anyEPY", m, xlevels = 2) |>
data.frame() |>
setDT()
# plot EPY effect
# merge data
e1[, type := '-5 to -1']
e2[, type := '0 to 2']
e = rbind(e1, e2)
# pairwise sample size
du = unique(dp[split == TRUE], by = c('pairID', 'nestID', 'date_rel_pair'))
du = du[date_rel_pair >= -5 & date_rel_pair <= 2]
du[, type := ifelse(date_rel_pair >= 0, '0 to 2', '-5 to -1')]
dsss = unique(du[anyEPY == FALSE], by = c('nestID', 'type'))
dsss = dsss[, .N, by = type]
dsss_epy = unique(du[anyEPY == TRUE], by = c('nestID', 'type'))
dsss_epy = dsss_epy[, .N, by = type]
# merge
dsss = merge(dsss, dsss_epy[, .(N_epy = N, type)], by = 'type', all.x = TRUE)
dsss[, N_epy_label := paste0(N, '/', N_epy)]
# p values
dsp = data.table(type = c( '-5 to -1', '0 to 2'),
p_value = c('p = 0.13', 'p = 0.72'))
p2 =
ggplot() +
geom_hline(yintercept = 0.5, color = 'black', linetype = 'dashed') +
geom_point(data = e, aes(type, fit, group = interaction(anyEPY, type), color = anyEPY),
position = position_dodge(width = 0.5)) +
geom_linerange(data = e, aes(x = type, ymin = upper, ymax = lower, color = anyEPY), size = 0.3,
position = position_dodge(width = 0.5)) +
scale_color_manual(values = c(col1, col2), name = '',
labels = c('No EPP', 'EPP'), drop = FALSE) +
geom_text(data = dsss, aes(type, Inf, label = N_epy_label), vjust = 1, size = sample_size_label) +
geom_text(data = dsp, aes(type, c(0.3, 0.35), label = p_value), vjust = 1, size = sample_size_label) +
scale_y_continuous(limits = c(-0.01, 1.01), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0, 0.05))) +
theme_classic(base_size = 10) +
theme(legend.position = 'none', legend.background = element_blank(),
plot.margin = unit(c(2, 2, 0, 2), 'pt'), axis.title.x = element_blank()) +
ylab('Proportion female moves away') +
xlab('Day relative to clutch initiation (= 0)')
p2
##
## Attache Paket: 'gridExtra'## Das folgende Objekt ist maskiert 'package:dplyr':
##
## combinep <- p1 + p2 +
plot_layout(ncol = 2) +
plot_annotation(tag_levels = 'a')
gt = patchworkGrob(p)
g = arrangeGrob(gt, bottom = textGrob(' Day relative to clutch initiation (= 0)', gp = gpar(fontsize = 10)))
# ggsave('./OUTPUTS/FIGURES/male_female_together_female_moving_epy.tiff', plot = g, width = 89, height = 89,
# units = c('mm'), dpi = 'print')Mate guarding intensity of polyandrous females
# Proportion of time together breeders
dps = dp[interaction == TRUE, .(N_int = .N), by = .(pairID, nestID, date_)]
du = unique(dp, by = c('pairID', 'nestID', 'date_'))
du = merge(du, dps, by = c('pairID', 'nestID', 'date_'), all.x = TRUE)
du[is.na(N_int), N_int := 0]
du[, int_prop := N_int / N]
# how many nests with polyandrous first clutch
du[f_polyandrous_first == TRUE, .N, by = .(nestID)]## nestID N
## 1: R211_19 8
## 2: R207_19 29
## 3: R405_19 13
## 4: R901_19 2du[, f_polyandrous_first_plot := ifelse(f_polyandrous_first == TRUE, '1st mate (polyandrous)', 'Other')]
# point sizes range
du[, .(min(N), max(N))]## V1 V2
## 1: 1 142# polyandrous females single plots
p1 =
ggplot() +
geom_rect(aes(xmin = as.Date('2019-06-11'), xmax = as.Date('2019-06-14'), ymin = -0.07, ymax = 1.07),
fill = egg_laying_color) +
geom_rect(aes(xmin = as.Date('2019-06-20'), xmax = as.Date('2019-06-23'), ymin = -0.07, ymax = 1.07),
fill = egg_laying_color) +
geom_path(data = du[ID2 == 273145121], aes(date_, int_prop, group = nestID,
color = f_polyandrous_first_plot), linewidth = 1) +
scale_color_manual(values = c(col1, col2), name = '',
labels = c(bquote("1"^st~ mate), bquote("2"^nd~ mate))) +
geom_point(data = du[ID2 == 273145121], aes(date_, int_prop, color = f_polyandrous_first_plot, size = N)) +
scale_size_area(max_size = 4, breaks=c(10, 50, 100)) +
scale_y_continuous(limits = c(-0.07, 1.07), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0, 0))) +
scale_x_date(date_breaks = '2 day', date_labels = '%d', # %b to add month
limits = c(as.Date('2019-06-11'), as.Date('2019-07-01')), guide = "axis_minor" ) +
theme_classic(base_size = 10) +
theme(legend.position = c(0.9, 0.88), legend.background = element_blank(), plot.margin = margin_,
legend.spacing.y = unit(-0.2, "cm"), legend.title = element_blank(), ggh4x.axis.ticks.length.minor = rel(1)) +
guides(size = "none") +
ylab('') +
xlab('')
p1## Warning: Removed 15 rows containing missing values (`geom_path()`).## Warning: Removed 15 rows containing missing values (`geom_point()`).
p2 =
ggplot() +
geom_rect(aes(xmin = as.Date('2019-06-09'), xmax = as.Date('2019-06-12'), ymin = -0.07, ymax = 1.07),
fill = egg_laying_color) +
geom_rect(aes(xmin = as.Date('2019-06-18'), xmax = as.Date('2019-06-21'), ymin = -0.07, ymax = 1.07),
fill = egg_laying_color) +
geom_path(data = du[ID2 == 270170935 & nestID == 'R405_19' | nestID == 'R406_19'],
aes(date_, int_prop, group = nestID, color = f_polyandrous_first_plot), linewidth = 1) +
scale_color_manual(values = c(col1, col2), name = '',
labels = c('1st mate', '2nd mate')) +
geom_point(data = du[ID2 == 270170935 & nestID == 'R405_19' | nestID == 'R406_19'],
aes(date_, int_prop, color = f_polyandrous_first_plot, size = N)) +
scale_size_area(max_size = 4, breaks=c(10, 50, 100)) +
scale_y_continuous(limits = c(-0.07, 1.07), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0, 0))) +
scale_x_date(date_breaks = '2 day', date_labels = '%d',guide = "axis_minor") +
theme_classic(base_size = 10) +
theme(legend.position = 'none', legend.background = element_blank(), plot.margin = margin_) +
ylab('') +
xlab('')
p2
p3 =
ggplot() +
geom_rect(aes(xmin = as.Date('2019-06-11'), xmax = as.Date('2019-06-13'), ymin = -0.07, ymax = 1.07),
fill = egg_laying_color) +
geom_rect(aes(xmin = as.Date('2019-06-15'), xmax = as.Date('2019-06-16'), ymin = -0.07, ymax = 1.07),
fill = egg_laying_color) +
geom_path(data = du[ID2 == 273145036],
aes(date_, int_prop, group = nestID, color = f_polyandrous_first_plot), linewidth = 1) +
scale_color_manual(values = c(col1, col2), name = '',
labels = c('1st mate', '2nd mate')) +
geom_point(data = du[ID2 == 273145036], aes(date_, int_prop, color = f_polyandrous_first_plot, size = N)) +
scale_size_area(max_size = 4, breaks=c(10, 50, 100)) +
scale_y_continuous(limits = c(-0.07, 1.07), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0, 0))) +
scale_x_date(date_breaks = '2 day', date_labels = '%d', guide = "axis_minor") +
theme_classic(base_size = 10) +
theme(legend.position = 'none', legend.background = element_blank(), plot.margin = margin_) +
ylab('') +
xlab('')
p3
p4 =
ggplot() +
geom_rect(aes(xmin = as.Date('2019-06-06'), xmax = as.Date('2019-06-09'), ymin = -0.07, ymax = 1.07),
fill = egg_laying_color) +
geom_rect(aes(xmin = as.Date('2019-06-15'), xmax = as.Date('2019-06-18'), ymin = -0.07, ymax = 1.07),
fill = egg_laying_color) +
geom_path(data = du[ID2 == 273145109],
aes(date_, int_prop, group = nestID, color = f_polyandrous_first_plot), linewidth = 1) +
scale_color_manual(values = c(col1, col2), name = '',
labels = c('1st mate', '2nd mate')) +
geom_point(data = du[ID2 == 273145109], aes(date_, int_prop, color = f_polyandrous_first_plot, size = N)) +
scale_size_area(max_size = 4, breaks=c(10, 50, 100)) +
scale_y_continuous(limits = c(-0.07, 1.07), breaks = seq(0, 1, 0.1),
labels = c('0.0', '', '0.2', '', '0.4', '', '0.6', '', '0.8', '', '1.0'),
expand = expansion(add = c(0, 0))) +
scale_x_date(date_breaks = '2 day', date_labels = '%d', guide = "axis_minor") +
theme_classic(base_size = 10) +
theme(legend.position = 'none', legend.background = element_blank(), plot.margin = margin_) +
ylab('Proportion of time together') +
xlab('Date (June)')
p4
# merge plots
p1 + p2 + p3 + p4 +
plot_layout(nrow = 4, ncol = 1) +
plot_annotation(tag_levels = 'a')## Warning: Removed 15 rows containing missing values (`geom_path()`).
## Removed 15 rows containing missing values (`geom_point()`).
Behavioural observations
# Data
d = fread('./DATA/OBSERVATIONS.txt', sep = '\t', header = TRUE) %>% data.table
### Following flights
# Total observations
unique(d, by = 'obs_id') |> nrow()## [1] 5524# How many ID at observations?
d[, N_ind_int := .N, by = obs_id]
# How many observations with at least two IDs and a flight?
d[, any_flights := any(flight %like% 'F'), by = obs_id]
# Total with at least 2 IDs
unique(d[N_ind_int > 1], by = 'obs_id') |> nrow()## [1] 2876# Total with a flight initiator scored
d[, any_F1 := any(flight == 'F1'), by = obs_id]
d[, any_F0 := any(flight == 'F0'), by = obs_id]
# subset relevant data
ds = d[!is.na(sex) & N_ind_int > 1 & any_F1 == TRUE]
# Flight imitator scored
ds[flight == 'F1'] |> nrow()## [1] 153## [1] 88# Any male following?
ds[, any_male_F0 := any(flight == 'F0' & sex == 'M'), by = obs_id]
unique(ds[N_ind_int > 1 & sex == 'F' & flight == 'F1' & any_male_F0 == TRUE], by = 'obs_id') |> nrow()## [1] 72## [1] 65# Any female following?
ds[, any_female_F0 := any(flight == 'F0' & sex == 'F'), by = obs_id]
unique(ds[N_ind_int > 1 & sex == 'M' & flight == 'F1' & any_female_F0 == TRUE], by = 'obs_id') |> nrow()## [1] 49## [1] 58## [1] 82## [1] 42## [1] 75### Aggressive interactions
# How many observations with at least two IDs and a aggression?
d[, any_aggres := any(!is.na(aggres)), by = obs_id]
# Total with at least 2 IDs
unique(d[N_ind_int > 1], by = 'obs_id') |> nrow()## [1] 2876## [1] 148## [1] 5.146036## [1] 72## [1] 2.503477# Total with a aggression initiator scored
d[, any_1_obs_id := any(aggres %like% '1'), by = obs_id]
d[, any_0_obs_id := any(aggres %like% '0'), by = obs_id]
# for each aggression
d[, any_1 := any(aggres %like% '1'), by = 1:nrow(d)]
d[, any_0 := any(aggres %like% '0'), by = 1:nrow(d)]
# both sexes aggressive in one observation
d[, pair_1 := any(any_1 == TRUE & sex == 'F') & any(any_1 == TRUE & sex == 'M'), by = obs_id]
# Which sex was target?
d[, any_female_0 := any(aggres %like% '0' & sex == 'F'), by = obs_id]
d[, any_male_0 := any(aggres %like% '0' & sex == 'M'), by = obs_id]
# subset relevant data
ds = d[!is.na(sex) & N_ind_int > 1 & any_1_obs_id == TRUE]
# Aggression imitator scored
ds[any_1 == TRUE] |> nrow()## [1] 144## [1] 62## [1] 82## [1] 9## [1] 40## [1] 25## [1] 46## [1] 34## [1] 43## [1] 66## [1] 41## [1] 57## [1] 57## [1] 42## [1] 6### Copulations
# subset data
ds = d[!is.na(cop)]
# initiators
ds[, initiator := like(cop, '1')]
ds[, .(cop, initiator)]## cop initiator
## 1: A0 FALSE
## 2: A1 TRUE
## 3: A0 FALSE
## 4: A1 TRUE
## 5: S1 TRUE
## ---
## 577: S FALSE
## 578: A FALSE
## 579: A FALSE
## 580: A0 FALSE
## 581: A1 TRUE## [1] 581## [1] 113## [1] 82## [1] 31## [1] 73## [1] 27## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=German_Germany.utf8 LC_CTYPE=German_Germany.utf8 LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C LC_TIME=German_Germany.utf8
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] gridExtra_2.3 DHARMa_0.4.6 ggh4x_0.2.3 performance_0.10.2 dplyr_1.1.0 officer_0.5.2 flextable_0.8.5 broomExtra_5.0.0
## [9] effects_4.2-2 carData_3.0-5 glmmTMB_1.1.5 activity_1.3.2 patchwork_1.1.2 doFuture_0.12.2 future_1.31.0 ggnewscale_0.4.8
## [17] stringr_1.5.0 knitr_1.42 sf_1.0-13 foreach_1.5.2 auksRuak_0.1.0 viridis_0.6.2 viridisLite_0.4.2 ggplot2_3.4.3
## [25] magrittr_2.0.3 data.table_1.14.8
##
## loaded via a namespace (and not attached):
## [1] backports_1.4.1 uuid_1.1-0 systemfonts_1.0.4 plyr_1.8.8 TMB_1.9.2 splines_4.2.2 listenv_0.9.0 gap.datasets_0.0.5
## [9] TH.data_1.1-2 digest_0.6.31 htmltools_0.5.4 fansi_1.0.4 memoise_2.0.1 doParallel_1.0.17 globals_0.16.2 sandwich_3.0-2
## [17] askpass_1.1 timechange_0.2.0 gfonts_0.2.0 colorspace_2.1-0 mitools_2.4 xfun_0.40 crayon_1.5.2 jsonlite_1.8.4
## [25] lme4_1.1-31 survival_3.4-0 zoo_1.8-11 iterators_1.0.14 glue_1.6.2 gtable_0.3.4 emmeans_1.8.4-1 scales_1.2.1
## [33] mvtnorm_1.1-3 DBI_1.1.3 Rcpp_1.0.10 xtable_1.8-4 units_0.8-2 proxy_0.4-27 survey_4.1-1 datawizard_0.6.5
## [41] ellipsis_0.3.2 farver_2.1.1 pkgconfig_2.0.3 qgam_1.3.4 sass_0.4.5 nnet_7.3-18 utf8_1.2.3 janitor_2.1.0
## [49] crul_1.3 labeling_0.4.3 tidyselect_1.2.0 rlang_1.1.1 ggspatial_1.1.7 later_1.3.0 munsell_0.5.0 tools_4.2.2
## [57] cachem_1.0.6 cli_3.6.0 generics_0.1.3 broom_1.0.3 evaluate_0.20 fastmap_1.1.0 yaml_2.3.7 zip_2.2.2
## [65] purrr_1.0.1 nlme_3.1-160 mime_0.12 xml2_1.3.3 gap_1.5-1 compiler_4.2.2 rstudioapi_0.14 curl_5.0.2
## [73] e1071_1.7-13 tibble_3.1.8 bslib_0.4.2 stringi_1.7.12 highr_0.10 parameters_0.20.2 forcats_1.0.0 gdtools_0.3.0
## [81] lattice_0.20-45 Matrix_1.5-1 classInt_0.4-9 nloptr_2.0.3 vctrs_0.5.2 furrr_0.3.1 pillar_1.9.0 lifecycle_1.0.3
## [89] jquerylib_0.1.4 sfext_0.1.0 estimability_1.4.1 insight_0.19.0 filenamr_0.1.0 httpuv_1.6.8 R6_2.5.1 promises_1.2.0.1
## [97] KernSmooth_2.23-20 parallelly_1.34.0 cliExtras_0.1.0 codetools_0.2-18 boot_1.3-28 MASS_7.3-58.1 openssl_2.0.5 rprojroot_2.0.3
## [105] withr_2.5.0 httpcode_0.3.0 broom.mixed_0.2.9.4 multcomp_1.4-25 mgcv_1.8-41 bayestestR_0.13.0 parallel_4.2.2 tidyr_1.3.0
## [113] coda_0.19-4 class_7.3-20 minqa_1.2.5 rmarkdown_2.20 snakecase_0.11.0 numDeriv_2016.8-1.1 shiny_1.7.4 lubridate_1.9.2
## [121] base64enc_0.1-3