Introduction
This report presents an analysis of the relationship between the heat index and the number of psychiatric admissions, using distributed lag non-linear models (DLNM). The heat index is a measure that combines air temperature and relative humidity to determine an apparent temperature.
Method
Data
The data used in this analysis includes daily records of the heat index, temperature, and the number of “passages” (psychiatric admissions ?) in CPOA. The heat index is calculated using the formula provided by the National Weather Service.
En vrai le Heat.index est tellement corrélé à la température, qu’utiliser l’un ou l’autre pour nous est égal.
##
## Pearson's product-moment correlation
##
## data: B$TNTXM and B$HEAT.INDEX
## t = 468.8, df = 1326, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9966559 0.9973030
## sample estimates:
## cor
## 0.9969968# plot de Heat_index (MA30) et températures (MA30) et températures
library(TTR)
# Calculate moving averages
B$TNTXM_MA7 <- SMA(B$TNTXM, n = 7)
B$HI_MA7 <- SMA(B$HEAT.INDEX, n = 7)
# Plot the original data and the moving average
plot(B$DATE, B$TNTXM, col = "lightblue", xlab = "Date", ylab = "Temperature",
main = "Temperature and Heat Index Over Time with 7-day Moving Averages")
lines(B$DATE, B$TNTXM_MA7, col = "#1611a5")
lines(B$DATE,B$HI_MA7, col ="green")
legend("topright",
legend = c("Raw Temperature", "Temperature MA(7)", "Heat Index MA(7)"),
col = c("lightblue", "#1611a5", "green"),
lty = c(1,1,1))
# Plot of the difference between Temperature and Heat Index
plot(B$DATE, SMA(B$HEAT.INDEX-B$TNTXM, n = 7), type="l",
main="Difference between Heat Index and Temperature (MA7)",
ylab="Difference (Heat Index - Temperature)",
xlab="Date")
abline(h=0, col="red")
#plot histogramme of HEAT.INDEX with values of median and mean.
hist(B$HEAT.INDEX, breaks = 100, main = "Histogram of Heat Index", xlab = "Heat Index", ylab = "Frequency")
abline(v = median(B$HEAT.INDEX), col = "red", lwd = 2, lty = 2)
abline(v = mean(B$HEAT.INDEX), col = "blue", lwd = 2, lty = 2)
text(10,70, labels = round(median(B$HEAT.INDEX), 2), col = "red", pos = 3, cex = 1.5)
text(17,70, labels = round(mean(B$HEAT.INDEX), 2), col = "blue", pos = 3, cex = 1.5)
legend("topright", legend = c("Median", "Mean"), col = c("red", "blue"), lwd = 2, lty = 2)
Model Specification
A distributed lag non-linear model (DLNM) was used to assess the effect of the heat index on the number of passages. The model includes a crossbasis function for the heat index with a maximum lag of 2 days, as recommended by previous studies. Natural splines were used for both the exposure-response and lag-response functions. After testing several specifications, we choose 2 degrees of freedom for both the exposure-response. For exposure-reponse, higher degrees seemed injustified as previous studies showed quasi linear relationship between temperatures and suicides (Lehmann) for examples. Higher degrees were tested in sensitivity analysis but did not improve the model in terms of AIC. Poisson regression model adjusted on time, day of week, bank holidays. According to Kwac et al.: > Based on previous studies, a maximum lag of 2 days was chosen for suicides and accidents, and the maximum lag was set to 21 days for other causes of death. This choice was examined in sensitivity analyses. The (time) lagged effects were modeled through a natural cubic spline with an intercept and 3 knots equally spaced on the logarithm scale (a total of 5 degrees of freedom), except for suicide, for which only 2 knots were required because of the short lag considered. Region-specific estimates were then summed over all lags to obtain region-specific lag-cumulative estimates.
cb_HI <- crossbasis(B$HEAT.INDEX,
lag = 2,
argvar = list(fun = "ns", df = 2),
arglag = list(fun = "ns", df = 2))
model_HI <- glm(nombre_passage ~ cb_HI
+weekday
+holiday
+ns(DATE,8*4)
+ns(B$Year,2),
,data=B,
family = poisson(link="log"))
aic_value <- AIC(model_HI)
print(aic_value)## [1] 7889.064Model Fit
The model’s AIC value is 7889.0637928.
Results : passages
Heat Index and Number of Passages
The following plot shows the cumulative relative risk (RR) of the number of passages as a function of the heat index.
pred.hi <- crosspred(cb_HI, model_HI, at= -5:35, cumul=TRUE, cen=-1)
plot(pred.hi, "overall",
col=4, cumul=TRUE,
xlab="Heat Index",
ylab="Cumulative RR",
main="Cumulative Relative Risk of Psychiatric Admissions by Heat Index")
Comment: The plot indicates that the relative risk increases with
higher heat index values, suggesting a potential impact of heat on the
number of passages. Time series analyses adjusted for time
(natural cubic spline with 8 degrees of freedom per year), day of week,
bank holidays
RRs are either expressed as the relative risk of suicide at the 99th temperature percentile versus the 1st percentile and for a +1°C increase above 1st percentile.
# 1. Define percentiles of interest
pct_1 <- quantile(B$HEAT.INDEX, probs = 0.01, na.rm = TRUE)
pct_99 <- quantile(B$HEAT.INDEX, probs = 0.99, na.rm = TRUE)
# 2. Create a "crosspred" object
# 'cen' is the reference (center) value for the exposure;
# if you want to interpret everything relative to the 1st percentile, set cen = pct_1.
pred_HI <- crosspred(cb_HI,
model_HI,
cen = pct_1, # or pick another reference, e.g. the median
by = 0.1, # step to evaluate the exposure
cumul = TRUE) # for cumulative effect across lags
# The "by=0.1" means the function evaluated the exposure at increments of 0.1.
# We want the row in pred_HI that is closest to pct_99:
idx_99 <- which.min(abs(pred_HI$predvar - pct_99))
RR_99 <- pred_HI$allRRfit [idx_99]
RR_99_low <- pred_HI$allRRlow [idx_99]
RR_99_hi <- pred_HI$allRRhigh[idx_99]
cat("RR(99th vs. 1st pct) =", RR_99, "(", RR_99_low, "-", RR_99_hi, ")\n")## RR(99th vs. 1st pct) = 1.157126 ( 1.02477 - 1.306577 )The results show that the relative risk (RR) at the 99th percentile (HI = 29°C) compared to the 1st percentile (HI = -1°C)is 1.157126 (1.0247701 - 1.3065767). This indicates that the relative risk of passages increases significantly with the increase in the heat index.
temp_x <- -1 # example reference temperature
temp_xp1 <- 0
idx_x <- which.min(abs(pred_HI$predvar - temp_x))
idx_xp1 <- which.min(abs(pred_HI$predvar - temp_xp1))
RR_x <- pred_HI$allRRfit[idx_x]
RR_xp1 <- pred_HI$allRRfit[idx_xp1]
RR_1C <- RR_xp1 / RR_x
RR_1C_low <- pred_HI$allRRlow[idx_xp1] / pred_HI$allRRhigh[idx_x]
RR_1C_hi <- pred_HI$allRRhigh[idx_xp1] / pred_HI$allRRlow[idx_x]
cat("RR for +1 HI from", temp_x, "to", temp_xp1, ":", RR_1C, "(", RR_1C_low, "-", RR_1C_hi, ")\n")## RR for +1 HI from -1 to 0 : 1.011402 ( 1.00336 - 1.019508 )For an increase of +1 in the heat index from the -1 value, the RR is 1.0114019 (1.0033603 - 1.0195079).
Analyse de sensibilité
L’analyse de sensibilité explore l’effet de la variation des degrés de liberté dans les fonctions splines pour l’exposition-réponse (var_df) et la réponse en décalage (lag_df), ainsi que l’effet du maximum lag. Les résultats sont présentés dans une grille de graphiques 3x4 et un tableau des AIC.
# Définition des valeurs de df pour var et lag
var_df_values <- c(1, 2, 3)
lag_df_values <- c(2, 3)
lag_max_values <- c(2, 3)
# Création d'une grille de graphiques 3x4
par(mfrow = c(3, 4), mar = c(4, 4, 3, 1))
# Stockage des AIC pour tous les modèles
aic_matrix <- matrix(NA, nrow = length(var_df_values), ncol = length(lag_df_values) * length(lag_max_values))
rownames(aic_matrix) <- paste("var_df =", var_df_values)
colnames(aic_matrix) <- paste("lag_df =", rep(lag_df_values, each = length(lag_max_values)), ", lag_max =", rep(lag_max_values, times = length(lag_df_values)))
# Boucle sur toutes les combinaisons
for (i in seq_along(var_df_values)) {
for (j in seq_along(lag_df_values)) {
for (k in seq_along(lag_max_values)) {
var_df <- var_df_values[i]
lag_df <- lag_df_values[j]
lag_max <- lag_max_values[k]
# Calcul de l'index pour la matrice AIC
col_idx <- (j - 1) * length(lag_max_values) + k
# Création du crossbasis
cb.heat_index <- crossbasis(B$HEAT.INDEX,
lag = lag_max,
argvar = list(fun = "ns", df = var_df),
arglag = list(fun = "ns", df = lag_df))
# Modèle
model_hi <- glm(nombre_hospitalisation ~ cb.heat_index
+ weekday
+ holiday
+ ns(DATE, 8*4)
+ ns(B$Year, 2),
data = B,
family = poisson(link = "log"))
# Stockage de l'AIC
aic_matrix[i, col_idx] <- AIC(model_hi)
# Prédiction et graphique
pred.hi <- crosspred(cb.heat_index, model_hi, at = seq(-5, 35, by = 0.1), cumul = TRUE, cen = -1)
plot(pred.hi$predvar, pred.hi$allRRfit, type = "l", col = "red", lwd = 2,
xlab = "Heat Index", ylab = "Cumulative RR",
main = paste("var_df =", var_df, ", lag_df =", lag_df, "\nlag_max =", lag_max),
ylim = c(0.8, 1.5), cex.main = 0.8)
lines(pred.hi$predvar, pred.hi$allRRlow, col = "red", lwd = 1, lty = 2)
lines(pred.hi$predvar, pred.hi$allRRhigh, col = "red", lwd = 1, lty = 2)
abline(h = 1, col = "black", lwd = 1)
}
}
}
# Retour à la disposition normale
par(mfrow = c(1, 1))
# Affichage du tableau des AIC
cat("Tableau des AIC pour l'analyse de sensibilité (Hospitalisations):\n")## Tableau des AIC pour l'analyse de sensibilité (Hospitalisations):## lag_df = 2 , lag_max = 2 lag_df = 2 , lag_max = 3
## var_df = 1 6876.42 6871.44
## var_df = 2 6879.74 6872.84
## var_df = 3 6883.60 6876.69
## lag_df = 3 , lag_max = 2 lag_df = 3 , lag_max = 3
## var_df = 1 6877.46 6873.44
## var_df = 2 6881.60 6875.58
## var_df = 3 6887.47 6880.96# Identification du meilleur modèle
best_idx <- which(aic_matrix == min(aic_matrix), arr.ind = TRUE)
best_var_df <- var_df_values[best_idx[1]]
best_col <- best_idx[2]
best_lag_df_idx <- ceiling(best_col / length(lag_max_values))
best_lag_max_idx <- best_col %% length(lag_max_values)
if (best_lag_max_idx == 0) best_lag_max_idx <- length(lag_max_values)
best_lag_df <- lag_df_values[best_lag_df_idx]
best_lag_max <- lag_max_values[best_lag_max_idx]
cat("\nMeilleur modèle: var_df =", best_var_df, ", lag_df =", best_lag_df, ", lag_max =", best_lag_max, "\n")##
## Meilleur modèle: var_df = 1 , lag_df = 2 , lag_max = 3## AIC minimum: 6871.44HOMMES / FEMMES
# Modèles pour hommes, femmes, et ensemble
cb_HI <- crossbasis(B$HEAT.INDEX, lag = 2, argvar = list(fun = "ns", df = 2), arglag = list(fun = "ns", df = 2))
model_HI_H <- glm(nombre_homme ~ cb_HI + weekday + holiday + ns(DATE, 8*4) + ns(B$Year, 2), data = B, family = poisson(link = "log"))
model_HI_F <- glm(nombre_femme ~ cb_HI + weekday + holiday + ns(DATE, 8*4) + ns(B$Year, 2), data = B, family = poisson(link = "log"))
model_HI_T <- glm(nombre_passage ~ cb_HI + weekday + holiday + ns(DATE, 8*4) + ns(B$Year, 2), data = B, family = poisson(link = "log"))
# Prédictions
pred_HI_H <- crosspred(cb_HI, model_HI_H, at = seq(-5, 35, by = 1),cen = -1, cumul = TRUE)
pred_HI_F <- crosspred(cb_HI, model_HI_F, at = seq(-5, 35, by = 1),cen = -1, cumul = TRUE)
pred_HI_T <- crosspred(cb_HI, model_HI_T, at = seq(-5, 35, by = 1),cen = -1, cumul = TRUE)# Reserve space in the bottom margin:
# Here I used oma = c(3, 0, 2, 0) so we have 3 lines for the bottom.
par(mfrow = c(1, 3), oma = c(3, 0, 2, 0))
plot(pred_HI_H, "overall", col = "#1ebc12", xlab = "Heat Index",
ylab = "Relative risk", main = "Men", ylim = c(0.8, 1.4))
points(29, pred_HI_H$allRRfit[which.min(abs(pred_HI_H$predvar - 29))],
col = "black", pch = 19)
text(29, pred_HI_H$allRRfit[which.min(abs(pred_HI_H$predvar - 29))] + 0.02,
labels = paste("RR =", round(pred_HI_H$allRRfit[which.min(abs(pred_HI_H$predvar - 29))], 2)),
col = "black")
plot(pred_HI_F, "overall", col = "#d58f0d", xlab = "Heat Index",
ylab = "Relative risk", main = "Women", ylim = c(0.8, 1.4))
points(29, pred_HI_F$allRRfit[which.min(abs(pred_HI_F$predvar - 29))],
col = "black", pch = 19)
text(29, pred_HI_F$allRRfit[which.min(abs(pred_HI_H$predvar - 29))],
labels = paste("RR =", round(pred_HI_F$allRRfit[which.min(abs(pred_HI_H$predvar - 29))], 2)),
col = "black")
plot(pred_HI_T, "overall", col = "#3f4040", xlab = "Heat Index",
ylab = "Relative risk", main = "Men and women together", ylim = c(0.8, 1.4))
points(29, pred_HI_T$allRRfit[which.min(abs(pred_HI_T$predvar - 29))],
col = "black", pch = 19)
text(29, pred_HI_T$allRRfit[which.min(abs(pred_HI_H$predvar - 29))],
labels = paste("RR =", round(pred_HI_T$allRRfit[which.min(abs(pred_HI_H$predvar - 29))], 2)),
col = "black")
# Overall title (top margin)
mtext("Overall Cumulative Risks for Psychiatric Admissions",
outer = TRUE, cex = 1.2)
# Text in the bottom margin (side=1, outer=TRUE for bottom outer margin)
mtext("Cumulative Risks at the 99th percentile (HI = 29) compared to 1st percentile (-1)",color="blue",
side = 1, outer = TRUE, line = 1, cex = 0.8, adj = 0.5)## Warning in mtext("Cumulative Risks at the 99th percentile (HI = 29) compared to
## 1st percentile (-1)", : "color" n'est pas un paramètre graphique
Autres Types de Graphes
pred.hi <- crosspred(cb_HI, model_HI, by=1, cen=5)
plot(pred.hi, xlab = "Heat index", theta = 240, phi = 40, ltheta = -185, zlab = "RR", main = "3D graph")
plot(pred.hi, "contour", plot.title = title(xlab = "Heat index", ylab = "Lag", main = "Contour graph"), key.title = title("RR"))
Analysis of Hospitalizations
Evolution of Hospitalizations and Ratio
Before diving into the model analysis, let’s examine the temporal evolution of hospitalizations and their ratio to total admissions.
# Création des moyennes mobiles pour lisser les données
B$MA7_passages <- SMA(B$nombre_passage, n = 7)
B$MA7_hospis <- SMA(B$nombre_hospitalisation, n = 7)
B$ratio_hospis_passages <- B$nombre_hospitalisation / B$nombre_passage
B$MA7_ratio <- SMA(B$ratio_hospis_passages, n = 7)
# Premier graphique : Évolution des passages et hospitalisations
p1 <- ggplot(B, aes(x = DATE)) +
geom_line(aes(y = MA7_passages, color = "Passages"), size = 1) +
geom_line(aes(y = MA7_hospis, color = "Hospitalisations"), size = 1) +
scale_color_manual(values = c("Passages" = "#1ebc12", "Hospitalisations" = "#d58f0d")) +
theme_minimal() +
labs(
title = "Évolution du nombre de passages et d'hospitalisations au CPOA",
subtitle = "Moyenne mobile sur 7 jours",
x = "Date",
y = "Nombre",
color = "Type"
) +
theme(
legend.position = "bottom",
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 12),
axis.title = element_text(size = 12),
legend.title = element_text(size = 12),
legend.text = element_text(size = 11)
)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.# Second graphique : Évolution du ratio hospitalisations/passages
p2 <- ggplot(B, aes(x = DATE)) +
geom_line(aes(y = MA7_ratio), color = "#3f4040", size = 1) +
theme_minimal() +
labs(
title = "Évolution du ratio hospitalisations/passages",
subtitle = "Moyenne mobile sur 7 jours",
x = "Date",
y = "Ratio (hospitalisations/passages)"
) +
theme(
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 12),
axis.title = element_text(size = 12)
)
# Création d'une grille avec les deux graphiques
grid.arrange(p1, p2, ncol = 1)## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).
Comment: The top graph shows the evolution of both admissions and hospitalizations over time, using a 7-day moving average to smooth the data. The bottom graph displays the ratio of hospitalizations to admissions, which can help identify periods where the proportion of patients requiring hospitalization changed.
Model Specification for Hospitalizations
We applied the same DLNM approach to analyze the relationship between heat index and the number of hospitalizations. The model specification remains consistent with the previous analysis, using a maximum lag of 2 days and natural splines for both exposure-response and lag-response functions.
# Création du crossbasis pour le Heat Index
cb_HI_hospis <- crossbasis(B$HEAT.INDEX,
lag = 2,
argvar = list(fun = "ns", df = 2),
arglag = list(fun = "ns", df = 2))
# Modèle pour nombre_hospitalisation
model_HI_hospis <- glm(nombre_hospitalisation ~ cb_HI_hospis
+ weekday
+ holiday
+ ns(DATE,8*4)
+ ns(B$Year,2),
data = B,
family = poisson(link = "log"))
# Calcul de l'AIC
aic_value_hospis <- AIC(model_HI_hospis)
print(paste("AIC du modèle hospitalisations:", aic_value_hospis))## [1] "AIC du modèle hospitalisations: 6879.74206400927"Results for Hospitalizations
The following plot shows the cumulative relative risk (RR) of hospitalizations as a function of the heat index.
pred.hi.hospis <- crosspred(cb_HI_hospis, model_HI_hospis, at= -5:35, cumul=TRUE, cen=-1)
plot(pred.hi.hospis, "overall",
col=4, cumul=TRUE,
xlab="Heat Index",
ylab="Cumulative RR",
main="Cumulative Relative Risk of Hospitalizations by Heat Index")
Comment: The plot shows the relationship between heat index and hospitalization risk, allowing for comparison with the previous analysis of total admissions.
Results : Hospitalizations
Heat Index and Number of Hospitalizations
The following plot shows the cumulative relative risk (RR) of hospitalizations as a function of the heat index.
# Création du crossbasis pour le Heat Index
cb_HI_hospis <- crossbasis(B$HEAT.INDEX,
lag = 2,
argvar = list(fun = "ns", df = 2),
arglag = list(fun = "ns", df = 2))
# Modèle pour nombre_hospitalisation
model_HI_hospis <- glm(nombre_hospitalisation ~ cb_HI_hospis
+ weekday
+ holiday
+ ns(DATE,8*4)
+ ns(B$Year,2),
data = B,
family = poisson(link = "log"))
pred.hi.hospis <- crosspred(cb_HI_hospis, model_HI_hospis, at= -5:35, by=0.1, cumul=TRUE, cen=-1)
plot(pred.hi.hospis, "overall",
col=4, cumul=TRUE,
xlab="Heat Index",
ylab="Cumulative RR",
main="Cumulative Relative Risk of Hospitalizations by Heat Index")
Comment: The plot shows the relationship between heat index and hospitalization risk, allowing for comparison with the previous analysis of total admissions.
# 1. Define percentiles of interest
pct_1 <- quantile(B$HEAT.INDEX, probs = 0.01, na.rm = TRUE)
pct_99 <- quantile(B$HEAT.INDEX, probs = 0.99, na.rm = TRUE)
# 2. Create a "crosspred" object
pred_HI_hospis <- crosspred(cb_HI_hospis,
model_HI_hospis,
cen = pct_1,
by = 0.1,
cumul = TRUE)
# We want the row in pred_HI_hospis that is closest to pct_99:
idx_99 <- which.min(abs(pred_HI_hospis$predvar - pct_99))
RR_99_hospis <- pred_HI_hospis$allRRfit [idx_99]
RR_99_low_hospis <- pred_HI_hospis$allRRlow [idx_99]
RR_99_hi_hospis <- pred_HI_hospis$allRRhigh[idx_99]
cat("RR(99th vs. 1st pct) =", RR_99_hospis, "(", RR_99_low_hospis, "-", RR_99_hi_hospis, ")\n")## RR(99th vs. 1st pct) = 1.172115 ( 0.9781486 - 1.404544 )# RR for +1 HI increase
temp_x <- -1
temp_xp1 <- 0
idx_x <- which.min(abs(pred_HI_hospis$predvar - temp_x))
idx_xp1 <- which.min(abs(pred_HI_hospis$predvar - temp_xp1))
RR_x_hospis <- pred_HI_hospis$allRRfit[idx_x]
RR_xp1_hospis <- pred_HI_hospis$allRRfit[idx_xp1]
RR_1C_hospis <- RR_xp1_hospis / RR_x_hospis
RR_1C_low_hospis <- pred_HI_hospis$allRRlow[idx_xp1] / pred_HI_hospis$allRRhigh[idx_x]
RR_1C_hi_hospis <- pred_HI_hospis$allRRhigh[idx_xp1] / pred_HI_hospis$allRRlow[idx_x]
cat("RR for +1 HI from", temp_x, "to", temp_xp1, ":", RR_1C_hospis, "(", RR_1C_low_hospis, "-", RR_1C_hi_hospis, ")\n")## RR for +1 HI from -1 to 0 : 1.005691 ( 0.9938976 - 1.017624 )The results show that the relative risk (RR) at the 99th percentile (HI = 29°C) compared to the 1st percentile (HI = -1°C) is 1.1721148 (0.9781486 - 1.4045445). This indicates that the relative risk of hospitalizations increases but not significantly with the increase in the heat index.
For an increase of +1 in the heat index from the -1 value, the RR is 1.0056909 (0.9938976 - 1.017624).
Analyse de sensibilité
# Définition des valeurs de df pour var et lag
var_df_values <- c(1, 2, 3)
lag_df_values <- c(2, 3)
lag_max_values <- c(2, 3)
# Création d'une grille de graphiques 3x4
par(mfrow = c(3, 4), mar = c(4, 4, 3, 1))
# Stockage des AIC pour tous les modèles
aic_matrix <- matrix(NA, nrow = length(var_df_values), ncol = length(lag_df_values) * length(lag_max_values))
rownames(aic_matrix) <- paste("var_df =", var_df_values)
colnames(aic_matrix) <- paste("lag_df =", rep(lag_df_values, each = length(lag_max_values)), ", lag_max =", rep(lag_max_values, times = length(lag_df_values)))
# Boucle sur toutes les combinaisons
for (i in seq_along(var_df_values)) {
for (j in seq_along(lag_df_values)) {
for (k in seq_along(lag_max_values)) {
var_df <- var_df_values[i]
lag_df <- lag_df_values[j]
lag_max <- lag_max_values[k]
# Calcul de l'index pour la matrice AIC
col_idx <- (j - 1) * length(lag_max_values) + k
# Création du crossbasis
cb.heat_index <- crossbasis(B$HEAT.INDEX,
lag = lag_max,
argvar = list(fun = "ns", df = var_df),
arglag = list(fun = "ns", df = lag_df))
# Modèle
model_hi <- glm(nombre_hospitalisation ~ cb.heat_index
+ weekday
+ holiday
+ ns(DATE, 8*4)
+ ns(B$Year, 2),
data = B,
family = poisson(link = "log"))
# Stockage de l'AIC
aic_matrix[i, col_idx] <- AIC(model_hi)
# Prédiction et graphique
pred.hi <- crosspred(cb.heat_index, model_hi, at = seq(-5, 35, by = 0.1), cumul = TRUE, cen = -1)
plot(pred.hi$predvar, pred.hi$allRRfit, type = "l", col = "darkred", lwd = 2,
xlab = "Heat Index", ylab = "Cumulative RR",
main = paste("var_df =", var_df, ", lag_df =", lag_df, "\nlag_max =", lag_max),
ylim = c(0.8, 1.5), cex.main = 0.8)
lines(pred.hi$predvar, pred.hi$allRRlow, col = "darkred", lwd = 1, lty = 2)
lines(pred.hi$predvar, pred.hi$allRRhigh, col = "darkred", lwd = 1, lty = 2)
abline(h = 1, col = "black", lwd = 1)
}
}
}
# Retour à la disposition normale
par(mfrow = c(1, 1))
# Affichage du tableau des AIC
cat("Tableau des AIC pour l'analyse de sensibilité (Hospitalisations):\n")## Tableau des AIC pour l'analyse de sensibilité (Hospitalisations):## lag_df = 2 , lag_max = 2 lag_df = 2 , lag_max = 3
## var_df = 1 6876.42 6871.44
## var_df = 2 6879.74 6872.84
## var_df = 3 6883.60 6876.69
## lag_df = 3 , lag_max = 2 lag_df = 3 , lag_max = 3
## var_df = 1 6877.46 6873.44
## var_df = 2 6881.60 6875.58
## var_df = 3 6887.47 6880.96# Identification du meilleur modèle
best_idx <- which(aic_matrix == min(aic_matrix), arr.ind = TRUE)
best_var_df <- var_df_values[best_idx[1]]
best_col <- best_idx[2]
best_lag_df_idx <- ceiling(best_col / length(lag_max_values))
best_lag_max_idx <- best_col %% length(lag_max_values)
if (best_lag_max_idx == 0) best_lag_max_idx <- length(lag_max_values)
best_lag_df <- lag_df_values[best_lag_df_idx]
best_lag_max <- lag_max_values[best_lag_max_idx]
cat("\nMeilleur modèle: var_df =", best_var_df, ", lag_df =", best_lag_df, ", lag_max =", best_lag_max, "\n")##
## Meilleur modèle: var_df = 1 , lag_df = 2 , lag_max = 3## AIC minimum: 6871.44Results : Constrained Hospitalizations
HI and Constrained Hospitalizations
The following plot shows the cumulative relative risk (RR) of constrained hospitalizations as a function of the heat index.
# Création du crossbasis pour contraintes
cb_HI_contrainte <- crossbasis(B$HEAT.INDEX,
lag = 2,
argvar = list(fun = "ns", df = 2),
arglag = list(fun = "ns", df = 2))
# Modèle pour contraintes
model_HI_contrainte <- glm(nombre_contrainte ~ cb_HI_contrainte
+ weekday
+ holiday
+ ns(DATE, 8*4)
+ ns(Year, 2),
data = B,
family = poisson(link = "log"))
pred.hi.contrainte <- crosspred(cb_HI_contrainte, model_HI_contrainte, at= -5:35, cumul=TRUE, cen=-1)
plot(pred.hi.contrainte, "overall",
col="red", cumul=TRUE,
xlab="Heat Index",
ylab="Cumulative RR",
main="Cumulative Relative Risk of Constrained Hospitalizations by Heat Index")
Comment: The plot shows the relationship between heat index and constrained hospitalization risk, which may differ from voluntary hospitalizations.
# 1. Define percentiles of interest
pct_1 <- quantile(B$HEAT.INDEX, probs = 0.01, na.rm = TRUE)
pct_99 <- quantile(B$HEAT.INDEX, probs = 0.99, na.rm = TRUE)
# 2. Create a "crosspred" object
pred_HI_contrainte <- crosspred(cb_HI_contrainte,
model_HI_contrainte,
cen = pct_1,
by = 0.1,
cumul = TRUE)
# We want the row in pred_HI_contrainte that is closest to pct_99:
idx_99 <- which.min(abs(pred_HI_contrainte$predvar - pct_99))
RR_99_contrainte <- pred_HI_contrainte$allRRfit [idx_99]
RR_99_low_contrainte <- pred_HI_contrainte$allRRlow [idx_99]
RR_99_hi_contrainte <- pred_HI_contrainte$allRRhigh[idx_99]
cat("RR(99th vs. 1st pct) =", RR_99_contrainte, "(", RR_99_low_contrainte, "-", RR_99_hi_contrainte, ")\n")## RR(99th vs. 1st pct) = 1.347602 ( 1.036204 - 1.752582 )# RR for +1 HI increase
temp_x <- -1
temp_xp1 <- 0
idx_x <- which.min(abs(pred_HI_contrainte$predvar - temp_x))
idx_xp1 <- which.min(abs(pred_HI_contrainte$predvar - temp_xp1))
RR_x_contrainte <- pred_HI_contrainte$allRRfit[idx_x]
RR_xp1_contrainte <- pred_HI_contrainte$allRRfit[idx_xp1]
RR_1C_contrainte <- RR_xp1_contrainte / RR_x_contrainte
RR_1C_low_contrainte <- pred_HI_contrainte$allRRlow[idx_xp1] / pred_HI_contrainte$allRRhigh[idx_x]
RR_1C_hi_contrainte <- pred_HI_contrainte$allRRhigh[idx_xp1] / pred_HI_contrainte$allRRlow[idx_x]
cat("RR for +1 HI from", temp_x, "to", temp_xp1, ":", RR_1C_contrainte, "(", RR_1C_low_contrainte, "-", RR_1C_hi_contrainte, ")\n")## RR for +1 HI from -1 to 0 : 1.011551 ( 0.9941967 - 1.029207 )The results show that the relative risk (RR) at the 99th percentile (HI = 29°C) compared to the 1st percentile (HI = -1°C) is 1.3476021 (1.0362036 - 1.7525816). This indicates that the relative risk of constrained hospitalizations increases but not significantly with the increase in the heat index.
For an increase of +1 in the heat index from the -1 value, the RR is 1.0115505 (0.9941967 - 1.0292073).
Analyse de sensibilité
# Définition des valeurs de df pour var et lag
var_df_values <- c(1, 2, 3)
lag_df_values <- c(2, 3)
lag_max_values <- c(2, 3)
# Création d'une grille de graphiques 3x4
par(mfrow = c(3, 4), mar = c(4, 4, 3, 1))
# Stockage des AIC pour tous les modèles
aic_matrix <- matrix(NA, nrow = length(var_df_values), ncol = length(lag_df_values) * length(lag_max_values))
rownames(aic_matrix) <- paste("var_df =", var_df_values)
colnames(aic_matrix) <- paste("lag_df =", rep(lag_df_values, each = length(lag_max_values)), ", lag_max =", rep(lag_max_values, times = length(lag_df_values)))
# Boucle sur toutes les combinaisons
for (i in seq_along(var_df_values)) {
for (j in seq_along(lag_df_values)) {
for (k in seq_along(lag_max_values)) {
var_df <- var_df_values[i]
lag_df <- lag_df_values[j]
lag_max <- lag_max_values[k]
# Calcul de l'index pour la matrice AIC
col_idx <- (j - 1) * length(lag_max_values) + k
# Création du crossbasis
cb.heat_index <- crossbasis(B$HEAT.INDEX,
lag = lag_max,
argvar = list(fun = "ns", df = var_df),
arglag = list(fun = "ns", df = lag_df))
# Modèle
model_hi <- glm(nombre_contrainte ~ cb.heat_index
+ weekday
+ holiday
+ ns(DATE, 8*4)
+ ns(B$Year, 2),
data = B,
family = poisson(link = "log"))
# Stockage de l'AIC
aic_matrix[i, col_idx] <- AIC(model_hi)
# Prédiction et graphique
pred.hi <- crosspred(cb.heat_index, model_hi, at = seq(-5, 35, by = 0.1), cumul = TRUE, cen = -1)
plot(pred.hi$predvar, pred.hi$allRRfit, type = "l", col = "darkred", lwd = 2,
xlab = "Heat Index", ylab = "Cumulative RR",
main = paste("var_df =", var_df, ", lag_df =", lag_df, "\nlag_max =", lag_max),
ylim = c(0.8, 1.5), cex.main = 0.8)
lines(pred.hi$predvar, pred.hi$allRRlow, col = "darkred", lwd = 1, lty = 2)
lines(pred.hi$predvar, pred.hi$allRRhigh, col = "darkred", lwd = 1, lty = 2)
abline(h = 1, col = "black", lwd = 1)
}
}
}
# Retour à la disposition normale
par(mfrow = c(1, 1))
# Affichage du tableau des AIC
cat("Tableau des AIC pour l'analyse de sensibilité (Hospitalisations contraintes):\n")## Tableau des AIC pour l'analyse de sensibilité (Hospitalisations contraintes):## lag_df = 2 , lag_max = 2 lag_df = 2 , lag_max = 3
## var_df = 1 5848.60 5844.83
## var_df = 2 5850.68 5845.54
## var_df = 3 5853.83 5848.44
## lag_df = 3 , lag_max = 2 lag_df = 3 , lag_max = 3
## var_df = 1 5846.64 5846.66
## var_df = 2 5849.09 5849.10
## var_df = 3 5854.11 5853.88# Identification du meilleur modèle
best_idx <- which(aic_matrix == min(aic_matrix), arr.ind = TRUE)
best_var_df <- var_df_values[best_idx[1]]
best_col <- best_idx[2]
best_lag_df_idx <- ceiling(best_col / length(lag_max_values))
best_lag_max_idx <- best_col %% length(lag_max_values)
if (best_lag_max_idx == 0) best_lag_max_idx <- length(lag_max_values)
best_lag_df <- lag_df_values[best_lag_df_idx]
best_lag_max <- lag_max_values[best_lag_max_idx]
cat("\nMeilleur modèle: var_df =", best_var_df, ", lag_df =", best_lag_df, ", lag_max =", best_lag_max, "\n")##
## Meilleur modèle: var_df = 1 , lag_df = 2 , lag_max = 3## AIC minimum: 5844.83Results : F1
HI and F1
The following plot shows the cumulative relative risk (RR) of F1 diagnoses (substance use disorders) as a function of the heat index.
# Création du crossbasis pour F1
cb_HI_F1 <- crossbasis(B$HEAT.INDEX,
lag = 2,
argvar = list(fun = "ns", df = 2),
arglag = list(fun = "ns", df = 2))
# Modèle pour F1
model_HI_F1 <- glm(nombre_F1 ~ cb_HI_F1
+ weekday
+ holiday
+ ns(DATE, 8*4)
+ ns(Year, 2),
data = B,
family = poisson(link = "log"))
pred.hi.F1 <- crosspred(cb_HI_F1, model_HI_F1, at= -5:35, cumul=TRUE, cen=-1)
plot(pred.hi.F1, "overall",
col="purple", cumul=TRUE,
xlab="Heat Index",
ylab="Cumulative RR",
main="Cumulative Relative Risk of F1 Diagnoses by Heat Index")
Comment: The plot shows the relationship between heat index and substance use disorder admissions.
# 1. Define percentiles of interest
pct_1 <- quantile(B$HEAT.INDEX, probs = 0.01, na.rm = TRUE)
pct_99 <- quantile(B$HEAT.INDEX, probs = 0.99, na.rm = TRUE)
# 2. Create a "crosspred" object
pred_HI_F1 <- crosspred(cb_HI_F1,
model_HI_F1,
cen = pct_1,
by = 0.1,
cumul = TRUE)
# We want the row in pred_HI_F1 that is closest to pct_99:
idx_99 <- which.min(abs(pred_HI_F1$predvar - pct_99))
RR_99_F1 <- pred_HI_F1$allRRfit [idx_99]
RR_99_low_F1 <- pred_HI_F1$allRRlow [idx_99]
RR_99_hi_F1 <- pred_HI_F1$allRRhigh[idx_99]
cat("RR(99th vs. 1st pct) =", RR_99_F1, "(", RR_99_low_F1, "-", RR_99_hi_F1, ")\n")## RR(99th vs. 1st pct) = 1.389044 ( 0.9192351 - 2.098966 )# RR for +1 HI increase
temp_x <- -1
temp_xp1 <- 0
idx_x <- which.min(abs(pred_HI_F1$predvar - temp_x))
idx_xp1 <- which.min(abs(pred_HI_F1$predvar - temp_xp1))
RR_x_F1 <- pred_HI_F1$allRRfit[idx_x]
RR_xp1_F1 <- pred_HI_F1$allRRfit[idx_xp1]
RR_1C_F1 <- RR_xp1_F1 / RR_x_F1
RR_1C_low_F1 <- pred_HI_F1$allRRlow[idx_xp1] / pred_HI_F1$allRRhigh[idx_x]
RR_1C_hi_F1 <- pred_HI_F1$allRRhigh[idx_xp1] / pred_HI_F1$allRRlow[idx_x]
cat("RR(99th vs. 1st pct) =", RR_99_F1, "(", RR_99_low_F1, "-", RR_99_hi_F1, ")\n")## RR(99th vs. 1st pct) = 1.389044 ( 0.9192351 - 2.098966 )cat("RR for +1 HI from", temp_x, "to", temp_xp1, ":", RR_1C_F1, "(", RR_1C_low_F1, "-", RR_1C_hi_F1, ")\n")## RR for +1 HI from -1 to 0 : 1.025229 ( 0.9980065 - 1.053194 )The results show that the relative risk (RR) at the 99th percentile (HI = 29°C) compared to the 1st percentile (HI = -1°C) is 1.3890439 (0.9192351 - 2.0989659). This indicates that the relative risk of F1 diagnoses increases but not significantly with the increase in the heat index.
For an increase of +1 in the heat index from the -1 value, the RR is 1.0252291 (0.9980065 - 1.0531943).
Analyse de sensibilité
# Définition des valeurs de df pour var et lag
var_df_values <- c(1, 2, 3)
lag_df_values <- c(2, 3)
lag_max_values <- c(2, 3)
# Création d'une grille de graphiques 3x4
par(mfrow = c(3, 4), mar = c(4, 4, 3, 1))
# Stockage des AIC pour tous les modèles
aic_matrix <- matrix(NA, nrow = length(var_df_values), ncol = length(lag_df_values) * length(lag_max_values))
rownames(aic_matrix) <- paste("var_df =", var_df_values)
colnames(aic_matrix) <- paste("lag_df =", rep(lag_df_values, each = length(lag_max_values)), ", lag_max =", rep(lag_max_values, times = length(lag_df_values)))
# Boucle sur toutes les combinaisons
for (i in seq_along(var_df_values)) {
for (j in seq_along(lag_df_values)) {
for (k in seq_along(lag_max_values)) {
var_df <- var_df_values[i]
lag_df <- lag_df_values[j]
lag_max <- lag_max_values[k]
# Calcul de l'index pour la matrice AIC
col_idx <- (j - 1) * length(lag_max_values) + k
# Création du crossbasis
cb.heat_index <- crossbasis(B$HEAT.INDEX,
lag = lag_max,
argvar = list(fun = "ns", df = var_df),
arglag = list(fun = "ns", df = lag_df))
# Modèle
model_hi <- glm(nombre_F1 ~ cb.heat_index
+ weekday
+ holiday
+ ns(DATE, 8*4)
+ ns(B$Year, 2),
data = B,
family = poisson(link = "log"))
# Stockage de l'AIC
aic_matrix[i, col_idx] <- AIC(model_hi)
# Prédiction et graphique
pred.hi <- crosspred(cb.heat_index, model_hi, at = seq(-5, 35, by = 0.1), cumul = TRUE, cen = -1)
plot(pred.hi$predvar, pred.hi$allRRfit, type = "l", col = "purple", lwd = 2,
xlab = "Heat Index", ylab = "Cumulative RR",
main = paste("var_df =", var_df, ", lag_df =", lag_df, "\nlag_max =", lag_max),
ylim = c(0.8, 1.5), cex.main = 0.8)
lines(pred.hi$predvar, pred.hi$allRRlow, col = "purple", lwd = 1, lty = 2)
lines(pred.hi$predvar, pred.hi$allRRhigh, col = "purple", lwd = 1, lty = 2)
abline(h = 1, col = "black", lwd = 1)
}
}
}
# Retour à la disposition normale
par(mfrow = c(1, 1))
# Affichage du tableau des AIC
cat("Tableau des AIC pour l'analyse de sensibilité (F1 - Troubles liés aux substances):\n")## Tableau des AIC pour l'analyse de sensibilité (F1 - Troubles liés aux substances):## lag_df = 2 , lag_max = 2 lag_df = 2 , lag_max = 3
## var_df = 1 4579.35 4575.53
## var_df = 2 4580.49 4577.21
## var_df = 3 4583.83 4580.82
## lag_df = 3 , lag_max = 2 lag_df = 3 , lag_max = 3
## var_df = 1 4579.85 4574.40
## var_df = 2 4577.15 4575.73
## var_df = 3 4579.33 4581.26# Identification du meilleur modèle
best_idx <- which(aic_matrix == min(aic_matrix), arr.ind = TRUE)
best_var_df <- var_df_values[best_idx[1]]
best_col <- best_idx[2]
best_lag_df_idx <- ceiling(best_col / length(lag_max_values))
best_lag_max_idx <- best_col %% length(lag_max_values)
if (best_lag_max_idx == 0) best_lag_max_idx <- length(lag_max_values)
best_lag_df <- lag_df_values[best_lag_df_idx]
best_lag_max <- lag_max_values[best_lag_max_idx]
cat("\nMeilleur modèle: var_df =", best_var_df, ", lag_df =", best_lag_df, ", lag_max =", best_lag_max, "\n")##
## Meilleur modèle: var_df = 1 , lag_df = 3 , lag_max = 3## AIC minimum: 4574.398Results : F2
HI and F2
The following plot shows the cumulative relative risk (RR) of F2 diagnoses (schizophrenia and related disorders) as a function of the heat index.
# Création du crossbasis pour F2
cb_HI_F2 <- crossbasis(B$HEAT.INDEX,
lag = 2,
argvar = list(fun = "ns", df = 2),
arglag = list(fun = "ns", df = 2))
# Modèle pour F2
model_HI_F2 <- glm(nombre_F2 ~ cb_HI_F2
+ weekday
+ holiday
+ ns(DATE, 8*4)
+ ns(Year, 2),
data = B,
family = poisson(link = "log"))
pred.hi.F2 <- crosspred(cb_HI_F2, model_HI_F2, at= -5:35, cumul=TRUE, cen=-1)
plot(pred.hi.F2, "overall",
col="orange", cumul=TRUE,
xlab="Heat Index",
ylab="Cumulative RR",
main="Cumulative Relative Risk of F2 Diagnoses by Heat Index")
Comment: The plot shows the relationship between heat index and schizophrenia-related disorder admissions.
# 1. Define percentiles of interest
pct_1 <- quantile(B$HEAT.INDEX, probs = 0.01, na.rm = TRUE)
pct_99 <- quantile(B$HEAT.INDEX, probs = 0.99, na.rm = TRUE)
# 2. Create a "crosspred" object
pred_HI_F2 <- crosspred(cb_HI_F2,
model_HI_F2,
cen = pct_1,
by = 0.1,
cumul = TRUE)
# We want the row in pred_HI_F2 that is closest to pct_99:
idx_99 <- which.min(abs(pred_HI_F2$predvar - pct_99))
RR_99_F2 <- pred_HI_F2$allRRfit [idx_99]
RR_99_low_F2 <- pred_HI_F2$allRRlow [idx_99]
RR_99_hi_F2 <- pred_HI_F2$allRRhigh[idx_99]
cat("RR(99th vs. 1st pct) =", RR_99_F2, "(", RR_99_low_F2, "-", RR_99_hi_F2, ")\n")## RR(99th vs. 1st pct) = 1.210726 ( 0.9421276 - 1.555901 )# RR for +1 HI increase
temp_x <- -1
temp_xp1 <- 0
idx_x <- which.min(abs(pred_HI_F2$predvar - temp_x))
idx_xp1 <- which.min(abs(pred_HI_F2$predvar - temp_xp1))
RR_x_F2 <- pred_HI_F2$allRRfit[idx_x]
RR_xp1_F2 <- pred_HI_F2$allRRfit[idx_xp1]
RR_1C_F2 <- RR_xp1_F2 / RR_x_F2
RR_1C_low_F2 <- pred_HI_F2$allRRlow[idx_xp1] / pred_HI_F2$allRRhigh[idx_x]
RR_1C_hi_F2 <- pred_HI_F2$allRRhigh[idx_xp1] / pred_HI_F2$allRRlow[idx_x]
cat("RR for +1 HI from", temp_x, "to", temp_xp1, ":", RR_1C_F2, "(", RR_1C_low_F2, "-", RR_1C_hi_F2, ")\n")## RR for +1 HI from -1 to 0 : 1.013901 ( 0.9971368 - 1.030947 )The results show that the relative risk (RR) at the 99th percentile (HI = 29°C) compared to the 1st percentile (HI = -1°C) is 1.2107261 (0.9421276 - 1.5559014). This indicates that the relative risk of F2 diagnoses increases but not significantly with the increase in the heat index.
For an increase of +1 in the heat index from the -1 value, the RR is 1.013901 (0.9971368 - 1.030947).
Analyse de sensibilité
# Définition des valeurs de df pour var et lag
var_df_values <- c(1, 2, 3)
lag_df_values <- c(2, 3)
lag_max_values <- c(2, 3)
# Création d'une grille de graphiques 3x4
par(mfrow = c(3, 4), mar = c(4, 4, 3, 1))
# Stockage des AIC pour tous les modèles
aic_matrix <- matrix(NA, nrow = length(var_df_values), ncol = length(lag_df_values) * length(lag_max_values))
rownames(aic_matrix) <- paste("var_df =", var_df_values)
colnames(aic_matrix) <- paste("lag_df =", rep(lag_df_values, each = length(lag_max_values)), ", lag_max =", rep(lag_max_values, times = length(lag_df_values)))
# Boucle sur toutes les combinaisons
for (i in seq_along(var_df_values)) {
for (j in seq_along(lag_df_values)) {
for (k in seq_along(lag_max_values)) {
var_df <- var_df_values[i]
lag_df <- lag_df_values[j]
lag_max <- lag_max_values[k]
# Calcul de l'index pour la matrice AIC
col_idx <- (j - 1) * length(lag_max_values) + k
# Création du crossbasis
cb.heat_index <- crossbasis(B$HEAT.INDEX,
lag = lag_max,
argvar = list(fun = "ns", df = var_df),
arglag = list(fun = "ns", df = lag_df))
# Modèle
model_hi <- glm(nombre_F2 ~ cb.heat_index
+ weekday
+ holiday
+ ns(DATE, 8*4)
+ ns(B$Year, 2),
data = B,
family = poisson(link = "log"))
# Stockage de l'AIC
aic_matrix[i, col_idx] <- AIC(model_hi)
# Prédiction et graphique
pred.hi <- crosspred(cb.heat_index, model_hi, at = seq(-5, 35, by = 0.1), cumul = TRUE, cen = -1)
plot(pred.hi$predvar, pred.hi$allRRfit, type = "l", col = "orange", lwd = 2,
xlab = "Heat Index", ylab = "Cumulative RR",
main = paste("var_df =", var_df, ", lag_df =", lag_df, "\nlag_max =", lag_max),
ylim = c(0.8, 1.5), cex.main = 0.8)
lines(pred.hi$predvar, pred.hi$allRRlow, col = "orange", lwd = 1, lty = 2)
lines(pred.hi$predvar, pred.hi$allRRhigh, col = "orange", lwd = 1, lty = 2)
abline(h = 1, col = "black", lwd = 1)
}
}
}
# Retour à la disposition normale
par(mfrow = c(1, 1))
# Affichage du tableau des AIC
cat("Tableau des AIC pour l'analyse de sensibilité (F2 - Schizophrénie et troubles délirants):\n")## Tableau des AIC pour l'analyse de sensibilité (F2 - Schizophrénie et troubles délirants):## lag_df = 2 , lag_max = 2 lag_df = 2 , lag_max = 3
## var_df = 1 6029.91 6025.19
## var_df = 2 6032.19 6027.79
## var_df = 3 6034.55 6031.09
## lag_df = 3 , lag_max = 2 lag_df = 3 , lag_max = 3
## var_df = 1 6031.87 6027.14
## var_df = 2 6035.31 6030.90
## var_df = 3 6039.65 6035.53# Identification du meilleur modèle
best_idx <- which(aic_matrix == min(aic_matrix), arr.ind = TRUE)
best_var_df <- var_df_values[best_idx[1]]
best_col <- best_idx[2]
best_lag_df_idx <- ceiling(best_col / length(lag_max_values))
best_lag_max_idx <- best_col %% length(lag_max_values)
if (best_lag_max_idx == 0) best_lag_max_idx <- length(lag_max_values)
best_lag_df <- lag_df_values[best_lag_df_idx]
best_lag_max <- lag_max_values[best_lag_max_idx]
cat("\nMeilleur modèle: var_df =", best_var_df, ", lag_df =", best_lag_df, ", lag_max =", best_lag_max, "\n")##
## Meilleur modèle: var_df = 1 , lag_df = 2 , lag_max = 3cat("\nMeilleur modèle: var_df =", var_df_values[best_idx[1]], ", lag_df =", lag_df_values[best_idx[2]], "\n")##
## Meilleur modèle: var_df = 1 , lag_df = 3## AIC minimum: 6025.191Results : F3
HI and F3
The following plot shows the cumulative relative risk (RR) of F3 diagnoses (mood disorders) as a function of the heat index.
# Création du crossbasis pour F3
cb_HI_F3 <- crossbasis(B$HEAT.INDEX,
lag = 2,
argvar = list(fun = "ns", df = 2),
arglag = list(fun = "ns", df = 2))
# Modèle pour F3
model_HI_F3 <- glm(nombre_F3 ~ cb_HI_F3
+ weekday
+ holiday
+ ns(DATE, 8*4)
+ ns(Year, 2),
data = B,
family = poisson(link = "log"))
pred.hi.F3 <- crosspred(cb_HI_F3, model_HI_F3, at= -5:35, cumul=TRUE, cen=-1)
plot(pred.hi.F3, "overall",
col="brown", cumul=TRUE,
xlab="Heat Index",
ylab="Cumulative RR",
main="Cumulative Relative Risk of F3 Diagnoses by Heat Index")
Comment: The plot shows the relationship between heat index and mood disorder admissions.
# 1. Define percentiles of interest
pct_1 <- quantile(B$HEAT.INDEX, probs = 0.01, na.rm = TRUE)
pct_99 <- quantile(B$HEAT.INDEX, probs = 0.99, na.rm = TRUE)
# 2. Create a "crosspred" object
pred_HI_F3 <- crosspred(cb_HI_F3,
model_HI_F3,
cen = pct_1,
by = 0.1,
cumul = TRUE)
# We want the row in pred_HI_F3 that is closest to pct_99:
idx_99 <- which.min(abs(pred_HI_F3$predvar - pct_99))
RR_99_F3 <- pred_HI_F3$allRRfit [idx_99]
RR_99_low_F3 <- pred_HI_F3$allRRlow [idx_99]
RR_99_hi_F3 <- pred_HI_F3$allRRhigh[idx_99]
cat("RR(99th vs. 1st pct) =", RR_99_F3, "(", RR_99_low_F3, "-", RR_99_hi_F3, ")\n")## RR(99th vs. 1st pct) = 1.072769 ( 0.8681733 - 1.325579 )# RR for +1 HI increase
temp_x <- -1
temp_xp1 <- 0
idx_x <- which.min(abs(pred_HI_F3$predvar - temp_x))
idx_xp1 <- which.min(abs(pred_HI_F3$predvar - temp_xp1))
RR_x_F3 <- pred_HI_F3$allRRfit[idx_x]
RR_xp1_F3 <- pred_HI_F3$allRRfit[idx_xp1]
RR_1C_F3 <- RR_xp1_F3 / RR_x_F3
RR_1C_low_F3 <- pred_HI_F3$allRRlow[idx_xp1] / pred_HI_F3$allRRhigh[idx_x]
RR_1C_hi_F3 <- pred_HI_F3$allRRhigh[idx_xp1] / pred_HI_F3$allRRlow[idx_x]
cat("RR for +1 HI from", temp_x, "to", temp_xp1, ":", RR_1C_F3, "(", RR_1C_low_F3, "-", RR_1C_hi_F3, ")\n")## RR for +1 HI from -1 to 0 : 1.009788 ( 0.9958611 - 1.023909 )The results show that the relative risk (RR) at the 99th percentile (HI = 29°C) compared to the 1st percentile (HI = -1°C) is 1.0727686 (0.8681733 - 1.3255792). This indicates that the relative risk of F3 diagnoses increasesbut not significantly with the increase in the heat index.
For an increase of +1 in the heat index from the -1 value, the RR is 1.0097876 (0.9958611 - 1.0239089).
Analyse de sensibilité
# Définition des valeurs de df pour var et lag
var_df_values <- c(1, 2, 3)
lag_df_values <- c(2, 3)
lag_max_values <- c(2, 3)
# Création d'une grille de graphiques 3x4
par(mfrow = c(3, 4), mar = c(4, 4, 3, 1))
# Stockage des AIC pour tous les modèles
aic_matrix <- matrix(NA, nrow = length(var_df_values), ncol = length(lag_df_values) * length(lag_max_values))
rownames(aic_matrix) <- paste("var_df =", var_df_values)
colnames(aic_matrix) <- paste("lag_df =", rep(lag_df_values, each = length(lag_max_values)), ", lag_max =", rep(lag_max_values, times = length(lag_df_values)))
# Boucle sur toutes les combinaisons
for (i in seq_along(var_df_values)) {
for (j in seq_along(lag_df_values)) {
for (k in seq_along(lag_max_values)) {
var_df <- var_df_values[i]
lag_df <- lag_df_values[j]
lag_max <- lag_max_values[k]
# Calcul de l'index pour la matrice AIC
col_idx <- (j - 1) * length(lag_max_values) + k
# Création du crossbasis
cb.heat_index <- crossbasis(B$HEAT.INDEX,
lag = lag_max,
argvar = list(fun = "ns", df = var_df),
arglag = list(fun = "ns", df = lag_df))
# Modèle
model_hi <- glm(nombre_F3 ~ cb.heat_index
+ weekday
+ holiday
+ ns(DATE, 8*4)
+ ns(B$Year, 2),
data = B,
family = poisson(link = "log"))
# Stockage de l'AIC
aic_matrix[i, col_idx] <- AIC(model_hi)
# Prédiction et graphique
pred.hi <- crosspred(cb.heat_index, model_hi, at = seq(-5, 35, by = 0.1), cumul = TRUE, cen = -1)
plot(pred.hi$predvar, pred.hi$allRRfit, type = "l", col = "brown", lwd = 2,
xlab = "Heat Index", ylab = "Cumulative RR",
main = paste("var_df =", var_df, ", lag_df =", lag_df, "\nlag_max =", lag_max),
ylim = c(0.8, 1.5), cex.main = 0.8)
lines(pred.hi$predvar, pred.hi$allRRlow, col = "brown", lwd = 1, lty = 2)
lines(pred.hi$predvar, pred.hi$allRRhigh, col = "brown", lwd = 1, lty = 2)
abline(h = 1, col = "black", lwd = 1)
}
}
}
# Retour à la disposition normale
par(mfrow = c(1, 1))
# Affichage du tableau des AIC
cat("Tableau des AIC pour l'analyse de sensibilité (F3 - Troubles de l'humeur):\n")## Tableau des AIC pour l'analyse de sensibilité (F3 - Troubles de l'humeur):## lag_df = 2 , lag_max = 2 lag_df = 2 , lag_max = 3
## var_df = 1 6499.76 6496.03
## var_df = 2 6502.32 6496.96
## var_df = 3 6506.11 6500.87
## lag_df = 3 , lag_max = 2 lag_df = 3 , lag_max = 3
## var_df = 1 6500.30 6497.94
## var_df = 2 6500.14 6497.38
## var_df = 3 6504.02 6501.67# Identification du meilleur modèle
best_idx <- which(aic_matrix == min(aic_matrix), arr.ind = TRUE)
cat("\nMeilleur modèle: var_df =", var_df_values[best_idx[1]], ", lag_df =", lag_df_values[best_idx[2]], "\n")##
## Meilleur modèle: var_df = 1 , lag_df = 3## AIC minimum: 6496.03Results : F4
HI and F4
The following plot shows the cumulative relative risk (RR) of F4 diagnoses (neurotic, stress-related and somatoform disorders) as a function of the heat index.
# Création du crossbasis pour F4
cb_HI_F4 <- crossbasis(B$HEAT.INDEX,
lag = 2,
argvar = list(fun = "ns", df = 2),
arglag = list(fun = "ns", df = 2))
# Modèle pour F4
model_HI_F4 <- glm(nombre_F4 ~ cb_HI_F4
+ weekday
+ holiday
+ ns(DATE, 8*4)
+ ns(Year, 2),
data = B,
family = poisson(link = "log"))
pred.hi.F4 <- crosspred(cb_HI_F4, model_HI_F4, at= -5:35, cumul=TRUE, cen=-1)
plot(pred.hi.F4, "overall",
col="pink", cumul=TRUE,
xlab="Heat Index",
ylab="Cumulative RR",
main="Cumulative Relative Risk of F4 Diagnoses by Heat Index")
Comment: The plot shows the relationship between heat index and neurotic disorder admissions.
# 1. Define percentiles of interest
pct_1 <- quantile(B$HEAT.INDEX, probs = 0.01, na.rm = TRUE)
pct_99 <- quantile(B$HEAT.INDEX, probs = 0.99, na.rm = TRUE)
# 2. Create a "crosspred" object
pred_HI_F4 <- crosspred(cb_HI_F4,
model_HI_F4,
cen = pct_1,
by = 0.1,
cumul = TRUE)
# We want the row in pred_HI_F4 that is closest to pct_99:
idx_99 <- which.min(abs(pred_HI_F4$predvar - pct_99))
RR_99_F4 <- pred_HI_F4$allRRfit [idx_99]
RR_99_low_F4 <- pred_HI_F4$allRRlow [idx_99]
RR_99_hi_F4 <- pred_HI_F4$allRRhigh[idx_99]
cat("RR(99th vs. 1st pct) =", RR_99_F4, "(", RR_99_low_F4, "-", RR_99_hi_F4, ")\n")## RR(99th vs. 1st pct) = 1.007203 ( 0.7560398 - 1.341806 )# RR for +1 HI increase
temp_x <- -1
temp_xp1 <- 0
idx_x <- which.min(abs(pred_HI_F4$predvar - temp_x))
idx_xp1 <- which.min(abs(pred_HI_F4$predvar - temp_xp1))
RR_x_F4 <- pred_HI_F4$allRRfit[idx_x]
RR_xp1_F4 <- pred_HI_F4$allRRfit[idx_xp1]
RR_1C_F4 <- RR_xp1_F4 / RR_x_F4
RR_1C_low_F4 <- pred_HI_F4$allRRlow[idx_xp1] / pred_HI_F4$allRRhigh[idx_x]
RR_1C_hi_F4 <- pred_HI_F4$allRRhigh[idx_xp1] / pred_HI_F4$allRRlow[idx_x]
cat("RR for +1 HI from", temp_x, "to", temp_xp1, ":", RR_1C_F4, "(", RR_1C_low_F4, "-", RR_1C_hi_F4, ")\n")## RR for +1 HI from -1 to 0 : 1.003171 ( 0.9848069 - 1.021878 )The results show that the relative risk (RR) at the 99th percentile (HI = 29°C) compared to the 1st percentile (HI = -1°C) is 1.0072035 (0.7560398 - 1.3418061). This indicates that the relative risk of F4 diagnoses increases but not significantly with the increase in the heat index.
For an increase of +1 in the heat index from the -1 value, the RR is 1.0031711 (0.9848069 - 1.0218778).
Analyse de sensibilité
# Définition des valeurs de df pour var et lag
var_df_values <- c(1, 2, 3)
lag_df_values <- c(2, 3)
lag_max_values <- c(2, 3)
# Création d'une grille de graphiques 3x4
par(mfrow = c(3, 4), mar = c(4, 4, 3, 1))
# Stockage des AIC pour tous les modèles
aic_matrix <- matrix(NA, nrow = length(var_df_values), ncol = length(lag_df_values) * length(lag_max_values))
rownames(aic_matrix) <- paste("var_df =", var_df_values)
colnames(aic_matrix) <- paste("lag_df =", rep(lag_df_values, each = length(lag_max_values)), ", lag_max =", rep(lag_max_values, times = length(lag_df_values)))
# Boucle sur toutes les combinaisons
for (i in seq_along(var_df_values)) {
for (j in seq_along(lag_df_values)) {
for (k in seq_along(lag_max_values)) {
var_df <- var_df_values[i]
lag_df <- lag_df_values[j]
lag_max <- lag_max_values[k]
# Calcul de l'index pour la matrice AIC
col_idx <- (j - 1) * length(lag_max_values) + k
# Création du crossbasis
cb.heat_index <- crossbasis(B$HEAT.INDEX,
lag = lag_max,
argvar = list(fun = "ns", df = var_df),
arglag = list(fun = "ns", df = lag_df))
# Modèle
model_hi <- glm(nombre_F4 ~ cb.heat_index
+ weekday
+ holiday
+ ns(DATE, 8*4)
+ ns(B$Year, 2),
data = B,
family = poisson(link = "log"))
# Stockage de l'AIC
aic_matrix[i, col_idx] <- AIC(model_hi)
# Prédiction et graphique
pred.hi <- crosspred(cb.heat_index, model_hi, at = seq(-5, 35, by = 0.1), cumul = TRUE, cen = -1)
plot(pred.hi$predvar, pred.hi$allRRfit, type = "l", col = "pink", lwd = 2,
xlab = "Heat Index", ylab = "Cumulative RR",
main = paste("var_df =", var_df, ", lag_df =", lag_df, "\nlag_max =", lag_max),
ylim = c(0.8, 1.5), cex.main = 0.8)
lines(pred.hi$predvar, pred.hi$allRRlow, col = "pink", lwd = 1, lty = 2)
lines(pred.hi$predvar, pred.hi$allRRhigh, col = "pink", lwd = 1, lty = 2)
abline(h = 1, col = "black", lwd = 1)
}
}
}
# Retour à la disposition normale
par(mfrow = c(1, 1))
# Affichage du tableau des AIC
cat("Tableau des AIC pour l'analyse de sensibilité (F4 - Troubles névrotiques et liés au stress):\n")## Tableau des AIC pour l'analyse de sensibilité (F4 - Troubles névrotiques et liés au stress):## lag_df = 2 , lag_max = 2 lag_df = 2 , lag_max = 3
## var_df = 1 5684.66 5676.97
## var_df = 2 5686.42 5678.96
## var_df = 3 5688.24 5680.95
## lag_df = 3 , lag_max = 2 lag_df = 3 , lag_max = 3
## var_df = 1 5686.15 5678.32
## var_df = 2 5689.58 5682.25
## var_df = 3 5693.36 5686.17# Identification du meilleur modèle
best_idx <- which(aic_matrix == min(aic_matrix), arr.ind = TRUE)
cat("\nMeilleur modèle: var_df =", var_df_values[best_idx[1]], ", lag_df =", lag_df_values[best_idx[2]], "\n")##
## Meilleur modèle: var_df = 1 , lag_df = 3## AIC minimum: 5676.967Guide méthodologique : Étude de sensibilité et choix des paramètres
Introduction
Jai fait faire à ChatGPT un guide méthodologique pour l’étude de sensibilité pour les modèles DLNM. J’ai vérifié et corrigé les erreurs. C’est ce qu’il y a en dessous. A mon sens, la litérature est assez claire sur le lag_max de 2-3 jours. Et nos résultats ne changent pas énormément entre 2 et 3. Dans un soucis de garder les modèles les plus simples, j’ai choisi de garder le lag_max de 2 et var_df de 2 aussi, j’ai aussi garder des relations quadratiques (var_df=2 ) pour être moins rigide que en linéaire. Et je ne vois pas de raison de faire des relations cubiques (var_df=3) ou plus complexes (var_df=4). Ce serait faire l’hypothèse que les effets du HEat.index sont positifs puis negatifs puis positifs… ca me parait improbable.
Paramètres à tester dans l’étude de sensibilité
1. Degrés de liberté des fonctions splines
Paramètre var_df (exposure-response) :
- df = 1 : Relation linéaire entre l’exposition et la réponse
- df = 2 : Relation quadratique, permet une courbe simple
- df = 3 : Relation cubique, permet des courbes plus complexes
- df ≥ 4 : Relations plus complexes, risque de surapprentissage
Paramètre lag_df (lag-response) :
- df = 1 : Effet de décalage constant
- df = 2 : Effet de décalage variable de manière simple
- df = 3 : Effet de décalage variable de manière complexe
2. Maximum lag
Le maximum lag définit la période pendant laquelle l’exposition peut avoir un effet. Pour les admissions psychiatriques, la littérature suggère :
- Lag court (1-3 jours) : Pour les effets immédiats
- Lag moyen (7-14 jours) : Pour les effets à court terme
- Lag long (21-30 jours) : Pour les effets à plus long terme
3. Fonctions de base
- Natural splines (ns) : Recommandées pour leur stabilité
- B-splines (bs) : Alternative flexible
- Polynomial splines : Moins utilisées en pratique
Indicateurs à évaluer
1. Critères d’information
AIC (Akaike Information Criterion) :
- Mesure l’équilibre entre la qualité d’ajustement et la complexité du modèle
- Plus l’AIC est faible, meilleur est le modèle
- Différence d’AIC > 2 considérée comme significative
BIC (Bayesian Information Criterion) :
- Pénalise plus fortement la complexité que l’AIC
- Utile pour éviter le surapprentissage
2. Qualité de l’ajustement
Deviance :
- Mesure la qualité de l’ajustement du modèle
- Plus la deviance est faible, meilleur est l’ajustement
Résidus :
- Analyse des résidus pour détecter les problèmes de spécification
- Tests de normalité et d’indépendance
3. Stabilité des résultats
Forme de la courbe exposition-réponse :
- Vérifier que la forme générale reste cohérente
- Identifier les changements majeurs dans la direction de l’effet
Intervalles de confiance :
- Évaluer la précision des estimations
- Vérifier la stabilité des intervalles de confiance
Interprétation des résultats
1. Sélection du modèle optimal
Approche par critères d’information :
- Choisir le modèle avec l’AIC le plus faible
- Considérer les modèles avec ΔAIC < 2 comme équivalents
Approche par stabilité :
- Privilégier les modèles produisant des résultats stables
- Éviter les modèles avec des courbes instables ou des intervalles de confiance très larges
2. Évaluation de la robustesse
Robustesse forte :
- Résultats similaires pour différentes spécifications
- Intervalles de confiance stables
- Forme de courbe cohérente
Robustesse faible :
- Résultats très variables selon les spécifications
- Intervalles de confiance instables
- Formes de courbes incohérentes
3. Signification pratique
Effet significatif et robuste :
- Association statistiquement significative
- Résultats stables dans les analyses de sensibilité
- Interprétation clinique claire
Effet non robuste :
- Association variable selon les spécifications
- Nécessite une interprétation prudente
- Peut indiquer un manque de puissance ou une spécification inappropriée
Recommandations pratiques
1. Stratégie d’analyse
- Commencer par des modèles simples (df = 1-2)
- Augmenter progressivement la complexité (df = 3-4)
- Évaluer la stabilité des résultats
- Choisir le modèle optimal basé sur AIC et stabilité
2. Présentation des résultats
- Tableau comparatif des AIC pour toutes les spécifications
- Graphiques de sensibilité montrant les courbes pour différentes spécifications
- Discussion de la robustesse des résultats principaux
3. Limitations
- Surapprentissage : Risque avec trop de degrés de liberté
- Puissance statistique : Diminue avec l’augmentation de la complexité
- Interprétabilité : Modèles complexes moins faciles à interpréter
Références méthodologiques
Les recommandations présentées ici sont basées sur :
- Gasparrini et al. (2010) : “Distributed lag linear and non-linear models in R”
- Armstrong (2006) : “Models for the relationship between ambient temperature and daily mortality”
- Gasparrini (2014) : “Modeling exposure-lag-response associations with distributed lag non-linear models”