ARTICLE 1 PHENOTYPIC VARIATION AND SELECTION
Pierre-Alexandre QUITTET
2026-04-02
DATA AND PACKAGES
Setting working directory:
## Current working directory: C:/Users/quittet/Documents/marmot-quantitative-genetics
Importing the required functions:
Importing the phenotypic data and the pedigree:
TEMPORAL VARIATION IN SPRING CONDITIONS
plot_index <- function(data, y, ylab, title) {
overall_mean <- mean(data[[y]], na.rm = TRUE)
data_temp <- data |>
mutate(anomaly = if_else(data[[y]] <= overall_mean, "#C44E52", "#4C72B0"))
plot <- data_temp |>
ggplot(aes(
x = as.numeric(as.character(year)),
y = !!(sym(y)),
group = 1
)) +
geom_hline(yintercept = overall_mean, linetype = 2, size = 1) +
geom_point(size = 3) +
geom_point(color = data_temp$anomaly, size = 2.5) +
geom_smooth(method = "lm",
color = "#474747",
fill = "#474747",
linewidth = 1.5,
alpha = 0.2,
se = TRUE) +
theme_classic(base_size = 20) +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5, face = "bold", size = 18),
panel.grid = element_line(color = "#CFCFCF", size = 0.2, linetype = 2)) +
ylab(ylab) +
ggtitle(title) +
theme(legend.position = "none", axis.title.x = element_blank())
return(plot)
}df_clim <- df_abiotic_short
df_abiotic_short$spring_index_1 <- df_abiotic_short$spring_index_1
ClimVars <- c("spring_index_1", "spring_index_2")
df_clim <- df_clim[ClimVars]
# Standardise
df_clim <- as.data.frame(apply(df_clim, 2, scale))
df_clim$year <- df_abiotic_short$year
# PCA Axis 1
plot1 <- plot_index(data = df_clim, y = "spring_index_1",
title = "Aridity & Temperature", ylab = "Spring index")
# PCA Axis 2
plot2 <- plot_index(data = df_clim, y = "spring_index_2",
ylab = "", title = "NDVI")
(plot1 | plot2)
PHENOTYPIC TEMPORAL VARIATION
RESIDUAL COMPUTATION AND PHENOTYPIC TEMPORAL TREND
Models were corrected for: - normalised capture time - sex - maternal identity
Partial residuals based on emergence date were plotted using the
compute_partial_residual() function available in
FUNCTIONS.R.
pheno <- pheno |>
arrange(capture_id)
pheno <- pheno |>
filter(age_annee == 0,
age_d == 0,
masse < 800,
annee_emergence %in% 1991:2024,
!is.na(masse),
!is.na(age_annee),
!is.na(mother)) |>
filter(!duplicated(individual_id))Mass
pheno |>
summarise(mean_mass = round(mean(masse, na.rm = TRUE), 2),
sd_mass = round(sd(masse, na.rm = TRUE), 2))## mean_mass sd_mass
## 1 401.24 113.12data_mass <- pheno[!is.na(pheno$masse), ]
data_mass$annee_emergence_reduced <- factor(
data_mass$annee_emergence,
levels = c(1991:2024),
labels = substr(c(1991:2024), start = 3, stop = 4)
)
mod_mass <- lme4::lmer(masse ~
scale(capture_time) +
sexe +
annee_emergence +
(1 | mother),
data = data_mass)
res_mass <- compute_partial_residual(data = data_mass, mod = mod_mass, marginal = TRUE)
mean_pheno <- tapply(X = res_mass$partial_res, data_mass$annee_emergence_reduced, mean)
compute_ci <- function(sd, n) qnorm(p = .975) * sd / sqrt(n)
ic_pheno <- tapply(X = res_mass$partial_res, data_mass$annee_emergence_reduced,
function(x) compute_ci(sd = sd(x), n = sum(!is.na(x))))
df_mass <- data.frame(y = mean_pheno, ic = ic_pheno, cohort = names(mean_pheno))
df_mass$cohort <- factor(df_mass$cohort, levels = substr(c(1991:2024), start = 3, stop = 4))
Size
pheno |>
summarise(mean_tibia = round(mean(tibia, na.rm = TRUE), 2),
sd_tibia = round(sd(tibia, na.rm = TRUE), 2))## mean_tibia sd_tibia
## 1 53.21 5.27data_tibia <- pheno[!is.na(pheno$tibia), ]
data_tibia <- data_tibia[data_tibia$annee_emergence %in% c(1991:2024), ]
data_tibia$annee_emergence_reduced <- factor(
data_tibia$annee_emergence,
levels = c(1991:2024),
labels = substr(c(1991:2024), start = 3, stop = 4)
)
mod_tibia <- lme4::lmer(tibia ~
scale(capture_time) +
sexe +
annee_emergence +
(1 | mother),
data = data_tibia)
res_tibia <- compute_partial_residual(data = data_tibia, mod = mod_tibia, marginal = TRUE)
mean_pheno <- tapply(X = res_tibia$partial_res, data_tibia$annee_emergence_reduced, mean)
ic_pheno <- tapply(X = res_tibia$partial_res, data_tibia$annee_emergence_reduced,
function(x) compute_ci(sd = sd(x), n = sum(!is.na(x))))
df_tibia <- data.frame(y = mean_pheno, ic = ic_pheno, cohort = names(mean_pheno))
df_tibia$cohort <- factor(df_tibia$cohort, levels = substr(c(1991:2024), start = 3, stop = 4))
# MASS
data_mass["res"] <- res_mass$partial_res
mean_data_mass <- data_mass |>
mutate(annee_emergence = as.numeric(as.character(annee_emergence))) |>
group_by(annee_emergence) |>
summarise(m = mean(masse),
ymin = m - compute_ci(sd = sd(masse), n = sum(!is.na(masse))),
ymax = m + compute_ci(sd = sd(masse), n = sum(!is.na(masse))),
precision = 1 / var(masse),
n = n())
plot_mass_gam <-
ggplot(data_mass,
aes(x = as.numeric(as.character(annee_emergence)), y = masse)) +
geom_jitter(alpha = 0.05, width = 0.1, size = 0.7) +
geom_errorbar(data = mean_data_mass,
aes(x = annee_emergence, y = m, ymin = ymin, ymax = ymax),
width = 0, color = "#616161", linewidth = 0.5) +
geom_point(data = mean_data_mass,
aes(x = annee_emergence, y = m, size = n),
col = "#383838", alpha = 1) +
geom_smooth(data = mean_data_mass,
aes(x = annee_emergence, y = m, weight = n),
method = "lm",
color = "#FFF5EE", fill = "#4682B4",
linewidth = 1, alpha = 0.6, se = TRUE) +
ylim(c(250, 700)) +
theme_classic(base_size = 20) +
theme(legend.position = "none",
panel.grid = element_line(color = "#CFCFCF", size = 0.2, linetype = 2)) +
xlab("") + ylab("Corrected mass (g)") +
theme(axis.title.x = element_blank())
# TIBIA
data_tibia["res"] <- res_tibia$partial_res
mean_data_tibia <- data_tibia |>
mutate(annee_emergence = as.numeric(as.character(annee_emergence))) |>
group_by(annee_emergence) |>
summarise(m = mean(tibia),
ymin = m - compute_ci(sd = sd(tibia), n = sum(!is.na(tibia))),
ymax = m + compute_ci(sd = sd(tibia), n = sum(!is.na(tibia))),
precision = 1 / var(tibia),
n = n())
plot_tibia_gam <-
ggplot(data_tibia, aes(x = as.numeric(as.character(annee_emergence)), y = tibia)) +
geom_jitter(alpha = 0.05, width = 0.1, size = 1) +
geom_errorbar(data = mean_data_tibia,
aes(x = annee_emergence, y = m, ymin = ymin, ymax = ymax),
width = 0, color = "#616161", linewidth = 0.5) +
geom_point(data = mean_data_tibia,
aes(x = annee_emergence, y = m, size = n),
col = "#383838", alpha = 1) +
geom_smooth(data = mean_data_tibia,
aes(x = annee_emergence, y = m, weight = n),
method = "gam",
formula = y ~ s(x, k = 4),
color = "#FFF5EE", fill = "#CD4F39",
linewidth = 1, alpha = 0.6, se = TRUE) +
theme_classic(base_size = 20) +
scale_x_continuous(limits = c(1990, 2024)) +
theme(legend.position = "none",
panel.grid = element_line(color = "#CFCFCF", size = 0.2, linetype = 2)) +
ylim(c(42, 62)) +
xlab("") + ylab("Corrected tibia length (mm)") +
theme(axis.title.x = element_blank())
gridExtra::grid.arrange(plot_mass_gam, plot_tibia_gam, ncol = 2)
Testing temporal trends
mod <- lmer(masse ~ sexe + scale(capture_time) +
scale(as.numeric(annee_emergence), scale = FALSE) +
(1 | mother),
data = pheno)
summary(mod)## Linear mixed model fit by REML ['lmerMod']
## Formula: masse ~ sexe + scale(capture_time) + scale(as.numeric(annee_emergence), scale = FALSE) + (1 | mother)
## Data: pheno
##
## REML criterion at convergence: 20195.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4390 -0.5816 -0.0891 0.4822 4.0149
##
## Random effects:
## Groups Name Variance Std.Dev.
## mother (Intercept) 5088 71.33
## Residual 7250 85.15
## Number of obs: 1697, groups: mother, 164
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 389.8348 6.6066 59.007
## sexem 15.1266 4.3400 3.485
## scale(capture_time) 32.7781 3.0904 10.606
## scale(as.numeric(annee_emergence), scale = FALSE) -2.0155 0.5835 -3.454
##
## Correlation of Fixed Effects:
## (Intr) sexem scl(_)
## sexem -0.363
## scl(cptr_t) -0.065 0.058
## s(.(_),s=FA 0.006 -0.001 0.180## 2.5 % 97.5 %
## .sig01 62.100664 81.0441996
## .sigma 82.166446 88.1972452
## (Intercept) 376.885104 402.7645556
## sexem 6.615348 23.6266780
## scale(capture_time) 26.641156 38.8488755
## scale(as.numeric(annee_emergence), scale = FALSE) -3.157175 -0.8723805mod <- lmer(tibia ~ sexe + scale(capture_time) +
scale(as.numeric(annee_emergence), scale = FALSE) +
(1 | mother),
data = pheno)
summary(mod)## Linear mixed model fit by REML ['lmerMod']
## Formula: tibia ~ sexe + scale(capture_time) + scale(as.numeric(annee_emergence), scale = FALSE) + (1 | mother)
## Data: pheno
##
## REML criterion at convergence: 8848.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0062 -0.5442 0.0207 0.5757 3.8756
##
## Random effects:
## Groups Name Variance Std.Dev.
## mother (Intercept) 11.60 3.406
## Residual 15.11 3.887
## Number of obs: 1539, groups: mother, 148
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 53.41555 0.33036 161.687
## sexem 1.08436 0.20839 5.203
## scale(capture_time) 0.64836 0.15542 4.172
## scale(as.numeric(annee_emergence), scale = FALSE) -0.21438 0.03199 -6.701
##
## Correlation of Fixed Effects:
## (Intr) sexem scl(_)
## sexem -0.346
## scl(cptr_t) -0.015 0.053
## s(.(_),s=FA -0.147 -0.011 0.066## 2.5 % 97.5 %
## .sig01 2.9507995 3.8890237
## .sigma 3.7440116 4.0331982
## (Intercept) 52.7692081 54.0640399
## sexem 0.6762483 1.4930696
## scale(capture_time) 0.3431340 0.9524159
## scale(as.numeric(annee_emergence), scale = FALSE) -0.2771740 -0.1517813
SELECTION ANALYSIS
Fitness proxies used were first-winter survival and LRS (lifetime reproductive success). Survival was estimated from a CMR (capture–mark–recapture) model. Individual life histories were extracted using the most probable path from the VITERBI algorithm. Confidence intervals were obtained by bootstrap (Mitchell-Olds & Shaw, 1987).
UTILITY FUNCTIONS
# Significance levels
p_values_function <- function(p) {
ifelse(p <= 0.001, "***",
ifelse(p <= 0.01, "**",
ifelse(p <= 0.05, "*",
ifelse(p < 0.1, "°", ""))))
}
################################################################################
# #
# Structured output: differentials and gradients #
# #
################################################################################
print_selection <- function(df,
title = "",
show_intercept = FALSE) {
cat("\n")
cat(strrep("═", 65), "\n")
cat(" ", title, "\n", sep = "")
cat(strrep("═", 65), "\n")
df_print <- df
if (!show_intercept) {
df_print <- df_print[!grepl("Intercept", df_print$params), ]
}
df_print <- df_print |>
transmute(
Type = method,
Parameter = params,
Estimate = formatC(y, format = "f", digits = 4),
`[2.5%` = formatC(ci_lower, format = "f", digits = 4),
`97.5%]` = formatC(ci_upper, format = "f", digits = 4),
`p-value` = formatC(p_value, format = "f", digits = 3),
Sig = p
)
# Split by linear / quadratic
for (m in unique(df_print$Type)) {
cat("\n ──", m, "──\n")
sub <- df_print[df_print$Type == m, -1]
print(sub, row.names = FALSE, right = FALSE)
}
cat("\n")
}
################################################################################
# #
# Structured output: correlational selection #
# #
################################################################################
print_selection_corr <- function(df, title = "") {
cat("\n")
cat(strrep("═", 65), "\n")
cat(" ", title, "\n", sep = "")
cat(strrep("═", 65), "\n\n")
df_print <- df |>
transmute(
Term = terme,
`γ_ij` = formatC(gamma_ij, format = "f", digits = 4),
`[2.5%` = formatC(ci_lower, format = "f", digits = 4),
`97.5%]` = formatC(ci_upper, format = "f", digits = 4),
`p-value` = formatC(p_value, format = "f", digits = 3),
Sig = p
)
print(df_print, row.names = FALSE, right = FALSE)
}
################################################################################
# #
# Structured output: climate selection (continuous covariate) #
# #
################################################################################
print_selection_climate <- function(df,
title = "",
show_intercept = FALSE) {
cat("\n")
cat(strrep("═", 65), "\n")
cat(" ", title, "\n", sep = "")
cat(strrep("═", 65), "\n")
df_print <- df
if (!show_intercept) {
df_print <- df_print[!grepl("Intercept", df_print$params), ]
}
df_print <- df_print |>
transmute(
Type = method,
Parameter = params,
Estimate = formatC(y, format = "f", digits = 4),
`[2.5%` = formatC(ci_lower, format = "f", digits = 4),
`97.5%]` = formatC(ci_upper, format = "f", digits = 4),
`p-value` = formatC(p_value, format = "f", digits = 3),
Sig = p
)
for (m in unique(df_print$Type)) {
sub <- df_print[df_print$Type == m, -1]
main_fx <- sub[!grepl(":", sub$Parameter), ]
inter <- sub[ grepl(":", sub$Parameter), ]
cat("\n ── ", m, " ──\n", sep = "")
if (nrow(main_fx) > 0) {
cat(" Main effects:\n")
print(main_fx, row.names = FALSE, right = FALSE)
}
if (nrow(inter) > 0) {
cat("\n Trait × climate interactions:\n")
print(inter, row.names = FALSE, right = FALSE)
}
}
cat("\n")
}COMPUTATION FUNCTIONS
################################################################################
# #
# Global differentials and gradients #
# #
################################################################################
compute_selection <- function(data = pheno,
trait = "masse_scaled",
fitness_proxy = "w1",
iteration = 10,
weights = NULL) {
for (trait.temp in trait) {
data <- data[!is.na(data[[trait.temp]]), ]
}
cov_boot_linear <- as.data.frame(matrix(
NA, ncol = length(trait) + 1, nrow = iteration))
cov_boot_non_linear <- as.data.frame(matrix(
NA, ncol = length(trait) * 2 + 1, nrow = iteration))
formula_mod_linear <- as.formula(paste0(
"fitness ~ ", paste0(trait, collapse = " + ")))
formula_mod_non_linear <- as.formula(paste0(
"fitness ~ ", paste0(trait, collapse = " + "),
paste0(" + I(0.5*", trait, "^2)", collapse = " + ")))
k <- gsub("w", "", fitness_proxy)
data$row_id <- seq_len(nrow(data))
for (i in seq_len(iteration)) {
index_random <- data |>
group_by(annee_emergence, territoire) |>
reframe(row_id = sample(row_id, size = n(), replace = TRUE)) |>
pull(row_id)
data_boot_i <- data[index_random, ]
fitness_abs <- if (fitness_proxy == "LRS_relative") sym("LRS") else sym(paste0("W", k))
if (!is.null(weights)) {
weight_var <- sym(weights)
data_boot_i <- data_boot_i |>
group_by(annee_emergence) |>
filter(!is.na(!!weight_var) & !is.na(!!fitness_abs)) |>
mutate(fitness = !!fitness_abs / weighted.mean(!!fitness_abs, w = !!weight_var)) |>
ungroup()
weights_mod <- proportions(data_boot_i[[weights]])
} else {
data_boot_i <- data_boot_i |>
group_by(annee_emergence) |>
filter(!is.na(!!fitness_abs)) |>
mutate(fitness = !!fitness_abs / mean(!!fitness_abs)) |>
ungroup()
weights_mod <- NULL
}
mod_linear <- lm(formula_mod_linear, data = data_boot_i, weights = weights_mod, na.action = "na.omit")
mod_non_linear <- lm(formula_mod_non_linear, data = data_boot_i, weights = weights_mod, na.action = "na.omit")
cov_boot_linear[i, ] <- coef(mod_linear)
cov_boot_non_linear[i, ] <- coef(mod_non_linear)
}
coef_linear_temp <- coef(mod_linear)
coef_non_linear_temp <- coef(mod_non_linear)
colnames(cov_boot_linear) <- names(coef_linear_temp)
colnames(cov_boot_non_linear) <- names(coef_non_linear_temp)
summarise_boot <- function(mat) {
list(
mean = apply(mat, 2, mean, na.rm = TRUE),
lower = apply(mat, 2, quantile, probs = 0.025, na.rm = TRUE),
upper = apply(mat, 2, quantile, probs = 0.975, na.rm = TRUE),
p = apply(mat, 2, function(x)
min(sum(x < 0, na.rm = TRUE), sum(x > 0, na.rm = TRUE)) / iteration)
)
}
sl <- summarise_boot(cov_boot_linear)
sn <- summarise_boot(cov_boot_non_linear)
data.frame(
params = c(names(coef_linear_temp), names(coef_non_linear_temp)),
y = round(c(sl$mean, sn$mean), 3),
ci_lower = round(c(sl$lower, sn$lower), 3),
ci_upper = round(c(sl$upper, sn$upper), 3),
trait = trait,
p_value = round(c(sl$p, sn$p), 3),
p = c(p_values_function(sl$p), p_values_function(sn$p)),
method = c(rep("Linear", length(coef_linear_temp)),
rep("Quadratic", length(coef_non_linear_temp)))
)
}
################################################################################
# #
# Correlational selection (γ_ij mass × tibia) #
# #
################################################################################
compute_correlational_selection <- function(data = pheno,
fitness_proxy = "w1",
iteration = 1000,
weights = NULL) {
k <- gsub("w", "", fitness_proxy)
data <- data[!is.na(data$masse_scaled) & !is.na(data$tibia_scaled), ]
data$row_id <- seq_len(nrow(data))
formula_mod <- fitness ~
masse_scaled + tibia_scaled +
I(0.5 * masse_scaled^2) +
I(0.5 * tibia_scaled^2) +
masse_scaled:tibia_scaled
boot_corr <- numeric(iteration)
for (i in seq_len(iteration)) {
index_random <- data |>
group_by(annee_emergence, territoire) |>
reframe(row_id = sample(row_id, size = n(), replace = TRUE)) |>
pull(row_id)
data_boot_i <- data[index_random, ]
fitness_abs <- if (fitness_proxy == "LRS_relative") sym("LRS") else sym(paste0("W", k))
if (!is.null(weights)) {
weight_var <- sym(weights)
data_boot_i <- data_boot_i |>
group_by(annee_emergence) |>
filter(!is.na(!!weight_var) & !is.na(!!fitness_abs)) |>
mutate(fitness = !!fitness_abs / weighted.mean(!!fitness_abs, w = !!weight_var)) |>
ungroup()
weights_mod <- proportions(data_boot_i[[weights]])
} else {
data_boot_i <- data_boot_i |>
group_by(annee_emergence) |>
filter(!is.na(!!fitness_abs)) |>
mutate(fitness = !!fitness_abs / mean(!!fitness_abs)) |>
ungroup()
weights_mod <- NULL
}
mod <- lm(formula_mod, data = data_boot_i, weights = weights_mod, na.action = "na.omit")
boot_corr[i] <- coef(mod)["masse_scaled:tibia_scaled"]
}
p_two_sided <- min(sum(boot_corr < 0, na.rm = TRUE),
sum(boot_corr > 0, na.rm = TRUE)) / iteration
data.frame(
terme = "masse_scaled:tibia_scaled",
gamma_ij = round(mean(boot_corr, na.rm = TRUE), 4),
ci_lower = round(quantile(boot_corr, 0.025, na.rm = TRUE), 4),
ci_upper = round(quantile(boot_corr, 0.975, na.rm = TRUE), 4),
p_value = round(p_two_sided, 4),
p = p_values_function(p_two_sided),
row.names = NULL
)
}
################################################################################
# #
# Temporal selection #
# #
################################################################################
compute_selection_across_time <- function(data,
trait,
fitness_proxy,
iteration = 1000,
weights = NULL,
minimum_sample_size = 30) {
data_boot <- data
for (i in seq_along(trait)) data_boot <- data_boot[!is.na(data_boot[[trait[i]]]), ]
year_n <- data_boot |> group_by(annee_emergence) |> mutate(n = n())
data_boot <- data_boot[year_n$n > minimum_sample_size, ]
data_boot$row_id <- seq_len(nrow(data_boot))
year <- as.character(sort(unique(data_boot$annee_emergence)))
colnames_temp <- unlist(lapply(trait, function(t) paste0(t, "_", year)))
cov_boot <- as.data.frame(matrix(NA, ncol = length(colnames_temp), nrow = iteration,
dimnames = list(NULL, colnames_temp)))
formula_mod <- as.formula(paste0(
fitness_proxy, "~",
paste0(trait, collapse = "+"), "+",
paste0(trait, ":annee_emergence", collapse = "+")))
k <- gsub("w", "", fitness_proxy)
data_boot$row_id <- seq_len(nrow(data_boot))
for (i in seq_len(iteration)) {
index_random <- data_boot |>
group_by(annee_emergence, territoire) |>
reframe(row_id = sample(row_id, size = n(), replace = TRUE)) |>
pull(row_id)
data_boot_i <- data_boot[index_random, ]
fitness_abs <- if (fitness_proxy == "LRS_relative") sym("LRS") else sym(paste0("W", k))
if (!is.null(weights)) {
weight_var <- sym(weights)
data_boot_i <- data_boot_i |>
group_by(annee_emergence) |>
filter(!is.na(!!weight_var) & !is.na(!!fitness_abs)) |>
mutate(fitness = !!fitness_abs / weighted.mean(!!fitness_abs, w = !!weight_var)) |>
ungroup()
weights_mod <- data_boot_i[[weights]]
weights_mod <- weights_mod / sum(weights_mod, na.rm = TRUE)
} else {
data_boot_i <- data_boot_i |>
group_by(annee_emergence) |>
filter(!is.na(!!fitness_abs)) |>
mutate(fitness = !!fitness_abs / mean(!!fitness_abs)) |>
ungroup()
weights_mod <- NULL
}
mod_temp <- lm(formula_mod, data = data_boot_i, weights = weights_mod, na.action = "na.omit")
params <- c()
for (j in seq_along(trait)) {
params <- c(params, as.data.frame(
emmeans::emtrends(mod_temp, ~ annee_emergence, var = trait[j]))[, 2])
}
cov_boot[i, ] <- params
}
mean_boot <- apply(cov_boot, 2, mean, na.rm = TRUE)
ci_lower <- apply(cov_boot, 2, quantile, probs = 0.025, na.rm = TRUE)
ci_upper <- apply(cov_boot, 2, quantile, probs = 0.975, na.rm = TRUE)
precision <- apply(cov_boot, 2, function(x) 1 / sd(x, na.rm = TRUE))
p <- apply(cov_boot, 2, function(x) p_values_function(sum(x < 0, na.rm = TRUE) / iteration))
df_gg <- data.frame(
y = round(mean_boot, 3),
ci_lower = round(ci_lower, 3),
ci_upper = round(ci_upper, 3),
trait = rep(trait, each = length(year)),
year = rep(year, length(trait)),
precision = round(precision, 3),
p = p
)
df_gg$trait <- factor(df_gg$trait, levels = trait)
df_gg$year <- factor(df_gg$year, levels = year)
df_gg
}
bind_selection_across_time <- function(trait = c("masse_scaled", "tibia_scaled"),
data = data,
fitness_proxy = fitness_proxy,
iteration = iteration,
weights = NULL,
minimum_sample_size = 5) {
output <- NULL
for (trait_temp in trait) {
output <- rbind(output,
compute_selection_across_time(data, fitness_proxy = fitness_proxy,
iteration = iteration, weights = weights,
trait = trait_temp,
minimum_sample_size = minimum_sample_size))
}
output
}
################################################################################
# #
# Climate selection (continuous covariate) #
# #
################################################################################
compute_selection_climate_covariable <- function(data = pheno,
climate_variable = "spring_index_2_y",
trait = "masse_scaled",
fitness_proxy = "w1",
iteration = 10,
weights = NULL) {
for (trait.temp in trait) data <- data[!is.na(data[[trait.temp]]), ]
formula_mod_linear <- as.formula(paste0(
"fitness ~ (", paste0(trait, collapse = " + "), ") * ", climate_variable))
formula_mod_non_linear <- as.formula(paste0(
"fitness ~ (", paste0(trait, collapse = " + "), ") * ", climate_variable,
" + (", paste0("I(0.5 * ", trait, "^2)", collapse = " + "), ")"))
k <- gsub("w", "", fitness_proxy)
data$row_id <- seq_len(nrow(data))
cov_boot_linear <- cov_boot_non_linear <- NULL
for (i in seq_len(iteration)) {
index_random <- data |>
group_by(annee_emergence, territoire) |>
reframe(row_id = sample(row_id, size = n(), replace = TRUE)) |>
pull(row_id)
data_boot_i <- data[index_random, ]
fitness_abs <- if (fitness_proxy == "LRS_relative") sym("LRS") else sym(paste0("W", k))
if (!is.null(weights)) {
weight_var <- sym(weights)
data_boot_i <- data_boot_i |>
group_by(annee_emergence) |>
filter(!is.na(!!weight_var) & !is.na(!!fitness_abs)) |>
mutate(fitness = !!fitness_abs / weighted.mean(!!fitness_abs, w = !!weight_var)) |>
ungroup()
weights_mod <- proportions(data_boot_i[[weights]])
} else {
data_boot_i <- data_boot_i |>
group_by(annee_emergence) |>
filter(!is.na(!!fitness_abs)) |>
mutate(fitness = !!fitness_abs / mean(!!fitness_abs)) |>
ungroup()
weights_mod <- NULL
}
mod_linear <- lm(formula_mod_linear, data = data_boot_i, weights = weights_mod, na.action = "na.omit")
mod_non_linear <- lm(formula_mod_non_linear, data = data_boot_i, weights = weights_mod, na.action = "na.omit")
coef_l <- coef(mod_linear)
coef_n <- coef(mod_non_linear)
if (is.null(cov_boot_linear)) {
cov_boot_linear <- as.data.frame(matrix(
NA, ncol = length(coef_l), nrow = iteration,
dimnames = list(NULL, names(coef_l))))
cov_boot_non_linear <- as.data.frame(matrix(
NA, ncol = length(coef_n), nrow = iteration,
dimnames = list(NULL, names(coef_n))))
}
cov_boot_linear[i, ] <- coef_l
cov_boot_non_linear[i, ] <- coef_n
}
summarise_boot <- function(mat) {
list(
mean = apply(mat, 2, mean, na.rm = TRUE),
lower = apply(mat, 2, quantile, probs = 0.025, na.rm = TRUE),
upper = apply(mat, 2, quantile, probs = 0.975, na.rm = TRUE),
p = apply(mat, 2, function(x)
min(sum(x < 0, na.rm = TRUE), sum(x > 0, na.rm = TRUE)) / iteration)
)
}
sl <- summarise_boot(cov_boot_linear)
sn <- summarise_boot(cov_boot_non_linear)
data.frame(
params = c(names(coef_l), names(coef_n)),
y = round(c(sl$mean, sn$mean), 3),
ci_lower = round(c(sl$lower, sn$lower), 3),
ci_upper = round(c(sl$upper, sn$upper), 3),
trait = trait,
p_value = round(c(sl$p, sn$p), 3),
p = c(p_values_function(sl$p), p_values_function(sn$p)),
method = c(rep("Linear", length(coef_l)),
rep("Quadratic", length(coef_n)))
)
}
################################################################################
# #
# Fitness surface visualisation #
# #
################################################################################
plot_selection <- function(mod_result = S_climatic_mass,
title = "",
xlab = "Trait",
ylab = "w",
type_selection = c("Linear", "Quadratic")) {
mod_result_filtered <- mod_result |> filter(method == type_selection)
b <- mod_result_filtered$y
b_l <- mod_result_filtered$ci_lower
b_u <- mod_result_filtered$ci_upper
trait_seq <- seq(-2, 2, length.out = 300)
envs <- c(-2, 0, 2)
env_labels <- c("Spring index = -2", "Spring index = 0", "Spring index = +2")
df_list <- lapply(seq_along(envs), function(i) {
e <- envs[i]
if (type_selection == "Linear") {
w <- b[1] + b[2] * trait_seq + b[3] * e + b[4] * trait_seq * e
w_l <- b_l[1] + b_l[2] * trait_seq + b_l[3] * e + b_l[4] * trait_seq * e
w_u <- b_u[1] + b_u[2] * trait_seq + b_u[3] * e + b_u[4] * trait_seq * e
} else {
w <- b[1] + b[2] * trait_seq + b[3] * e + b[4] * (0.5 * trait_seq^2) + b[5] * trait_seq * e
w_l <- b_l[1] + b_l[2] * trait_seq + b_l[3] * e + b_l[4] * (0.5 * trait_seq^2) + b_l[5] * trait_seq * e
w_u <- b_u[1] + b_u[2] * trait_seq + b_u[3] * e + b_u[4] * (0.5 * trait_seq^2) + b_u[5] * trait_seq * e
}
data.frame(trait = trait_seq, w = w, w_lower = w_l, w_upper = w_u,
env = factor(env_labels[i], levels = env_labels))
})
df_fit <- do.call(rbind, df_list)
pal <- c("Spring index = -2" = "#C44E52",
"Spring index = 0" = "#474747",
"Spring index = +2" = "#4C72B0")
ggplot(df_fit, aes(x = trait, y = w, colour = env, fill = env)) +
geom_hline(yintercept = 1, colour = "#454545", linetype = 2, linewidth = 0.8) +
geom_ribbon(aes(ymin = w_lower, ymax = w_upper),
alpha = 0.2, colour = NA, show.legend = FALSE) +
geom_line(linewidth = 1.5) +
scale_colour_manual(name = "Spring index", values = pal, labels = c("-2", " 0", "+2")) +
scale_fill_manual (name = "Spring index", values = pal, labels = c("-2", " 0", "+2")) +
labs(title = title, x = xlab, y = ylab) +
theme_classic(base_size = 20) +
theme(
legend.position = c(0.04, 1.04),
legend.justification = c(0, 1),
legend.background = element_rect(fill = alpha("white", 0.4),
colour = "#212121", linewidth = 0.5),
legend.key.size = unit(0.8, "lines"),
legend.text = element_text(size = 13),
legend.title = element_text(size = 13, face = "bold"),
plot.title = element_text(size = 20),
panel.grid = element_line(color = "#CFCFCF", linewidth = 0.2, linetype = 2)
) +
guides(colour = guide_legend(override.aes = list(linewidth = 1), keywidth = unit(1, "cm")))
}DATA CLEANING
# Binarise survival
pheno$W1 <- ifelse(pheno$W1 < 0.5, 0, 1)
# Relative fitness per year
pheno <- pheno |>
group_by(annee_emergence) |>
mutate(w1 = W1 / mean(W1, na.rm = TRUE)) |>
as.data.frame()
# Standardise traits
pheno <- pheno |>
mutate(
masse_scaled = scale(masse)[, 1],
tibia_scaled = scale(tibia)[, 1],
longueur_corps_scaled = scale(longueur_corps)
)GLOBAL SELECTION
Survival
# ── Selection differentials (S) — one trait at a time ────────────────────────
S_mass_w1 <- compute_selection(
trait = "masse_scaled",
data = pheno,
fitness_proxy = "w1",
iteration = iteration
)
S_tibia_w1 <- compute_selection(
trait = "tibia_scaled",
data = pheno,
fitness_proxy = "w1",
iteration = iteration
)
# ── Selection gradients (β) — multivariate model ─────────────────────────────
beta_w1 <- compute_selection(
trait = c("masse_scaled", "tibia_scaled"),
data = pheno,
fitness_proxy = "w1",
iteration = iteration
)
# ── Correlational selection (γ_ij) ───────────────────────────────────────────
corr_w1 <- compute_correlational_selection(
data = pheno,
fitness_proxy = "w1",
iteration = iteration
)Selection differential — Mass
##
## ═════════════════════════════════════════════════════════════════
## S (Mass) — Proxy: First-winter survival
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.1180 0.0840 0.1530 0.000 ***
##
## ── Quadratic ──
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.1280 0.0930 0.1640 0.000 ***
## I(0.5 * masse_scaled^2) -0.0390 -0.0850 0.0080 0.053 °Selection differential — Tibia
##
## ═════════════════════════════════════════════════════════════════
## S (Tibia) — Proxy: First-winter survival
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled 0.0960 0.0660 0.1270 0.000 ***
##
## ── Quadratic ──
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled 0.0980 0.0660 0.1340 0.000 ***
## I(0.5 * tibia_scaled^2) -0.0100 -0.0550 0.0340 0.312Selection gradient — Mass
print_selection(
beta_w1 |> filter(grepl("masse", params) | grepl("Intercept", params)),
title = "β gradient (Mass) — Proxy: First-winter survival"
)##
## ═════════════════════════════════════════════════════════════════
## β gradient (Mass) — Proxy: First-winter survival
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.1060 0.0650 0.1510 0.000 ***
##
## ── Quadratic ──
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.1220 0.0780 0.1700 0.000 ***
## I(0.5 * masse_scaled^2) -0.0410 -0.1010 0.0240 0.093 °Selection gradient — Tibia
print_selection(
beta_w1 |> filter(grepl("tibia", params) | grepl("Intercept", params)),
title = "β gradient (Tibia) — Proxy: First-winter survival"
)##
## ═════════════════════════════════════════════════════════════════
## β gradient (Tibia) — Proxy: First-winter survival
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled 0.0260 -0.0100 0.0600 0.066 °
##
## ── Quadratic ──
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled 0.0160 -0.0270 0.0580 0.202
## I(0.5 * tibia_scaled^2) 0.0090 -0.0350 0.0550 0.366Correlational selection
##
## ═════════════════════════════════════════════════════════════════
## γ_ij (Mass × Tibia) — Proxy: First-winter survival
## ═════════════════════════════════════════════════════════════════
##
## Term γ_ij [2.5% 97.5%] p-value Sig
## masse_scaled:tibia_scaled 0.0072 -0.0638 0.0714 0.413
LRS
data_lrs <- pheno[pheno$annee_emergence %in% 1990:2018, ]
# ── Selection differentials (S) ───────────────────────────────────────────────
S_mass_lrs <- compute_selection(
trait = "masse_scaled",
data = data_lrs,
fitness_proxy = "LRS_relative",
iteration = iteration
)
S_tibia_lrs <- compute_selection(
trait = "tibia_scaled",
data = data_lrs,
fitness_proxy = "LRS_relative",
iteration = iteration
)
# ── Selection gradients (β) ───────────────────────────────────────────────────
beta_lrs <- compute_selection(
trait = c("masse_scaled", "tibia_scaled"),
data = data_lrs,
fitness_proxy = "LRS_relative",
iteration = iteration
)
# ── Correlational selection (γ_ij) ───────────────────────────────────────────
corr_lrs <- compute_correlational_selection(
data = data_lrs,
fitness_proxy = "LRS_relative",
iteration = iteration
)Selection differential — Mass
##
## ═════════════════════════════════════════════════════════════════
## S (Mass) — Proxy: LRS
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.2620 0.1400 0.3850 0.000 ***
##
## ── Quadratic ──
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.2880 0.1440 0.4250 0.000 ***
## I(0.5 * masse_scaled^2) -0.0770 -0.2880 0.1250 0.242Selection differential — Tibia
##
## ═════════════════════════════════════════════════════════════════
## S (Tibia) — Proxy: LRS
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled 0.0720 -0.0160 0.1640 0.060 °
##
## ── Quadratic ──
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled 0.0410 -0.0760 0.1620 0.254
## I(0.5 * tibia_scaled^2) 0.1080 -0.0630 0.2870 0.124Selection gradient — Mass
print_selection(
beta_lrs |> filter(grepl("masse", params) | grepl("Intercept", params)),
title = "β gradient (Mass) — Proxy: LRS"
)##
## ═════════════════════════════════════════════════════════════════
## β gradient (Mass) — Proxy: LRS
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.2130 0.0670 0.3660 0.002 **
##
## ── Quadratic ──
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.3410 0.1720 0.5130 0.000 ***
## I(0.5 * masse_scaled^2) -0.3060 -0.5090 -0.0580 0.005 **Selection gradient — Tibia
print_selection(
beta_lrs |> filter(grepl("tibia", params) | grepl("Intercept", params)),
title = "β gradient (Tibia) — Proxy: LRS"
)##
## ═════════════════════════════════════════════════════════════════
## β gradient (Tibia) — Proxy: LRS
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled -0.0450 -0.1560 0.0600 0.215
##
## ── Quadratic ──
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled -0.1440 -0.2970 0.0010 0.027 *
## I(0.5 * tibia_scaled^2) 0.1990 0.0050 0.3850 0.024 *Correlational selection
##
## ═════════════════════════════════════════════════════════════════
## γ_ij (Mass × Tibia) — Proxy: LRS
## ═════════════════════════════════════════════════════════════════
##
## Term γ_ij [2.5% 97.5%] p-value Sig
## masse_scaled:tibia_scaled -0.0565 -0.4350 0.3113 0.373
TEMPORAL VARIATION OF SELECTION
Survival
S_temporal_w1 <- bind_selection_across_time(
trait = c("masse_scaled", "tibia_scaled"),
data = pheno[pheno$annee_emergence %in% 1991:2024, ],
minimum_sample_size = 5,
fitness_proxy = "w1",
iteration = iteration
)Temporal S — Mass
## y ci_lower ci_upper trait year precision p
## masse_scaled_1991 0.138 0.078 0.213 masse_scaled 1991 26.682 ***
## masse_scaled_1992 0.064 -0.108 0.253 masse_scaled 1992 10.792
## masse_scaled_1993 0.006 -0.014 0.025 masse_scaled 1993 97.647
## masse_scaled_1994 0.239 0.095 0.374 masse_scaled 1994 13.870 ***
## masse_scaled_1995 0.003 -0.208 0.204 masse_scaled 1995 9.705
## masse_scaled_1996 -0.011 -0.286 0.240 masse_scaled 1996 7.050
## masse_scaled_1997 0.197 -0.047 0.432 masse_scaled 1997 8.397 °
## masse_scaled_1998 0.270 0.106 0.438 masse_scaled 1998 11.608 **
## masse_scaled_1999 0.002 -0.005 0.009 masse_scaled 1999 302.428
## masse_scaled_2000 0.612 0.435 0.771 masse_scaled 2000 11.434 ***
## masse_scaled_2001 0.211 -0.357 0.646 masse_scaled 2001 3.901
## masse_scaled_2002 0.219 -0.080 0.516 masse_scaled 2002 6.528 °
## masse_scaled_2003 0.453 0.364 0.542 masse_scaled 2003 21.484 ***
## masse_scaled_2004 -0.008 -0.417 0.422 masse_scaled 2004 4.741
## masse_scaled_2005 0.004 -0.008 0.015 masse_scaled 2005 167.062
## masse_scaled_2006 0.318 0.130 0.516 masse_scaled 2006 10.178 ***
## masse_scaled_2007 0.006 -0.013 0.024 masse_scaled 2007 104.902
## masse_scaled_2008 0.107 0.005 0.199 masse_scaled 2008 19.392 *
## masse_scaled_2009 0.195 -0.156 0.583 masse_scaled 2009 5.271
## masse_scaled_2010 0.166 -0.487 0.554 masse_scaled 2010 3.822
## masse_scaled_2011 0.517 0.373 0.673 masse_scaled 2011 13.375 ***
## masse_scaled_2012 0.136 -0.203 0.475 masse_scaled 2012 5.748
## masse_scaled_2013 -0.004 -0.150 0.126 masse_scaled 2013 14.496
## masse_scaled_2014 0.233 0.052 0.414 masse_scaled 2014 10.443 **
## masse_scaled_2015 -0.033 -0.250 0.210 masse_scaled 2015 8.218
## masse_scaled_2016 0.193 -0.004 0.390 masse_scaled 2016 9.823 *
## masse_scaled_2017 -0.155 -0.477 0.296 masse_scaled 2017 5.080
## masse_scaled_2018 -0.083 -0.196 0.035 masse_scaled 2018 16.754
## masse_scaled_2019 -0.008 -0.033 0.018 masse_scaled 2019 73.740
## masse_scaled_2020 0.071 -0.038 0.161 masse_scaled 2020 19.728 °
## masse_scaled_2021 -0.040 -0.190 0.120 masse_scaled 2021 12.458
## masse_scaled_2022 0.199 0.049 0.348 masse_scaled 2022 13.171 **
## masse_scaled_2023 0.055 -0.106 0.240 masse_scaled 2023 10.935
## masse_scaled_2024 0.058 -0.288 0.403 masse_scaled 2024 5.779Temporal S — Tibia
## y ci_lower ci_upper trait year precision p
## tibia_scaled_1997 0.075 -0.053 0.204 tibia_scaled 1997 15.142
## tibia_scaled_1998 0.102 -0.005 0.204 tibia_scaled 1998 18.880 *
## tibia_scaled_1999 0.006 -0.018 0.032 tibia_scaled 1999 77.250
## tibia_scaled_2000 0.348 0.192 0.489 tibia_scaled 2000 13.339 ***
## tibia_scaled_2001 0.294 0.027 0.550 tibia_scaled 2001 7.466 *
## tibia_scaled_2002 0.385 0.104 0.725 tibia_scaled 2002 6.708 ***
## tibia_scaled_2003 0.379 0.179 0.564 tibia_scaled 2003 9.870 ***
## tibia_scaled_2004 -0.041 -0.251 0.189 tibia_scaled 2004 8.673
## tibia_scaled_2005 -0.007 -0.041 0.023 tibia_scaled 2005 59.099
## tibia_scaled_2006 0.119 -0.026 0.254 tibia_scaled 2006 13.986 °
## tibia_scaled_2007 -0.006 -0.035 0.019 tibia_scaled 2007 71.423
## tibia_scaled_2008 -0.008 -0.138 0.122 tibia_scaled 2008 14.466
## tibia_scaled_2009 0.062 -0.289 0.431 tibia_scaled 2009 5.547
## tibia_scaled_2010 0.256 -0.177 0.623 tibia_scaled 2010 4.544
## tibia_scaled_2011 0.482 0.223 0.770 tibia_scaled 2011 7.016 ***
## tibia_scaled_2012 0.090 -0.251 0.496 tibia_scaled 2012 5.145
## tibia_scaled_2013 0.016 -0.151 0.163 tibia_scaled 2013 12.086
## tibia_scaled_2014 0.302 0.078 0.512 tibia_scaled 2014 8.948 **
## tibia_scaled_2015 0.003 -0.276 0.259 tibia_scaled 2015 7.407
## tibia_scaled_2016 0.263 0.049 0.475 tibia_scaled 2016 8.922 **
## tibia_scaled_2017 -0.206 -0.677 0.439 tibia_scaled 2017 3.443
## tibia_scaled_2018 -0.174 -0.379 0.054 tibia_scaled 2018 9.014
## tibia_scaled_2019 -0.006 -0.035 0.020 tibia_scaled 2019 69.441
## tibia_scaled_2020 0.098 -0.017 0.213 tibia_scaled 2020 16.990 *
## tibia_scaled_2021 0.038 -0.129 0.213 tibia_scaled 2021 11.431
## tibia_scaled_2022 0.256 0.048 0.445 tibia_scaled 2022 9.765 **
## tibia_scaled_2023 -0.141 -0.332 0.113 tibia_scaled 2023 8.821
## tibia_scaled_2024 0.109 -0.256 0.440 tibia_scaled 2024 5.660Test for temporal trend
mod_trend_masse_w1 <- lm(y ~ as.numeric(year), weights = precision,
data = S_temporal_w1[S_temporal_w1$trait == "masse_scaled", ])
cat("── Temporal trend S (Mass) — Survival ──\n")## ── Temporal trend S (Mass) — Survival ──##
## Call:
## lm(formula = y ~ as.numeric(year), data = S_temporal_w1[S_temporal_w1$trait ==
## "masse_scaled", ], weights = precision)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -0.8654 -0.2267 0.1549 0.4487 1.8950
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0533877 0.0422178 1.265 0.215
## as.numeric(year) -0.0001806 0.0025628 -0.070 0.944
##
## Residual standard error: 0.7085 on 32 degrees of freedom
## Multiple R-squared: 0.0001552, Adjusted R-squared: -0.03109
## F-statistic: 0.004966 on 1 and 32 DF, p-value: 0.9443mod_trend_tibia_w1 <- lm(y ~ as.numeric(year), weights = precision,
data = S_temporal_w1[S_temporal_w1$trait == "tibia_scaled", ])
cat("── Temporal trend S (Tibia) — Survival ──\n")## ── Temporal trend S (Tibia) — Survival ──##
## Call:
## lm(formula = y ~ as.numeric(year), data = S_temporal_w1[S_temporal_w1$trait ==
## "tibia_scaled", ], weights = precision)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -0.63918 -0.32424 0.05319 0.62179 1.13881
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.091566 0.062698 1.460 0.156
## as.numeric(year) -0.001881 0.003099 -0.607 0.549
##
## Residual standard error: 0.566 on 26 degrees of freedom
## Multiple R-squared: 0.01398, Adjusted R-squared: -0.02395
## F-statistic: 0.3686 on 1 and 26 DF, p-value: 0.5491No significant temporal trend in selection.
LRS
S_temporal_lrs <- bind_selection_across_time(
trait = c("masse_scaled", "tibia_scaled"),
data = pheno[pheno$annee_emergence %in% 1990:2018, ],
minimum_sample_size = 5,
fitness_proxy = "LRS_relative",
iteration = iteration
)Temporal S — Mass
## y ci_lower ci_upper trait year precision p
## masse_scaled_1991 -0.404 -0.488 -0.327 masse_scaled 1991 23.877
## masse_scaled_1992 0.164 -0.459 0.743 masse_scaled 1992 3.318
## masse_scaled_1993 1.425 -0.559 3.872 masse_scaled 1993 0.869 °
## masse_scaled_1994 0.949 -0.476 3.702 masse_scaled 1994 0.847
## masse_scaled_1995 0.492 -0.376 1.269 masse_scaled 1995 2.353 °
## masse_scaled_1996 1.227 -0.308 3.569 masse_scaled 1996 0.883
## masse_scaled_1997 0.556 -0.176 1.323 masse_scaled 1997 2.565 °
## masse_scaled_1998 0.756 0.069 1.530 masse_scaled 1998 2.591 *
## masse_scaled_1999 0.835 -0.066 2.079 masse_scaled 1999 1.774 °
## masse_scaled_2000 0.671 -0.065 1.823 masse_scaled 2000 2.081 *
## masse_scaled_2001 1.019 -0.190 2.475 masse_scaled 2001 1.473 °
## masse_scaled_2002 0.149 -0.241 0.508 masse_scaled 2002 5.195
## masse_scaled_2003 0.488 -0.350 1.527 masse_scaled 2003 2.013
## masse_scaled_2004 -0.775 -1.979 -0.079 masse_scaled 2004 1.856
## masse_scaled_2005 -0.358 -0.988 0.125 masse_scaled 2005 3.498
## masse_scaled_2006 0.113 -0.459 0.684 masse_scaled 2006 3.420
## masse_scaled_2007 -0.296 -0.606 0.022 masse_scaled 2007 6.222
## masse_scaled_2008 0.186 -0.011 0.393 masse_scaled 2008 9.279 *
## masse_scaled_2009 -0.760 -1.942 0.110 masse_scaled 2009 1.842
## masse_scaled_2010 -0.661 -2.001 0.300 masse_scaled 2010 1.569
## masse_scaled_2011 0.636 -0.067 1.373 masse_scaled 2011 2.685 *
## masse_scaled_2012 0.420 -0.056 0.845 masse_scaled 2012 4.165 *
## masse_scaled_2013 0.344 0.024 0.677 masse_scaled 2013 5.951 *
## masse_scaled_2014 0.473 -0.100 1.121 masse_scaled 2014 3.280 °
## masse_scaled_2015 0.088 -0.495 0.616 masse_scaled 2015 3.516
## masse_scaled_2016 -0.001 -0.730 0.618 masse_scaled 2016 2.834
## masse_scaled_2017 0.361 -0.680 2.086 masse_scaled 2017 1.455
## masse_scaled_2018 -0.243 -0.913 0.419 masse_scaled 2018 2.893Temporal S — Tibia
## y ci_lower ci_upper trait year precision p
## tibia_scaled_1997 0.226 -0.145 0.624 tibia_scaled 1997 4.994
## tibia_scaled_1998 0.280 -0.104 0.746 tibia_scaled 1998 4.414
## tibia_scaled_1999 -0.032 -0.512 0.509 tibia_scaled 1999 3.908
## tibia_scaled_2000 0.626 -0.137 1.656 tibia_scaled 2000 2.168 °
## tibia_scaled_2001 0.591 -0.311 1.548 tibia_scaled 2001 2.109
## tibia_scaled_2002 0.667 0.000 1.402 tibia_scaled 2002 2.800 *
## tibia_scaled_2003 -0.480 -2.219 0.607 tibia_scaled 2003 1.320
## tibia_scaled_2004 -0.351 -1.387 0.436 tibia_scaled 2004 2.000
## tibia_scaled_2005 -0.233 -0.955 0.367 tibia_scaled 2005 2.949
## tibia_scaled_2006 0.005 -0.796 0.579 tibia_scaled 2006 2.649
## tibia_scaled_2007 0.028 -0.701 0.602 tibia_scaled 2007 2.885
## tibia_scaled_2008 -0.272 -0.611 0.018 tibia_scaled 2008 6.223
## tibia_scaled_2009 -1.117 -2.660 0.178 tibia_scaled 2009 1.340
## tibia_scaled_2010 -0.212 -0.823 0.349 tibia_scaled 2010 3.287
## tibia_scaled_2011 0.676 -0.115 1.556 tibia_scaled 2011 2.270 °
## tibia_scaled_2012 0.265 -0.122 0.676 tibia_scaled 2012 4.926 °
## tibia_scaled_2013 0.417 -0.027 0.915 tibia_scaled 2013 4.034 *
## tibia_scaled_2014 0.640 0.063 1.279 tibia_scaled 2014 3.139 *
## tibia_scaled_2015 0.002 -0.716 0.715 tibia_scaled 2015 2.645
## tibia_scaled_2016 -0.529 -1.858 0.494 tibia_scaled 2016 1.553
## tibia_scaled_2017 -0.456 -1.521 0.685 tibia_scaled 2017 1.754
## tibia_scaled_2018 0.067 -0.488 0.639 tibia_scaled 2018 3.449Test for temporal trend
mod_trend_masse_lrs <- lm(y ~ as.numeric(year), weights = precision,
data = S_temporal_lrs[S_temporal_lrs$trait == "masse_scaled", ])
cat("── Temporal trend S (Mass) — LRS ──\n")## ── Temporal trend S (Mass) — LRS ──##
## Call:
## lm(formula = y ~ as.numeric(year), data = S_temporal_lrs[S_temporal_lrs$trait ==
## "masse_scaled", ], weights = precision)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -1.6622 -0.4692 0.4005 0.8909 1.3674
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.074810 0.156259 -0.479 0.636
## as.numeric(year) 0.010973 0.009924 1.106 0.279
##
## Residual standard error: 0.9067 on 26 degrees of freedom
## Multiple R-squared: 0.04492, Adjusted R-squared: 0.008184
## F-statistic: 1.223 on 1 and 26 DF, p-value: 0.2789mod_trend_tibia_lrs <- lm(y ~ as.numeric(year), weights = precision,
data = S_temporal_lrs[S_temporal_lrs$trait == "tibia_scaled", ])
cat("── Temporal trend S (Tibia) — LRS ──\n")## ── Temporal trend S (Tibia) — LRS ──##
## Call:
## lm(formula = y ~ as.numeric(year), data = S_temporal_lrs[S_temporal_lrs$trait ==
## "tibia_scaled", ], weights = precision)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -1.3780 -0.6016 -0.0681 0.6069 1.0833
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.243690 0.242438 1.005 0.327
## as.numeric(year) -0.008963 0.013259 -0.676 0.507
##
## Residual standard error: 0.7016 on 20 degrees of freedom
## Multiple R-squared: 0.02234, Adjusted R-squared: -0.02655
## F-statistic: 0.4569 on 1 and 20 DF, p-value: 0.5068No significant temporal trend in selection.
SELECTION AND CLIMATE CHANGE
CLIMATE VARIABLE AS COVARIATE
PCA AXIS 1 — Aridity & Temperature
Survival
# Differentials
S_clim1_mass_w1 <- compute_selection_climate_covariable(
data = pheno, trait = "masse_scaled",
fitness_proxy = "w1", climate_variable = "spring_index_1_y", iteration = iteration)
S_clim1_tibia_w1 <- compute_selection_climate_covariable(
data = pheno, trait = "tibia_scaled",
fitness_proxy = "w1", climate_variable = "spring_index_1_y", iteration = iteration)
# Gradients (multivariate model)
S_clim1_gradient_w1 <- compute_selection_climate_covariable(
data = pheno, trait = c("masse_scaled", "tibia_scaled"),
fitness_proxy = "w1", climate_variable = "spring_index_1_y", iteration = iteration)Differential S — Mass
##
## ═════════════════════════════════════════════════════════════════
## S (Mass) × PCA Axis 1 — Proxy: Survival
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.1160 0.0820 0.1500 0.000 ***
## spring_index_1_y -0.0080 -0.0110 -0.0060 0.000 ***
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled:spring_index_1_y 0.0240 -0.0010 0.0480 0.028 *
##
## ── Quadratic ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.1290 0.0930 0.1640 0.000 ***
## spring_index_1_y -0.0100 -0.0130 -0.0070 0.000 ***
## I(0.5 * masse_scaled^2) -0.0470 -0.0950 0.0040 0.035 *
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled:spring_index_1_y 0.0270 0.0030 0.0510 0.020 *Differential S — Tibia
##
## ═════════════════════════════════════════════════════════════════
## S (Tibia) × PCA Axis 1 — Proxy: Survival
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled 0.0960 0.0660 0.1270 0.000 ***
## spring_index_1_y 0.0030 0.0000 0.0070 0.020 *
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled:spring_index_1_y 0.0240 -0.0080 0.0570 0.068 °
##
## ── Quadratic ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled 0.0980 0.0650 0.1320 0.000 ***
## spring_index_1_y 0.0020 -0.0020 0.0070 0.148
## I(0.5 * tibia_scaled^2) -0.0090 -0.0490 0.0340 0.333
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled:spring_index_1_y 0.0240 -0.0080 0.0570 0.073 °Gradient β — Mass
print_selection_climate(
S_clim1_gradient_w1 |> filter(grepl("masse|Intercept", params)),
title = "β gradient (Mass) × PCA Axis 1 — Proxy: Survival"
)##
## ═════════════════════════════════════════════════════════════════
## β gradient (Mass) × PCA Axis 1 — Proxy: Survival
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.0920 0.0460 0.1400 0.000 ***
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled:spring_index_1_y 0.0580 0.0230 0.0930 0.001 ***
##
## ── Quadratic ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.1120 0.0620 0.1620 0.000 ***
## I(0.5 * masse_scaled^2) -0.0530 -0.1180 0.0110 0.054 °
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled:spring_index_1_y 0.0600 0.0250 0.0960 0.001 ***Gradient β — Tibia
print_selection_climate(
S_clim1_gradient_w1 |> filter(grepl("tibia|Intercept", params)),
title = "β gradient (Tibia) × PCA Axis 1 — Proxy: Survival"
)##
## ═════════════════════════════════════════════════════════════════
## β gradient (Tibia) × PCA Axis 1 — Proxy: Survival
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled 0.0310 -0.0040 0.0660 0.041 *
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled:spring_index_1_y -0.0290 -0.0700 0.0120 0.088 °
##
## ── Quadratic ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled 0.0190 -0.0220 0.0610 0.194
## I(0.5 * tibia_scaled^2) 0.0100 -0.0370 0.0530 0.340
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled:spring_index_1_y -0.0270 -0.0680 0.0160 0.105plot_survival_s_mass_1 <- plot_selection(
mod_result = S_clim1_mass_w1, type_selection = "Quadratic",
ylab = "w (survival)", xlab = "Scaled mass")
plot_survival_gradient_mass_1 <- plot_selection(
mod_result = S_clim1_gradient_w1 |> filter(!grepl("tibia", params)),
type_selection = "Quadratic",
ylab = "w (survival)", xlab = "Scaled mass",
title = "Aridity & Temperature")
LRS
# Differentials
S_clim1_mass_lrs <- compute_selection_climate_covariable(
data = data_lrs, trait = "masse_scaled",
fitness_proxy = "LRS_relative", climate_variable = "spring_index_1_y", iteration = iteration)
S_clim1_tibia_lrs <- compute_selection_climate_covariable(
data = data_lrs, trait = "tibia_scaled",
fitness_proxy = "LRS_relative", climate_variable = "spring_index_1_y", iteration = iteration)
# Gradients
S_clim1_gradient_lrs <- compute_selection_climate_covariable(
data = data_lrs, trait = c("masse_scaled", "tibia_scaled"),
fitness_proxy = "LRS_relative", climate_variable = "spring_index_1_y", iteration = iteration)Differential S — Mass
##
## ═════════════════════════════════════════════════════════════════
## S (Mass) × PCA Axis 1 — Proxy: LRS
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.2740 0.1560 0.4040 0.000 ***
## spring_index_1_y 0.0030 -0.0170 0.0210 0.376
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled:spring_index_1_y -0.1290 -0.2170 -0.0290 0.007 **
##
## ── Quadratic ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.2980 0.1620 0.4250 0.000 ***
## spring_index_1_y -0.0010 -0.0240 0.0210 0.471
## I(0.5 * masse_scaled^2) -0.0690 -0.2730 0.1470 0.236
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled:spring_index_1_y -0.1270 -0.2160 -0.0240 0.008 **Differential S — Tibia
##
## ═════════════════════════════════════════════════════════════════
## S (Tibia) × PCA Axis 1 — Proxy: LRS
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled 0.0710 -0.0230 0.1700 0.061 °
## spring_index_1_y 0.0080 -0.0010 0.0180 0.028 *
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled:spring_index_1_y -0.1100 -0.2130 0.0150 0.044 *
##
## ── Quadratic ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled 0.0440 -0.0730 0.1710 0.221
## spring_index_1_y 0.0190 -0.0000 0.0370 0.027 *
## I(0.5 * tibia_scaled^2) 0.0960 -0.0720 0.2690 0.131
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled:spring_index_1_y -0.1000 -0.1980 0.0160 0.045 *Gradient β — Mass
print_selection_climate(
S_clim1_gradient_lrs |> filter(grepl("masse|Intercept", params)),
title = "β gradient (Mass) × PCA Axis 1 — Proxy: LRS"
)##
## ═════════════════════════════════════════════════════════════════
## β gradient (Mass) × PCA Axis 1 — Proxy: LRS
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.2580 0.1080 0.4140 0.001 ***
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled:spring_index_1_y -0.0610 -0.1880 0.0520 0.152
##
## ── Quadratic ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.3720 0.2060 0.5280 0.000 ***
## I(0.5 * masse_scaled^2) -0.2810 -0.4840 -0.0270 0.022 *
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled:spring_index_1_y -0.0650 -0.1880 0.0460 0.132Gradient β — Tibia
print_selection_climate(
S_clim1_gradient_lrs |> filter(grepl("tibia|Intercept", params)),
title = "β gradient (Tibia) × PCA Axis 1 — Proxy: LRS"
)##
## ═════════════════════════════════════════════════════════════════
## β gradient (Tibia) × PCA Axis 1 — Proxy: LRS
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled -0.0660 -0.1830 0.0470 0.130
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled:spring_index_1_y -0.0840 -0.2100 0.0600 0.102
##
## ── Quadratic ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled -0.1530 -0.3060 -0.0030 0.022 *
## I(0.5 * tibia_scaled^2) 0.1760 -0.0100 0.3530 0.031 *
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled:spring_index_1_y -0.0600 -0.1840 0.0770 0.170plot_lrs_s_mass_1 <- plot_selection(
mod_result = S_clim1_mass_lrs, type_selection = "Linear",
ylab = "w (LRS)", xlab = "Scaled mass")
plot_lrs_gradient_mass_1 <- plot_selection(
mod_result = S_clim1_gradient_lrs |> filter(!grepl("tibia", params)),
type_selection = "Quadratic",
ylab = "w (LRS)", xlab = "Scaled mass")
PCA AXIS 2 — NDVI
Survival
# Differentials
S_clim2_mass_w1 <- compute_selection_climate_covariable(
data = pheno, trait = "masse_scaled",
fitness_proxy = "w1", climate_variable = "spring_index_2_y", iteration = iteration)
S_clim2_tibia_w1 <- compute_selection_climate_covariable(
data = pheno, trait = "tibia_scaled",
fitness_proxy = "w1", climate_variable = "spring_index_2_y", iteration = iteration)
# Gradients
S_clim2_gradient_w1 <- compute_selection_climate_covariable(
data = pheno, trait = c("masse_scaled", "tibia_scaled"),
fitness_proxy = "w1", climate_variable = "spring_index_2_y", iteration = iteration)Differential S — Mass
##
## ═════════════════════════════════════════════════════════════════
## S (Mass) × PCA Axis 2 — Proxy: Survival
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.1300 0.0960 0.1650 0.000 ***
## spring_index_2_y -0.0360 -0.0460 -0.0270 0.000 ***
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled:spring_index_2_y 0.0430 0.0050 0.0810 0.012 *
##
## ── Quadratic ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.1480 0.1130 0.1830 0.000 ***
## spring_index_2_y -0.0340 -0.0450 -0.0250 0.000 ***
## I(0.5 * masse_scaled^2) -0.0610 -0.1130 -0.0090 0.014 *
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled:spring_index_2_y 0.0580 0.0160 0.0980 0.002 **Differential S — Tibia
##
## ═════════════════════════════════════════════════════════════════
## S (Tibia) × PCA Axis 2 — Proxy: Survival
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled 0.1160 0.0800 0.1520 0.000 ***
## spring_index_2_y -0.0620 -0.0800 -0.0460 0.000 ***
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled:spring_index_2_y 0.0250 -0.0040 0.0530 0.050 *
##
## ── Quadratic ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled 0.1170 0.0810 0.1550 0.000 ***
## spring_index_2_y -0.0610 -0.0780 -0.0420 0.000 ***
## I(0.5 * tibia_scaled^2) -0.0190 -0.0940 0.0530 0.304
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled:spring_index_2_y 0.0330 -0.0130 0.0820 0.080 °Gradient β — Mass
print_selection_climate(
S_clim2_gradient_w1 |> filter(grepl("masse|Intercept", params)),
title = "β gradient (Mass) × PCA Axis 2 — Proxy: Survival"
)##
## ═════════════════════════════════════════════════════════════════
## β gradient (Mass) × PCA Axis 2 — Proxy: Survival
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.1240 0.0780 0.1690 0.000 ***
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled:spring_index_2_y 0.0650 -0.0040 0.1360 0.037 *
##
## ── Quadratic ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.1520 0.1040 0.2010 0.000 ***
## I(0.5 * masse_scaled^2) -0.0550 -0.1210 0.0110 0.051 °
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled:spring_index_2_y 0.0770 0.0050 0.1500 0.018 *Gradient β — Tibia
print_selection_climate(
S_clim2_gradient_w1 |> filter(grepl("tibia|Intercept", params)),
title = "β gradient (Tibia) × PCA Axis 2 — Proxy: Survival"
)##
## ═════════════════════════════════════════════════════════════════
## β gradient (Tibia) × PCA Axis 2 — Proxy: Survival
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled 0.0360 -0.0080 0.0770 0.061 °
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled:spring_index_2_y 0.0010 -0.0500 0.0510 0.493
##
## ── Quadratic ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled 0.0200 -0.0260 0.0660 0.198
## I(0.5 * tibia_scaled^2) -0.0150 -0.0850 0.0520 0.328
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled:spring_index_2_y 0.0090 -0.0500 0.0720 0.387plot_survival_s_mass_2 <- plot_selection(
mod_result = S_clim2_mass_w1, type_selection = "Quadratic",
ylab = "w (survival)", xlab = "Scaled mass")
plot_survival_gradient_mass_2 <- plot_selection(
mod_result = S_clim2_gradient_w1 |> filter(!grepl("tibia", params)),
type_selection = "Quadratic",
ylab = "w (survival)", xlab = "Scaled mass", title = "NDVI")
LRS
# Differentials
S_clim2_mass_lrs <- compute_selection_climate_covariable(
data = data_lrs, trait = "masse_scaled",
fitness_proxy = "LRS_relative", climate_variable = "spring_index_2_y", iteration = iteration)
S_clim2_tibia_lrs <- compute_selection_climate_covariable(
data = data_lrs, trait = "tibia_scaled",
fitness_proxy = "LRS_relative", climate_variable = "spring_index_2_y", iteration = iteration)
# Gradients
S_clim2_gradient_lrs <- compute_selection_climate_covariable(
data = data_lrs, trait = c("masse_scaled", "tibia_scaled"),
fitness_proxy = "LRS_relative", climate_variable = "spring_index_2_y", iteration = iteration)Differential S — Mass
##
## ═════════════════════════════════════════════════════════════════
## S (Mass) × PCA Axis 2 — Proxy: LRS
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.2430 0.1220 0.3630 0.001 ***
## spring_index_2_y -0.0770 -0.1140 -0.0420 0.000 ***
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled:spring_index_2_y 0.2500 0.1320 0.3650 0.000 ***
##
## ── Quadratic ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.2960 0.1650 0.4240 0.001 ***
## spring_index_2_y -0.0730 -0.1110 -0.0370 0.000 ***
## I(0.5 * masse_scaled^2) -0.1680 -0.3520 0.0420 0.048 *
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled:spring_index_2_y 0.2780 0.1560 0.3810 0.000 ***Differential S — Tibia
##
## ═════════════════════════════════════════════════════════════════
## S (Tibia) × PCA Axis 2 — Proxy: LRS
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled 0.0660 -0.0740 0.2060 0.183
## spring_index_2_y -0.0990 -0.1640 -0.0360 0.000 ***
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled:spring_index_2_y 0.1100 0.0200 0.2030 0.006 **
##
## ── Quadratic ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled 0.0650 -0.0750 0.2030 0.182
## spring_index_2_y -0.1090 -0.1830 -0.0360 0.002 **
## I(0.5 * tibia_scaled^2) 0.0790 -0.1850 0.3470 0.279
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled:spring_index_2_y 0.0710 -0.0870 0.2210 0.179Gradient β — Mass
print_selection_climate(
S_clim2_gradient_lrs |> filter(grepl("masse|Intercept", params)),
title = "β gradient (Mass) × PCA Axis 2 — Proxy: LRS"
)##
## ═════════════════════════════════════════════════════════════════
## β gradient (Mass) × PCA Axis 2 — Proxy: LRS
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.2420 0.0770 0.4000 0.001 ***
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled:spring_index_2_y 0.2200 -0.0240 0.4660 0.040 *
##
## ── Quadratic ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled 0.3690 0.1870 0.5430 0.000 ***
## I(0.5 * masse_scaled^2) -0.3350 -0.5190 -0.0760 0.008 **
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## masse_scaled:spring_index_2_y 0.2640 0.0240 0.4950 0.018 *Gradient β — Tibia
print_selection_climate(
S_clim2_gradient_lrs |> filter(grepl("tibia|Intercept", params)),
title = "β gradient (Tibia) × PCA Axis 2 — Proxy: LRS"
)##
## ═════════════════════════════════════════════════════════════════
## β gradient (Tibia) × PCA Axis 2 — Proxy: LRS
## ═════════════════════════════════════════════════════════════════
##
## ── Linear ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled -0.0820 -0.2600 0.0840 0.176
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled:spring_index_2_y 0.0070 -0.1530 0.1840 0.479
##
## ── Quadratic ──
## Main effects:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled -0.1420 -0.3280 0.0330 0.063 °
## I(0.5 * tibia_scaled^2) 0.1300 -0.1740 0.3880 0.172
##
## Trait × climate interactions:
## Parameter Estimate [2.5% 97.5%] p-value Sig
## tibia_scaled:spring_index_2_y -0.0500 -0.2680 0.1790 0.325plot_lrs_s_mass_2 <- plot_selection(
mod_result = S_clim2_mass_lrs, type_selection = "Linear",
ylab = "w (LRS)", xlab = "Scaled mass")
plot_lrs_gradient_mass_2 <- plot_selection(
mod_result = S_clim2_gradient_lrs |> filter(!grepl("tibia", params)),
type_selection = "Quadratic",
ylab = "w (LRS)", xlab = "Scaled mass")
FINAL FIGURES
Selection differentials
gridExtra::grid.arrange(
plot_survival_s_mass_1,
plot_lrs_s_mass_1,
plot_survival_s_mass_2,
plot_lrs_s_mass_2,
ncol = 2
)
