Ce document reproduit uniquement les analyses et les figures présentées dans le manuscrit soumis. Chaque section correspond à une figure ou à un tableau de l’article. Les libellés, titres et bornes des graphiques reprennent ceux de l’article.

1 Méthode (rappel)

Un modèle non-linéaire à décalage distribué (DLNM) est utilisé pour estimer l’effet du Heat Index (HI) sur le nombre de passages aux urgences psychiatriques.

  • Cross-basis sur le HI : décalage (lag) maximal de 2 jours.
  • Splines naturelles (ns) à 2 degrés de liberté pour la relation exposition-réponse et pour la relation lag-réponse.
  • Régression de Poisson ajustée sur le temps (ns(DATE, 8 df/an)), la tendance annuelle (ns(Year, 2)), le jour de la semaine et les jours fériés.
  • Valeur de référence : HI = -1 °C (1er percentile de la distribution).
  • Résultats exprimés en risques relatifs (RR) avec intervalles de confiance à 95 %, au 99e percentile (HI = 28.4 °C) vs 1er percentile (HI = -1 °C), et par +1 °C au-dessus de -1 °C.
# Ajuste un DLNM pour une variable de comptage et renvoie la crossbasis, le modèle
# et la prédiction cumulée centrée sur -1 °C.
# NB : la crossbasis locale s'appelle littéralement `cb` afin que crosspred()
# retrouve correctement les coefficients du modèle.
fit_hi <- function(outcome, data = B, at = seq(-5, 35, by = 0.1), cen = -1) {
  cb <- crossbasis(data$HEAT.INDEX,
                   lag = 2,
                   argvar = list(fun = "ns", df = 2),
                   arglag = list(fun = "ns", df = 2))
  f <- reformulate(c("cb", "weekday", "holiday", "ns(DATE, 8*4)", "ns(Year, 2)"),
                   response = outcome)
  model <- glm(f, data = data, family = poisson(link = "log"))
  pred  <- crosspred(cb, model, at = at, cen = cen, cumul = TRUE)
  list(cb = cb, model = model, pred = pred, aic = AIC(model))
}

# RR (et IC 95 %) à une valeur donnée de HI, à partir d'une prédiction centrée sur -1.
rr_at <- function(pred, x) {
  i <- which.min(abs(pred$predvar - x))
  # unname() est indispensable : pred$allRRfit[i] est un élément nommé
  # (nom = valeur de predvar). Sans unname(), c(fit = ...) produirait des noms
  # composés ("fit.29"), et fmt_rr() ne retrouverait plus v["fit"] -> NA.
  c(fit  = unname(pred$allRRfit[i]),
    low  = unname(pred$allRRlow[i]),
    high = unname(pred$allRRhigh[i]))
}

# Mise en forme "RR (IC bas–IC haut)".
fmt_rr <- function(v, d = 2) sprintf(paste0("%.", d, "f (%.", d, "f–%.", d, "f)"),
                                     v["fit"], v["low"], v["high"])

# Percentiles de référence (pour information)
pct_1  <- quantile(B$HEAT.INDEX, 0.01, na.rm = TRUE)
pct_99 <- quantile(B$HEAT.INDEX, 0.99, na.rm = TRUE)

Percentiles observés du Heat Index : 1er percentile = -0.3 °C, 99e percentile = 28.4 °C.

2 Description de la population (Table 1)

Les caractéristiques sociodémographiques et diagnostiques détaillées (sexe, âge, provenance, mode d’arrivée, diagnostic principal, orientation) ne figurent pas dans la base agrégée par jour que j’ai utilisée pour les modèles DLNM de ce rapport. La Table 1 ci-dessous reproduit les valeurs publiées dans le manuscrit (n = 69 764 consultations aux urgences psychiatriques, 2016–2023, portant sur 50 014 patients).

# Reproduction de la Table 1 du manuscrit à partir des données individuelles publiées
# (non recalculable depuis la base agrégée B). Indentation via espaces insécables
# pour que la hiérarchie soit conservée dans le rendu HTML.
ind <- function(x) paste0("\u00a0\u00a0", x)

table1 <- data.frame(
  `Characteristic` = c(
    "Total psychiatric ED consultations",
    "Gender",
    ind("Female"), ind("Male"),
    "Age (years)",
    "Provenance before consultation",
    ind("Home"), ind("General hospital"), ind("Street"),
    ind("Institution"), ind("Work"), ind("School"),
    "Mode of arrival at the ED",
    ind("Personal means"), ind("Ambulance"), ind("Police or fire brigade"),
    ind("Social services"), ind("Other"),
    "Main diagnosis (ICD-10)",
    ind("F1 – Substance use disorders"), ind("F2 – Psychotic disorders"),
    ind("F3 – Mood disorders"), ind("F4 – Anxiety disorders"),
    ind("F5 – Eating disorders"), ind("F6 – Personality disorders"),
    "Orientation after ED",
    ind("Hospitalization"), ind("Outpatient")
  ),
  n = c(
    "69,764",
    "",
    "37,138", "32,626",
    "Mean: 35.9 (15.0–98.2)",
    "",
    "40,936", "11,160", "9,460", "3,482", "1,210", "521",
    "",
    "52,657", "11,317", "2,331", "469", "1,172",
    "",
    "5,914", "15,181", "24,055", "13,781", "1,258", "5,589",
    "",
    "28,665", "41,099"
  ),
  `%` = c(
    "",
    "",
    "53.2", "46.8",
    "–",
    "",
    "60.0", "16.4", "13.9", "5.1", "1.8", "0.8",
    "",
    "77.2", "16.6", "3.4", "0.6", "1.7",
    "",
    "8.5", "21.8", "34.5", "19.8", "1.8", "8.0",
    "",
    "41.1", "58.9"
  ),
  check.names = FALSE
)

kable(table1, row.names = FALSE, align = c("l", "r", "r"),
      caption = "Table 1: Sociodemographic characteristics of patients presenting to the psychiatric emergency department during the study period. Data presented as n (%), unless otherwise specified. ICD-10: F1 substance use disorders; F2 psychotic disorders; F3 mood disorders; F4 anxiety disorders; F5 eating disorders; F6 personality disorders.")
Table 1: Sociodemographic characteristics of patients presenting to the psychiatric emergency department during the study period. Data presented as n (%), unless otherwise specified. ICD-10: F1 substance use disorders; F2 psychotic disorders; F3 mood disorders; F4 anxiety disorders; F5 eating disorders; F6 personality disorders.
Characteristic n %
Total psychiatric ED consultations 69,764
Gender
  Female 37,138 53.2
  Male 32,626 46.8
Age (years) Mean: 35.9 (15.0–98.2) –
Provenance before consultation
  Home 40,936 60.0
  General hospital 11,160 16.4
  Street 9,460 13.9
  Institution 3,482 5.1
  Work 1,210 1.8
  School 521 0.8
Mode of arrival at the ED
  Personal means 52,657 77.2
  Ambulance 11,317 16.6
  Police or fire brigade 2,331 3.4
  Social services 469 0.6
  Other 1,172 1.7
Main diagnosis (ICD-10)
  F1 – Substance use disorders 5,914 8.5
  F2 – Psychotic disorders 15,181 21.8
  F3 – Mood disorders 24,055 34.5
  F4 – Anxiety disorders 13,781 19.8
  F5 – Eating disorders 1,258 1.8
  F6 – Personality disorders 5,589 8.0
Orientation after ED
  Hospitalization 28,665 41.1
  Outpatient 41,099 58.9

3 Figure 1 — Total number of visits

Association entre le Heat Index et l’ensemble des passages aux urgences psychiatriques.

res_total <- fit_hi("nombre_passage")

plot(res_total$pred, "overall",
     col = 4, cumul = TRUE,
     xlab = "Heat Index (°C)",
     ylab = "Cumulative RR")

# Statistiques rapportées dans l'article
rr99_total <- rr_at(res_total$pred, 28.4)   # 99e pct vs référence -1
rr1_total  <- rr_at(res_total$pred, 0)    # +1 °C au-dessus de -1

Figure 1 — Cumulative Relative Risk of visits in psychiatric ED, by Heat Index, over a lag period of 2 days (centered on -1°C).

Cumulative Risks at the 99th percentile (HI = 28.4) compared to 1st percentile (-1°C).

Résultats : au 99e percentile (HI = 28.4 °C) vs 1er percentile (HI = -1 °C), RR = 1.12 (1.04–1.21). Pour chaque +1 °C au-dessus de -1 °C, RR = 1.007 (1.003–1.012). (AIC = 17 621.9.)

4 Figure 2 — Gender differences

Effet du Heat Index sur les passages selon le sexe.

res_men   <- fit_hi("nombre_homme")
res_women <- fit_hi("nombre_femme")

rr99_men   <- rr_at(res_men$pred, 28.4)
rr99_women <- rr_at(res_women$pred, 28.4)

pred_HI_H <- res_men$pred
pred_HI_F <- res_women$pred

par(mfrow = c(1, 2))

plot(pred_HI_H, "overall", col = "#1ebc12", xlab = "Heat Index (°C)",
     ylab = "Relative risk", main = "Men", ylim = c(0.8, 1.4))
points(28.4, pred_HI_H$allRRfit[which.min(abs(pred_HI_H$predvar - 28.4))],
       col = "black", pch = 19)
text(28.4, pred_HI_H$allRRfit[which.min(abs(pred_HI_H$predvar - 28.4))] + 0.02,
     labels = paste("RR =", round(pred_HI_H$allRRfit[which.min(abs(pred_HI_H$predvar - 28.4))], 2)),
     col = "black")

plot(pred_HI_F, "overall", col = "#d58f0d", xlab = "Heat Index (°C)",
     ylab = "Relative risk", main = "Women", ylim = c(0.8, 1.4))
points(28.4, pred_HI_F$allRRfit[which.min(abs(pred_HI_F$predvar - 28.4))],
       col = "black", pch = 19)
text(28.4, pred_HI_F$allRRfit[which.min(abs(pred_HI_F$predvar - 28.4))] + 0.02,
     labels = paste("RR =", round(pred_HI_F$allRRfit[which.min(abs(pred_HI_F$predvar - 28.4))], 2)),
     col = "black")

par(mfrow = c(1, 1))

Figure 2 — Cumulative Relative Risks of visits in psychiatric ED, by gender, over a lag period of 2 days (centered on -1°C).

Cumulative Risks at the 99th percentile (HI = 28.4) compared to 1st percentile (-1°C).

# Test formel de l'interaction Sexe × Heat Index (format long)
B_long <- bind_rows(
  transmute(B, DATE, HEAT.INDEX, weekday, holiday, Year,
            sexe = "Homme", y = nombre_homme),
  transmute(B, DATE, HEAT.INDEX, weekday, holiday, Year,
            sexe = "Femme", y = nombre_femme)
) %>% mutate(sexe = factor(sexe, levels = c("Homme", "Femme")))

cb <- crossbasis(B_long$HEAT.INDEX, lag = 2,
                 argvar = list(fun = "ns", df = 2),
                 arglag = list(fun = "ns", df = 2))

m_inter <- glm(y ~ cb * sexe + weekday + holiday + ns(DATE, 8*4) + ns(Year, 2),
               data = B_long, family = poisson(link = "log"))
m_noint <- glm(y ~ cb + sexe + weekday + holiday + ns(DATE, 8*4) + ns(Year, 2),
               data = B_long, family = poisson(link = "log"))

test_int <- anova(m_noint, m_inter, test = "Chisq")
p_int <- test_int$`Pr(>Chi)`[2]

Résultats : chez les femmes, RR = 1.17 (1.05–1.29) au 99e percentile ; chez les hommes, RR = 1.08 (0.97–1.20). L’interaction Sexe × Heat Index n’est pas statistiquement significative (test du rapport de vraisemblance, p = 0.343).

5 Figure 3 — Hospitalizations

RR des hospitalisations (A) et des hospitalisations sous contrainte / involuntary commitment (B).

res_hospi     <- fit_hi("nombre_hospitalisation")
res_contrainte <- fit_hi("nombre_contrainte")

rr99_hospi      <- rr_at(res_hospi$pred, 28.4)
rr99_contrainte <- rr_at(res_contrainte$pred, 28.4)

par(mfrow = c(1, 2))

plot(res_hospi$pred, "overall", col = 4,
     xlab = "Heat Index (°C)", ylab = "Cumulative RR")
mtext("A.", side = 1, line = 3, adj = 0, font = 2, cex = 1.1)

plot(res_contrainte$pred, "overall", col = "red",
     xlab = "Heat Index (°C)", ylab = "Cumulative RR")
mtext("B.", side = 1, line = 3, adj = 0, font = 2, cex = 1.1)

par(mfrow = c(1, 1))

Figure 3 - Cumulative Relative Risks of hospitalization (A) and involuntary commitment (B), by Heat Index, over a lag period of 2 days.

Cumulative Risks at the 99th percentile (HI = 28.4) compared to 1st percentile (-1°C).

Résultats : hospitalisations (toutes), RR = 1.15 (1.02–1.29) ; hospitalisations sous contrainte, RR = 1.23 (1.04–1.47).

6 Figure 4 — Diagnostic categories

RR des passages par catégorie diagnostique : (A) F1 troubles liés aux substances, (B) F2 troubles psychotiques, (C) F3 troubles de l’humeur, (D) F4 troubles anxieux.

res_F1 <- fit_hi("nombre_F1")
res_F2 <- fit_hi("nombre_F2")
res_F3 <- fit_hi("nombre_F3")
res_F4 <- fit_hi("nombre_F4")

par(mfrow = c(2, 2))

plot(res_F1$pred, "overall", col = "purple",
     xlab = "Heat Index (°C)", ylab = "Cumulative RR")
mtext("A.", side = 1, line = 3, adj = 0, font = 2, cex = 1.1)

plot(res_F2$pred, "overall", col = "orange",
     xlab = "Heat Index (°C)", ylab = "Cumulative RR")
mtext("B.", side = 1, line = 3, adj = 0, font = 2, cex = 1.1)

plot(res_F3$pred, "overall", col = "darkred",
     xlab = "Heat Index (°C)", ylab = "Cumulative RR")
mtext("C.", side = 1, line = 3, adj = 0, font = 2, cex = 1.1)

plot(res_F4$pred, "overall", col = "pink",
     xlab = "Heat Index (°C)", ylab = "Cumulative RR")
mtext("D.", side = 1, line = 3, adj = 0, font = 2, cex = 1.1)

par(mfrow = c(1, 1))

Figure 4 - Cumulative Relative Risks of emergency visits for substance use disorders (A), psychotic disorders (B), Mood Disorders (C) and anxiety disorders (D), by Heat Index, over a lag period of 2 days.

Cumulative Risks at the 99th percentile (HI = 28.4) compared to 1st percentile (-1°C).

7 Table 2 — Relative Risks by diagnosis

RR des passages par diagnostic, au 99e percentile (HI = 28.4 °C) vs 1er percentile (HI = -1 °C), et par +1 °C au-dessus de -1 °C.

diag_res <- list(
  "Substance use disorders (F1)" = res_F1,
  "Psychotic disorders (F2)"     = res_F2,
  "Mood Disorders (F3)"          = res_F3,
  "Anxiety disorders (F4)"       = res_F4
)

table2 <- do.call(rbind, lapply(names(diag_res), function(nm) {
  p <- diag_res[[nm]]$pred
  data.frame(
    `Diagnostic category` = nm,
    `Relative Risk (RR) at the 99th vs. 1st pct, and 95% Confidence Interval (IC)` = fmt_rr(rr_at(p, 28.4)),
    `RR for +1°C HI from -1°C and 95% Confidence Interval (IC)` = fmt_rr(rr_at(p, 0), 2),
    check.names = FALSE
  )
}))

kable(table2, row.names = FALSE,
      caption = "Table 2 - Relative Risks of emergency visits by diagnosis, at the 99th percentile (HI = 28.4°C) compared to the 1st percentile (HI = -1°C).")
Table 2 - Relative Risks of emergency visits by diagnosis, at the 99th percentile (HI = 28.4°C) compared to the 1st percentile (HI = -1°C).
Diagnostic category Relative Risk (RR) at the 99th vs. 1st pct, and 95% Confidence Interval (IC) RR for +1°C HI from -1°C and 95% Confidence Interval (IC)
Substance use disorders (F1) 1.27 (0.98–1.65) 1.02 (1.01–1.04)
Psychotic disorders (F2) 1.17 (1.00–1.38) 1.01 (1.00–1.02)
Mood Disorders (F3) 1.11 (0.97–1.26) 1.00 (1.00–1.01)
Anxiety disorders (F4) 1.10 (0.93–1.30) 1.00 (0.99–1.01)

8 Figure S1 — Sensitivity analysis (supplementary)

Analyse de sensibilité de l’association Heat Index — passages aux urgences psychiatriques, en faisant varier les degrés de liberté et le décalage maximal.

var_df_values  <- c(1, 2, 3)
lag_df_values  <- c(2, 3)
lag_max_values <- c(2, 3)

par(mfrow = c(3, 4), mar = c(4, 4, 3, 1))

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)))

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]
      col_idx <- (j - 1) * length(lag_max_values) + k

      cb <- crossbasis(B$HEAT.INDEX, lag = lag_max,
                       argvar = list(fun = "ns", df = var_df),
                       arglag = list(fun = "ns", df = lag_df))
      model_hi <- glm(nombre_passage ~ cb + weekday + holiday +
                        ns(DATE, 8*4) + ns(Year, 2),
                      data = B, family = poisson(link = "log"))
      aic_matrix[i, col_idx] <- AIC(model_hi)

      pred.hi <- crosspred(cb, model_hi, at = seq(-5, 35, by = 0.1),
                           cumul = TRUE, cen = -1)
      plot(pred.hi$predvar, pred.hi$allRRfit, type = "l", col = 4, lwd = 2,
           xlab = "Heat Index (°C)", ylab = "Cumulative RR",
           main = paste0("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 = 4, lty = 2)
      lines(pred.hi$predvar, pred.hi$allRRhigh, col = 4, lty = 2)
      abline(h = 1, col = "black")
    }
  }
}

par(mfrow = c(1, 1))

Figure S1: Sensitivity analysis of the association between the heat index and psychiatric emergency visits using different degrees of freedom.

kable(round(aic_matrix, 1),
      caption = "AIC des différentes spécifications (analyse de sensibilité, total des passages).")
AIC des différentes spécifications (analyse de sensibilité, total des passages).
lag_df = 2 , lag_max = 2 lag_df = 2 , lag_max = 3 lag_df = 3 , lag_max = 2 lag_df = 3 , lag_max = 3
var_df = 1 17620.6 17614.5 17622.6 17615.3
var_df = 2 17621.9 17615.8 17625.8 17618.5
var_df = 3 17625.7 17619.5 17631.5 17624.0 
best_idx <- which(aic_matrix == min(aic_matrix), arr.ind = TRUE)

Le modèle retenu (var_df = 2, lag_df = 2, lag_max = 2) offre des estimations robustes et interprétables ; les spécifications plus complexes (df = 3, lag_max = 3) produisent des courbes de forme similaire sans amélioration substantielle de l’ajustement.

9 Réponses aux relecteurs

9.1 Commentaire 1 — Description des températures

Commentaire 1. When describing the data, please go into more detail regarding temperatures, and visualize this information in an appropriate manner (this can also be done in the supplement). You state „caution given widened confidence intervals at extreme values, reflecting sparse data”. Yet this is relevant for transdisciplinary science communication as you may be able to demonstrate a rising temperature trend in your data. According to data available from iied, during 2014–2023, 6 out of 10 years experienced 20 or more days of at least 30°C (https://www.iied.org/sites/default/files/uploads/2024/06/Hot_cities_Paris.pdf). Please illustrate the actual temperature and the number of days above 30 degrees for example, using scatter plots or violin plots.

9.1.1 Statistiques descriptives des indicateurs de température et d’humidité

# Résumé des variables thermiques et hygrométriques sur toute la période d'étude.
temp_summary <- data.frame(
  Variable = c("Daily mean temperature (°C)",
               "Daily maximum temperature (°C)",
               "Daily minimum temperature (°C)",
               "Heat Index (°C)",
               "Relative humidity (%)"),
  Min    = c(min(B$TNTXM),    min(B$TX),    min(B$TN),    min(B$HEAT.INDEX, na.rm = TRUE),    min(B$UM)),
  Mean   = c(mean(B$TNTXM),   mean(B$TX),   mean(B$TN),   mean(B$HEAT.INDEX, na.rm = TRUE),   mean(B$UM)),
  Median = c(median(B$TNTXM), median(B$TX), median(B$TN), median(B$HEAT.INDEX, na.rm = TRUE), median(B$UM)),
  P99    = c(quantile(B$TNTXM, .99), quantile(B$TX, .99), quantile(B$TN, .99),
             quantile(B$HEAT.INDEX, .99, na.rm = TRUE), quantile(B$UM, .99)),
  Max    = c(max(B$TNTXM),    max(B$TX),    max(B$TN),    max(B$HEAT.INDEX, na.rm = TRUE),    max(B$UM)),
  check.names = FALSE
)
temp_summary[, -1] <- round(temp_summary[, -1], 1)

kable(temp_summary, row.names = FALSE, align = c("l", rep("r", 5)),
      caption = "Table S1 — Descriptive statistics of daily temperature, Heat Index and relative humidity (Paris-Montsouris, 2016–2023).")
Table S1 — Descriptive statistics of daily temperature, Heat Index and relative humidity (Paris-Montsouris, 2016–2023).
Variable Min Mean Median P99 Max
Daily mean temperature (°C) -3.8 13.7 13.3 27.9 33.8
Daily maximum temperature (°C) -2.4 17.5 17.2 34.9 42.6
Daily minimum temperature (°C) -7.0 9.9 9.8 21.2 25.0
Heat Index (°C) -3.8 13.1 12.6 28.4 33.6
Relative humidity (%) 31.0 70.3 71.5 94.0 98.0

Sur la période 2016–2023, la température moyenne quotidienne était de 13.7 °C (min -3.8 °C ; max 33.8 °C). Le Heat Index avait une médiane de 12.6 °C et un 99e percentile de 28.4 °C, qui correspond à la valeur de comparaison « haute » (HI = 28.4 °C) utilisée dans les modèles (vs. le 1er percentile, -1 °C).

9.1.2 Nombre de jours chauds par année

days30 <- B %>%
  filter(!is.na(TX)) %>%
  group_by(Year) %>%
  summarise(`Days Tmax ≥ 30°C`   = sum(TX >= 30),
            `Days HI ≥ 30°C`     = sum(HEAT.INDEX >= 30, na.rm = TRUE),
            `Max Tmax (°C)`      = round(max(TX), 1),
            .groups = "drop") %>%
  rename(Year = Year)

n_years_20 <- sum(days30$`Days Tmax ≥ 30°C` >= 20)

kable(days30, row.names = FALSE, align = c("l", "r", "r", "r"),
      caption = "Table S2 — Number of hot days per year (maximum temperature and Heat Index ≥ 30°C), Paris-Montsouris.")
Table S2 — Number of hot days per year (maximum temperature and Heat Index ≥ 30°C), Paris-Montsouris.
Year Days Tmax ≥ 30°C Days HI ≥ 30°C Max Tmax (°C)
2016 14 0 36.6
2017 17 2 36.9
2018 26 1 37.4
2019 20 2 42.6
2020 23 3 39.3
2021 12 0 33.3
2022 23 1 40.5
2023 26 0 35.5

Sur nos 8 années d’étude (2016–2023), 5 années ont connu ≥ 20 jours avec une température maximale ≥ 30 °C. Ce résultat corrobore, sur notre propre jeu de données, la tendance rapportée par l’IIED (6 années sur 10 pour 2014–2023), et illustre la fréquence élevée et récurrente des jours de forte chaleur à Paris.

9.1.3 Figure S2 — Days above 30°C per year

ggplot(days30, aes(x = factor(Year), y = `Days Tmax ≥ 30°C`)) +
  geom_col(fill = "#c0392b", width = 0.7) +
  geom_hline(yintercept = 20, linetype = "dashed", colour = "grey30") +
  geom_text(aes(label = `Days Tmax ≥ 30°C`), vjust = -0.4, size = 3.6) +
  annotate("text", x = 0.7, y = 21.3, label = "20 days (IIED reference)",
           hjust = 0, size = 3, colour = "grey30") +
  labs(x = "Year", y = "Number of days with Tmax ≥ 30°C",
       title = "Figure S2 — Days per year with maximum temperature ≥ 30°C",
       subtitle = "Paris-Montsouris weather station, 2016–2023") +
  expand_limits(y = max(days30$`Days Tmax ≥ 30°C`) + 3)

Figure S2 — Number of days per year with a daily maximum temperature ≥ 30 °C (Paris-Montsouris, 2016–2023). The dashed line marks the 20-day threshold reported by the IIED for Paris.

9.1.4 Figure S3 — Distribution of temperature by year (violin plots)

ggplot(B, aes(x = factor(Year), y = TNTXM)) +
  geom_violin(fill = "#3498db", alpha = 0.45, colour = NA) +
  geom_boxplot(width = 0.12, outlier.size = 0.5, alpha = 0.8) +
  labs(x = "Year", y = "Daily mean temperature (°C)",
       title = "Figure S3 — Distribution of daily mean temperature by year",
       subtitle = "Paris-Montsouris weather station, 2016–2023")

Figure S3 — Yearly distribution of daily mean temperature (violin + box plots). The width of each violin reflects the density of days at a given temperature; the upper tails correspond to the (relatively rare) extreme-heat days.

9.1.5 Figure S4 — Heat Index vs. temperature (scatter plot)

ggplot(B, aes(x = TNTXM, y = HEAT.INDEX, colour = UM)) +
  geom_point(alpha = 0.5, size = 1) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = "grey40") +
  scale_colour_viridis_c(name = "Relative\nhumidity (%)", option = "C") +
  labs(x = "Daily mean temperature (°C)", y = "Heat Index (°C)",
       title = "Figure S4 — Heat Index vs. daily mean temperature",
       subtitle = "Coloured by relative humidity; dashed line = identity (HI = T)")

Figure S4 — Relationship between the Heat Index and daily mean temperature, coloured by relative humidity. At higher temperatures, the Heat Index increasingly exceeds air temperature when humidity is high (points above the identity line), illustrating the added value of the Heat Index over temperature alone.

9.2 Commentaire 2 — Données brutes de température et d’humidité

Commentaire 2. “It would be helpful to include the extracted raw data on humidity and temperature in the supplement, especially since the HI is yet less well known.” On peut donner, dans les suppléments, l’ensemble des données météorologiques quotidiennes ayant servi à calculer le Heat Index (température moyenne, maximale, minimale, humidité relative et HI dérivé), sous forme de fichiers (CSV et Excel) dont voici un aperçu ci-dessous :

# Jeu de données brut quotidien (T° + humidité + HI) exporté pour le supplément.
raw_weather <- data.frame(
  Date                  = B$DATE,
  Year                  = B$Year,
  Mean_temperature_C    = round(B$TNTXM, 1),
  Max_temperature_C     = round(B$TX, 1),
  Min_temperature_C     = round(B$TN, 1),
  Relative_humidity_pct = round(B$UM, 0),
  Heat_index_C          = round(B$HEAT.INDEX, 1)
)
raw_weather <- raw_weather[order(raw_weather$Date), ]

# Dossier de destination (à côté du manuscrit)
supp_dir <- "/home/pablo/Dropbox/CPOA - Saint Anne/Marine AKKAOUI/article soumis"
csv_path  <- file.path(supp_dir, "Supplementary_Data_Weather_2016-2023.csv")
xlsx_path <- file.path(supp_dir, "Supplementary_Data_Weather_2016-2023.xlsx")

write.csv(raw_weather, csv_path, row.names = FALSE, fileEncoding = "UTF-8")
if (requireNamespace("openxlsx", quietly = TRUE)) {
  openxlsx::write.xlsx(raw_weather, xlsx_path)
}

Le jeu de données supplémentaire contient 2910 jours couvrant la période du 12/01/2016 au 30/12/2023. Il a été exporté en fichier CSV.

Aperçu des 10 premières lignes :

kable(head(raw_weather, 10), row.names = FALSE,
      caption = "Table S3 — Extract of the supplementary raw weather dataset (daily temperature, relative humidity and derived Heat Index; Paris-Montsouris, 2016–2023).")
Table S3 — Extract of the supplementary raw weather dataset (daily temperature, relative humidity and derived Heat Index; Paris-Montsouris, 2016–2023).
Date Year Mean_temperature_C Max_temperature_C Min_temperature_C Relative_humidity_pct Heat_index_C
2016-01-12 2016 7.1 8.4 5.7 74 5.8
2016-01-13 2016 5.5 7.4 3.6 73 4.0
2016-01-14 2016 5.7 8.1 3.2 77 4.3
2016-01-15 2016 4.1 5.8 2.4 68 4.1
2016-01-16 2016 3.6 5.9 1.2 70 3.6
2016-01-17 2016 3.8 4.9 2.6 72 3.8
2016-01-18 2016 0.1 1.3 -1.2 68 0.1
2016-01-19 2016 -0.4 1.6 -2.3 71 -0.4
2016-01-20 2016 -2.4 -0.4 -4.4 81 -2.4
2016-01-21 2016 -0.2 4.0 -4.4 78 -0.2

9.2.1 Figure S5 — Time series of temperature, Heat Index and humidity

long_weather <- rbind(
  data.frame(Date = B$DATE, variable = "Mean temperature (°C)", value = B$TNTXM),
  data.frame(Date = B$DATE, variable = "Heat Index (°C)",       value = B$HEAT.INDEX),
  data.frame(Date = B$DATE, variable = "Relative humidity (%)", value = B$UM)
)
long_weather$variable <- factor(long_weather$variable,
  levels = c("Mean temperature (°C)", "Heat Index (°C)", "Relative humidity (%)"))

ggplot(long_weather, aes(x = Date, y = value, colour = variable)) +
  geom_line(linewidth = 0.3) +
  facet_wrap(~ variable, ncol = 1, scales = "free_y") +
  scale_colour_manual(values = c("#e67e22", "#c0392b", "#2980b9"), guide = "none") +
  labs(x = "Date", y = NULL,
       title = "Figure S5 — Daily temperature, Heat Index and relative humidity",
       subtitle = "Paris-Montsouris weather station, 2016–2023")

Figure S5 — Daily time series of mean temperature, Heat Index and relative humidity (Paris-Montsouris, 2016–2023). This illustrates the two raw inputs (temperature and humidity) used to derive the Heat Index and their seasonal dynamics.

9.3 Commentaire 3 — Homogénéité inter-annuelle du risque

Commentaire 3. It is unclear if this is driven by a homogeneous annual temperature distribution over the course of the years, or if more high HI data points from recent years are driving this effect (as temperatures reports might suggest). Please explore this, e.g. comparing 8 yearly risk analyses on HI distributions and cumulative risk.

9.3.1 Analyses de risque annuelles (8 modèles DLNM séparés)

# Un DLNM par année : ns(DATE, 8) pour l'ajustement temporel intra-annuel,
# ns(Year, 2) est retiré (constant au sein d'une année).
years_vec <- sort(unique(B$Year))

yearly_fits <- lapply(years_vec, function(y) {
  d  <- B[B$Year == y, ]
  cb <- crossbasis(d$HEAT.INDEX, lag = 2,
                   argvar = list(fun = "ns", df = 2),
                   arglag = list(fun = "ns", df = 2))
  m  <- glm(nombre_passage ~ cb + weekday + holiday + ns(DATE, 8),
            data = d, family = poisson(link = "log"))
  crosspred(cb, m, at = seq(-5, 35, by = 0.1), cen = -1, cumul = TRUE)
})
names(yearly_fits) <- years_vec

# RR cumulé pour chaque année, selon DEUX définitions de la valeur "haute" :
#  (A) contraste fixe commun : HI = 28.4 °C (99e pct groupé) vs -1 °C
#  (B) 99e percentile PROPRE à chaque année vs -1 °C (analyse de robustesse)
rr29_year <- do.call(rbind, lapply(years_vec, function(y) {
  p    <- yearly_fits[[as.character(y)]]
  p99y <- as.numeric(quantile(B$HEAT.INDEX[B$Year == y], .99, na.rm = TRUE))
  iA   <- which.min(abs(p$predvar - 28.4))
  iB   <- which.min(abs(p$predvar - p99y))
  data.frame(Year = y,
             # (A) contraste fixe (utilisé aussi par le forest plot Figure S8)
             RR   = round(unname(p$allRRfit[iA]),  2),
             low  = round(unname(p$allRRlow[iA]),  2),
             high = round(unname(p$allRRhigh[iA]), 2),
             # (B) 99e percentile annuel
             p99_year = round(p99y, 1),
             RR_own   = round(unname(p$allRRfit[iB]),  2),
             low_own  = round(unname(p$allRRlow[iB]),  2),
             high_own = round(unname(p$allRRhigh[iB]), 2))
}))

# Table S5 : présentation des deux définitions côte à côte
tableS5 <- data.frame(
  Year = rr29_year$Year,
  `RR at fixed HI = 28.4 °C (95% CI)` =
    sprintf("%.2f (%.2f-%.2f)", rr29_year$RR, rr29_year$low, rr29_year$high),
  `Yearly 99th pct (°C)` = rr29_year$p99_year,
  `RR at yearly 99th pct (95% CI)` =
    sprintf("%.2f (%.2f-%.2f)", rr29_year$RR_own, rr29_year$low_own, rr29_year$high_own),
  check.names = FALSE
)

kable(tableS5, row.names = FALSE, align = c("l", "r", "r", "r"),
      caption = "Table S5 — Yearly cumulative RR vs. -1 °C, using two definitions of the high-exposure value: (A) a fixed common contrast at HI = 28.4 °C (pooled 99th percentile), and (B) each year's own 99th percentile of HI (robustness).")
Table S5 — Yearly cumulative RR vs. -1 °C, using two definitions of the high-exposure value: (A) a fixed common contrast at HI = 28.4 °C (pooled 99th percentile), and (B) each year’s own 99th percentile of HI (robustness).
Year RR at fixed HI = 28.4 °C (95% CI) Yearly 99th pct (°C) RR at yearly 99th pct (95% CI)
2016 1.42 (1.09-1.85) 28.0 1.41 (1.09-1.84)
2017 1.11 (0.88-1.41) 28.0 1.11 (0.88-1.40)
2018 1.26 (1.03-1.55) 28.3 1.26 (1.03-1.55)
2019 1.04 (0.81-1.32) 28.8 1.03 (0.80-1.31)
2020 1.36 (1.06-1.74) 29.8 1.40 (1.08-1.80)
2021 1.06 (0.84-1.34) 24.3 1.04 (0.86-1.27)
2022 0.98 (0.79-1.20) 28.3 0.97 (0.79-1.20)
2023 1.58 (1.25-2.00) 28.1 1.57 (1.24-1.98)

Deux définitions de la valeur « haute » sont présentées. (A) Le contraste fixe (28,4 °C pour toutes les années) maintient le même niveau d’exposition d’une année à l’autre : il isole l’effet de la chaleur des différences de distribution entre années, et constitue donc l’analyse la plus interprétable pour la question posée (homogénéité de l’effet). (B) Le 99e percentile propre à chaque année reflète l’extrême réellement observé cette année-là, mais mélange l’effet et la variation du percentile lui-même (p. ex. 2021, année fraîche, dont le 99e percentile n’est que de 24.3 °C). Les deux définitions donnent des résultats quasi identiques, ce qui confirme la robustesse : le RR annuel ne dépend pas du choix de la valeur de comparaison.

9.3.2 Figure S8 — Cumulative RR at HI = 28.4 °C, by year (forest plot)

ggplot(rr29_year, aes(x = RR, y = factor(Year))) +
  geom_vline(xintercept = 1, linetype = "dashed", colour = "grey50") +
  geom_errorbarh(aes(xmin = low, xmax = high), height = 0.25, colour = "#2c3e50") +
  geom_point(size = 2.6, colour = "#c0392b") +
  labs(x = "Cumulative RR at HI = 28.4 °C (vs. -1 °C)", y = "Year",
       title = "Figure S8 — Yearly cumulative risk at HI = 28.4 °C",
       subtitle = "Separate DLNM per year; wide CIs reflect sparse extreme-HI days")

Figure S8 — Risques relatifs cumulés à HI = 28.4 °C, estimés année par année (8 DLNM séparés). Les estimations fluctuent autour de la valeur groupée sans tendance temporelle monotone, et les intervalles de confiance sont larges, reflétant la rareté des jours de HI extrême au sein d’une seule année. Cela justifie l’analyse groupée sur l’ensemble de la période, plus stable et mieux alimentée aux valeurs extrêmes.

9.3.3 Figure S9 — Cumulative RR curves by year (overlay)

curves <- do.call(rbind, lapply(years_vec, function(y) {
  p <- yearly_fits[[as.character(y)]]
  data.frame(Year = factor(y), HI = p$predvar, RR = p$allRRfit)
}))

ggplot(curves, aes(x = HI, y = RR, colour = Year)) +
  geom_hline(yintercept = 1, linetype = "dashed", colour = "grey60") +
  geom_line(linewidth = 0.7) +
  coord_cartesian(ylim = c(0.7, 1.8)) +
  scale_colour_viridis_d() +
  labs(x = "Heat Index (°C)", y = "Cumulative RR (centered on -1 °C)",
       title = "Figure S9 — Exposure–response curves by year",
       subtitle = "Separate DLNM per year, 2016–2023")

Figure S9 — Courbes exposition–réponse (RR cumulé) estimées séparément pour chaque année. Les courbes annuelles ont une forme globalement cohérente (augmentation du risque avec le HI), avec une variabilité attendue liée à la taille d’échantillon annuelle.

9.4 Commentaire 4 — Distribution annuelle du Heat Index

Commentaire 4. Please report the yearly HI distributions (does the 99th percentile increase over time?)

9.4.1 Table S4 — Distribution annuelle du Heat Index

hi_year <- B %>%
  filter(!is.na(HEAT.INDEX)) %>%
  group_by(Year) %>%
  summarise(N = n(),
            Mean   = round(mean(HEAT.INDEX), 1),
            Median = round(median(HEAT.INDEX), 1),
            P90    = round(quantile(HEAT.INDEX, .90), 1),
            P95    = round(quantile(HEAT.INDEX, .95), 1),
            P99    = round(quantile(HEAT.INDEX, .99), 1),
            Max    = round(max(HEAT.INDEX), 1),
            .groups = "drop")

global_p99 <- round(quantile(B$HEAT.INDEX, .99, na.rm = TRUE), 1)

kable(hi_year, row.names = FALSE, align = c("l", rep("r", 7)),
      caption = paste0("Table S4 — Yearly distribution of the Heat Index (°C), Paris-Montsouris. Pooled 99th percentile = ", global_p99, " °C."))
Table S4 — Yearly distribution of the Heat Index (°C), Paris-Montsouris. Pooled 99th percentile = 28.4 °C.
Year N Mean Median P90 P95 P99 Max
2016 355 12.3 11.7 21.8 23.1 28.0 29.0
2017 365 12.9 12.6 22.1 24.1 28.0 30.8
2018 365 13.4 12.9 23.1 25.1 28.3 30.8
2019 365 13.0 12.1 22.1 24.4 28.8 33.6
2020 366 13.6 13.4 22.0 24.1 29.8 31.9
2021 365 12.2 11.6 21.2 23.3 24.3 27.4
2022 365 13.7 13.9 23.2 24.3 28.3 31.5
2023 364 13.8 13.0 23.2 24.9 28.1 29.1

Le 99e percentile global du HI est de 28.4 °C, qui est la valeur de comparaison « haute » utilisée dans les modèles de l’article. Les 99e percentiles annuels sont remarquablement stables (28, 28, 28.3, 28.8, 29.8, 24.3, 28.3, 28.1 °C), à l’exception de 2021, année plus fraîche (P99 = 24.3 °C).

9.4.2 Tendance temporelle du 99e percentile

lm_p99 <- lm(P99 ~ Year, data = hi_year)
slope_p99 <- round(coef(lm_p99)[2], 3)
p_p99     <- round(summary(lm_p99)$coefficients[2, 4], 3)
rho_p99   <- round(cor(hi_year$Year, hi_year$P99, method = "spearman"), 2)

Il n’y a pas de tendance croissante du 99e percentile annuel sur 2016–2023 : la pente estimée est de -0.105 °C/an (régression linéaire, p = 0.703) et la corrélation de rang de Spearman est faible (ρ = 0.16). Autrement dit, le 99e percentile global n’est pas tiré par les seules années récentes : les jours de HI élevé se répartissent de façon relativement homogène sur l’ensemble de la période.

9.4.3 Figure S6 — Distribution annuelle du Heat Index (violin plots)

ggplot(B[!is.na(B$HEAT.INDEX), ], aes(x = factor(Year), y = HEAT.INDEX)) +
  geom_violin(fill = "#8e44ad", alpha = 0.4, colour = NA) +
  geom_boxplot(width = 0.12, outlier.size = 0.5, alpha = 0.8) +
  geom_hline(yintercept = global_p99, linetype = "dashed", colour = "#c0392b") +
  annotate("text", x = 0.7, y = global_p99 + 1, hjust = 0, size = 3, colour = "#c0392b",
           label = paste0("Pooled 99th pct = ", global_p99, "°C")) +
  labs(x = "Year", y = "Heat Index (°C)",
       title = "Figure S6 — Yearly distribution of the Heat Index",
       subtitle = "Paris-Montsouris, 2016–2023; dashed line = pooled 99th percentile")

Figure S6 — Distribution annuelle du Heat Index (violin + box plots). La ligne rouge marque le 99e percentile groupé (28.4 °C). Les distributions annuelles sont comparables d’une année à l’autre, sans dérive nette du haut de la distribution.

9.4.4 Figure S7 — 99e percentile annuel du HI dans le temps

ggplot(hi_year, aes(x = Year, y = P99)) +
  geom_hline(yintercept = global_p99, linetype = "dashed", colour = "grey50") +
  geom_smooth(method = "lm", se = TRUE, colour = "#c0392b", fill = "#c0392b", alpha = 0.12) +
  geom_point(size = 2.5, colour = "#2c3e50") +
  geom_line(colour = "#2c3e50", alpha = 0.5) +
  scale_x_continuous(breaks = hi_year$Year) +
  labs(x = "Year", y = "Yearly 99th percentile of HI (°C)",
       title = "Figure S7 — Yearly 99th percentile of the Heat Index over time",
       subtitle = paste0("Linear trend: ", slope_p99, " °C/year, p = ", p_p99))

Figure S7 — 99e percentile annuel du Heat Index (2016–2023) avec tendance linéaire. La pente n’est pas significative (-0.105 °C/an, p = 0.703), indiquant l’absence d’augmentation du 99e percentile sur la période étudiée.