################################################################################
##
## Supplementary R script
##
## Manuscript: "Whispers in the night: An experimental approach to understanding
## low-amplitude signalling in tropical screech-owls" by LFPeixoto, PCPaiva and LPGonzaga
##
## This script reproduces all statistical analyses and figures reported in
## the manuscript and in the Supplementary Material. It is divided into:
##
##   1. Data loading and preparation
##   2. Experiment I — sender experiment (primary analyses)
##        2.1  Predictor coding
##        2.2  Random-effect structure justification (LRT)
##        2.3  Three linear mixed-effects models (LMMs)
##        2.4  Summary table for the manuscript
##        2.5  Correlation matrix among response variables
##        2.6  Density-ridge plots (Fig. 2)
##   3. Supplementary analyses for Experiment I
##        3.1  Principal Component Analysis (PCA) on the three response
##             variables (concordance check)
##        3.2  LMM on PC1
##   4. Experiment II — receiver experiment
##        4.1  Contingency table
##        4.2  Fisher's exact test
##        4.3  Effect size: Newcombe 95% confidence interval
##        4.4  Figure 4
##
## Data file: "Whispers in the night DATASETS.xlsx" (two sheets: "Exp. I" and "Exp. II"),
## available as Supplementary Material.
##
################################################################################


# ============================================================================
# 1. PACKAGES AND DATA LOADING
# ============================================================================

library(readxl)
library(lme4)
library(lmerTest)
library(ggplot2)
library(dplyr)
library(tidyr)

path <- "Data whispers.xlsx"

expI  <- read_excel(path, sheet = "Exp. I")
expII <- read_excel(path, sheet = "Exp. II")

# Trim whitespace in column names
names(expI)  <- trimws(names(expI))
names(expII) <- trimws(names(expII))

# check
str(expI)
str(expII)


# ============================================================================
# 2. EXPERIMENT I — SENDER EXPERIMENT
# ============================================================================

# ----------------------------------------------------------------------------
# 2.1  Predictor coding
# ----------------------------------------------------------------------------
# Model: 0 = sound only (no taxidermic mount), 1 = sound+model

expI$Trial_type <- factor(
  ifelse(expI$Model == 0, "sound_only", "sound_model"),
  levels = c("sound_model", "sound_only")
)

# Sampling point as factor
expI$Sampling_point <- factor(expI$`Sampling point`)

# Sanity checks
table(expI$Trial_type)
table(expI$Sampling_point, expI$Trial_type)


# ----------------------------------------------------------------------------
# 2.2  Random-effect structure justification (LRT)
# ----------------------------------------------------------------------------

response_vars <- c("Num. Phrases", "Mean Phrase", "Mean Int.")

# Helper to extract a clean column name to use as variable inside the formula
get_var <- function(df, name) {
  df[[name]]
}

cat("\n--- LRT: does (1|Distance) improve fit beyond (1|Sampling_point)? ---\n")
for (v in response_vars) {
  expI$y <- expI[[v]]
  m_simple <- lmer(y ~ Trial_type + (1 | Sampling_point),
                   data = expI, REML = FALSE)
  m_full   <- lmer(y ~ Trial_type + (1 | Sampling_point) + (1 | Distance),
                   data = expI, REML = FALSE)
  test <- anova(m_simple, m_full)
  cat("\n", v, "\n")
  print(test)
  cat("Singular full model? ", isSingular(m_full), "\n")
}
expI$y <- NULL

# Distance does not improve the models and is dropped.


# ----------------------------------------------------------------------------
# 2.3  Three linear mixed-effects models (primary analyses)
# ----------------------------------------------------------------------------

m_num <- lmer(`Num. Phrases` ~ Trial_type + (1 | Sampling_point), data = expI)
m_phr <- lmer(`Mean Phrase`  ~ Trial_type + (1 | Sampling_point), data = expI)
m_int <- lmer(`Mean Int.`    ~ Trial_type + (1 | Sampling_point), data = expI)

cat("\n\n=== Number of phrases ===\n")
print(summary(m_num))
cat("\n\n=== Mean phrase duration ===\n")
print(summary(m_phr))
cat("\n\n=== Mean interval duration ===\n")
print(summary(m_int))


# ----------------------------------------------------------------------------
# 2.4  Summary table for the manuscript
# ----------------------------------------------------------------------------

results_table <- function(model, label) {
  s <- summary(model)
  coef_row <- s$coefficients["Trial_typesound_only", ]
  vc <- as.data.frame(VarCorr(model))
  data.frame(
    Variable = label,
    Estimate = round(coef_row["Estimate"], 4),
    SE       = round(coef_row["Std. Error"], 4),
    df       = round(coef_row["df"], 2),
    t_value  = round(coef_row["t value"], 3),
    p_value  = signif(coef_row["Pr(>|t|)"], 3),
    Var_Sampling_point = round(vc$vcov[vc$grp == "Sampling_point"], 4),
    Var_Residual       = round(vc$vcov[vc$grp == "Residual"], 4),
    ICC      = round(
      vc$vcov[vc$grp == "Sampling_point"] /
        (vc$vcov[vc$grp == "Sampling_point"] + vc$vcov[vc$grp == "Residual"]),
      3
    )
  )
}

results_summary <- rbind(
  results_table(m_num, "Number of phrases"),
  results_table(m_phr, "Mean phrase duration (s)"),
  results_table(m_int, "Mean interval duration (s)")
)
cat("\n\n=== Summary table for the manuscript ===\n")
print(results_summary, row.names = FALSE)


# ----------------------------------------------------------------------------
# 2.5  Correlation matrix among the three response variables
# ----------------------------------------------------------------------------

cor_mat <- cor(expI[, response_vars], use = "complete.obs")
cat("\n\n=== Correlation matrix among response variables ===\n")
print(round(cor_mat, 3))


# ----------------------------------------------------------------------------
# 2.6  Density-ridge plots (Figure 2)
# ----------------------------------------------------------------------------

library(ggridges)

plot_ridge <- function(var_name, label_x) {
  ggplot(expI, aes(x = .data[[var_name]], y = Trial_type, fill = Trial_type)) +
    geom_density_ridges(alpha = 0.7, scale = 1.2) +
    scale_fill_manual(values = c("sound_model" = "#7CB07C",
                                 "sound_only"  = "#2C5985")) +
    theme_minimal() +
    theme(legend.position = "none") +
    labs(x = label_x, y = NULL)
}

p1 <- plot_ridge("Num. Phrases", "Number of phrases")
p2 <- plot_ridge("Mean Phrase",  "Mean phrase duration (s)")
p3 <- plot_ridge("Mean Int.",    "Mean interval duration (s)")
print(p1); print(p2); print(p3)


# ============================================================================
# 3. SUPPLEMENTARY ANALYSES FOR EXPERIMENT I
# ============================================================================

# ----------------------------------------------------------------------------
# 3.1  Principal Component Analysis (concordance check)
# ----------------------------------------------------------------------------

pca_data <- expI[, response_vars]
pca <- prcomp(pca_data, scale. = TRUE)

cat("\n\n=== Supplementary PCA on the three response variables ===\n")
print(summary(pca))
cat("\nLoadings:\n")
print(round(pca$rotation, 3))


# ----------------------------------------------------------------------------
# 3.2  LMM on PC1
# ----------------------------------------------------------------------------

pca_scores$Sampling_point <- expI$Sampling_point
pca_scores$PC1_inv <- -pca_scores$PC1

m_pc1 <- lmer(PC1_inv ~ Trial_type + (1 | Sampling_point), data = pca_scores)
cat("\n=== Supplementary LMM on PC1 (inverted: high = more signalling) ===\n")
print(summary(m_pc1))


# ============================================================================
# 4. EXPERIMENT II — RECEIVER EXPERIMENT
# ============================================================================

# ----------------------------------------------------------------------------
# 4.1  Contingency table
# ----------------------------------------------------------------------------

expII$Stimulus <- factor(expII$Stimulus,
                         levels = c("Control", "Low amp. PB"))
expII$Outcome  <- factor(expII$Outcome,
                         levels = c("Singing", "Silence"))

tab2x2 <- table(expII$Stimulus, expII$Outcome)
cat("\n\n=== Experiment II — 2x2 table ===\n")
print(tab2x2)


# ----------------------------------------------------------------------------
# 4.2  Fisher's exact test
# ----------------------------------------------------------------------------

fisher_result <- fisher.test(tab2x2)
cat("\n=== Fisher's exact test ===\n")
print(fisher_result)


# ----------------------------------------------------------------------------
# 4.3  Effect size: Newcombe 95% confidence interval
# ----------------------------------------------------------------------------
# Using DescTools::BinomDiffCI with method = "scorecc" (Newcombe's score with cc)

library(DescTools)

# Order: low-amp silenced, low-amp total, broadcast silenced, broadcast total
n_lowamp_sil   <- tab2x2["Low amp. PB", "Silence"]
n_lowamp_tot   <- sum(tab2x2["Low amp. PB", ])
n_control_sil  <- tab2x2["Control", "Silence"]
n_control_tot  <- sum(tab2x2["Control", ])

diff_ci <- BinomDiffCI(x1 = n_lowamp_sil, n1 = n_lowamp_tot,
                       x2 = n_control_sil, n2 = n_control_tot,
                       method = "scorecc")
cat("\n=== Difference in proportions (low-amp − control), 95% CI (Newcombe) ===\n")
print(diff_ci)


# ----------------------------------------------------------------------------
# 4.4  Figure 4
# ----------------------------------------------------------------------------

library(ggplot2)

plot_data <- data.frame(
  Stimulus = factor(c("Broadcast", "Broadcast",
                      "Low-amplitude", "Low-amplitude"),
                    levels = c("Broadcast", "Low-amplitude")),
  Outcome  = factor(c("Continued singing", "Silenced",
                      "Continued singing", "Silenced"),
                    levels = c("Continued singing", "Silenced")),
  Count    = c(5, 1, 0, 6)
)

p_fig4 <- ggplot(plot_data, aes(x = Outcome, y = Count, fill = Outcome)) +
  geom_bar(stat = "identity", width = 0.65,
           colour = "black", linewidth = 0.3) +
  geom_text(aes(label = Count), vjust = -0.4, size = 4) +
  facet_wrap(~ Stimulus, strip.position = "top") +
  scale_fill_manual(values = c("Continued singing" = "#233d4d",
                               "Silenced"          = "#fe7f2d")) +
  scale_y_continuous(breaks = 0:6, limits = c(0, 6.5), expand = c(0, 0)) +
  labs(x = NULL, y = "Number of trials") +
  theme_minimal(base_size = 12) +
  theme(
    legend.position    = "none",
    strip.text         = element_text(size = 12, face = "bold"),
    strip.background   = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.minor   = element_blank(),
    panel.spacing      = unit(1.5, "lines")
  )

print(p_fig4)

# Para salvar:
ggsave("Fig4_ExpII.png", p_fig4, width = 7, height = 4.5, dpi = 300)
ggsave("Fig4_ExpII.pdf", p_fig4, width = 7, height = 4.5)


################################################################################
## End of script
################################################################################