RSF for sleeping sites and Moon illumination
1 Libraries
set.seed(533)
library(tidyverse)
library(terra)
library(suncalc)
library(purrr)
library(hexbin)
library(viridis)
library(geosphere)
library(scales)
library(spdep)
library(performance)
library(flextable)
library(amt)
library(showtext)
library(patchwork)
library(lme4)
library(mapview)
library(sf)
library(sjPlot)
library(usdm)
library(visreg)
library(tmap)
library(broom)
library(broom.mixed)
library(AICcmodavg)2 Resource Selection Functions (RSF)
Resource Selection Functions (RSF) are statistical models used to quantify the probability of use of specific resources or habitats by animals, based on covariates measured at both used and available sites. In the context of bat roost-site selection, RSFs allow us to evaluate the relative importance of environmental variables—such as canopy height and moonlight, on the likelihood that a given site is selected for roosting.
How does an RSF work? - RSF models compare environmental covariates at observed (used) sites with a set of randomly generated (available) points. - By including both fixed effects (e.g., canopy height, moonlight, and their interaction) and random effects (e.g., individual, group, or night), we can account for both general patterns and non-independence in the data. - The RSF typically uses a logistic regression framework to estimate the relative probability of site selection as a function of covariates.
Application to Canopy Height and Moonlight: - In this study, we used an RSF to specifically test whether bats preferentially select roosts under taller canopies, and whether this preference changes with nocturnal light conditions (moonlight). - The model includes an interaction between canopy height and moonlight to test whether the effect of tall canopies on selection probability depends on lunar illumination. - By interpreting the model coefficients, we can determine… not only whether each factor matters, but also how their combination shapes roost-site selection strategies under different environmental conditions.
In summary: - RSFs provide a robust, quantitative approach to test ecological hypotheses about resource selection. - In this case, they allow us to rigorously assess how both canopy structure and moonlight influence bat roosting behavior, and whether these effects interact.
2.1 Load Pre-processed Data
We start by loading a list of all data frames representing different design combinations (number of pseudo-absence points and buffer size), each already processed to include the relevant spatial and environmental covariates.
2.2 Select the Optimal Sampling Scheme
From this list, we select the design previously identified as optimal (here, 10 pseudo-absence points per used location and a 300 m buffer), balancing estimator stability and minimizing spatial autocorrelation.
2.3 Prepare Covariates for Modeling
For comparability and model convergence, continuous covariates are standardized (mean-centered and scaled by standard deviation).
3 RSF Model Fitting
After preparing the data, we fit a series of Resource Selection Function (RSF) models to quantify how covariates affect roost-site selection probability.
- The models are mixed-effects logistic regressions
(
glmer), with a random intercept for each group and night (group/fecha) to account for repeated measures and hierarchical structure. - Covariates include standardized (scaled) canopy height, moonlight, NDVI, precipitation, and maximum temperature.
- We fit three models for model comparison:
1. Null Model:
Baseline model with only the random effects, representing no effect of
the environmental covariates. 2. Main Effects & Interaction
Model: Model testing the main effects of canopy height and
moonlight, as well as their interaction. 3. Full Model:
Model including the interaction plus additional covariates (NDVI,
precipitation, and max temperature) to account for further environmental
influences.
rsf_null <- glmer(
used ~ 1 +
(1 | group/fecha),
data = rsf,
family = binomial(link = "logit"))
rsf_model <- glmer(
used ~ moon_light_sc * canopy_height_sc +
(1 | group/fecha),
data = rsf,
family = binomial(link = "logit"))
rsf_model_full <- glmer(
used ~ moon_light_sc * canopy_height_sc + moon_light_sc * ndvi_sc + pp_sc + t_max_sc +
(1 | group/fecha),
data = rsf,
family = binomial(link = "logit"))
rsf_model_quad_int <- glmer(
used ~ moon_light_sc * canopy_height_sc +
moon_light_sc * I(canopy_height_sc^2) +
moon_light_sc * ndvi_sc + pp_sc + t_max_sc +
(1 | group/fecha),
data = rsf,
family = binomial(link = "logit")
)
rsf_model_quad <- glmer(
used ~ moon_light_sc * canopy_height_sc +
I(canopy_height_sc^2) +
(1 | group/fecha),
data = rsf,
family = binomial(link = "logit")
)
rsf_model_full_inter <- glmer(
used ~ moon_light_sc * canopy_height_sc + ndvi_sc* canopy_height_sc + pp_sc* canopy_height_sc + t_max_sc* canopy_height_sc + moon_light_sc * ndvi_sc +
(1 | group/fecha),
data = rsf,
family = binomial(link = "logit"))
summary(rsf_model_full_inter)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## used ~ moon_light_sc * canopy_height_sc + ndvi_sc * canopy_height_sc +
## pp_sc * canopy_height_sc + t_max_sc * canopy_height_sc +
## moon_light_sc * ndvi_sc + (1 | group/fecha)
## Data: rsf
##
## AIC BIC logLik deviance df.resid
## 2402.8 2485.1 -1188.4 2376.8 4113
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5288 -0.3287 -0.2768 -0.2290 11.3982
##
## Random effects:
## Groups Name Variance Std.Dev.
## fecha:group (Intercept) 0 0
## group (Intercept) 0 0
## Number of obs: 4126, groups: fecha:group, 384; group, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.61082 0.07302 -35.755 < 2e-16 ***
## moon_light_sc 0.04511 0.06667 0.677 0.498659
## canopy_height_sc 0.58775 0.07651 7.683 1.56e-14 ***
## ndvi_sc -0.37989 0.11268 -3.371 0.000748 ***
## pp_sc -0.03429 0.06605 -0.519 0.603667
## t_max_sc 0.01093 0.05747 0.190 0.849115
## moon_light_sc:canopy_height_sc -0.65485 0.05926 -11.051 < 2e-16 ***
## canopy_height_sc:ndvi_sc 0.01378 0.05310 0.259 0.795269
## canopy_height_sc:pp_sc 0.09965 0.06782 1.469 0.141761
## canopy_height_sc:t_max_sc 0.01815 0.05385 0.337 0.736129
## moon_light_sc:ndvi_sc 0.30571 0.09591 3.187 0.001436 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mn_lg_ cnpy__ ndv_sc pp_sc t_mx_s m__:__
## mon_lght_sc -0.065
## cnpy_hght_s -0.344 0.537
## ndvi_sc -0.204 -0.054 -0.242
## pp_sc 0.041 0.054 -0.029 -0.050
## t_max_sc 0.011 -0.067 0.008 -0.108 0.074
## mn_lght_:__ 0.543 0.122 -0.259 -0.048 -0.039 -0.006
## cnpy_hght_sc:n_ -0.291 0.069 0.156 0.465 0.005 0.007 -0.069
## cnpy_hght_sc:p_ -0.017 0.008 0.057 -0.042 -0.514 -0.022 0.089
## cnpy_hg_:__ 0.043 0.005 -0.021 -0.073 -0.018 -0.218 -0.040
## mn_lght_s:_ -0.205 -0.100 0.052 0.061 0.075 0.046 -0.460
## cnpy_hght_sc:n_ cnpy_hght_sc:p_ c__:__
## mon_lght_sc
## cnpy_hght_s
## ndvi_sc
## pp_sc
## t_max_sc
## mn_lght_:__
## cnpy_hght_sc:n_
## cnpy_hght_sc:p_ -0.072
## cnpy_hg_:__ -0.198 0.076
## mn_lght_s:_ 0.417 -0.080 -0.007
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')models <- list(
"Null" = rsf_null,
"Main Effects" = rsf_model,
"Full Model" = rsf_model_full,
"Quadratic Interaction" = rsf_model_quad_int
)
get_nparams <- function(mod) {
length(fixef(mod)) + length(VarCorr(mod))
}
model_stats <- lapply(models, function(mod) broom.mixed::glance(mod) %>% dplyr::select(AIC, nobs))
n_params <- sapply(models, get_nparams)
get_r2safe <- function(mod) {
r2out <- tryCatch(performance::r2(mod), error = function(e) NULL)
if (is.null(r2out)) {
return(list(R2_marginal = NA, R2_conditional = NA))
} else if (!("R2_conditional" %in% names(r2out))) {
return(list(R2_marginal = unname(r2out$R2_marginal), R2_conditional = NA))
} else {
return(list(R2_marginal = unname(r2out$R2_marginal), R2_conditional = unname(r2out$R2_conditional)))
}
}
r2_tab <- lapply(models, get_r2safe)
r2_df <- do.call(rbind, lapply(r2_tab, as.data.frame))
summary_tab <- tibble(
Model = names(models),
AIC = sapply(model_stats, function(x) x$AIC),
n_params = n_params,
n_obs = sapply(model_stats, function(x) x$nobs),
R2_marginal = round(r2_df$R2_marginal, 3),
R2_conditional = round(r2_df$R2_conditional, 3)
) %>%
mutate(
delta_AIC = round(AIC - min(AIC), 2)
) %>%
arrange(AIC)
ft <- flextable(summary_tab) %>%
autofit() %>%
set_header_labels(
Model = "Model",
AIC = "AIC",
delta_AIC = "ΔAIC",
n_params = "Parameters",
n_obs = "N",
R2_marginal = "R² Marginal",
R2_conditional = "R² Conditional"
)
ft## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## used ~ moon_light_sc * canopy_height_sc + moon_light_sc * I(canopy_height_sc^2) +
## moon_light_sc * ndvi_sc + pp_sc + t_max_sc + (1 | group/fecha)
## Data: rsf
##
## AIC BIC logLik deviance df.resid
## 2310.1 2386.0 -1143.0 2286.1 4114
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1804 -0.3059 -0.2632 -0.2099 7.1984
##
## Random effects:
## Groups Name Variance Std.Dev.
## fecha:group (Intercept) 1.766e-09 4.203e-05
## group (Intercept) 1.493e-03 3.864e-02
## Number of obs: 4126, groups: fecha:group, 384; group, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.94978 0.08889 -33.186 < 2e-16 ***
## moon_light_sc 0.28631 0.07279 3.934 8.37e-05 ***
## canopy_height_sc 0.55402 0.07742 7.156 8.32e-13 ***
## I(canopy_height_sc^2) 0.25236 0.03372 7.485 7.17e-14 ***
## ndvi_sc 0.05182 0.12339 0.420 0.675
## pp_sc -0.01078 0.05900 -0.183 0.855
## t_max_sc -0.01641 0.05877 -0.279 0.780
## moon_light_sc:canopy_height_sc -0.68590 0.07223 -9.496 < 2e-16 ***
## moon_light_sc:I(canopy_height_sc^2) -0.18342 0.02959 -6.199 5.69e-10 ***
## moon_light_sc:ndvi_sc -0.02443 0.09652 -0.253 0.800
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mn_lg_ cnpy__ I(__^2 ndv_sc pp_sc t_mx_s m__:__ m__:I(
## mon_lght_sc -0.249
## cnpy_hght_s -0.266 0.368
## I(cnpy__^2) -0.566 0.321 0.387
## ndvi_sc -0.380 0.156 -0.214 0.548
## pp_sc 0.035 0.060 0.035 -0.027 -0.062
## t_max_sc -0.002 -0.083 0.055 0.014 -0.076 0.056
## mn_lght_:__ 0.290 0.049 0.112 0.246 0.116 0.022 -0.021
## m__:I(__^2) 0.277 -0.349 0.269 -0.074 -0.202 0.051 0.021 0.532
## mn_lght_s:_ 0.144 -0.293 0.140 -0.224 -0.378 0.051 0.060 -0.219 0.428
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')3.1 Evaluating Model Assumptions with DHARMa
Before interpreting the results of the RSF models, it is essential to
evaluate the validity of model assumptions.
Generalized linear mixed models (GLMMs), such as our RSF, require
checking for violations such as residual non-uniformity, overdispersion,
and outliers, which can influence inference and lead to misleading
conclusions.
We use the DHARMa package to simulate scaled residuals
from the fitted model and visually inspect the diagnostics:
library(car)
vif(lm(used ~ moon_light_sc * canopy_height_sc + ndvi_sc + pp_sc + t_max_sc, data = rsf))## moon_light_sc canopy_height_sc
## 1.144704 1.409962
## ndvi_sc pp_sc
## 1.302608 1.018436
## t_max_sc moon_light_sc:canopy_height_sc
## 1.019237 1.057280
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.98866, p-value = 0.78
## alternative hypothesis: two.sided
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.012854, p-value = 0.503
## alternative hypothesis: two-sided
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: sim_res
## outliers at both margin(s) = 9, observations = 4126, p-value = 0.7257
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.0009978944 0.0041367105
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.0021812894 Visualizing Fixed Effects from the RSF Model
This script extracts and visualizes the estimated coefficients (“betas”) from the fitted full RSF model. The goal is to display the effect size and associated uncertainty (standard error) for each fixed effect, making interpretation of the model’s main covariates clear.
Key steps in the script: - Extract model
estimates:
The broom::tidy() function pulls the fixed-effect
coefficients from rsf_model_full. - Remove the
intercept:
We filter out the intercept term to focus on meaningful covariates. -
Relabel covariate names:
For clarity and publication-quality plots, we recode the technical
variable names (like moon_light_sc) to more descriptive
labels (e.g., "Moonlight (scaled)"). -
Plotting:
We use ggplot2 to create a coefficient plot: - Each point
represents a covariate’s estimated effect (log odds). - Error bars show
±1 standard error. - A horizontal dashed line at 0 highlights which
effects are positive, negative, or near zero.
This visualization helps to quickly assess which covariates have strong, weak, or uncertain effects on roost-site selection, as estimated by the RSF model.
coefs_rsf <- tidy(rsf_model_full, effects = "fixed") %>%
filter(term != "(Intercept)") %>%
mutate(term_label = case_when(
term == "moon_light_sc" ~ "Moonlight",
term == "canopy_height_sc" ~ "Canopy Height",
term == "moon_light_sc:canopy_height_sc" ~ "Moon × Canopy",
term == "moon_light_sc:ndvi_sc" ~ "Moon × NDVI",
term == "ndvi_sc" ~ "NDVI",
term == "pp_sc" ~ "Precipitation",
term == "t_max_sc" ~ "Max Temp",
TRUE ~ term))
est_plot <-
ggplot(coefs_rsf, aes(x = term_label, y = estimate, color = estimate)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray60", linewidth = 0.7) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate + std.error),
width = 0.15, linewidth = 1) +
scale_color_gradient2(low = "firebrick", mid = "grey36", high = "blue4", midpoint = 0) +
coord_cartesian(ylim = c(-1, 1)) +
labs(
x = "",
y = "Estimate (log odds, ± SE)",
color = "Estimate",
title = "RSF Model Fixed Effects") +
theme_minimal(base_size = 15) +
theme(
axis.text.x = element_text(size = 13, angle = 45, hjust = 1))
est_plot
Curves
library(ggeffects)
rsf$used <- factor(rsf$used, levels = c(FALSE, TRUE))
rsf_model_full <- glmer(
used ~ moon_light_sc * canopy_height_sc + ndvi_sc + pp_sc + t_max_sc +
(1 | group/fecha),
data = rsf ,
family = binomial(link = "logit"))
moon_quantiles <- quantile(rsf$moon_light_effective_avg, c(0.1 ,0.2, 0.4, 0.5, 0.6), na.rm=TRUE)
mean_moon <- mean(rsf$moon_light_effective_avg, na.rm=TRUE)
sd_moon <- sd(rsf$moon_light_effective_avg, na.rm=TRUE)
moon_vals_sc <- (moon_quantiles - mean_moon) / sd_moon
pred_df <- ggpredict(
rsf_model_full,
terms = c(
"canopy_height_sc",
sprintf("moon_light_sc [%.3f,%.3f,%.3f,%.3f,%.3f]",
moon_vals_sc[1], moon_vals_sc[2], moon_vals_sc[3], moon_vals_sc[4], moon_vals_sc[5])
),
type = "fixed"
)
mean_can <- mean(rsf$canopy_height, na.rm=TRUE)
sd_can <- sd(rsf$canopy_height, na.rm=TRUE)
pred_df$x_orig <- pred_df$x * sd_can + mean_can
moon_labs <- round(moon_quantiles, 2)
pred_df$group <- factor(pred_df$group, labels = paste0("Moonlight: ", moon_labs))moon_range <- rsf$moon_light_effective_avg
moon_range <- moon_range[moon_range >= 0 & moon_range <= 0.3]
moon_breaks <- quantile(
rsf$moon_light_effective_avg,
probs = seq(0, 1, length.out = 8),
na.rm = TRUE
)
moon_breaks <- unique(moon_breaks)
labels <- paste0(
sprintf("%.2f", head(moon_breaks, -1)),
"–",
sprintf("%.2f", tail(moon_breaks, -1))
)
rsf$moon_cat <- cut(
rsf$moon_light_effective_avg,
breaks = moon_breaks,
labels = labels,
include.lowest = TRUE,
right = TRUE
)
canopy_seq <- seq(8, 30, length.out = 80)
preds_list <- list()
for(mcat in levels(rsf$moon_cat)) {
moon_val <- median(rsf$moon_light_effective_avg[rsf$moon_cat == mcat], na.rm=TRUE)
for(g in unique(rsf$group)) {
df_pred <- data.frame(
canopy_height = canopy_seq,
moon_light_effective_avg = moon_val,
group = g,
ndvi_sentinel = mean(rsf$ndvi_sentinel, na.rm=TRUE),
pp = mean(rsf$pp, na.rm=TRUE),
t_max = mean(rsf$t_max, na.rm=TRUE)
)
# Escalar
df_pred$canopy_height_sc <- (df_pred$canopy_height - mean(rsf$canopy_height, na.rm=TRUE)) / sd(rsf$canopy_height, na.rm=TRUE)
df_pred$moon_light_sc <- (df_pred$moon_light_effective_avg - mean(rsf$moon_light_effective_avg, na.rm=TRUE)) / sd(rsf$moon_light_effective_avg, na.rm=TRUE)
df_pred$ndvi_sc <- (df_pred$ndvi_sentinel - mean(rsf$ndvi_sentinel, na.rm=TRUE)) / sd(rsf$ndvi_sentinel, na.rm=TRUE)
df_pred$pp_sc <- (df_pred$pp - mean(rsf$pp, na.rm=TRUE)) / sd(rsf$pp, na.rm=TRUE)
df_pred$t_max_sc <- (df_pred$t_max - mean(rsf$t_max, na.rm=TRUE)) / sd(rsf$t_max, na.rm=TRUE)
# Predicción con IC
pred <- predict(
rsf_model_full,
newdata = df_pred,
type = "link",
se.fit = TRUE,
re.form = NA
)
df_pred$fit <- plogis(pred$fit)
df_pred$lo <- plogis(pred$fit - 1.96 * pred$se.fit)
df_pred$hi <- plogis(pred$fit + 1.96 * pred$se.fit)
df_pred$moon_cat <- mcat
preds_list[[paste(mcat, g, sep = "_")]] <- df_pred
}
}
pred_all <- bind_rows(preds_list)
plot <- ggplot(pred_all, aes(x = canopy_height, y = fit, color = moon_cat, fill = moon_cat)) +
geom_line(size = 1.4) +
geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.14, color = NA) +
scale_color_viridis_d(name = "Moonlight", option = "plasma", end = 0.88) +
scale_fill_viridis_d(name = "Moonlight", option = "plasma", end = 0.88) +
labs(
x = "Canopy Height (m)",
y = "Relative Selection Strength"
) +
theme_bw(base_size = 15) +
theme(
legend.position = "top",
strip.text = element_text(face = "bold", size = 15),
axis.text.x = element_text(angle = 45, hjust = 1) )p_main <- ggplot(pred_all, aes(x = canopy_height, y = fit, color = moon_cat, fill = moon_cat)) +
geom_line(size = 1.9) +
geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.16, color = NA) +
scale_color_viridis_d(option = "inferno", end = 0.95, guide = "none") +
scale_fill_viridis_d(option = "inferno", end = 0.95, guide = "none") +
labs(
x = "Canopy Height (m)",
y = "Relative Selection Strength"
) +
theme_bw(base_size = 16) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
plot.margin = margin(8, 30, 8, 0)
)
grad_df <- data.frame(
y = seq(0, 0.6, length.out = 400),
x = 1
)
breaks_vals <- seq(0, 0.6, by = 0.1)
breaks_labs <- format(breaks_vals, nsmall = 1)
grad_plot <- ggplot(grad_df, aes(x = x, y = y, fill = y)) +
geom_tile(width = 0.33) +
scale_fill_viridis_b(option = "inferno", limits = c(0, 0.6), guide = "none") +
scale_y_continuous(
expand = c(0.03, 0.03),
breaks = breaks_vals,
labels = NULL,
limits = c(0, 0.6)
) +
theme_void() +
annotate("text", x = 0.5, y = breaks_vals, label = breaks_labs, hjust = 0, size = 4.3, fontface = "plain") +
annotate("text", x = 1.36, y = 0.3, label = "Moon Effective Light", hjust = 0.5, vjust = 0.5, angle = 270, size = 4.8, fontface = "plain") +
theme(
plot.margin = margin(8, 10, 8, 0) # top, right, bottom, left
)
final_plot <- p_main + grad_plot + plot_layout(widths = c(6,1.2))
final_plot
ndvi_seq <- seq(min(rsf$ndvi_sentinel, na.rm=TRUE),
max(rsf$ndvi_sentinel, na.rm=TRUE),
length.out = 80)
preds_list_ndvi <- list()
for(mcat in levels(rsf$moon_cat)) {
moon_val <- median(rsf$moon_light_effective_avg[rsf$moon_cat == mcat], na.rm=TRUE)
for(g in unique(rsf$group)) {
df_pred <- data.frame(
ndvi_sentinel = ndvi_seq,
moon_light_effective_avg = moon_val,
group = g,
canopy_height = mean(rsf$canopy_height, na.rm=TRUE),
pp = mean(rsf$pp, na.rm=TRUE),
t_max = mean(rsf$t_max, na.rm=TRUE)
)
# Escalado de predictores
df_pred$canopy_height_sc <- (df_pred$canopy_height - mean(rsf$canopy_height, na.rm=TRUE)) / sd(rsf$canopy_height, na.rm=TRUE)
df_pred$moon_light_sc <- (df_pred$moon_light_effective_avg - mean(rsf$moon_light_effective_avg, na.rm=TRUE)) / sd(rsf$moon_light_effective_avg, na.rm=TRUE)
df_pred$ndvi_sc <- (df_pred$ndvi_sentinel - mean(rsf$ndvi_sentinel, na.rm=TRUE)) / sd(rsf$ndvi_sentinel, na.rm=TRUE)
df_pred$pp_sc <- (df_pred$pp - mean(rsf$pp, na.rm=TRUE)) / sd(rsf$pp, na.rm=TRUE)
df_pred$t_max_sc <- (df_pred$t_max - mean(rsf$t_max, na.rm=TRUE)) / sd(rsf$t_max, na.rm=TRUE)
# Predicción
pred <- predict(
rsf_model_full,
newdata = df_pred,
type = "link",
se.fit = TRUE,
re.form = NA
)
df_pred$fit <- plogis(pred$fit)
df_pred$lo <- plogis(pred$fit - 1.96 * pred$se.fit)
df_pred$hi <- plogis(pred$fit + 1.96 * pred$se.fit)
df_pred$moon_cat <- mcat
preds_list_ndvi[[paste(mcat, g, sep = "_")]] <- df_pred
}
}
pred_all_ndvi <- bind_rows(preds_list_ndvi)
pred_all_ndvi$moon_cat <- factor(pred_all_ndvi$moon_cat)p_ndvi <- ggplot(pred_all_ndvi, aes(x = ndvi_sentinel, y = fit, color = moon_cat, fill = moon_cat)) +
geom_line(size = 1.9) +
geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.16, color = NA) +
scale_color_viridis_d(option = "green", end = 0.95, guide = "none") +
scale_fill_viridis_d(option = "green", end = 0.95, guide = "none") +
labs(
x = "NDVI",
y = "Relative Selection Strength"
) +
theme_bw(base_size = 16) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
plot.margin = margin(8, 30, 8, 0)
)
greens <- colorRampPalette(c("#e5f5e0", "#a1d99b", "#31a354", "#006d2c"))
n_colors <- length(levels(pred_all_ndvi$moon_cat))
green_colors <- greens(n_colors)
p_ndvi <- ggplot(pred_all_ndvi, aes(x = ndvi_sentinel, y = fit, color = moon_cat, fill = moon_cat)) +
geom_line(size = 1.9) +
geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.16, color = NA) +
scale_color_manual(values = green_colors, guide = "none") +
scale_fill_manual(values = green_colors, guide = "none") +
labs(
x = "NDVI",
y = "Relative Selection Strength"
) +
theme_bw(base_size = 16) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
plot.margin = margin(8, 30, 8, 0)
)
grad_df <- data.frame(
y = seq(0, 0.6, length.out = 400),
x = 1
)
breaks_vals <- seq(0, 0.6, by = 0.1)
breaks_labs <- format(breaks_vals, nsmall = 1)
grad_plot <- ggplot(grad_df, aes(x = x, y = y, fill = y)) +
geom_tile(width = 0.33) +
scale_fill_gradientn(colours = greens(100), limits = c(0, 0.6), guide = "none") +
scale_y_continuous(
expand = c(0.03, 0.03),
breaks = breaks_vals,
labels = NULL,
limits = c(0, 0.6)
) +
theme_void() +
annotate("text", x = 0.5, y = breaks_vals, label = breaks_labs, hjust = 0, size = 4.3, fontface = "plain") +
annotate("text", x = 1.36, y = 0.3, label = "Moon Effective Light", hjust = 0.5, vjust = 0.5, angle = 270, size = 4.8, fontface = "plain") +
theme(
plot.margin = margin(8, 10, 8, 0)
)
final_plot_ndvi <- p_ndvi + grad_plot + plot_layout(widths = c(6, 1.2))
final_plot_ndvi
5 Looping for Model Fitting and Metric Extraction
This script systematically fits Resource Selection Function (RSF) models for each pseudo-absence design (combination of buffer size and number of pseudo-absence points) and extracts key metrics for model comparison and interpretation.
5.1 What the function does
Input:
Each element inall_dfs_procis a data.frame with different combinations of pseudo-absence points (n_fake) and buffer size (buffer).For each dataset:
- Prepare covariates:
- Scales/standardizes all covariates (moonlight, canopy height, NDVI, precipitation, max temperature) for comparability.
- Ensures grouping and response variables are factors.
- Fit two models:
- Null model:
Only random effects (1 | group/fecha), no predictors. - Full model:
All fixed effects and interaction (moon_light_sc * canopy_height_sc + ndvi_sc + pp_sc + t_max_sc) plus random effects. - Error handling:
If any model fails to converge, returns NAs to avoid crashing the loop.
- Null model:
- Extract model metrics:
- AIC: Akaike Information Criterion for both models.
- ΔAIC: Difference in AIC between full and null models.
- AUC: Area Under the Curve for the full model, reflecting model discriminative ability.
- Model coefficients: Extracts all fixed-effect betas and their standard errors.
- Combine results:
- Outputs a tidy data.frame for each combo, with design info (buffer, n_fake), model type, metrics (AIC, delta AIC, AUC), and coefficient values.
- Prepare covariates:
5.2 Purpose and why this is important
- This approach allows direct comparison across all 30 design combinations.
- You can see how model fit and estimated effects vary with pseudo-absence design.
- It enables rigorous sensitivity analysis to ensure inference is robust to choices in spatial sampling and buffer size.
- Key for selecting an “optimal” design that balances bias, variance, and spatial autocorrelation.
combos <- names(all_dfs_proc)
fit_models_and_stats <- function(df) {
if (nrow(df) < 50) return(tibble(
AIC_null = NA, AIC_full = NA, delta_AIC = NA,
R2_null = NA, R2_full = NA, n = nrow(df)
))
df <- df %>%
mutate(
moon_light_sc = scale(moon_light_effective_avg),
canopy_height_sc = scale(canopy_height),
ndvi_sc = scale(ndvi_sentinel),
pp_sc = scale(pp),
t_max_sc = scale(t_max),
group = as.factor(group),
fecha = as.factor(fecha)
)
df <- df %>% filter(
!is.na(moon_light_sc), !is.na(canopy_height_sc),
!is.na(ndvi_sc), !is.na(pp_sc), !is.na(t_max_sc),
!is.na(used), !is.na(group), !is.na(fecha)
)
# Nulo
mod_null <- tryCatch(
glmer(used ~ 1 + (1 | group/fecha), data = df, family = binomial),
error = function(e) NULL)
# Full
mod_full <- tryCatch(
glmer(used ~ moon_light_sc * canopy_height_sc + ndvi_sc + pp_sc + t_max_sc +
(1 | group/fecha),
data = df, family = binomial),
error = function(e) NULL)
# Extraer métricas
AIC_null <- if (!is.null(mod_null)) AIC(mod_null) else NA
AIC_full <- if (!is.null(mod_full)) AIC(mod_full) else NA
r2_null <- tryCatch(performance::r2(mod_null)$R2_marginal, error = function(e) NA)
r2_full <- tryCatch(performance::r2(mod_full)$R2_marginal, error = function(e) NA)
tibble(
AIC_null = AIC_null,
AIC_full = AIC_full,
delta_AIC = AIC_null - AIC_full,
R2_null = r2_null,
R2_full = r2_full,
n = nrow(df)
)
}
results_list <- map(all_dfs_proc, fit_models_and_stats)
rsf_compare_table <- bind_rows(results_list, .id = "combo") %>%
mutate(
n_fake = as.numeric(sub("n([0-9]+)_b[0-9]+", "\\1", combo)),
buffer = as.numeric(sub("n[0-9]+_b([0-9]+)", "\\1", combo))
)
rsf_compare_table <- rsf_compare_table %>%
select(combo, n_fake, buffer, n, everything())
rsf_compare_tableggplot(rsf_compare_table, aes(x = n_fake, y = delta_AIC, color = factor(buffer))) +
geom_point(size = 2) +
geom_line(aes(group = buffer), size = 2) +
scale_color_viridis_d(name = "Buffer (m)") +
labs(x = "Pseudo-absence Points per Location", y = "ΔAIC (Null - Full)", title = "Model Improvement Across Designs") +
theme_minimal(base_size = 15)
extract_metrics <- function(df, nm) {
n_fake <- as.numeric(sub("n([0-9]+)_b[0-9]+", "\\1", nm))
buffer <- as.numeric(sub("n[0-9]+_b([0-9]+)", "\\1", nm))
df <- df %>%
mutate(
moon_light_sc = scale(moon_light_effective_avg),
canopy_height_sc = scale(canopy_height),
ndvi_sc = scale(ndvi_sentinel),
pp_sc = scale(pp),
t_max_sc = scale(t_max)
)
df$group <- as.factor(df$group)
df$fecha <- as.factor(df$fecha)
df$used <- as.factor(df$used)
mod_null <- tryCatch(
glmer(
used ~ 1 + (1 | group/fecha),
data = df,
family = binomial(link = "logit"),
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e5))
),
error = function(e) NULL
)
mod_full <- tryCatch(
glmer(
used ~ moon_light_sc * canopy_height_sc + moon_light_sc * ndvi_sc + pp_sc + t_max_sc +
(1 | group/fecha),
data = df,
family = binomial(link = "logit"),
control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e5))
),
error = function(e) NULL
)
if (is.null(mod_null) | is.null(mod_full)) {
return(tibble(
combo = nm, n_fake = n_fake, buffer = buffer,
model = c("null", "full"),
aic = NA, delta_aic = NA, auc = NA,
term = NA, estimate = NA, std.error = NA
))
}
aic_null <- AIC(mod_null)
aic_full <- AIC(mod_full)
delta_aic <- aic_full - aic_null
pred_full <- predict(mod_full, type = "response", re.form = NA)
auc_full <- tryCatch(
performance::performance_auc(prediction = pred_full, response = df$used)$AUC,
error = function(e) NA_real_
)
est <- tidy(mod_full, effects = "fixed") %>%
select(term, estimate, std.error)
est <- est %>%
mutate(
combo = nm, n_fake = n_fake, buffer = buffer,
model = "full", aic = aic_full, delta_aic = delta_aic, auc = auc_full
) %>%
select(combo, n_fake, buffer, model, aic, delta_aic, auc, term, estimate, std.error)
est_null <- tidy(mod_null, effects = "fixed") %>%
select(term, estimate, std.error) %>%
mutate(
combo = nm, n_fake = n_fake, buffer = buffer,
model = "null", aic = aic_null, delta_aic = 0, auc = NA
) %>%
select(combo, n_fake, buffer, model, aic, delta_aic, auc, term, estimate, std.error)
bind_rows(est_null, est)
}
results_table2 <- imap_dfr(all_dfs_proc, extract_metrics)plot_estimates <-
results_table %>%
filter(model == "full", !term %in% "(Intercept)") %>%
mutate(
term_label = case_when(
term == "moon_light_sc" ~ "Moonlight",
term == "canopy_height_sc" ~ "Canopy Height",
term == "moon_light_sc:canopy_height_sc" ~ "Moon × Canopy",
term == "moon_light_sc:ndvi_sc" ~ "Moon × NDVI",
term == "ndvi_sc" ~ "NDVI",
term == "pp_sc" ~ "Precipitation",
term == "t_max_sc" ~ "Max Temperature",
TRUE ~ as.character(term)
),
buffer = as.factor(buffer)
) %>%
ggplot(aes(x = n_fake, y = estimate, color = buffer, group = buffer)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray65", linewidth = 0.7) +
geom_line(size = 1.15, alpha = 0.89) +
geom_point(size = 2.3, alpha = 0.94) +
geom_errorbar(
aes(ymin = estimate - std.error, ymax = estimate + std.error),
width = 0.13, alpha = 0.60, size = 0.7
) +
facet_wrap(~ term_label, scales = "free_y") +
scale_color_viridis_d(name = "Buffer (m)", option = "viridis", end = 0.92) +
scale_y_continuous(
expand = expansion(mult = c(0.60, 1.02))
) +
labs(
x = "Pseudo-absence Points per Location",
y = "Model Coefficient (β ± SE)"
) +
theme_classic(base_size = 16) +
theme(
strip.background = element_blank(),
strip.text = element_text( size = 15),
legend.position = "top",
axis.title = element_text(),
axis.text.x = element_text(size = 13, angle = 45, hjust = 1))
plot_estimates
est_plot <- results_table %>%
filter(model == "full", term != "(Intercept)") %>%
mutate(term_label = case_when(
term == "moon_light_sc" ~ "Moonlight",
term == "canopy_height_sc" ~ "Canopy",
term == "moon_light_sc:canopy_height_sc" ~ "Moon × Canopy",
term == "moon_light_sc:ndvi_sc" ~ "Moon × NDVI",
term == "ndvi_sc" ~ "NDVI",
term == "pp_sc" ~ "Precip.",
term == "t_max_sc" ~ "Max Temp",
TRUE ~ term
))
est_plot <- ggplot(est_plot, aes(x = term_label, y = estimate, fill = as.factor(buffer))) +
geom_boxplot(outlier.shape = NA, alpha = 0.52, width = 0.95, color = NA, position = position_dodge(width = 0.7)) +
geom_jitter(
aes(color = as.factor(buffer)),
position = position_jitterdodge(jitter.width = 0.23, dodge.width = 0.8),
size = 2.9, alpha = 0.55
) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray60", linewidth = 0.7) +
scale_fill_viridis_d(name = "Buffer (m)", option = "inferno", end = 0.89) +
scale_color_viridis_d(name = "Buffer (m)", option = "inferno", end = 0.89) +
labs(
x = "",
y = "Estimate (log odds, ± SE)"
) +
theme_minimal(base_size = 18) +
theme(
axis.text.x = element_text(size = 14, angle = 45, color = "black", hjust = 1),
legend.position = "top" )
est_plot