Vocal interactions in Spix’s disc-winged bats (Thyroptera tricolor)
Data analysis for the paper:
G Chaverri, M Sagot, JL Stynoski, M Araya-Salas, Y Araya-Ajoy, M Nagy, M Knörnschild, S Chaves-Ramírez, N Rose, M Sánchez-Chavarría, Y Jiménez-Torres, D Ulloa-Sanabria, H Solís-Hernández, GG Carter. In review. Contact-calling rates within groups of disc-winged bats vary by caller but not by the receiver’s identity, kinship, or association. Philosophical Transactions of the Royal Society B.
# load function from sketchy
source("https://raw.githubusercontent.com/maRce10/sketchy/main/R/load_packages.R")
# install/ load packages
sketchy::load_packages(packages = c("igraph", "bipartite", "asnipe",
"assortnet", "ggplot2", "ggmap", "ecodist", "igraphdata", "statnet",
"RColorBrewer", "tidyverse", "geosphere", "mapproj", "sna", "stringr",
"lme4", "glmmTMB", "broom.mixed", "boot", "patchwork", "performance",
"MASS", "viridis"))# set ggplot theme
theme_set(theme_bw(base_size = 14))
# function to get mean and 95% CI of values x via bootstrapping
boot_ci <- function(x, perms = 5000, bca = F) {
get_mean <- function(x, d) {
return(mean(x[d]))
}
x <- as.vector(na.omit(x))
mean <- mean(x)
if (bca) {
boot <- boot.ci(boot(data = x, statistic = get_mean, R = perms,
parallel = "multicore", ncpus = 4), type = "bca")
low <- boot$bca[1, 4]
high <- boot$bca[1, 5]
} else {
boot <- boot.ci(boot(data = x, statistic = get_mean, R = perms,
parallel = "multicore", ncpus = 4), type = "perc")
low <- boot$perc[1, 4]
high <- boot$perc[1, 5]
}
c(low = low, mean = mean, high = high, N = round(length(x)))
}
# function to get mean and 95% CI via bootstrapping of values y
# within grouping variable x
boot_ci2 <- function(d = d, y = d$y, x = d$x, perms = 5000, bca = F) {
df <- data.frame(effect = unique(x))
df$low <- NA
df$mean <- NA
df$high <- NA
df$n.obs <- NA
for (i in 1:nrow(df)) {
ys <- y[which(x == df$effect[i])]
if (length(ys) > 1 & var(ys) > 0) {
b <- boot_ci(y[which(x == df$effect[i])], perms = perms,
bca = bca)
df$low[i] <- b[1]
df$mean[i] <- b[2]
df$high[i] <- b[3]
df$n.obs[i] <- b[4]
} else {
df$low[i] <- min(ys)
df$mean[i] <- mean(ys)
df$high[i] <- max(ys)
df$n.obs[i] <- length(ys)
}
}
df
}
# create function to convert matrix to data-frame
matrix_to_df <- function(m1) {
data.frame(dyad = paste(rownames(m1)[col(m1)], colnames(m1)[row(m1)],
sep = "_"), value = c(t(m1)), stringsAsFactors = FALSE)
}
# create function to plot variable by kinship with means and 95%
# CIs
plot_kinship_effect <- function(d = d, y = d$sri, label = "label") {
# plot association by kinship
set.seed(123)
means <- d %>%
mutate(y = y) %>%
filter(!is.na(kinship)) %>%
group_by(dyad) %>%
summarize(kinship = mean(kinship), y = mean(y), udyad.sex = first(udyad.sex)) %>%
mutate(kin = ifelse(kinship > 0, "kin", "nonkin")) %>%
dplyr::select(kin, y) %>%
boot_ci2(y = .$y, x = .$kin, bca = T) %>%
rename(kin = effect) %>%
mutate(kin = factor(kin, levels = c("nonkin", "kin")))
points <- d %>%
mutate(y = y) %>%
filter(!is.na(kinship)) %>%
group_by(dyad) %>%
summarize(kinship = mean(kinship), y = mean(y), udyad.sex = first(udyad.sex)) %>%
mutate(kin = ifelse(kinship > 0, "kin", "nonkin")) %>%
mutate(kin = factor(kin, levels = c("nonkin", "kin")))
# plot means and 95% CI
means %>%
ggplot(aes(x = kin, y = mean, color = kin)) + geom_jitter(data = points,
aes(y = y), size = 2, alpha = 0.5, width = 0.1) + geom_point(position = position_nudge(x = 0.25),
size = 3) + geom_errorbar(aes(ymin = low, ymax = high, width = 0.1),
position = position_nudge(x = 0.25), size = 1) + xlab("") +
ylab(label) + scale_color_viridis_d(alpha = 0.8, begin = 0.4,
end = 0.9, direction = -1, option = "G") + coord_flip() +
theme_bw(base_size = 20) + theme(legend.position = "none")
}# Using
# https://dshizuka.github.io/networkanalysis/example_make_sparrow_network.html
# Input data
batdat21 = read.csv("./data/raw/associations_2021.csv", header = TRUE,
sep = ";", colClasses = c("numeric", "character", "Date", "character",
"character", "character", "character", "character", "character",
"character", "character", "character", "character"))
batdat22 = read.csv("./data/raw/associations_2022.csv", header = TRUE,
sep = ";", colClasses = c("numeric", "character", "Date", "character",
"character", "character", "character", "character", "character",
"character", "character"))
# Add attribute data
attrib_21 <- read.csv("./data/raw/supp_21.csv", header = TRUE, sep = ";",
colClasses = c("character", "character"))
attrib_22 <- read.csv("./data/raw/supp_22.csv", header = TRUE, sep = ";",
colClasses = c("character", "character"))
# Prepare association data 2021
batcols = grep("Bat", colnames(batdat21))
# batcols batdat21[,batcols] gather(batdat21[,batcols])
bat.ids = unique(gather(batdat21[, batcols])$value)
# bat.ids
bat.ids = bat.ids[!is.na(bat.ids)]
# bat.ids
bat.ids_df_21 <- data.frame(bat.ids)
# csv_path <- 'bats_21.csv' write.csv(bat.ids_df_21, csv_path,
# row.names = FALSE)
m1_21 = apply(batdat21[, batcols], 1, function(x) match(bat.ids, x))
# m1_21
m1_21[is.na(m1_21)] = 0
m1_21[m1_21 > 0] = 1
# m1_21
rownames(m1_21) = bat.ids #rows are bat ids
colnames(m1_21) = paste("group", 1:ncol(m1_21), sep = "_") #columns are group IDs (just 'group_#')
# m1_21
# rowSums(m1_21) set.seed(2) png(file='./output/Histogram
# captures 2021.png', width=600, height=400) par(oma=c(1,1,1,1))
# # all sides have 3 lines of space par(mar=c(5,4,4,2) + 0.8)
# hist(rowSums(m1_21), main='', xlab='Number of observations',
# ylab='Frequency', cex.lab = 2.5, cex.axis = 1.5) dev.off()
average_n_21 <- mean(rowSums(m1_21)) #average number of times a single bat was observed
# average_n_21
sd_n_21 <- sd(rowSums(m1_21)) #sd number of times a single bat was observed
# sd_n_21
average_gs_21 <- mean(colSums(m1_21)) #average group size
# average_gs_21
sd_gs_21 <- sd(colSums(m1_21)) #sd group size
# sd_gs_21
# Prepare association data 2022
batcols = grep("Bat", colnames(batdat22))
# batcols batdat22[,batcols] gather(batdat22[,batcols])
bat.ids = unique(gather(batdat22[, batcols])$value)
# bat.ids
bat.ids = bat.ids[is.na(bat.ids) == F]
# bat.ids
bat.ids_df_22 <- data.frame(bat.ids)
# csv_path = 'bats_22.csv' write.csv(bat.ids_df_22, csv_path,
# row.names = FALSE)
m1_22 = apply(batdat22[, batcols], 1, function(x) match(bat.ids, x))
# m1_22
m1_22[is.na(m1_22)] = 0
m1_22[m1_22 > 0] = 1
# m1_22
rownames(m1_22) = bat.ids #rows are bat ids
colnames(m1_22) = paste("group", 1:ncol(m1_22), sep = "_") #columns are group IDs (just 'group_#')
# m1_22
# rowSums(m1_22) set.seed(2) png(file='./output/Histogram
# captures 2022.png', width=600, height=400) par(oma=c(1,1,1,1))
# # all sides have 3 lines of space par(mar=c(5,4,4,2) + 0.8)
# hist (rowSums(m1_22), main='', xlab='Number of observations',
# ylab='Frequency', cex.lab = 2.5, cex.axis = 1.5) dev.off()
average_n_22 <- mean(rowSums(m1_22)) #average number of times a single bat was observed
# average_n_22
sd_n_22 <- sd(rowSums(m1_22)) #sd number of times a single bat was observed
# sd_n_22
average_gs_22 <- mean(colSums(m1_22)) #average group size
# average_gs_22
sd_gs_22 <- sd(colSums(m1_22)) #average group size
# sd_gs_22Generating 76 x 76 matrixg_21=graph_from_adjacency_matrix(adj_21, "undirected", weighted=T) #the igraph object
V(g_21)$sex <- factor(attrib_21[match(V(g_21)$name, attrib_21$bat_id), "sex"])
mismatched_names <- setdiff(V(g_21)$name, attrib_21$bat_id) #check for issues
#Create adjacency data for associations (2022)
adj_22=get_network(t(m1_22), data_format="GBI", association_index = "SRI") # the adjacency matrixGenerating 76 x 76 matrixg_22=graph_from_adjacency_matrix(adj_22, "undirected", weighted=T) #the igraph object
V(g_22)$sex <- factor(attrib_22[match(V(g_22)$name, attrib_22$bat_id), "sex"])
mismatched_names <- setdiff(V(g_22)$name, attrib_22$bat_id) #check for issues
#Select individuals for which more than 10 observations were made (not incorporated into adjacency data for associations, above)
#m2_21=m1_21[which(rowSums(m1_21)>10),]
#adj_21=get_network(t(m2_21), data_format="GBI","SRI")
#write.csv(adj_21, "assoc_21.csv", row.names=TRUE)
#rowSums(m2_21)
#m2_22=m1_22[which(rowSums(m1_22)>10),]
#adj_22=get_network(t(m2_22), data_format="GBI","SRI")
#rowSums(m2_22)
#write.csv(adj_22, "assoc_22.csv", row.names=TRUE)
#rowSums(m2_22)
#Create communities (2021)
com_21=fastgreedy.community(g_21) #community detection method
# com_21[[11]]
# Assign community names
mapping <- c(
11, # Old membership 1 should be corrected to 11
8, # Old membership 2 should be corrected to 8
1, # Old membership 3 should be corrected to 1
2, # Old membership 4 should be corrected to 2
6, # Old membership 5 should be corrected to 6
9, # Old membership 6 should be corrected to 9
4, # Old membership 7 should be corrected to 4
10, # Old membership 8 should be corrected to 10
7, # Old membership 9 should be corrected to 7
5, # Old membership 10 should be corrected to 5
3 # Old membership 11 should be corrected to 3
)
corrected_memberships <- mapping[com_21$membership]
com_21$membership <- corrected_memberships
community_assignments <- com_21$membership
color_mapping <- c(mako(n = 2, alpha = 0.7, begin = 0.4, end = 0.9, direction = -1), "gray")
names(color_mapping) <- c("female", "male", "unknown")
node_sex <- V(g_21)$sex
node_colors <- sapply(node_sex, function(sex) {
color_mapping[sex]
})
V(g_21)$color <- node_colors
#Prepare to save figure
# png(file="./output/Association networks.png", width=1200 * 3, height=600 * 3, res = 300)
# pdf(file="./output/Association networks.pdf", width=1200 * 1, height=800 * 1)
par(mfrow = c(1, 2), xpd = TRUE)
# Iterate through com_21 and assign community names to nodes in g_21
set.seed(123)
plot(
g_21,
vertex.color = node_colors, # Use the calculated node colors
vertex.size = 10, # Adjust the size of the nodes as needed
vertex.label = community_assignments, # Use the extracted community assignments as labels
edge.width=E(g_21)$weight*10
)
title(main = "2021", cex.main = 2)
# legend("topleft", legend = names(color_mapping), fill = color_mapping, bty = "n", inset=c(0.66, 0.6), x.intersp = 0.2, y.intersp = 0.55, cex = 1.1)
points(x = c(-1.12, -0.32, 0.34) + 0.2, y = rep(-1.4, 3), pch = 21, cex = 4, bg = color_mapping, col = "black")
text(x = c(-1.12, -0.32, 0.34) + 0.26, y = rep(-1.4, 3), labels = names(color_mapping), pos = 4, cex = 2)
#Create communities (2022)
com_22=fastgreedy.community(g_22) #community detection method
# Assign community names
mapping <- c(
3,
4,
1,
9,
10,
11,
7,
5,
6,
2,
12,
8
)
corrected_memberships <- mapping[com_22$membership]
com_22$membership <- corrected_memberships
community_assignments <- com_22$membership
node_sex <- V(g_22)$sex
node_colors <- sapply(node_sex, function(sex) {
color_mapping[sex]
})
V(g_22)$color <- node_colors
#Prepare to save figure
# png(file="./output/Association networks 2022.png",width=600, height=600)
# Iterate through com_21 and assign community names to nodes in g_21
set.seed(123)
plot(
g_22,
vertex.color = node_colors, # Use the calculated node colors
vertex.size = 10, # Adjust the size of the nodes as needed
vertex.label = community_assignments, # Use the extracted community assignments as labels
edge.width=E(g_22)$weight*10
)
title(main = "2022", cex.main = 2)
Modularity:
# Data manipulation Joining the two data sets
d1 <- read.csv("./data/raw/vocal_interactions_2021.csv", sep = ";")
d2 <- read.csv("./data/raw/vocal_interactions_2022.csv", sep = ";")
colnames(d2)[6] <- "Dyad"
d <- rbind(d1, d2)
## Adding observation level random effect
d$ob <- 1:nrow(d)
# Stats models Model for inquiry
mod_inquiry <- glmer.nb(n_inquiry ~ 1 + (1 | Bat_roosting) + (1 |
Bat_flying) + (1 | Group), data = d)
summary(mod_inquiry)Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: Negative Binomial(2.7219) ( log )
Formula: n_inquiry ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) + (1 |
Group)
Data: d
AIC BIC logLik deviance df.resid
4596.7 4617.2 -2293.3 4586.7 446
Scaled residuals:
Min 1Q Median 3Q Max
-1.62529 -0.48778 0.07217 0.47606 2.82516
Random effects:
Groups Name Variance Std.Dev.
Bat_roosting (Intercept) 8.029e-11 8.961e-06
Bat_flying (Intercept) 2.348e-01 4.845e-01
Group (Intercept) 2.141e-02 1.463e-01
Number of obs: 451, groups: Bat_roosting, 74; Bat_flying, 73; Group, 14
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.12062 0.07845 52.53 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')### Extracting variance components
mu_l <- fixef(mod_inquiry)[1]
mu <- exp(mu_l)
theta <- getME(mod_inquiry, "glmer.nb.theta")
VBF_I_l <- VarCorr(mod_inquiry)$Bat_flying[1, 1] #Variance Bats Flying
VBF_I <- (exp(VBF_I_l) - 1) * exp(2 * mu_l + VBF_I_l)
VBR_I_l <- VarCorr(mod_inquiry)$Bat_roosting[1, 1] #Variance Bats Roosting
VBR_I <- (exp(VBR_I_l) - 1) * exp(2 * mu_l + VBR_I_l)
VBG_I_l <- VarCorr(mod_inquiry)$Group[1, 1] #Variance GRoup
VBG_I <- (exp(VBG_I_l) - 1) * exp(2 * mu_l + VBG_I_l)
NBV_I <- mu + mu^2/theta
#### Repeatabilty roosting individual
print("Repeatabilty roosting individual")[1] "Repeatabilty roosting individual" (Intercept)
1.084486e-10 [1] "Repeatabilty flying individual"(Intercept)
0.465914 Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: Negative Binomial(4.214) ( log )
Formula: n_inquiry ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) + (1 |
Group)
Data: d2
AIC BIC logLik deviance df.resid
2606.7 2624.4 -1298.4 2596.7 250
Scaled residuals:
Min 1Q Median 3Q Max
-1.93767 -0.52024 0.04805 0.49889 2.73912
Random effects:
Groups Name Variance Std.Dev.
Bat_roosting (Intercept) 7.930e-13 8.905e-07
Bat_flying (Intercept) 2.221e-01 4.712e-01
Group (Intercept) 3.879e-02 1.970e-01
Number of obs: 255, groups: Bat_roosting, 71; Bat_flying, 69; Group, 14
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.21595 0.08999 46.85 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')### Extracting variance components
mu_l <- fixef(mod_inquiry)[1]
mu <- exp(mu_l)
theta <- getME(mod_inquiry, "glmer.nb.theta")
VBF_I_l <- VarCorr(mod_inquiry)$Bat_flying[1, 1] #Variance Bats Flying
VBF_I <- (exp(VBF_I_l) - 1) * exp(2 * mu_l + VBF_I_l)
VBR_I_l <- VarCorr(mod_inquiry)$Bat_roosting[1, 1] #Variance Bats Roosting
VBR_I <- (exp(VBR_I_l) - 1) * exp(2 * mu_l + VBR_I_l)
VBG_I_l <- VarCorr(mod_inquiry)$Group[1, 1] #Variance GRoup
VBG_I <- (exp(VBG_I_l) - 1) * exp(2 * mu_l + VBG_I_l)
NBV_I <- mu + mu^2/theta
#### Repeatabilty roosting individual
print("Repeatabilty roosting individual")[1] "Repeatabilty roosting individual" (Intercept)
1.409571e-12 [1] "Repeatabilty flying individual"(Intercept)
0.5519081 Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: Negative Binomial(0.1622) ( log )
Formula: n_response ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) + (1 |
Group)
Data: d
AIC BIC logLik deviance df.resid
3742.3 3762.8 -1866.1 3732.3 446
Scaled residuals:
Min 1Q Median 3Q Max
-0.4026 -0.4017 -0.3616 0.2471 4.1603
Random effects:
Groups Name Variance Std.Dev.
Bat_roosting (Intercept) 1.118e+00 1.057e+00
Bat_flying (Intercept) 5.331e-13 7.301e-07
Group (Intercept) 1.489e-12 1.220e-06
Number of obs: 451, groups: Bat_roosting, 74; Bat_flying, 73; Group, 14
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.074 0.213 19.12 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')### Extracting variance components
mu_l <- fixef(mod_response)[1] #+ fixef(mod_response)[2]*mean(d$n_inquiry, na.rm=TRUE)
mu <- exp(mu_l)
theta <- getME(mod_response, "glmer.nb.theta")
VBF_R_l <- VarCorr(mod_response)$Bat_flying[1, 1] #Variance Bats Flying
VBF_R <- (exp(VBF_R_l) - 1) * exp(2 * mu_l + VBF_R_l)
VBR_R_l <- VarCorr(mod_response)$Bat_roosting[1, 1] #Variance Bats Roosting
VBR_R <- (exp(VBR_R_l) - 1) * exp(2 * mu_l + VBR_R_l)
VBG_R_l <- VarCorr(mod_response)$Group[1, 1] #Variance GRoup
VBG_R <- (exp(VBG_R_l) - 1) * exp(2 * mu_l + VBG_R_l)
NBV_R <- mu + mu^2/theta
#### Repeatabilty roosting individual
print("Repeatabilty roosting individual")[1] "Repeatabilty roosting individual"(Intercept)
0.5044 [1] "Repeatabilty flying individual" (Intercept)
4.274015e-14 # #Data manipulation ## Joining the two data sets
# d1<-read.csv('vocal_interactions_2021.csv', sep=';')
# d2<-read.csv('vocal_interactions_2022.csv', sep=';')
# colnames(d2)[6]<-'Dyad' d<-rbind(d1,d2)
## Adding observation level random effect
d$ob <- 1:nrow(d)
# Stats models Model for inquiry with all observations
mod_inquiry <- glmer(n_inquiry ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) +
(1 | Group) + (1 | ob), data = d, family = "poisson")
### Extracting variance components
VBF_I <- VarCorr(mod_inquiry)$Bat_flying[1, 1] #Variance Bats Flying
VBR_I <- VarCorr(mod_inquiry)$Bat_roosting[1, 1] #Variance Bats Roosting
VBG_I <- VarCorr(mod_inquiry)$Group[1, 1] #Variance Group
VOD_I <- VarCorr(mod_inquiry)$ob[1, 1] #Variance Overddispersion
PV_I <- log(1/exp(fixef(mod_inquiry)[1]) + 1) #Variance poisson processs
#### Repeatabilty roosting individual
VBR_I/(VBR_I + VOD_I + VBF_I + PV_I + VBG_I)(Intercept)
0.07582562 (Intercept)
0.3913205 ## Model for inquiry only when someone responded
d2 <- d[d$n_response > 0, ]
mod_inquiry2 <- glmer(n_inquiry ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) +
(1 | ob), data = d2, family = "poisson")
### Extracting variance components
VBF_I <- VarCorr(mod_inquiry2)$Bat_flying[1, 1] #Variance Bats Flying
VBR_I <- VarCorr(mod_inquiry2)$Bat_roosting[1, 1] #Variance Bats Roosting
VBG_I <- VarCorr(mod_inquiry)$Group[1, 1]
VOD_I <- VarCorr(mod_inquiry2)$ob[1, 1] #Variance Overddispersion
PV_I <- log(1/exp(fixef(mod_inquiry2)[1]) + 1) #Variance poisson processs
#### Repeatabilty roosting individual
print("Repeatabilty roosting individual")[1] "Repeatabilty roosting individual"(Intercept)
0.02516122 [1] "Repeatabilty flying individual"(Intercept)
0.5240666 Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: poisson ( log )
Formula: n_response ~ 1 + (1 | Bat_roosting) + (1 | Bat_flying) + (1 |
Group) + (1 | ob)
Data: d
AIC BIC logLik deviance df.resid
3816.4 3837.0 -1903.2 3806.4 446
Scaled residuals:
Min 1Q Median 3Q Max
-0.75743 -0.39416 -0.00652 0.01611 0.17248
Random effects:
Groups Name Variance Std.Dev.
ob (Intercept) 8.3383 2.8876
Bat_roosting (Intercept) 7.8244 2.7972
Bat_flying (Intercept) 0.4551 0.6746
Group (Intercept) 1.0819 1.0401
Number of obs: 451, groups:
ob, 451; Bat_roosting, 74; Bat_flying, 73; Group, 14
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.1595061 0.0004604 2518 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (Nelder_Mead) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.15353 (tol = 0.002, component 1)
Model is nearly unidentifiable: very large eigenvalue
- Rescale variables?### Extracting variance components
VBF_R <- VarCorr(mod_response)$Bat_flying[1, 1] #Variance Bats Flying
VBR_R <- VarCorr(mod_response)$Bat_roosting[1, 1] #Variance Bats Roosting
VBG_R <- VarCorr(mod_response)$Group[1, 1]
VOD_R <- VarCorr(mod_response)$ob[1, 1] #Variance Overddispersion
PV_R <- log(1/exp(fixef(mod_response)[1]) + 1) #Variance poisson processs
## Repeatabilty roosting individual
print("Repeatabilty roosting individual")[1] "Repeatabilty roosting individual"(Intercept)
0.4353549 [1] "Repeatabilty flying individual"(Intercept)
0.02532077 # Code by Gerald Carter, gcarter1640@gmail.com
### clean and wrangle data
# get kinship data
k <- as.matrix(read.csv('./data/raw/KinshipMatrix.csv', sep= ",", row.names = 1))
# check symmetry
mean(k==t(k), na.rm=T)
# get sex data
sex1 <- read.csv('supp_21.csv', sep= ";")
sex2 <- read.csv('supp_22.csv', sep= ";")
# get association data
a21 <- read.csv('associations_2021.csv', sep= ";")
a22 <- read.csv('associations_2022.csv', sep= ";")
# get calling data
i21 <- read.csv('vocal_interactions_2021.csv', sep= ";")
i22 <- read.csv('vocal_interactions_2022.csv', sep= ";")
# tidy sex data
sex <-
rbind(sex1,sex2) %>%
mutate(bat= paste0('bat', bat_id)) %>%
group_by(bat, sex) %>%
summarize(n= n()) %>%
dplyr::select(bat, sex) %>%
ungroup()
# tidy 2021 association data
assoc21 <-
a21 %>%
pivot_longer(Bat1:Bat10, names_to = 'names', values_to = "bat") %>%
filter(!is.na(bat)) %>%
mutate(bat= paste0('bat', bat)) %>%
mutate(group= paste0('group', Group)) %>%
mutate(date= dmy(Date)) %>%
dplyr::select(obs= data_id, group, date, bat)
# tidy 2022 association data
assoc22 <-
a22 %>%
pivot_longer(Bat1:Bat8, names_to = 'names', values_to = "bat") %>%
filter(!is.na(bat)) %>%
mutate(bat= paste0('bat', bat)) %>%
mutate(group= paste0('group', Group)) %>%
mutate(date= dmy(Date)) %>%
dplyr::select(obs= data_id, group, date, bat)
# combine association data
assoc <-
rbind(assoc21, assoc22) %>%
mutate(obs= paste(obs,group,date, sep= "_")) %>%
dplyr::select(bat, obs) %>%
as.data.frame()
# get number of observations per bat
obs_per_bat <-
assoc %>%
group_by(bat) %>%
summarize(n.obs=n()) %>%
arrange(desc(n.obs))
# range is 1 to 51 times
range(obs_per_bat$n.obs)
# pick threshold of observations for including bats in analysis
# if they have fewer observations than we make their association rates NA
threshold <- 4
# plot number of observations per bat
assoc %>%
group_by(bat) %>%
summarize(n.obs=n()) %>%
ggplot(aes(x=n.obs))+
geom_histogram(fill="light blue", color="black")+
geom_vline(xintercept= threshold, color= 'red', size=1)+
xlab("number of observations of bat")+
ylab("count of bats")
# get bats seen fewer than threshold number
low.n.bats <-
assoc %>%
group_by(bat) %>%
summarize(n.obs=n()) %>%
filter(n.obs < threshold) %>%
pull(bat)
low.n.bats
length(low.n.bats)
# 14 bats were seen fewer than 4 times
### get SRI from asnipe
# get association rates as group-by-individual
gbi <- get_group_by_individual(assoc, data_format="individuals")
# get association network of SRI
assoc.net<- get_network(association_data=gbi, data_format = "GBI", association_index = "SRI")
# check symmetry
mean(assoc.net==t(assoc.net))
# get SRI values for every undirected pair as dataframe
SRI <- matrix_to_df(assoc.net)
# tidy kinship
kinship <-
k %>%
matrix_to_df() %>%
separate(dyad, into= c("bat1", "bat2")) %>%
mutate(bat2= str_remove(bat2, "X")) %>%
mutate(bat1= paste0('bat', bat1), bat2= paste0('bat', bat2)) %>%
filter(bat1!=bat2) %>%
# label undirected pair
mutate(dyad= ifelse(bat1<bat2, paste(bat1,bat2, sep="_"), paste(bat2,bat1, sep="_"))) %>%
group_by(dyad) %>%
summarize(kinship= first(value))
# count trials
nrow(i21)
nrow(i22)
# trials with missing data
rbind(i21,i22) %>%
as_tibble() %>%
filter(is.na(n_response)) %>%
nrow()
#2
# trials where neither bat called
rbind(i21,i22) %>%
as_tibble() %>%
filter(n_inquiry==0) %>%
nrow()
# 5
# fix and tidy calling data
calling <-
rbind(i21,i22) %>%
mutate(Date= as.character(mdy(Date))) %>%
mutate(group= paste(substr(Date, start=1, stop=4),Group, sep="_")) %>%
# need to recreate these labels because messed up by excel in raw data
separate(Dyad, into=c("bat_flying", "bat_roosting"), remove=F) %>%
mutate(bat_flying= paste0('bat', bat_flying), bat_roosting= paste0('bat', bat_roosting)) %>%
# label undirected dyads
mutate(dyad= ifelse(bat_flying<bat_roosting, paste(bat_flying, bat_roosting, sep= "_"), paste(bat_roosting, bat_flying, sep= "_"))) %>%
dplyr::select(date= Date, group, bat_flying, bat_roosting, dyad, flight_time, n_inquiry, n_response)
# combine all data
d <-
# start with calling data
calling %>%
# add SRI
mutate(sri= SRI$value[match(.$dyad, SRI$dyad)]) %>%
# add kinship
mutate(kinship= kinship$kinship[match(.$dyad, kinship$dyad)]) %>%
# add sexes of flying and roosting bats
mutate(flying.sex= sex$sex[match(.$bat_flying, sex$bat)]) %>%
mutate(roosting.sex= sex$sex[match(.$bat_roosting, sex$bat)]) %>%
# label sex combination for flying-->roosting dyad
mutate(dyad.sex= paste(flying.sex, roosting.sex, sep= "-->")) %>%
# label sex combination for dyad (male, female, both)
mutate(udyad.sex = case_when(
flying.sex == 'female' & roosting.sex == 'female' ~ "female-female",
flying.sex == 'male' & roosting.sex == 'male' ~ "male-male",
TRUE ~ "mixed-sex")) %>%
# replace zero sri (never previously seen together) with NA because association is not 0
mutate(sri = if_else(sri==0, NA, sri)) %>%
# if flying bat seen fewer than 4 times, then SRI is NA
mutate(sri = if_else(bat_flying %in% low.n.bats, NA, sri)) %>%
# if roosting bat seen fewer than 4 times, then SRI is NA
mutate(sri = if_else(bat_roosting %in% low.n.bats, NA, sri))
# inspect and fix cases where flying bats in more than 2 groups
t <-
d %>%
group_by(bat_flying, date, group) %>%
summarize(n=n()) %>%
arrange(date) %>%
group_by(bat_flying, group) %>%
summarize(n=n()) %>%
group_by(bat_flying) %>%
summarize(n=n()) %>%
filter(n>2) %>%
pull(bat_flying)
d %>%
dplyr::select(date, group, bat_flying) %>%
filter(bat_flying %in% t) %>%
arrange(bat_flying) %>%
as.data.frame()
# some of these seem impossible and must be errors
# fix impossible group labels
d$group[which(d$date=="2021-07-03" & d$bat_flying=="bat982126051278521")] <- "2021_9"
d$group[which(d$date=="2021-07-03" & d$bat_flying=="bat982126058484300")] <- "2021_9"
# fix cases where roosting bats in more than 2 groups
t <-
d %>%
group_by(bat_roosting, date, group) %>%
summarize(n=n()) %>%
arrange(date) %>%
group_by(bat_roosting, group) %>%
summarize(n=n()) %>%
group_by(bat_roosting) %>%
summarize(n=n()) %>%
filter(n>2) %>%
pull(bat_roosting)
d %>%
dplyr::select(date, group, bat_roosting) %>%
filter(bat_roosting %in% t) %>%
arrange(bat_roosting) %>%
as.data.frame()
# fix group labels
d$group[which(d$date=="2021-07-03" & d$bat_roosting=="bat982126058484300")] <- "2021_9"
# group 9 split into groups 9,1 and 9,2 during 2022
# There are interesting movements between group 9 and 10 and between groups 9,1 and 9,2
# For this analysis, I'm going to consider group 9,1 and 9,2 as the same group
d$group[which(d$group=="2022_9,1" | d$group=="2022_9,2")] <- "2021_9"
# remove trials without inquiry or response calls
d <-
d %>% # 453 rows
as_tibble() %>%
filter(! (is.na(n_inquiry) & is.na(n_response))) %>% # 451 rows
filter(n_inquiry>0) %>% # 446 rows
mutate(year= substr(date, 1,4)) # add year
# save data as csv
write.csv(d, file = paste0('./processed/calling-association-kinship-data.csv'), row.names= F)# get data
raw <- read.csv("./data/processed/calling-association-kinship-data.csv")
# look at data head(raw) group is the year and social group dyad
# is the undirected pair (flying and roosting bat) in
# alphanumeric order sri is simple ratio index of association
# flying.sex is sex of flying bat roosting.sex is sex of
# roosting bat dyad.sex is the sex of the flying and roosting
# bat udyad.sex is females, males, or mixed
# add a few variables to the data
d <- raw %>%
# get log count of inquiry calls
mutate(log_inquiry = log(n_inquiry)) %>%
# rescale kinship so 1 unit is 0.5
mutate(kinship2 = kinship * 2) %>%
# if the roosting bat leaves the roost (escapes) then flight
# time is < 300 s
mutate(escape = as.numeric(flight_time < 300)) %>%
# convert flight time from seconds to minutes (for easier
# interpretation)
mutate(flight_time2 = flight_time/60)
###### Describe the data-------------
# how many trials? d %>% nrow() 446
# # how many undirected pairs have vocal data? d %>% pull(dyad)
# %>% unique() %>% length 139 undirected pairs
# d %>% group_by(bat_flying, bat_roosting) %>% summarize(n=n())
# 254 directed pairs
# how many groups? n_distinct(d$group) 23
# how many of these have association data d %>% group_by(dyad)
# %>% summarize(sri= mean(sri)) %>% filter(!is.na(sri)) 133
# pairs have association
# how many of these have kinship data d %>% group_by(dyad) %>%
# summarize(kinship= mean(kinship)) %>% filter(!is.na(kinship))
# 82 pairs have kinship data
# how many related undirected pairs d %>% group_by(dyad) %>%
# summarize(kinship= mean(kinship)) %>% filter(kinship>0) 32 kin
# pairs
# how many unrelated undirected pairs d %>% group_by(dyad) %>%
# summarize(kinship= mean(kinship)) %>% filter(kinship==0) 50
# nonkin pairs
# what is mean kinship in group?
grp.kin <- d %>%
group_by(dyad, group) %>%
summarize(kinship = mean(kinship, na.rm = T)) %>%
group_by(group) %>%
summarize(kinship = mean(kinship, na.rm = T), n = n()) %>%
as.data.frame()
set.seed(123)
# grp.kin %>% pull(kinship) %>% boot_ci(bca=T) 0.23, 95% CI =
# [0.16, 0.31]
# how many groups have kin? grp.kin %>% filter(kinship>0) 17
# how many groups have no kin? grp.kin %>% filter(kinship==0) 4
# plot distribution of kinship
d %>%
ggplot(aes(x = kinship)) + facet_wrap(~udyad.sex, ncol = 1, scales = "free_y") +
geom_histogram(fill = viridis(10, alpha = 0.6)[3], color = "black") +
ggtitle("pairwise kinship")
# # check categories d %>% filter(!is.na(kinship)) %>%
# group_by(dyad) %>% summarize(kinship= mean(kinship)) %>%
# group_by(kinship) %>% summarize(n=n())
##### Effect of kinship on association-------------
# plot distribution of assoc
d %>%
ggplot(aes(x = sri)) + facet_wrap(~udyad.sex, ncol = 1, scale = "free_y") +
geom_histogram(fill = viridis(10, alpha = 0.6)[3], color = "black") +
xlab("association rate") + ggtitle("pairwise SRI values")
# looks ok
# what is mean and 95% CI of assoc among flying and roosting
# bats? set.seed(123) d %>% group_by(dyad) %>% summarize(sri =
# mean(sri)) %>% pull(sri) %>% boot_ci(bca=T) 0.51, 95% CI =
# [0.47, 0.55] a bit low, expected from past work is 0.76
# now only nonkin
set.seed(123)
k1 <- d %>%
filter(kinship == 0) %>%
group_by(dyad) %>%
summarize(sri = mean(sri)) %>%
pull(sri) %>%
boot_ci(bca = T) %>%
c(kinship = 0)
# k1 0.572, 95% CI = [0.500, 0.636]
# now only kin set.seed(123) d %>% filter(kinship>0) %>%
# group_by(dyad) %>% summarize(sri = mean(sri)) %>% pull(sri)
# %>% boot_ci(bca=T) 0.686, 95% CI = [0.631, 0.741]
# now only close kin
set.seed(123)
k2 <- d %>%
filter(kinship == 0.5) %>%
group_by(dyad) %>%
summarize(sri = mean(sri)) %>%
pull(sri) %>%
boot_ci(bca = T) %>%
c(kinship = 0.5)
# k2 0.687, 95% CI = [0.628, 0.750]
# now only 0.25 kin
set.seed(123)
k3 <- d %>%
filter(kinship == 0.25) %>%
group_by(dyad) %>%
summarize(sri = mean(sri)) %>%
pull(sri) %>%
boot_ci(bca = T) %>%
c(kinship = 0.25)
# k3 0.684, 95% CI = [0.542, 0.783]
# compile means
k <- rbind(k1, k2, k3) %>%
data.frame() %>%
mutate(kin = as.character(kinship)) %>%
mutate(kin2 = kinship > 0) %>%
mutate(assoc = mean)
# plot association by kinship
(kin.dyads.plot <- d %>%
filter(!is.na(kinship)) %>%
group_by(dyad) %>%
summarize(kinship = mean(kinship), assoc = mean(sri), udyad.sex = first(udyad.sex)) %>%
mutate(kin2 = kinship > 0) %>%
mutate(kin = as.character(kinship)) %>%
ggplot(aes(x = kin, y = assoc, group = kin, color = kin2)) + geom_jitter(size = 2,
width = 0.1, height = 0, aes(shape = udyad.sex)) + geom_point(data = k,
position = position_nudge(x = 0.25), size = 3) + geom_errorbar(data = k,
aes(ymin = low, ymax = high, width = 0.1), position = position_nudge(x = 0.25),
size = 1) + xlab("kinship") + ylab("association rate (simple ratio index)") +
scale_color_viridis_d(alpha = 0.8, begin = 0.4, end = 0.9, direction = -1,
option = "G") + theme_bw(base_size = 20) + theme(legend.position = "none",
legend.title = element_blank()))
# save as PDF ggsave( 'kin_association.pdf', plot =
# kin.dyads.plot, scale = 1, width = 3, height = 5, units =
# 'in', dpi = 300)
##### Effect of flight time on calling ----------
### How does flight time and count of inquiry calls predict
### count of response calls? how many trials ended early
length(d$flight_time)[1] 446[1] 70[1] 4 date group bat_flying bat_roosting
1 2021-07-02 2021_7 bat982126051278504 bat982126051278564
2 2021-07-03 2021_10 bat982126051278540 bat982126052945921
3 2021-07-04 2021_2 bat982126057845162 bat982126057845205
4 2021-06-28 2021_3 bat900200000279490 bat982126058484334
dyad flight_time n_inquiry n_response
1 bat982126051278504_bat982126051278564 360 203 0
2 bat982126051278540_bat982126052945921 320 37 22
3 bat982126057845162_bat982126057845205 301 166 283
4 bat900200000279490_bat982126058484334 315 30 0
sri kinship flying.sex roosting.sex dyad.sex udyad.sex year
1 0.8636364 0.5 female female female-->female female-female 2021
2 0.4603175 0.5 female female female-->female female-female 2021
3 0.8437500 0.5 male female male-->female mixed-sex 2021
4 0.8709677 0.5 female male female-->male mixed-sex 2021
log_inquiry kinship2 escape flight_time2
1 5.313206 1 0 6.000000
2 3.610918 1 0 5.333333
3 5.111988 1 0 5.016667
4 3.401197 1 0 5.250000# plot number of response calls by flight time
t1 <- d %>%
mutate(kinship = kinship > 0.1) %>%
ggplot(aes(x = flight_time, y = n_response)) + geom_point(size = 2,
alpha = 0.7, aes(color = kinship, shape = kinship)) + geom_smooth(method = "glm.nb",
color = mako(10)[9]) + scale_color_viridis_d(alpha = 1, begin = 0.4,
end = 0.9, direction = -1, option = "G") + xlab("flight time (seconds)") +
ylab("response call count") + theme_bw(base_size = 20) + theme(legend.position = "none")
# plot number of inquiry calls by flight time
t2 <- d %>%
mutate(kinship = kinship > 0.1) %>%
ggplot(aes(x = flight_time, y = n_inquiry)) + geom_point(size = 2,
alpha = 0.7, aes(color = kinship, shape = kinship)) + geom_smooth(method = "glm.nb",
color = mako(10)[9]) + scale_color_viridis_d(alpha = 1, begin = 0.4,
end = 0.9, direction = -1, option = "G") + xlab("flight time (seconds)") +
ylab("inquiry call count") + theme_bw(base_size = 20) + theme(legend.position = "none")
# plot together
t <- t1 + t2
t
# Overdispersion test
dispersion ratio = 6.288
Pearson's Chi-Squared = 2773.032
p-value = < 0.001# Overdispersion test
dispersion ratio = 0.786
Pearson's Chi-Squared = 345.749
p-value = 1# A tibble: 2 × 9
effect component term estimate std.error statistic p.value conf.low
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 fixed cond (Intercept) 11.9 1.51 19.5 7.77e-85 9.30
2 fixed cond flight_time2 1.43 0.0336 15.1 2.60e-51 1.36
# ℹ 1 more variable: conf.high <dbl># for every minute of flight time, the inquiry call count
# increases by a factor of 1.43 [1.36, 1.49]
# try poisson model for effect of flight time on response
# calling
t <- glmer(n_response ~ flight_time2 + n_inquiry + (1 | bat_flying) +
(1 | bat_roosting) + (1 | group), data = d, family = poisson)
check_overdispersion(t)# Overdispersion test
dispersion ratio = 62.090
Pearson's Chi-Squared = 27319.658
p-value = < 0.001# Overdispersion test
dispersion ratio = 0.564
Pearson's Chi-Squared = 246.887
p-value = 1# Check for zero-inflation
Observed zeros: 191
Predicted zeros: 29
Ratio: 0.15# A tibble: 4 × 9
effect component term estimate std.error statistic p.value conf.low conf.high
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 fixed cond (Int… 6.71 3.13 4.08 4.49e-5 2.79 17.7
2 fixed cond flig… 1.67 0.176 4.84 1.28e-6 1.34 2.03
3 fixed cond n_in… 1.01 0.00256 2.32 2.02e-2 1.00 1.01
4 fixed zi (Int… 0.647 0.0723 -3.89 1.00e-4 2.79 17.7 # for every minute of flight time, the response call count
# increases by a factor of 1.67 [1.34, 2.03]
##### Effect of inquiry calling on response calling------------
# plot with log counts
d %>%
ggplot(aes(x = n_inquiry, y = n_response)) + geom_point(size = 2) +
geom_smooth(method = "lm") + xlab("inquiry call count") + ylab("response call count")
# plot negative binomial curve plot all bats
d %>%
mutate(label = "all pairs") %>%
mutate(escape = as.logical(escape)) %>%
ggplot(aes(x = n_inquiry, y = n_response)) + facet_wrap(~label) +
geom_point(size = 2, alpha = 0.5) + geom_smooth(method = "glm.nb") +
xlab("inquiry call count") + ylab("response call count") + theme(legend.position = "none")
# plot with kinship
d3 <- d %>%
filter(!is.na(kinship)) %>%
mutate(escape = as.logical(escape)) %>%
mutate(kin = if_else(kinship > 0.1, "kin", "nonkin")) %>%
mutate(kin = factor(kin, levels = c("nonkin", "kin")))
d2$kin <- "all pairs"
shared_cols <- intersect(names(d3), names(d2))
dall <- rbind(d3[, shared_cols], d2[, shared_cols])
dall$kin <- factor(dall$kin, levels = c("all pairs", "nonkin", "kin"))
t2 <- ggplot(dall, aes(x = n_inquiry, y = (n_response), color = kin,
group = kin)) + facet_wrap(~kin) + geom_point(size = 2, alpha = 0.5) +
geom_smooth(method = "glm.nb", color = mako(10)[9]) + xlab("inquiry call count") +
ylab("response call count") + scale_color_viridis_d(alpha = 0.8,
begin = 0.4, end = 0.9, direction = -1, option = "G") + theme(legend.position = "none") +
theme_bw(base_size = 20)
t2
# save as PDF
ggsave("./output/fig_.pdf", plot = t2, width = 14, height = 7, units = "in",
dpi = 300)
# plot all bats log transform x-axis for better visualization
(t1 <- d %>%
mutate(label = "pairs with known kinship") %>%
filter(!is.na(kinship)) %>%
mutate(escape = as.logical(escape)) %>%
ggplot(aes(x = log_inquiry, y = n_response)) + facet_wrap(~label) +
geom_point(size = 2, alpha = 0.5, aes(shape = escape), color = mako(10)[9]) +
geom_smooth(method = "glm.nb", color = mako(10)[9]) + xlab("") +
ylab("response call count") + coord_cartesian(ylim = c(0, 800)) +
theme_bw(base_size = 20) + theme(legend.position = "none"))
# plot with kinship
(t2 <- d %>%
filter(!is.na(kinship)) %>%
mutate(escape = as.logical(escape)) %>%
mutate(kin = if_else(kinship > 0.1, "kin", "nonkin")) %>%
mutate(kin = factor(kin, levels = c("nonkin", "kin"))) %>%
ggplot(aes(x = log_inquiry, y = (n_response), color = kin, group = kin)) +
facet_wrap(~kin) + geom_point(size = 2, alpha = 0.5, aes(shape = escape)) +
geom_smooth(method = "glm.nb", color = mako(10)[9]) + xlab("log of inquiry call count ") +
ylab("response call count") + scale_color_viridis_d(alpha = 0.8,
begin = 0.4, end = 0.9, direction = -1, option = "G") + coord_cartesian(ylim = c(0,
800)) + theme_bw(base_size = 20) + theme(legend.position = "none") +
theme(axis.text.y = element_blank(), axis.title.y = element_blank(),
legend.position = "none"))
# save as PDF ggsave( './output/fig_2_v2.pdf', plot =
# inquiry.response, width = 10, height = 5, units = 'in', dpi =
# 300)
# for plots, get residuals from negative binomial mixed effect
# model predicting number of responses by number of inquiry
# calls
fit <- glmmTMB(n_response ~ log_inquiry + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group), data = d,
ziformula = ~1, family = nbinom2)
# get difference between observed response count and predicted
# response count add this response score to data for plotting
# this value represents response calling more than expected from
# inquiry calls
d$resid_response <- d$n_response - predict(fit, type = "response",
newdata = d, allow.new.levels = T)
# plot number of inquiry calls by kinship with partner
d %>%
filter(!is.na(kinship)) %>%
mutate(kin = as.character(kinship)) %>%
ggplot(aes(x = kin, y = n_inquiry, group = kin)) + geom_violin(fill = NA,
width = 1) + scale_color_viridis_d(alpha = 0.8, begin = 0.4, end = 0.9,
direction = -1, option = "G") + geom_jitter(size = 2, width = 0.1,
height = 0, aes(color = udyad.sex)) + xlab("kinship") + ylab("number of inquiry calls") +
theme(legend.position = "top", legend.title = element_blank())
# plot number of response calls by kinship with partner
d %>%
filter(!is.na(kinship)) %>%
mutate(kin = as.character(kinship)) %>%
ggplot(aes(x = kin, y = n_response, group = kin)) + geom_violin(fill = NA,
width = 1) + geom_jitter(size = 2, width = 0.1, height = 0, aes(color = udyad.sex)) +
scale_color_viridis_d(alpha = 0.8, begin = 0.4, end = 0.9, direction = -1,
option = "G") + xlab("kinship") + ylab("number of response calls") +
theme(legend.position = "top", legend.title = element_blank())
# plot residual response call variation
d %>%
filter(!is.na(kinship)) %>%
mutate(kin = as.character(kinship)) %>%
ggplot(aes(x = kin, y = resid_response, group = kin)) + geom_violin(fill = NA,
width = 1) + scale_color_viridis_d(alpha = 0.8, begin = 0.4, end = 0.9,
direction = -1, option = "G") + geom_jitter(size = 2, width = 0.1,
height = 0, aes(color = udyad.sex)) + xlab("kinship") + ylab("residual response call variation") +
theme(legend.position = "top", legend.title = element_blank())
# plot number of inquiry calls by association and kinship (by
# dyad type)
t1 <- d %>%
filter(!is.na(kinship)) %>%
mutate(kin = ifelse(kinship > 0.1, "kin", "nonkin")) %>%
ggplot(aes(x = sri, y = log_inquiry)) + facet_wrap(~dyad.sex,
scales = "free_y") + geom_point(aes(color = kin), size = 2, alpha = 0.6) +
geom_smooth(method = "lm", color = mako(10)[9]) + xlab("association rate (simple ratio index)") +
ylab("inquiry call count (log transformed)") + scale_color_viridis_d(alpha = 0.8,
begin = 0.4, end = 0.9, direction = -1, option = "G") + theme_bw() +
theme(legend.position = "none", legend.title = element_blank())
# plot number of response calls by association and kinship (by
# dyad type)
t2 <- d %>%
filter(!is.na(kinship)) %>%
mutate(kin = ifelse(kinship > 0.1, "kin", "nonkin")) %>%
ggplot(aes(x = sri, y = log(n_response + 1))) + facet_wrap(~dyad.sex,
scales = "free_y") + geom_point(aes(color = kin), size = 2, alpha = 0.6) +
geom_smooth(method = "lm", color = mako(10)[9]) + xlab("association rate (simple ratio index)") +
ylab("response call count (log transformed)") + scale_color_viridis_d(alpha = 0.8,
begin = 0.4, end = 0.9, direction = -1, option = "G") + theme_bw() +
theme(legend.position = "none", legend.title = element_blank())
# plot residual response calling by association (by dyad type)-
# controlling for inquiry calls
t3 <- d %>%
filter(!is.na(kinship)) %>%
mutate(kin = ifelse(kinship > 0.1, "kin", "nonkin")) %>%
ggplot(aes(x = sri, y = resid_response)) + facet_wrap(~dyad.sex,
scales = "free_y") + geom_point(aes(color = kin), size = 2, alpha = 0.6) +
geom_smooth(method = "lm", color = mako(10)[9]) + xlab("association rate (simple ratio index)") +
ylab("residual response call variation") + scale_color_viridis_d(alpha = 0.8,
begin = 0.4, end = 0.9, direction = -1, option = "G") + theme_bw() +
theme(legend.position = "none", legend.title = element_blank())
# combine plots
(by.sex.plot <- t1 + t2 + t3 + plot_layout(ncol = 1) + plot_annotation(tag_levels = "A"))
# scale all predictors
d$n_inquiry2 <- as.numeric(scale(d$n_inquiry))
d$n_response2 <- as.numeric(scale(d$n_response))
d$kinship2 <- as.numeric(scale(d$kinship))
d$sri2 <- as.numeric(scale(d$sri))
# poisson model is overdispersed
t <- glmer(n_response ~ kinship2*sri2 + n_inquiry2 + offset(log(flight_time)) + (1|bat_flying) + (1|bat_roosting) + (1|group), family= poisson, data = d)
check_overdispersion(t)# Overdispersion test
dispersion ratio = 68.790
Pearson's Chi-Squared = 21462.465
p-value = < 0.001# Overdispersion test
dispersion ratio = 0.330
Pearson's Chi-Squared = 102.647
p-value = 1# Check for zero-inflation
Observed zeros: 133
Predicted zeros: 116
Ratio: 0.87[1] "AIC"[1] 2698.025[1] "BIC"[1] 2731.94[1] "AIC"[1] 2646.976[1] "BIC"[1] 2684.66 Family: nbinom2 ( log )
Formula:
n_response ~ kinship2 * sri2 + n_inquiry2 + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Zero inflation: ~1
Data: d
AIC BIC logLik deviance df.resid
2647.0 2684.7 -1313.5 2627.0 310
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 5.031e-08 0.0002243
bat_roosting (Intercept) 2.374e-01 0.4872649
group (Intercept) 1.792e-01 0.4233283
Number of obs: 320, groups: bat_flying, 50; bat_roosting, 50; group, 21
Dispersion parameter for nbinom2 family (): 0.717
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.96566 0.18854 -5.122 3.03e-07 ***
kinship2 -0.03469 0.12451 -0.279 0.7805
sri2 0.06313 0.15705 0.402 0.6877
n_inquiry2 0.24756 0.11835 2.092 0.0365 *
kinship2:sri2 0.28108 0.15151 1.855 0.0636 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4359 0.1241 -3.512 0.000444 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1[1] "AIC"[1] 2648.433[1] "BIC"[1] 2682.348 Family: nbinom2 ( log )
Formula:
n_response ~ kinship2 + sri2 + n_inquiry2 + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Zero inflation: ~1
Data: d
AIC BIC logLik deviance df.resid
2648.4 2682.3 -1315.2 2630.4 311
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 0.0005841 0.02417
bat_roosting (Intercept) 0.2157783 0.46452
group (Intercept) 0.1779363 0.42183
Number of obs: 320, groups: bat_flying, 50; bat_roosting, 50; group, 21
Dispersion parameter for nbinom2 family (): 0.691
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.86322 0.18200 -4.743 2.11e-06 ***
kinship2 0.05982 0.12065 0.496 0.6200
sri2 -0.06581 0.15247 -0.432 0.6660
n_inquiry2 0.28596 0.12515 2.285 0.0223 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4445 0.1254 -3.544 0.000394 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# get predicted values
d$n_response_predicted <- predict(fit, type= "response", newdata= d, allow.new.levels = T)
# plot model performance
d %>%
filter(!is.na(kinship)) %>%
mutate(kinship= kinship>0) %>%
ggplot(aes(x=n_response, y=n_response_predicted))+
geom_point(size=2, aes(color=kinship))+
geom_smooth(method= "lm", color = mako(10)[9])+
xlab("observed count of response calls")+
ylab("predicted count of response calls")+
scale_color_viridis_d(alpha = 0.8, begin = 0.4, end = 0.9, direction = -1, option = "G") +
theme(legend.position= "top")
Call:
lm(formula = n_response_predicted ~ n_response, data = d)
Residuals:
Min 1Q Median 3Q Max
-105.112 -20.733 -3.588 17.534 170.622
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 56.43475 2.18989 25.77 <2e-16 ***
n_response 0.20018 0.01186 16.87 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 33.65 on 318 degrees of freedom
(126 observations deleted due to missingness)
Multiple R-squared: 0.4724, Adjusted R-squared: 0.4707
F-statistic: 284.7 on 1 and 318 DF, p-value: < 2.2e-16# get model coefficients with 95% CIs
coefs <-
tidy(fit,conf.int=TRUE,exponentiate=F,effects="fixed", conf.method="profile") %>%
filter(component=="cond") %>%
dplyr::select(-effect, -component) %>%
mutate(type= "full model")
# coefs
# save model results
# write.csv(coefs, file="response_model_results.csv")
# plot model results
theme_set(theme_bw(base_size = 12))
plot1 <-
coefs %>%
filter(term != "(Intercept)") %>%
mutate(term= case_when(
term == 'sri2' ~ "association",
term == 'kinship2' ~ "kinship",
term == 'n_inquiry2' ~ "inquiry calls")) %>%
mutate(term= factor(term, levels =c("inquiry calls",
"association",
"kinship"))) %>%
mutate(term = fct_rev(term)) %>%
ggplot(aes(x=estimate, y=term))+
geom_point(size=3)+
geom_errorbarh(aes(xmin=conf.low, xmax=conf.high, height=0.2), size=1.5)+
geom_vline(xintercept = 0, linetype= "dashed")+
ylab("")+
xlab("Effect size")+
coord_cartesian(xlim= c(-1.2, 1.2))+
theme(axis.text=element_text(size=12), strip.text = element_text(size=12)) +
theme_bw(base_size = 20)
plot1
# save as PDF
# ggsave(
# "response_model_results.pdf",
# plot = plot1,
# scale = 1,
# width = 5,
# height = 2,
# units = "in",
# dpi = 300)
# check that sri and kinship do not predict response calls as single predictors
# fit zero-inflated negative binomial model with only sri
t <- glmmTMB(n_response ~ sri2 + n_inquiry2 + offset(log(flight_time)) + (1|bat_flying) + (1|bat_roosting) + (1|group),
data=d,
ziformula=~1,
family=nbinom2)
summary(t) Family: nbinom2 ( log )
Formula:
n_response ~ sri2 + n_inquiry2 + offset(log(flight_time)) + (1 |
bat_flying) + (1 | bat_roosting) + (1 | group)
Zero inflation: ~1
Data: d
AIC BIC logLik deviance df.resid
3578.6 3611.2 -1781.3 3562.6 428
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 2.041e-08 0.0001429
bat_roosting (Intercept) 6.709e-01 0.8190675
group (Intercept) 6.850e-02 0.2617247
Number of obs: 436, groups: bat_flying, 72; bat_roosting, 73; group, 23
Dispersion parameter for nbinom2 family (): 0.657
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.94322 0.16923 -5.574 2.5e-08 ***
sri2 -0.03882 0.11711 -0.332 0.74024
n_inquiry2 0.30604 0.10779 2.839 0.00452 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4642 0.1155 -4.02 5.81e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Family: nbinom2 ( log )
Formula:
n_response ~ kinship2 + n_inquiry2 + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Zero inflation: ~1
Data: d
AIC BIC logLik deviance df.resid
2646.6 2676.8 -1315.3 2630.6 312
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 1.025e-07 0.0003202
bat_roosting (Intercept) 2.188e-01 0.4677805
group (Intercept) 1.883e-01 0.4338782
Number of obs: 320, groups: bat_flying, 50; bat_roosting, 50; group, 21
Dispersion parameter for nbinom2 family (): 0.693
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.88195 0.17682 -4.988 6.1e-07 ***
kinship2 0.03995 0.10960 0.364 0.7155
n_inquiry2 0.27751 0.11647 2.383 0.0172 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4437 0.1253 -3.542 0.000397 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Overdispersion test
dispersion ratio = 6.822
Pearson's Chi-Squared = 2128.461
p-value = < 0.001# Overdispersion test
dispersion ratio = 0.813
Pearson's Chi-Squared = 252.949
p-value = 0.993[1] 3130.339[1] 3164.254# Overdispersion test
dispersion ratio = 0.813
Pearson's Chi-Squared = 252.949
p-value = 0.993[1] 3130.339[1] 3164.254 Family: nbinom2 ( log )
Formula:
n_inquiry ~ kinship2 * sri2 + n_response2 + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d
AIC BIC logLik deviance df.resid
3130.3 3164.3 -1556.2 3112.3 311
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 1.833e-01 4.282e-01
bat_roosting (Intercept) 2.049e-09 4.526e-05
group (Intercept) 3.400e-02 1.844e-01
Number of obs: 320, groups: bat_flying, 50; bat_roosting, 50; group, 21
Dispersion parameter for nbinom2 family (): 4.68
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.49241 0.08338 -17.898 <2e-16 ***
kinship2 -0.06398 0.03962 -1.615 0.106
sri2 -0.01259 0.04996 -0.252 0.801
n_response2 0.01701 0.03447 0.494 0.622
kinship2:sri2 0.02409 0.05057 0.476 0.634
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1[1] 3128.567[1] 3158.713 Family: nbinom2 ( log )
Formula:
n_inquiry ~ kinship + sri + n_response + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d
AIC BIC logLik deviance df.resid
3128.6 3158.7 -1556.3 3112.6 312
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 1.861e-01 0.4314410
bat_roosting (Intercept) 2.571e-09 0.0000507
group (Intercept) 3.441e-02 0.1855110
Number of obs: 320, groups: bat_flying, 50; bat_roosting, 50; group, 21
Dispersion parameter for nbinom2 family (): 4.69
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.3851696 0.1446519 -9.576 <2e-16 ***
kinship -0.2524912 0.1637266 -1.542 0.123
sri -0.1002280 0.1789141 -0.560 0.575
n_response 0.0001123 0.0002065 0.544 0.587
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# get predicted values
d$n_inquiry_predicted <- predict(fit2, type= "response", newdata= d, allow.new.levels = T)
# plot model performance
d %>%
filter(!is.na(kinship)) %>%
mutate(kinship= kinship>0) %>%
ggplot(aes(x=n_inquiry, y=n_inquiry_predicted))+
geom_point(size=2, aes(color=kinship))+
geom_smooth(method= "lm", color = mako(10)[9])+
xlab("observed count of inquiry calls")+
ylab("predicted count of inquiry calls")+
scale_color_viridis_d(alpha = 0.8, begin = 0.4, end = 0.9, direction = -1, option = "G") +
theme(legend.position= "top")
Call:
lm(formula = n_inquiry_predicted ~ n_inquiry, data = d)
Residuals:
Min 1Q Median 3Q Max
-92.258 -12.967 -0.546 11.381 87.999
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.40700 2.06189 9.412 <2e-16 ***
n_inquiry 0.71686 0.02554 28.063 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 19.59 on 318 degrees of freedom
(126 observations deleted due to missingness)
Multiple R-squared: 0.7124, Adjusted R-squared: 0.7114
F-statistic: 787.5 on 1 and 318 DF, p-value: < 2.2e-16 Family: nbinom2 ( log )
Formula:
n_inquiry ~ kinship + sri + n_response + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d
AIC BIC logLik deviance df.resid
3128.6 3158.7 -1556.3 3112.6 312
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 1.861e-01 0.4314410
bat_roosting (Intercept) 2.571e-09 0.0000507
group (Intercept) 3.441e-02 0.1855110
Number of obs: 320, groups: bat_flying, 50; bat_roosting, 50; group, 21
Dispersion parameter for nbinom2 family (): 4.69
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.3851696 0.1446519 -9.576 <2e-16 ***
kinship -0.2524912 0.1637266 -1.542 0.123
sri -0.1002280 0.1789141 -0.560 0.575
n_response 0.0001123 0.0002065 0.544 0.587
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coefs2 <-
tidy(fit2,conf.int=TRUE,exponentiate=F,effects="fixed", conf.method="profile") %>%
filter(component=="cond") %>%
dplyr::select(-effect, -component) %>%
mutate(type= "full model")
# coefs2
# plot model results
theme_set(theme_bw(base_size = 12))
plot2 <-
coefs2 %>%
filter(term != "(Intercept)") %>%
mutate(term= case_when(
term == 'sri2' ~ "association",
term == 'kinship2' ~ "kinship",
term == 'n_response2' ~ "response calls")) %>%
mutate(term= factor(term, levels =c("response calls",
"association",
"kinship"))) %>%
mutate(term = fct_rev(term)) %>%
ggplot(aes(x=estimate, y=term))+
geom_point(size=2)+
geom_errorbarh(aes(xmin=conf.low, xmax=conf.high, height=0.2), size=1)+
geom_vline(xintercept = 0, linetype= "dashed")+
#coord_cartesian(xlim =c(-0.2,0.2))+
coord_cartesian(xlim =c(-1,1))+
ylab("")+
xlab("coefficient")+
ggtitle("including trials without response calls")+
theme(axis.text=element_text(size=12), strip.text = element_text(size=12))
plot2
# A tibble: 3 × 5
# Groups: effect [1]
effect group variance sum prop
<chr> <chr> <dbl> <dbl> <dbl>
1 ran_pars bat_flying 0.186 0.221 0.844
2 ran_pars bat_roosting 0.00000000257 0.221 0.0000000117
3 ran_pars group 0.0344 0.221 0.156 # Overdispersion test
dispersion ratio = 4.731
Pearson's Chi-Squared = 846.888
p-value = < 0.001# Overdispersion test
dispersion ratio = 0.721
Pearson's Chi-Squared = 128.308
p-value = 0.998[1] 1841.646[1] 1870.726 Family: nbinom2 ( log )
Formula:
n_inquiry ~ kinship2 * sri2 + n_response2 + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d2
AIC BIC logLik deviance df.resid
1841.6 1870.7 -911.8 1823.6 178
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 1.799e-01 4.242e-01
bat_roosting (Intercept) 9.667e-10 3.109e-05
group (Intercept) 5.945e-02 2.438e-01
Number of obs: 187, groups: bat_flying, 46; bat_roosting, 50; group, 20
Dispersion parameter for nbinom2 family (): 5.49
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.54031 0.09929 -15.514 <2e-16 ***
kinship2 -0.05930 0.05254 -1.129 0.259
sri2 0.10011 0.06959 1.439 0.150
n_response2 0.04173 0.04091 1.020 0.308
kinship2:sri2 0.07321 0.06709 1.091 0.275
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1[1] 1840.832[1] 1866.68 Family: nbinom2 ( log )
Formula:
n_inquiry ~ kinship2 + sri2 + n_response2 + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d2
AIC BIC logLik deviance df.resid
1840.8 1866.7 -912.4 1824.8 179
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 1.886e-01 4.343e-01
bat_roosting (Intercept) 1.555e-09 3.943e-05
group (Intercept) 6.285e-02 2.507e-01
Number of obs: 187, groups: bat_flying, 46; bat_roosting, 50; group, 20
Dispersion parameter for nbinom2 family (): 5.53
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.51837 0.09916 -15.313 <2e-16 ***
kinship2 -0.04529 0.05108 -0.887 0.375
sri2 0.06428 0.06197 1.037 0.300
n_response2 0.04913 0.04045 1.215 0.224
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# get predicted values
d2$n_inquiry_predicted <- predict(fit2, type= "response", newdata= d2, allow.new.levels = T)
# plot model performance
d2 %>%
filter(!is.na(kinship)) %>%
mutate(kinship= kinship>0) %>%
ggplot(aes(x=n_inquiry, y=n_inquiry_predicted))+
geom_point(size=2, aes(color=kinship))+
geom_smooth(method= "lm")+
xlab("observed count of inquiry calls")+
ylab("predicted count of inquiry calls")+
scale_color_viridis_d(alpha = 0.8, begin = 0.4, end = 0.9, direction = -1, option = "G") +
theme(legend.position= "top")
Call:
lm(formula = n_inquiry_predicted ~ n_inquiry, data = d2)
Residuals:
Min 1Q Median 3Q Max
-47.186 -9.972 -0.709 10.055 48.381
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 18.44086 2.16188 8.53 5.16e-15 ***
n_inquiry 0.72092 0.02607 27.66 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 15.51 on 185 degrees of freedom
(68 observations deleted due to missingness)
Multiple R-squared: 0.8052, Adjusted R-squared: 0.8042
F-statistic: 764.8 on 1 and 185 DF, p-value: < 2.2e-16 Family: nbinom2 ( log )
Formula:
n_inquiry ~ kinship2 + sri2 + n_response2 + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d2
AIC BIC logLik deviance df.resid
1840.8 1866.7 -912.4 1824.8 179
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 1.886e-01 4.343e-01
bat_roosting (Intercept) 1.555e-09 3.943e-05
group (Intercept) 6.285e-02 2.507e-01
Number of obs: 187, groups: bat_flying, 46; bat_roosting, 50; group, 20
Dispersion parameter for nbinom2 family (): 5.53
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.51837 0.09916 -15.313 <2e-16 ***
kinship2 -0.04529 0.05108 -0.887 0.375
sri2 0.06428 0.06197 1.037 0.300
n_response2 0.04913 0.04045 1.215 0.224
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# A tibble: 4 × 8
term estimate std.error statistic p.value conf.low conf.high type
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
1 (Intercept) -1.52 0.0992 -15.3 6.26e-53 -1.73 -1.32 full mod…
2 kinship2 -0.0453 0.0511 -0.887 3.75e- 1 -0.147 0.0551 full mod…
3 sri2 0.0643 0.0620 1.04 3.00e- 1 -0.0560 0.188 full mod…
4 n_response2 0.0491 0.0404 1.21 2.24e- 1 -0.0301 0.130 full mod…# plot model results
theme_set(theme_bw(base_size = 12))
plot3 <-
coefs2 %>%
filter(term != "(Intercept)") %>%
mutate(term= case_when(
term == 'sri2' ~ "association",
term == 'kinship2' ~ "kinship",
term == 'n_response2' ~ "response calls")) %>%
mutate(term= factor(term, levels =c("response calls",
"association",
"kinship"))) %>%
mutate(term = fct_rev(term)) %>%
ggplot(aes(x=estimate, y=term))+
geom_point(size=2)+
geom_errorbarh(aes(xmin=conf.low, xmax=conf.high, height=0.2), size=1)+
geom_vline(xintercept = 0, linetype= "dashed")+
ylab("")+
xlab("coefficient")+
coord_cartesian(xlim =c(-1,1))+
#coord_cartesian(xlim =c(-0.2,0.2))+
ggtitle("excluding trials without response calls")+
theme(axis.text=element_text(size=12),
# remove y-axis labels for combining
axis.text.y=element_blank(),
strip.text = element_text(size=12))
plot3
Family: nbinom2 ( log )
Formula:
n_inquiry ~ sri2 + n_response2 + offset(log(flight_time)) + (1 |
bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d
AIC BIC logLik deviance df.resid
4257.0 4285.5 -2121.5 4243.0 429
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 0.230975 0.48060
bat_roosting (Intercept) 0.004481 0.06694
group (Intercept) 0.012463 0.11164
Number of obs: 436, groups: bat_flying, 72; bat_roosting, 73; group, 23
Dispersion parameter for nbinom2 family (): 5.1
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.45077 0.06911 -20.993 <2e-16 ***
sri2 -0.04204 0.03200 -1.314 0.1889
n_response2 0.05036 0.02823 1.784 0.0745 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Family: nbinom2 ( log )
Formula:
n_inquiry ~ kinship2 + n_response2 + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d
AIC BIC logLik deviance df.resid
3126.9 3153.3 -1556.4 3112.9 313
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 1.825e-01 4.272e-01
bat_roosting (Intercept) 4.167e-09 6.456e-05
group (Intercept) 3.491e-02 1.868e-01
Number of obs: 320, groups: bat_flying, 50; bat_roosting, 50; group, 21
Dispersion parameter for nbinom2 family (): 4.67
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.49114 0.08130 -18.341 <2e-16 ***
kinship2 -0.06301 0.03668 -1.718 0.0858 .
n_response2 0.02134 0.03409 0.626 0.5314
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Family: nbinom2 ( log )
Formula:
n_inquiry ~ sri2 + n_response + offset(log(flight_time)) + (1 |
bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d2
AIC BIC logLik deviance df.resid
2491.4 2516.1 -1238.7 2477.4 244
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 2.249e-01 0.4741932
bat_roosting (Intercept) 9.361e-10 0.0000306
group (Intercept) 3.377e-02 0.1837609
Number of obs: 251, groups: bat_flying, 68; bat_roosting, 71; group, 23
Dispersion parameter for nbinom2 family (): 6.55
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.4406584 0.0827752 -17.404 <2e-16 ***
sri2 -0.0127403 0.0394305 -0.323 0.7466
n_response 0.0003509 0.0001850 1.897 0.0579 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Family: nbinom2 ( log )
Formula:
n_inquiry ~ kinship2 + n_response + offset(log(flight_time)) +
(1 | bat_flying) + (1 | bat_roosting) + (1 | group)
Data: d2
AIC BIC logLik deviance df.resid
1839.9 1862.5 -913.0 1825.9 180
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
bat_flying (Intercept) 2.113e-01 0.4596623
bat_roosting (Intercept) 1.204e-09 0.0000347
group (Intercept) 5.273e-02 0.2296381
Number of obs: 187, groups: bat_flying, 46; bat_roosting, 50; group, 20
Dispersion parameter for nbinom2 family (): 5.57
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.5274704 0.1026826 -14.876 <2e-16 ***
kinship2 -0.0343894 0.0498479 -0.690 0.490
n_response 0.0002562 0.0002399 1.068 0.286
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1[1] 280
R version 4.3.1 (2023-06-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3; LAPACK version 3.9.0
locale:
[1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C
[3] LC_TIME=es_CR.UTF-8 LC_COLLATE=pt_BR.UTF-8
[5] LC_MONETARY=es_CR.UTF-8 LC_MESSAGES=pt_BR.UTF-8
[7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
time zone: America/Costa_Rica
tzcode source: system (glibc)
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] viridis_0.6.4 viridisLite_0.4.2 MASS_7.3-60
[4] performance_0.10.8 patchwork_1.1.3 boot_1.3-28.1
[7] broom.mixed_0.2.9.4 glmmTMB_1.1.8 lme4_1.1-35.1
[10] Matrix_1.6-2 mapproj_1.2.11 maps_3.4.1.1
[13] geosphere_1.5-18 lubridate_1.9.3 forcats_1.0.0
[16] stringr_1.5.0 dplyr_1.1.3 purrr_1.0.2
[19] readr_2.1.4 tidyr_1.3.0 tibble_3.2.1
[22] tidyverse_2.0.0 RColorBrewer_1.1-3 statnet_2019.6
[25] tsna_0.3.5 ergm.count_4.1.1 tergm_4.2.0
[28] networkDynamic_0.11.4 ergm_4.5.0 igraphdata_1.0.1
[31] ecodist_2.1.3 ggmap_4.0.0 ggplot2_3.4.4
[34] assortnet_0.20 asnipe_1.1.17 bipartite_2.19
[37] sna_2.7-2 network_1.18.2 statnet.common_4.9.0
[40] vegan_2.6-4 lattice_0.22-5 permute_0.9-7
[43] igraph_1.5.1
loaded via a namespace (and not attached):
[1] rstudioapi_0.15.0 jsonlite_1.8.7 magrittr_2.0.3
[4] sketchy_1.0.2 estimability_1.4.1 farver_2.1.1
[7] nloptr_2.0.3 rmarkdown_2.25 ragg_1.2.6
[10] fields_15.2 vctrs_0.6.4 memoise_2.0.1
[13] minqa_1.2.6 htmltools_0.5.7 broom_1.0.5
[16] parallelly_1.36.0 htmlwidgets_1.6.2 plyr_1.8.9
[19] emmeans_1.8.9 cachem_1.0.8 TMB_1.9.9
[22] lpSolveAPI_5.5.2.0-17.11 lifecycle_1.0.4 pkgconfig_2.0.3
[25] R6_2.5.1 fastmap_1.1.1 rbibutils_2.2.16
[28] future_1.33.0 digest_0.6.33 numDeriv_2016.8-1.1
[31] colorspace_2.1-0 ergm.multi_0.2.0 furrr_0.3.1
[34] textshaping_0.3.7 labeling_0.4.3 fansi_1.0.5
[37] timechange_0.2.0 httr_1.4.7 mgcv_1.9-0
[40] compiler_4.3.1 remotes_2.4.2.1 withr_2.5.2
[43] backports_1.4.1 tools_4.3.1 trust_0.1-8
[46] rle_0.9.2 glue_1.6.2 networkLite_1.0.5
[49] nlme_3.1-163 grid_4.3.1 cluster_2.1.4
[52] generics_0.1.3 gtable_0.3.4 tzdb_0.4.0
[55] hms_1.1.3 sp_2.1-1 utf8_1.2.4
[58] pillar_1.9.0 spam_2.10-0 robustbase_0.99-1
[61] splines_4.3.1 tidyselect_1.2.0 knitr_1.45
[64] gridExtra_2.3 xfun_0.41 DEoptimR_1.1-3
[67] stringi_1.7.12 yaml_2.3.7 xaringanExtra_0.7.0
[70] evaluate_0.23 codetools_0.2-19 cli_3.6.1
[73] systemfonts_1.0.5 xtable_1.8-4 Rdpack_2.6
[76] munsell_0.5.0 Rcpp_1.0.11 globals_0.16.2
[79] coda_0.19-4 png_0.1-8 parallel_4.3.1
[82] dotCall64_1.1-0 jpeg_0.1-10 bitops_1.0-7
[85] listenv_0.9.0 mvtnorm_1.2-3 scales_1.2.1
[88] insight_0.19.6 packrat_0.9.2 crayon_1.5.2
[91] rlang_1.1.2 formatR_1.14