# EWEs and STSH - Big Green Island

### Identifying and quantifying the effects of extreme weather events on short-tailed shearwater colony at Big Green Island

#### Milan Sojitra

#### 2025-11-18

    # Load libraries for data processing, modelling, and visualisation
    library(tidyverse)
    library(openxlsx2)
    library(MASS)
    library(corrplot)
    library(DT)
    library(climwin)
    library(jtools)
    library(DHARMa)
    library(lavaan)
    library(DiagrammeR)
    library(lmtest)
    library(glmmTMB)
    library(scales)

# 1 Load and Process Data

    # Load daily extreme weather event (EWE) data
    extreme_weather <- wb_to_df("Breeding_colony_ewes/Big_Green_ewes.xlsx") %>%
                       dplyr::mutate(across(4:8, ~ ifelse(is.na(.), 0, .)))  # Replace NA with 0 for analysis

    # Create binary version: 1 = event occurred, 0 = no event, NA = missing
    extreme_weather_binary <- extreme_weather %>%
                              dplyr::mutate(across(4:8, ~ ifelse(!is.na(.) & . != 0, 1, ifelse(is.na(.), NA, 0))))

    # Note:
    # Missing values in extreme weather data are replaced with zero. This is critical because the slidingwin method ("method1" and "method2") internally calculates means when NA values are present, which is not suitable when assessing extreme values. We are specifically interested in whether an extreme event occurred, not in average conditions.
    # A small number of missing values are present in the dataset, and replacing them with zero ensures consistency without introducing bias in this context.

    # Load Short-tailed Shearwater breeding data
    breeding_data <- wb_to_df("Breeding_data/STSH.xlsx", sheet = "Big_Green_Island") %>%
                     dplyr::filter(!is.na(bs))  # Remove seasons without breeding success data

    # Record sample size
    sample_size <- nrow(breeding_data)

    # Assess normality of response variable
    # If p > 0.05, the data does not significantly deviate from normality.
    shapiro.test(breeding_data$bs)

    ##
    ##  Shapiro-Wilk normality test
    ##
    ## data:  breeding_data$bs
    ## W = 0.95047, p-value = 0.2569

    # Histogram with density curve
    hist(breeding_data$bs,
         main = "Histogram of Breeding Success",
         xlab = "Breeding Success",
         col = "#a6d6fa",
         border = "white",
         prob = TRUE
    )

    # Overlay kernel density estimate
    lines(density(breeding_data$bs, na.rm = TRUE), col = "#0D92F4", lwd = 2)

    # Q-Q plot
    ggplot(breeding_data, aes(sample = bs)) +
      stat_qq() +
      stat_qq_line(colour = "red") +
      labs(title = "Q-Q Plot of Breeding Success",
           x = "Theoretical Quantiles",
           y = "Sample Quantiles") +
      theme_classic()

The Shapiro–Wilk test did not indicate a significant deviation from normality (W = 0.95047, p = 0.2569); therefore, we failed to reject the null hypothesis that the data are normally distributed. Visual assessments, including the histogram and Q–Q plot, also support the conclusion that the data approximate a normal distribution.

# 2 Sliding window analysis

## 2.1 Actual above threshold values

    # Run the sliding window analysis using actual (non-binary) values
    output1 <- slidingwin(xvar = list(wet_day             = extreme_weather$wet_day,
                                      heavy_rain_day      = extreme_weather$heavy_rain_day,
                                      very_heavy_rain_day = extreme_weather$very_heavy_rain_day,
                                      ewdp                = extreme_weather$ewdp,
                                      vwdp                = extreme_weather$vwdp),
                          cdate     = extreme_weather$date,                       # Climate date
                          bdate     = breeding_data$date,                         # Biological event date
                          baseline  = lm(bs ~ 1,
                                         data = breeding_data),                   # Baseline model
                          cohort    = breeding_data$season,                       # Group by season
                          cinterval = "day",                                      # Daily resolution
                          range     = c(96, 0),                                   # Check windows from 1 December to 7 March
                          refday    = c(07, 03),                                  # Reference date: 7 March
                          type      = "absolute",                                 # Absolute window type
                          stat      = "sum",                                      # Sum values within each window
                          func      = "lin"                                       # For linear relationship
                          )

## 2.2 Binary above threshold values

    # Run the sliding window analysis using binary event indicators
    output2 <- slidingwin(xvar = list(wet_day_bi             = extreme_weather_binary$wet_day,
                                      heavy_rain_day_bi      = extreme_weather_binary$heavy_rain_day,
                                      very_heavy_rain_day_bi = extreme_weather_binary$very_heavy_rain_day,
                                      ewdp_bi                = extreme_weather_binary$ewdp,
                                      vwdp_bi                = extreme_weather_binary$vwdp),
                          cdate     = extreme_weather_binary$date,                # Climate date
                          bdate     = breeding_data$date,                         # Biological event date
                          baseline  = lm(bs ~ 1,
                                         data = breeding_data),                   # Baseline model
                          cohort    = breeding_data$season,                       # Group by season
                          cinterval = "day",                                      # Daily resolution
                          range     = c(96, 0),                                   # Check windows from 1 December to 7 March
                          refday    = c(07, 03),                                  # Reference date: 7 March
                          type      = "absolute",                                 # Absolute window type
                          stat      = "sum",                                      # Sum values within each window
                          func      = "lin"                                       # For linear relationship
                          )

## 2.3 Merge the results

    # Combine output from actual and binary sliding window analyses
    output <- merge_results(output1, output2)

    # View merged model combinations with calculated window duration
    datatable(output$combos %>%
              dplyr::mutate(WindowDuration = WindowOpen - WindowClose + 1),
              options = list(pageLength = 10, orderClasses = TRUE))

## 2.4 Check best model for each variable

Before running the randomisation process, we need to identify the best-performing model for each extreme weather variable. This ensures that we are testing the most likely biologically relevant window against random expectation.

What we are doing here: For each weather variable (e.g., heavy rain, wet days), we extract the model with:

The lowest AIC value, and

A window duration longer than 14 days, to focus on ecologically meaningful timeframes.

These best models represent the strongest climate–breeding success relationships, and will be used for the randomisation test to assess whether the relationship is likely to have occurred by chance.

No models fit the criteria, so they are excluded from further analysis.

# 3 Check for temporal trends

## 3.1 Check for autocorrelation

    ############################## Write function to run Durbin-Watson test ##############################

    # Run Durbin-Watson test for autocorrelation
    run_dw_test <- function(data, vars) {
      results <- lapply(vars, function(var) {
        formula_obj <- as.formula(paste(var, "~ season"))
        model <- lm(formula_obj, data = data)
        dw <- dwtest(model)

        data.frame(
          variable = var,
          DW_statistic = round(dw$statistic[[1]], 3),
          p_value = round(dw$p.value, 4),
          autocorrelation = ifelse(dw$p.value < 0.05,
                                   ifelse(dw$statistic < 2, "positive", "negative"),
                                   "none"),
          row.names = NULL
        )
      })

      do.call(rbind, results)
    }

    # Check for autocorrelation in breeding success
    run_dw_test(breeding_data, c("bs"))

    ##   variable DW_statistic p_value autocorrelation
    ## 1       bs        1.647   0.129            none

## 3.2 Trend in response variable

    # Function to fit glmmTMB without autocorrelation
    fit_glmmTMB_no_ar1 <- function(data, variables_families) {
      for (entry in variables_families) {
        var <- entry$var
        fam <- entry$family

        model <- tryCatch({
          glmmTMB::glmmTMB(as.formula(paste(var, "~ season")), data = data, family = fam)
        }, error = function(e) {
          message("Could not fit model for ", var, ": ", e$message)
          return(NULL)
        })

        if (is.null(model)) next

        # Summary
        cat("\n============================\n")
        cat("Model summary for:", var, "\n")
        print(summary(model))

        # Extract p-value safely
        coefs <- summary(model)$coefficients$cond
        season_row <- grep("^season$", rownames(coefs))
        p <- coefs[season_row, "Pr(>|z|)"]

        # Base plot
        p_plot <- ggplot(data, aes(x = season, y = .data[[var]])) +
          geom_point(shape = 19, size = 2.5, color = "#7B8FA1") +
          geom_line(linewidth = 0.5, color = "#7B8FA1", na.rm = TRUE)

        # Add model trend and annotation only if p < 0.05
        if (!is.na(p) && p < 0.05) {
          # Predictions
          newdata <- data.frame(season = seq(min(data$season), max(data$season)))
          preds <- predict(model, newdata = newdata, se.fit = TRUE, type = "response")

          # Annotation
          p_label <- if (p < 0.001) "p < 0.001" else paste0("p = ", round(p, 3))
          n <- nrow(data)
          ann_text <- paste0("n = ", n, "\n", p_label)

          p_plot <- p_plot +
            geom_ribbon(data = data.frame(
              season = newdata$season,
              ymin = preds$fit - 1.96 * preds$se.fit,
              ymax = preds$fit + 1.96 * preds$se.fit
            ), aes(x = season, ymin = ymin, ymax = ymax),
            fill = "#B6D0E2", alpha = 0.4, inherit.aes = FALSE) +
            geom_line(data = data.frame(
              season = newdata$season,
              fit = preds$fit
            ), aes(x = season, y = fit),
            color = "#0D92F4", linewidth = 0.75, inherit.aes = FALSE) +
            annotate("text",
                     x = Inf, y = Inf,
                     hjust = 1.1, vjust = 1.2,
                     size = 3, color = "#1a1a1a",
                     label = ann_text)
        }

        # Final styling and print
        p_plot <- p_plot +
          scale_x_continuous(breaks = scales::pretty_breaks(n = 5),
                             expand = ggplot2::expansion(mult = c(0.02, 0.02))) +
          scale_y_continuous(breaks = scales::pretty_breaks(n = 5),
                             expand = ggplot2::expansion(mult = c(0.02, 0.02))) +
          labs(
            title = paste("Trend in", gsub("_", " ", var)),
            x = "Season", y = var
          ) +
          theme_classic(base_family = "Helvetica") +
          theme(
            plot.title = element_text(size = 10, margin = margin(b = 1), colour = "#0d0d0d"),
            axis.title = element_text(size = 9, colour = "#0d0d0d"),
            axis.text = element_text(size = 8),
            axis.ticks.length = unit(1, "pt"),
            axis.ticks = element_line(linewidth = 0.5),
            axis.line = element_blank(),
            panel.border = element_rect(color = "#1a1a1a", fill = NA, size = 0.25)
          )

        print(p_plot)
      }
    }

    # Fit glmmTMB and plot
    fit_glmmTMB_no_ar1(breeding_data,
                         variables_families = list(list(var = "bs",
                                                        family = gaussian())
                                                   )
                       )

    ##
    ## ============================
    ## Model summary for: bs
    ##  Family: gaussian  ( identity )
    ## Formula:          bs ~ season
    ## Data: data
    ##
    ##       AIC       BIC    logLik -2*log(L)  df.resid
    ##     -20.6     -17.0      13.3     -26.6        22
    ##
    ##
    ## Dispersion estimate for gaussian family (sigma^2): 0.0202
    ##
    ## Conditional model:
    ##               Estimate Std. Error z value Pr(>|z|)
    ## (Intercept) -12.830562   6.720053  -1.909   0.0562 .
    ## season        0.006694   0.003342   2.003   0.0452 *
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1