R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

Preliminaries

Load libraries

library(knitr)
library(readxl)
library(dplyr)
library(broom)
library(ggplot2)
library(forcats)
library(ggpubr)
library(grid)  # for text_grob()

Set 1: Within-individual comparison with separation near pond

Data Processing

dataSET1 <- read_excel("CallContext_DataSubmission2.xlsx", 
    sheet = "Within-Male (Natural)", range = "A1:O54")

Individual Reaction Norms - Separated Recordings Occured in Natural Environment

DurNormSet1 <- ggplot(dataSET1, aes(x = fct_rev(Treatment), y = `TempCor Call duration (ms)`, group = `Male ID`)) +
  geom_line(alpha = 0.5, color = "dodgerblue2", linewidth=1) +  # one line per individual
  theme_classic() +
  geom_point(data = dataSET1, color = "dodgerblue2", size = 2) +
  labs(x = "Amplexus Status", y = "Call Duration (ms)") +
  scale_x_discrete(labels = c("POST-AMPLEX" = "Separated", "PRE-AMPLEX" = "Baseline"))

RateNormSet1 <- ggplot(dataSET1, aes(x = fct_rev(Treatment), y = `TempCor Call rate (call/min)`, group = `Male ID`)) +
  geom_line(alpha = 0.5, color = "darkolivegreen3", linewidth=1) +  # one line per individual
  theme_classic() +
  geom_point(data = dataSET1, color = "darkolivegreen3", size = 2) +
  labs(x = "Amplexus Status", y = "Call Rate (calls/min)") +
  scale_x_discrete(labels = c("POST-AMPLEX" = "Separated", "PRE-AMPLEX" = "Baseline"))

DCNormSet1 <- ggplot(dataSET1, aes(x = fct_rev(Treatment), y = `TempCor Duty cycle`, group = `Male ID`)) +
  geom_line(alpha = 0.5, color = "coral2", linewidth=1) +  # one line per individual
  theme_classic() +
  geom_point(data = dataSET1, color = "coral2", size = 2) +
  labs(x = "Amplexus Status", y = "Duty Cycle") +
  scale_x_discrete(labels = c("POST-AMPLEX" = "Separated", "PRE-AMPLEX" = "Baseline"))


ggarrange(DurNormSet1, RateNormSet1, DCNormSet1, 
          labels = c("A", "B", "C"),
          ncol = 3, nrow = 1,
          common.legend = TRUE,
          legend = "bottom")

Individual Box Plots - Separated Recordings Occured in Natural Environment

This prep allows me to add in the significance asteriks marks on the plots

stat.test.DurSet1 <- data.frame(
  group1 = "POST-AMPLEX",
  group2 = "PRE-AMPLEX",
  y.position = 1900,  # adjust if needed
  label = "***",
  x = 1.5,
  xend = 1.5
)

stat.test.RateSet1 <- data.frame(
  group1 = "POST-AMPLEX",
  group2 = "PRE-AMPLEX",
  y.position = 48,  # adjust if needed
  label = "***",
  x = 1.5,
  xend = 1.5
)

Generating box plots

BoxDurSet1 <- ggplot(dataSET1, aes(x=fct_rev(Treatment), y=`TempCor Call duration (ms)`) ) +
  geom_boxplot(fill=c("dodgerblue", "dodgerblue4"), 
               color=c("dodgerblue", "dodgerblue4"), 
               alpha=0.7,
               outlier.shape = NA) +
              coord_cartesian(ylim = c(0, 2000)) +
    geom_jitter(width = 0.15,  # horizontal jitter
              size = 1.2,
              color = "dodgerblue4",
              alpha=0.8) +
  scale_x_discrete(labels = c("POST-AMPLEX" = "Separated", "PRE-AMPLEX" = "Baseline")) +
  theme_classic() +
  xlab("Calling Context") + ylab("Call Duration (ms)") +
  theme(axis.text=element_text(size=10), axis.title=element_text(size=12)) +
    stat_pvalue_manual(stat.test.DurSet1, 
                     xmin = "group1", 
                     xend = "group2",
                     y.position = "y.position",
                     label = "label",
                     tip.length = 0.03)

BoxRateSet1 <- ggplot(dataSET1, aes(x=fct_rev(Treatment), y=`TempCor Call rate (call/min)`) ) +
  geom_boxplot(fill=c("darkolivegreen3", "darkolivegreen"), 
               color=c("darkolivegreen3", "darkolivegreen"), 
               alpha=0.7,
               outlier.shape = NA) +
              coord_cartesian(ylim = c(0, 50)) +
    geom_jitter(width = 0.15,  # horizontal jitter
              size = 1.2,
              color = "darkolivegreen",
              alpha=0.8) +
  scale_x_discrete(labels = c("POST-AMPLEX" = "Separated", "PRE-AMPLEX" = "Baseline")) +
  theme_classic() +
  xlab("Calling Context") + ylab("Call Rate (calls/min)") +
  theme(axis.text=element_text(size=10), axis.title=element_text(size=12)) +
    stat_pvalue_manual(stat.test.RateSet1, 
                     xmin = "group1", 
                     xend = "group2",
                     y.position = "y.position",
                     label = "label",
                     tip.length = 0.03)

BoxDCSet1 <- ggplot(dataSET1, aes(x=fct_rev(Treatment), y=`TempCor Duty cycle`) ) +
  geom_boxplot(fill=c("coral", "coral4"), 
               color=c("coral", "coral4"), 
               alpha=0.7,
               outlier.shape = NA) +
              coord_cartesian(ylim = c(0, 0.3)) +
    geom_jitter(width = 0.15,  # horizontal jitter
              size = 1.2,
              color = "coral4",
              alpha=0.8) +
  scale_x_discrete(labels = c("POST-AMPLEX" = "Separated", "PRE-AMPLEX" = "Baseline")) +
  theme_classic() +
  xlab("Calling Context") + ylab("Duty Cycle") +
  theme(axis.text=element_text(size=10), axis.title=element_text(size=12))

Combining the box plots of the call parameters

withinfieldfig <- ggarrange(BoxDurSet1, BoxRateSet1, BoxDCSet1, 
            labels = c("a", "b", "c"),
            label.x = c(0.6, 0.55, 0.55),
            label.y = 1.02,
            align = "none",
          ncol = 3, nrow = 1,
          common.legend = TRUE,
          legend = "bottom")

withinfieldfig <- annotate_figure(
  withinfieldfig,
  left = text_grob("Within-Male Field", rot = 90, size = 16, face = "bold")
)

Within-individual comparison with separation in lab

Data Processing

dataSET2 <- read_excel("CallContext_DataSubmission2.xlsx", 
    sheet = "Within-Male (Lab)", range = "A1:O45")

Individual Reaction Norms - Separated Recordings Occured in Lab Environment

DurNormSet2 <- ggplot(dataSET2, aes(x = fct_rev(Treatment), y = `TempCor Call duration (ms)`, group = `Male ID`)) +
  geom_line(alpha = 0.5, color = "dodgerblue2", linewidth=1) +  # one line per individual
  theme_classic() +
  geom_point(data = dataSET2, color = "dodgerblue2", size = 2) +
  labs(x = "Amplexus Status", y = "Call Duration (ms)") +
  scale_x_discrete(labels = c("POST-AMPLEX" = "Separated", "PRE-AMPLEX" = "Baseline"))

RateNormSet2 <- ggplot(dataSET2, aes(x = fct_rev(Treatment), y = `TempCor Call rate (call/min)`, group = `Male ID`)) +
  geom_line(alpha = 0.5, color = "darkolivegreen3", linewidth=1) +  # one line per individual
  theme_classic() +
  geom_point(data = dataSET2, color = "darkolivegreen3", size = 2) +
  labs(x = "Amplexus Status", y = "Call Rate (calls/min)") +
  scale_x_discrete(labels = c("POST-AMPLEX" = "Separated", "PRE-AMPLEX" = "Baseline"))

DCNormSet2 <- ggplot(dataSET2, aes(x = fct_rev(Treatment), y = `TempCor Duty cycle`, group = `Male ID`)) +
  geom_line(alpha = 0.5, color = "coral2", linewidth=1) +  # one line per individual
  theme_classic() +
  geom_point(data = dataSET2, color = "coral2", size = 2) +
  labs(x = "Amplexus Status", y = "Duty Cycle") +
  scale_x_discrete(labels = c("POST-AMPLEX" = "Separated", "PRE-AMPLEX" = "Baseline"))


ggarrange(DurNormSet2, RateNormSet2, DCNormSet2, 
          labels = c("a", "b", "c"),
          ncol = 3, nrow = 1,
          common.legend = TRUE,
          legend = "bottom")

Individual Box Plots - Separated Recordings Occured in Lab Environment

To add the asteriks needed to denote significance

stat.test.DurSet2 <- data.frame(
  group1 = "POST-AMPLEX",
  group2 = "PRE-AMPLEX",
  y.position = 1900,  # adjust if needed
  label = "***",
  x = 1.5,
  xend = 1.5
)

stat.test.RateSet2 <- data.frame(
  group1 = "POST-AMPLEX",
  group2 = "PRE-AMPLEX",
  y.position = 48,  # adjust if needed
  label = "***",
  x = 1.5,
  xend = 1.5
)

Generating the box plots

BoxDurSet2 <- ggplot(dataSET2, aes(x=fct_rev(Treatment), y=`TempCor Call duration (ms)`) ) +
  geom_boxplot(fill=c("dodgerblue", "dodgerblue4"), 
               color=c("dodgerblue", "dodgerblue4"), 
               alpha=0.7,
               outlier.shape = NA) +
              coord_cartesian(ylim = c(0, 2000)) +
    geom_jitter(width = 0.15,  # horizontal jitter
              size = 1.2,
              color = "dodgerblue4",
              alpha=0.8) +
  scale_x_discrete(labels = c("POST-AMPLEX" = "Separated", "PRE-AMPLEX" = "Baseline")) +
  theme_classic() +
  xlab("Calling Context") + ylab("Call Duration (ms)") +
  theme(axis.text=element_text(size=10), axis.title=element_text(size=12)) +
    stat_pvalue_manual(stat.test.DurSet2, 
                     xmin = "group1", 
                     xend = "group2",
                     y.position = "y.position",
                     label = "label",
                     tip.length = 0.03)

BoxRateSet2 <- ggplot(dataSET2, aes(x=fct_rev(Treatment), y=`TempCor Call rate (call/min)`) ) +
  geom_boxplot(fill=c("darkolivegreen3", "darkolivegreen"), 
               color=c("darkolivegreen3", "darkolivegreen"), 
               alpha=0.7,
               outlier.shape = NA) +
              coord_cartesian(ylim = c(0, 50)) +
    geom_jitter(width = 0.15,  # horizontal jitter
              size = 1.2,
              color = "darkolivegreen",
              alpha=0.8) +
  scale_x_discrete(labels = c("POST-AMPLEX" = "Separated", "PRE-AMPLEX" = "Baseline")) +
  theme_classic() +
  xlab("Calling Context") + ylab("Call Rate (calls/min)") +
  theme(axis.text=element_text(size=10), axis.title=element_text(size=12)) +
    stat_pvalue_manual(stat.test.RateSet2, 
                     xmin = "group1", 
                     xend = "group2",
                     y.position = "y.position",
                     label = "label",
                     tip.length = 0.03)

BoxDCSet2 <- ggplot(dataSET2, aes(x=fct_rev(Treatment), y=`TempCor Duty cycle`) ) +
  geom_boxplot(fill=c("coral", "coral4"), 
               color=c("coral", "coral4"), 
               alpha=0.7,
               outlier.shape = NA) +
              coord_cartesian(ylim = c(0, 0.3)) +
    geom_jitter(width = 0.15,  # horizontal jitter
              size = 1.2,
              color = "coral4",
              alpha=0.8) +
  scale_x_discrete(labels = c("POST-AMPLEX" = "Separated", "PRE-AMPLEX" = "Baseline")) +
  theme_classic() +
  xlab("Calling Context") + ylab("Duty Cycle") +
  theme(axis.text=element_text(size=10), axis.title=element_text(size=12))

Combining the box plots

withinmalelabfig <- ggarrange(BoxDurSet2, BoxRateSet2, BoxDCSet2, 
            labels = c("a", "b", "c"),
            label.x = c(0.6, 0.55, 0.55),
            label.y = 1.02,
            align = "none",
          ncol = 3, nrow = 1,
          common.legend = TRUE,
          legend = "bottom")

withinmalelabfig <- annotate_figure(
  withinmalelabfig,
  left = text_grob("Within-Male Lab", rot = 90, size = 16, face = "bold")
)

Between-Individual Comparison Across Sampling Years

Data Processing

dataYEAR <- read_excel("CallContext_DataSubmission2.xlsx", 
    sheet = "Between-Males", range = "A1:G238")

year_B2018 <- subset(dataYEAR, Year == "B-2018")
year_B2021 <- subset(dataYEAR, Year == "B-2021")
year_B2017 <- subset(dataYEAR, Year == "B-2017")
year_C2019 <- subset(dataYEAR, Year =="C-2019")
year_C2021 <- subset(dataYEAR, Year =="C-2021")
stat.test.DurYearly1 <- data.frame(
  group1 = "B-2017",
  group2 = "B-2021",
  y.position = 1700,  # adjust if needed
  label = c("NS"),
  x = 1.5,
  xend = 1.5
)

stat.test.DurYearly2 <- data.frame(
  group1 = "C-2019",
  group2 = "C-2021",
  y.position = 1700,  # adjust if needed
  label = c("NS"),
  x = 1.5,
  xend = 1.5
)

stat.test.DurYearly3 <- data.frame(
  group1 = "B-2017",
  group2 = "C-2021",
  y.position = 1900,  # adjust if needed
  label = c("***"),
  x = 1.5,
  xend = 1.5
)

stat.test.RateYearly1 <- data.frame(
  group1 = "B-2017",
  group2 = "B-2021",
  y.position = 43,  # adjust if needed
  label = c("NS"),
  x = 1.5,
  xend = 1.5
)

stat.test.RateYearly2 <- data.frame(
  group1 = "C-2019",
  group2 = "C-2021",
  y.position = 43,  # adjust if needed
  label = c("NS"),
  x = 1.5,
  xend = 1.5
)

stat.test.RateYearly3 <- data.frame(
  group1 = "B-2017",
  group2 = "C-2021",
  y.position = 48,  # adjust if needed
  label = c("***"),
  x = 1.5,
  xend = 1.5
)

Plotting

Generating the box plots

yearDUR <- ggplot(dataYEAR, aes(x=Year, y=`Call duration (ms)`) ) +
  geom_boxplot(fill=c("dodgerblue", "dodgerblue", "dodgerblue", "dodgerblue4", "dodgerblue4"), 
               color=c("dodgerblue", "dodgerblue", "dodgerblue", "dodgerblue4", "dodgerblue4"), 
               alpha=0.7,
               outlier.shape = NA) +
              coord_cartesian(ylim = c(0, 2000)) +
  stat_pvalue_manual(stat.test.DurYearly1,
                     label = "label", 
                     y.position = "y.position",
                     tip.length = 0.01) +
  stat_pvalue_manual(stat.test.DurYearly2,
                     label = "label", 
                     y.position = "y.position",
                     tip.length = 0.01) +
  stat_pvalue_manual(stat.test.DurYearly3,
                     label = "label", 
                     y.position = "y.position",
                     tip.length = 0.01) +
  scale_x_discrete(labels = c("B-2017" = "2017    ",
                              "B-2018" = "2018    ",
                              "B-2021" = "2021    ",
                              "C-2019" = "2019    ",
                              "C-2021" = "2021    ")) +
  theme_classic() +
  xlab("Sampling Year") + ylab("Call Duration (ms)") +
  theme(axis.text=element_text(size=9), axis.title=element_text(size=12)) 

yearRATE <- ggplot(dataYEAR, aes(x=Year, y=`Call rate (calls/min)`) ) +
  geom_boxplot(fill=c("darkolivegreen3", "darkolivegreen3", "darkolivegreen3", "darkolivegreen", "darkolivegreen"), 
               color=c("darkolivegreen3", "darkolivegreen3", "darkolivegreen3", "darkolivegreen", "darkolivegreen"), 
               alpha=0.7,
               outlier.shape = NA) +
              coord_cartesian(ylim = c(0, 50)) +
    stat_pvalue_manual(stat.test.RateYearly1,
                     label = "label", 
                     y.position = "y.position",
                     tip.length = 0.01) +
  stat_pvalue_manual(stat.test.RateYearly2,
                     label = "label", 
                     y.position = "y.position",
                     tip.length = 0.01) +
  stat_pvalue_manual(stat.test.RateYearly3,
                     label = "label", 
                     y.position = "y.position",
                     tip.length = 0.01) +
  scale_x_discrete(labels = c("B-2017" = "2017     ",
                              "B-2018" = "2018     ",
                              "B-2021" = "2021     ",
                              "C-2019" = "2019     ",
                              "C-2021" = "2021     ")) +
  theme_classic() +
  xlab("Sampling Year") + ylab("Call Rate (calls/min)") +
  theme(axis.text=element_text(size=9), axis.title=element_text(size=12))

yearDC <- ggplot(dataYEAR, aes(x=Year, y=`Duty cycle`) ) +
  geom_boxplot(fill=c("coral2", "coral2", "coral2", "coral4", "coral4"), 
               color=c("coral2", "coral2", "coral2", "coral4", "coral4"), 
               alpha=0.7,
               outlier.shape = NA) +
              coord_cartesian(ylim = c(0, 0.3)) + scale_y_continuous(breaks = seq(0, 0.3, by = 0.1)) +
  scale_x_discrete(labels = c("B-2017" = "2017     ",
                              "B-2018" = "2018     ",
                              "B-2021" = "2021     ",
                              "C-2019" = "2019     ",
                              "C-2021" = "2021     ")) +
  theme_classic() +
  xlab("Sampling Year") + ylab("Duty Cycle") +
  theme(axis.text=element_text(size=9), axis.title=element_text(size=12))

Combines the plots for different call measures

yearlyfig <- ggarrange(yearDUR, yearRATE, yearDC, 
          labels = c("", "", ""),
          ncol = 3, nrow = 1,
          common.legend = TRUE,
          legend = "bottom")
yearlyfig <- annotate_figure(
  yearlyfig,
  left = text_grob("Between-Male Field", rot = 90, size = 16, face = "bold")
)

FINAL FIGURE 1 - Big combination of all box plots with experimental data

row1D <- annotate_figure(
  ggarrange(BoxDurSet1, BoxRateSet1, BoxDCSet1, 
            ncol = 3, 
            labels = c("a", "b", "c"),
            label.x = c(0.55, 0.55, 0.55),
            label.y = 1.02,
            align = "none"),
  left = text_grob("Within-Male Field", rot = 90, face = "bold", size = 14)
)


row2D <- annotate_figure(
  ggarrange(BoxDurSet2, BoxRateSet2, BoxDCSet2, 
            ncol = 3, labels = c("", "", "")),
  left = text_grob("Within-Male Lab", rot = 90, face = "bold", size = 14)
)

row3D <- annotate_figure(
  ggarrange(yearDUR, yearRATE, yearDC, 
            ncol = 3, labels = c("", "", "")),
  left = text_grob("Between-Male Field", rot = 90, face = "bold", size = 14)
)

data_combined <- ggarrange(row1D, row2D, row3D, 
                            ncol = 1,
                            heights = c(1, 1, 1))


print(data_combined)

ggsave("data_combined.png", data_combined, width = 9, height = 9)

CONTROL for handling time

Data Processing

dataHANDLING <- read_excel("CallContext_DataSubmission2.xlsx", 
    sheet = "Handling Control", range = "A1:M25")

Individual Reaction Norms - Handling Time

DurNormCONT <- ggplot(dataHANDLING, aes(x = fct_rev(Treatment), y = `TempCor Call duration (ms)`, group = `Male ID`)) +
  geom_line(alpha = 0.5, color = "grey20", linewidth=1) +  # one line per individual
  theme_classic() +
  geom_point(data = dataHANDLING, color = "grey20", size = 2) +
  labs(x = "Relative Handling", y = "Call Duration (ms)") +
  scale_x_discrete(labels = c("POST-TRANSFER" = "Control", "PRE-TRANSFER" = "Baseline"))

RateNormCONT <- ggplot(dataHANDLING, aes(x = fct_rev(Treatment), y = `TempCor Call rate (call/min)`, group = `Male ID`)) +
  geom_line(alpha = 0.5, color = "grey20", linewidth=1) +  # one line per individual
  theme_classic() +
  geom_point(data = dataHANDLING, color = "grey20", size = 2) +
  labs(x = "Relative Handling", y = "Call Rate (calls/min)") +
  scale_x_discrete(labels = c("POST-TRANSFER" = "Handling", "PRE-TRANSFER" = "Baseline"))

DCNormCONT <- ggplot(dataHANDLING, aes(x = fct_rev(Treatment), y = `TempCor Duty cycle`, group = `Male ID`)) +
  geom_line(alpha = 0.5, color = "grey20", linewidth=1) +  # one line per individual
  theme_classic() +
  geom_point(data = dataHANDLING, color = "grey20", size = 2) +
  labs(x = "Relative Handling", y = "Duty Cycle") +
  scale_x_discrete(labels = c("POST-TRANSFER" = "Handling", "PRE-TRANSFER" = "Baseline"))

Combining the reaction norms for handling control

ggarrange(DurNormCONT, RateNormCONT, DCNormCONT,
          ncol = 3, nrow = 1,
          common.legend = TRUE,
          legend = "bottom")

Individual Box Plots - Handling time

stat.test.Ratehandle <- data.frame(
  group1 = "PRE-TRANSFER",
  group2 = "POST-TRANSFER",
  y.position = 48,  # adjust if needed
  label = "*",
  x = 1.5,
  xend = 1.5
)

Generating the box plots

BoxDurCONT <- ggplot(dataHANDLING, aes(x=fct_rev(Treatment), y=`TempCor Call duration (ms)`)) +
  geom_boxplot(fill="grey45", 
               color="dodgerblue4", 
               alpha=0.7,
               outlier.shape = NA) +
              coord_cartesian(ylim = c(0, 2000)) +
    geom_jitter(width = 0.15,  # horizontal jitter
              size = 1.2,
              color = "dodgerblue4",
              alpha=0.6) +
  scale_x_discrete(labels = c("POST-TRANSFER" = "Handled", "PRE-TRANSFER" = "Baseline")) +
  theme_classic() +
  xlab("Relative Handling") + ylab("Call Duration (ms)") +
  theme(axis.text=element_text(size=10), axis.title=element_text(size=12)) 

BoxRateCONT <- ggplot(dataHANDLING, aes(x=fct_rev(Treatment), y=`TempCor Call rate (call/min)`) ) +
  geom_boxplot(fill="grey45", 
               color="darkolivegreen", 
               alpha=0.7,
               outlier.shape = NA) +
              coord_cartesian(ylim = c(0, 50)) +
    geom_jitter(width = 0.15,  # horizontal jitter
              size = 1.2,
              color = "darkolivegreen",
              alpha=0.6) +
  scale_x_discrete(labels = c("POST-TRANSFER" = "Handled", "PRE-TRANSFER" = "Baseline")) +
  theme_classic() +
  xlab("Relative Handling") + ylab("Call Rate (calls/min)") +
  theme(axis.text=element_text(size=10), axis.title=element_text(size=12)) +
    stat_pvalue_manual(stat.test.Ratehandle, 
                     xmin = "group1", 
                     xend = "group2",
                     y.position = "y.position",
                     label = "label",
                     tip.length = 0.03)

BoxDCCONT <- ggplot(dataHANDLING, aes(x=fct_rev(Treatment), y=`TempCor Duty cycle`) ) +
  geom_boxplot(fill="grey45", 
               color="coral4", 
               alpha=0.7,
               outlier.shape = NA) +
              coord_cartesian(ylim = c(0, 0.3)) +
    geom_jitter(width = 0.15,  # horizontal jitter
              size = 1.2,
              color = "coral4",
              alpha=0.6) +
  scale_x_discrete(labels = c("POST-TRANSFER" = "Handled", "PRE-TRANSFER" = "Baseline")) +
  theme_classic() +
  xlab("Relative Handling") + ylab("Duty Cycle") +
  theme(axis.text=element_text(size=10), axis.title=element_text(size=12)) 

Combining the box plots

handlingcontfig <- ggarrange(BoxDurCONT, BoxRateCONT, BoxDCCONT,
          ncol = 3, nrow = 1,
          common.legend = TRUE,
          legend = "bottom")

handlingcontfig <- annotate_figure(
  handlingcontfig,
  left = text_grob("Handling Control", rot = 90, size = 16, face = "bold"),
)

CONTROL for recording time

Data Processing

dataTIME <- read_excel("CallContext_DataSubmission2.xlsx", 
    sheet = "Time Control", range = "A1:F57")

Individual Reaction Norms - control for recording time

DurNormTIME <- ggplot(dataTIME, aes(x = factor(`time`), y = `Call duration (ms)`, group = `Male ID`)) +
  geom_line(alpha = 0.5, color = "grey20", linewidth=1) +  # one line per individual
  theme_classic() +
  geom_point(data = dataTIME, color = "grey20", size = 2) +
  labs(x = "Recording Order", y = "Call Duration (ms)") +
  scale_x_discrete(labels = c("1" = "Early", "2" = "Late")) 

RateNormTIME <- ggplot(dataTIME, aes(x = factor(`time`), y = `Call rate (call/min)`, group = `Male ID`)) +
  geom_line(alpha = 0.5, color = "grey20", linewidth=1) +  # one line per individual
  theme_classic() +
  geom_point(data = dataTIME, color = "grey20", size = 2) +
  labs(x = "Recording Order", y = "Call Rate (calls/min)") +
  scale_x_discrete(labels = c("1" = "Early", "2" = "Late"))

DCNormTIME <- ggplot(dataTIME, aes(x = factor(`time`), y = `Duty cycle`, group = `Male ID`)) +
  geom_line(alpha = 0.5, color = "grey20", linewidth=1) +  # one line per individual
  theme_classic() +
  geom_point(data = dataTIME, color = "grey20", size = 2) +
  labs(x = "Recording Order", y = "Duty Cycle") +
  scale_x_discrete(labels = c("1" = "Early", "2" = "Late"))

ggarrange(DurNormTIME, RateNormTIME, DCNormTIME,
          ncol = 3, nrow = 1, dpi=300,
          common.legend = TRUE,
          legend = "bottom")
## $`1`

## 
## $`2`

## 
## attr(,"class")
## [1] "list"      "ggarrange"

Individual Box Plots - Time control recordings

stat.test.Durtime <- data.frame(
  group1 = "1",
  group2 = "2",
  y.position = 1900,  # adjust if needed
  label = "*",
  x = 1.5,
  xend = 1.5
)

stat.test.Ratetime <- data.frame(
  group1 = "1",
  group2 = "2",
  y.position = 48,  # adjust if needed
  label = "**",
  x = 1.5,
  xend = 1.5
)

stat.test.DCtime <- data.frame(
  group1 = "1",
  group2 = "2",
  y.position = 0.285,  # adjust if needed
  label = "***",
  x = 1.5,
  xend = 1.5
)

Generating the box plots

BoxDurTIME <- ggplot(dataTIME, aes(x=factor(`time`), y=`Call duration (ms)`)) +
  geom_boxplot(fill="grey45", 
               color="dodgerblue4", 
               alpha=0.7,
               outlier.shape = NA) +
              coord_cartesian(ylim = c(0, 2000)) +
    geom_jitter(width = 0.15,  # horizontal jitter
              size = 1.2,
              color = "dodgerblue4",
              alpha=0.6) +
  scale_x_discrete(labels = c("1" = "Early", "2" = "Late")) +
  theme_classic() +
  xlab("Recording Order") + ylab("Call Duration (ms)") +
  theme(axis.text=element_text(size=10), axis.title=element_text(size=12)) +
    stat_pvalue_manual(stat.test.Durtime, 
                     xmin = "group1", 
                     xend = "group2",
                     y.position = "y.position",
                     label = "label",
                     tip.length = 0.03)

BoxRateTIME <- ggplot(dataTIME, aes(x=factor(`time`), y=`Call rate (call/min)`) ) +
  geom_boxplot(fill="grey45", 
               color="darkolivegreen", 
               alpha=0.7,
               outlier.shape = NA) +
              coord_cartesian(ylim = c(0, 50)) +
    geom_jitter(width = 0.15,  # horizontal jitter
              size = 1.2,
              color = "darkolivegreen",
              alpha=0.6) +
  scale_x_discrete(labels = c("1" = "Early", "2" = "Late")) +
  theme_classic() +
  xlab("Recording Order") + ylab("Call Rate (calls/min)") +
  theme(axis.text=element_text(size=10), axis.title=element_text(size=12))+
    stat_pvalue_manual(stat.test.Ratetime, 
                     xmin = "group1", 
                     xend = "group2",
                     y.position = "y.position",
                     label = "label",
                     tip.length = 0.03)

BoxDCTIME <- ggplot(dataTIME, aes(x=factor(`time`), y=`Duty cycle`) ) +
  geom_boxplot(fill="grey45", 
               color="coral4", 
               alpha=0.7,
               outlier.shape = NA) +
              coord_cartesian(ylim = c(0, 0.3)) +
    geom_jitter(width = 0.15, 
              size = 1.2,
              color = "coral4",
              alpha=0.6) +
  scale_x_discrete(labels = c("1" = "Early", "2" = "Late")) +
  theme_classic() +
  xlab("Recording Order") + ylab("Duty Cycle") +
  theme(axis.text=element_text(size=10), axis.title=element_text(size=12))+
    stat_pvalue_manual(stat.test.DCtime, 
                     xmin = "group1", 
                     xend = "group2",
                     y.position = "y.position",
                     label = "label",
                     tip.length = 0.03)

Combining the box plots

timecontfig <- ggarrange(BoxDurTIME, BoxRateTIME, BoxDCTIME,
          ncol = 3, nrow = 1,
          common.legend = TRUE,
          legend = "bottom")
timecontfig <- annotate_figure(
  timecontfig,
  left = text_grob("Time Control", rot = 90, size = 16, face = "bold")
)

FINAL FIGURE 2 - Combining the control plots

row1C <- annotate_figure(
  ggarrange(BoxDurCONT, BoxRateCONT, BoxDCCONT, 
            ncol = 3,
            labels = c("a", "b", "c"),
            label.x = c(0.55, 0.55, 0.55),
            label.y = 1.02,
            align = "none"),
  left = text_grob("Handling Control", rot = 90, face = "bold", size = 14),
)

row2C <- annotate_figure(
  ggarrange(BoxDurTIME, BoxRateTIME, BoxDCTIME, 
            ncol = 3),
  left = text_grob("Time Control", rot = 90, face = "bold", size = 14)
)

control_combined <- ggarrange(row1C, row2C, 
                            ncol = 1,  # stack rows
                            heights = c(1, 1, 1))  # equal height, or adjust

print(control_combined)

ggsave("control_combined.png", control_combined, width = 9, height = 6, dpi=300)