Ce document a été généré avec l’outil R Markdown. Le code R et les données qui ont été utilisées sont ainsi mis à disposition et permettent donc la reproductibilité des résultats obtenus.

Par ailleurs, le document est mis à jour automatiquement chaque jour. Pour consulter les archives, cliquer ici.

Source de données utilisées:

Questions de recherche:

En fonction des réponses aux questions précédentes :

Packages et fonctions locales à charger:

library(cartogram)
library(cartography)
library(forecast)
library(kableExtra)
library(rgdal)
library(tidyverse)
library(vistime)
library(zoo)
source("fonctions.R")

Version de R utilisée:

R.version
##                _                           
## platform       x86_64-pc-linux-gnu         
## arch           x86_64                      
## os             linux-gnu                   
## system         x86_64, linux-gnu           
## status                                     
## major          4                           
## minor          0.3                         
## year           2020                        
## month          10                          
## day            10                          
## svn rev        79318                       
## language       R                           
## version.string R version 4.0.3 (2020-10-10)
## nickname       Bunny-Wunnies Freak Out

1 Données hospitalières relatives à l’épidémie de COVID-19

1.1 Présentation des données

1.1.1 Données par département

Dans un premier temps, on met à jour les données tous les jours de façon automatique.

Comment se présente les données du Ministère de la Santé ?

dep jour incid_hosp incid_rea incid_dc incid_rad nom_dep region
17847 69 2020-11-25 97 29 13 106 Rhône Auvergne-Rhône-Alpes
17848 69 2020-11-26 75 11 13 96 Rhône Auvergne-Rhône-Alpes
17849 69 2020-11-27 100 8 15 105 Rhône Auvergne-Rhône-Alpes
17850 69 2020-11-28 49 9 11 57 Rhône Auvergne-Rhône-Alpes

Ici il s’agit des données qui donnent chaque jour par département :

  • le nombre de nouvelles entrées en hospitalisations
  • le nombre de nouvelles entrées en réanimations
  • le nombre de décés
  • le nombre de sorties

On créé des fenêtres de 7 jour à partir du dernier jour observé. Par exemple si nous avons les données d’hospitalisation jusqu’au 28 novembre 2020 la semaine qui correspond à la semaine t0 correspond à la fenêtre [22 novembre 2020; 28 novembre 2020]. Dans chaque fenêtre, on calcule le nombre de nouvelles hospitalisations, réanimations et décès par département.

1.1.2 Données par région

On dispose égaglement des données d’hospitalisations/réanimations/décès par classe d’âge à la différence qu’il s’agit de données régionales et qu’il s’agit des données de stock (nombre d’hospitalisations et réanimations en cours) et pas du nombre de nouvelles hospitalisations. On peut toutefois estimer le nombre de nouvelles hospitalisations ou décès en faisant les différences des valeurs d’un jour sur l’autre.

my_url <- "https://www.data.gouv.fr/fr/datasets/r/08c18e08-6780-452d-9b8c-ae244ad529b3"
if (!file.exists(paste0(getwd(), "/data/age", to_day, ".csv"))) {
  download.file(my_url, destfile = paste0(getwd(), "/data/age", to_day, ".csv"))
}
hospital_age <- read.csv(paste0(getwd(), "/data/age", to_day, ".csv"), sep = ";")
# On ajoute le nom des régions:
hospital_age <- merge(hospital_age, code_region, by.x = "reg", by.y = "code")
# On utilise le format date pour coder le jour:
hospital_age$jour <- as.Date(hospital_age$jour)
# on affecte la semaine
hospital_age$semaine <- num_semaine(hospital_age$jour)

# on calcule les nouvelles hospitalisations/réanimations/décès
hospital_age$new_hosp <- 0
hospital_age$new_rea <- 0
hospital_age$new_dc <- 0

for (k in nrow(hospital_age):1) {
  age_k <- hospital_age$cl_age90[k] 
  jour_k <- hospital_age$jour[k] 
  reg_k <- hospital_age$reg[k] 
  rad_k <- hospital_age$rad[k] 
  dc_k <- hospital_age$dc[k] 
  
  ind_k <- which(hospital_age$reg == reg_k & hospital_age$cl_age90 == age_k & hospital_age$jour == jour_k - 1)
  if (length(ind_k) == 1) {
    hospital_age$new_hosp[k] <- max((hospital_age$hosp[k] - hospital_age$hosp[ind_k]) + 
                                (hospital_age$dc[k] - hospital_age$dc[ind_k]) +     
                              (hospital_age$rad[k] - hospital_age$rad[ind_k]) +
                              (hospital_age$rea[k] - hospital_age$rea[ind_k]), 0)
    hospital_age$new_rea[k] <- max((hospital_age$dc[k] - hospital_age$dc[ind_k]) +     
                              (hospital_age$rea[k] - hospital_age$rea[ind_k]), 0) 
    hospital_age$new_dc[k] <- max((hospital_age$dc[k] - hospital_age$dc[ind_k]), 0) 
  }
}
# on aggrege par semaine
my_basis_age <- hospital_age %>%
  group_by(region, semaine, cl_age90) %>%
  dplyr::summarize(hosp = sum(new_hosp),
                   rea = sum(new_rea),
                   dc = sum(new_dc),
            jour = max(jour),
            region = unique(region))
# On met au format wide
my_basis_age_wide <- tidyr::pivot_wider(my_basis_age,
                           id_cols = c("semaine", "region", "jour", "hosp", "rea", "dc", "cl_age90"),
                           names_from = "cl_age90",
                           values_from = c("hosp", "rea", "dc"))

1.2 Quelle est la situation cette semaine ?

On va calculer quelques chiffres clés pour mesure la situation des régions sur les 7 derniers jours qui viennent de s’écouler : [22 novembre 2020; 28 novembre 2020].

1.2.1 Résumé des hospitalisations

On représente par région:

  • le nombre total de nouvelles hospitalisations (semaine [22 novembre 2020; 28 novembre 2020]).

  • le nombre moyen journalier de nouvelles hospitalisations (semaine [22 novembre 2020; 28 novembre 2020]).

  • l’évolution (en pourcentage) entre la semaine [15 novembre 2020; 21 novembre 2020] et la semaine [22 novembre 2020; 28 novembre 2020].

region total semaine moyenne jour evolution en %
Auvergne-Rhône-Alpes 1937 277 -34.7
Ile-de-France 1637 234 -34.2
Grand Est 1047 150 -23.4
Hauts-de-France 992 142 -31.1
Provence-Alpes-Côte d’Azur 959 137 -33.1
Bourgogne-Franche-Comté 690 99 -28.1
Nouvelle-Aquitaine 569 81 -29.7
Occitanie 513 73 -42.2
Pays de la Loire 394 56 -27.2
Normandie 392 56 -29.6
Centre-Val de Loire 310 44 -33.5
Bretagne 136 19 -44.3
DOM-TOM 73 10 -31.1
Corse 7 1 -77.4
France entière 9656 1379 -32.4

On représente par département la carte des nouvelles hospitalisations sur la dernière semaine observée ([22 novembre 2020; 28 novembre 2020])

## OGR data source with driver: ESRI Shapefile 
## Source: "/media/thibault/My Passport/confinement/covid/departements 2015/DEPARTEMENT", layer: "DEPARTEMENT"
## with 96 features
## It has 11 fields

1.2.2 Résumé des réanimations

On représente par région:

  • le nombre total de nouvelles réanimations (semaine [22 novembre 2020; 28 novembre 2020]).

  • le nombre moyen journalier de nouvelles réanimations (semaine [22 novembre 2020; 28 novembre 2020]).

  • l’évolution (en pourcentage) entre la semaine [15 novembre 2020; 21 novembre 2020] et la semaine [22 novembre 2020; 28 novembre 2020].

region total semaine moyenne jour evolution en %
Auvergne-Rhône-Alpes 314 45 -21.9
Ile-de-France 275 39 -35
Hauts-de-France 159 23 -38.1
Provence-Alpes-Côte d’Azur 129 18 -33.2
Grand Est 125 18 -32.8
Bourgogne-Franche-Comté 120 17 -13
Occitanie 99 14 -43.1
Nouvelle-Aquitaine 58 8 -52.5
Centre-Val de Loire 56 8 -28.2
Pays de la Loire 49 7 -19.7
Normandie 38 5 -47.9
DOM-TOM 20 3 -39.4
Bretagne 18 3 -33.3
Corse 3 0 -57.1
France entière 1463 209 -32.7

On représente par département la carte des nouvelles réanimations sur la dernière semaine observée ([22 novembre 2020; 28 novembre 2020])

1.2.3 Résumé des décès

On représente par région:

  • le nombre total de nouveaux décès (semaine [22 novembre 2020; 28 novembre 2020]).

  • le nombre moyen journalier de nouveaux décès (semaine [22 novembre 2020; 28 novembre 2020]).

  • l’évolution (en pourcentage) entre la semaine [15 novembre 2020; 21 novembre 2020] et la semaine [22 novembre 2020; 28 novembre 2020].

region total semaine moyenne jour evolution en %
Auvergne-Rhône-Alpes 556 79 -13.3
Ile-de-France 366 52 -19.6
Grand Est 259 37 7.9
Hauts-de-France 248 35 1.2
Provence-Alpes-Côte d’Azur 239 34 -19.5
Bourgogne-Franche-Comté 181 26 -0.5
Nouvelle-Aquitaine 155 22 2.6
Occitanie 151 22 -21.8
Normandie 107 15 -14.4
Pays de la Loire 96 14 17.1
Centre-Val de Loire 86 12 -2.3
Bretagne 52 7 13
DOM-TOM 10 1 -33.3
Corse 0 0 -100
France entière 2506 358 -9.6

On représente par département la carte des nouveaux décès sur la dernière semaine observée ([22 novembre 2020; 28 novembre 2020])

1.3 Comment a évolué la situation depuis le début de l’épidémie ?

1.3.1 Graphique d’évolution

1.3.1.1 Hospitalisations

Ici, on représente le nombre d’entrée en hospitalisations par semaine en fonction du temps sur la France entière.

On représente la même figure mais en mettant en relief la répartition des valeurs par région :

On représente la même figure mais en mettant en relief la répartition des valeurs par classe d’âge :

On met les valeurs en pourcentages pour que le graphique soit plus visible

1.3.1.2 Réanimations

On représente le nombre cummulé d’entrée en réanimations par semaine en fonction du temps sur la France entière.

On représente la même figure mais en mettant en relief la répartition des valeurs par région :

On représente la même figure mais en mettant en relief la répartition des valeurs par classe d’âge :

On met les valeurs en pourcentages pour que le graphique soit plus visibles

Enfin, on représente le ratio réanimations / hospitalisations :

1.3.1.3 Décès

On représente le nombre cummulé de nouveaux décès par semaine en fonction du temps sur la France entière.

On représente la même figure mais en mettant en relief la répartition des valeurs par région :

On représente la même figure mais en mettant en relief la répartition des valeurs par classe d’âge :

On met les valeurs en pourcentages pour que le graphique soit plus visibles

Enfin, on représente le ratio décès / réanimations :

1.3.2 Graphique d’évolution du nombre d’hospitalisations par départements groupés par région

On va s’intéresser au nombre d’hospitalisations. On peut représenter cette information département par département. Ici, on représente le nombre cummulé d’entrée par semaine en fonction du temps.

On représente d’abord les 4 régions actuellement les plus touchées et pour lesquelles l’axe des ordonnées va de 0 à 1200.

On représente ensuite les 8 régions suivantes les plus touchées mais avec une échelle différente sur l’axe des ordonnées (0 à 400):

Enfin, on représente les 2 régions les moins touchées et avec une échelle différente sur l’axe des ordonnées (0 à 200):

1.3.3 Cartes d’évolution sur les 6 dernières semaines

On représente l’évolution des hospitalisations sur les 6 dernières semaines:

On représente l’évolution des réanimations sur les 6 dernières semaines:

On représente l’évolution des décès sur les 6 dernières semaines:

1.4 Départements avec les plus fortes évolutions en valeurs absolues par rapport à la semaine précédente

On calcule la différence entre le nombre de nouveaux patients hospitalisés sur la période [22 novembre 2020; 28 novembre 2020] et sur la période [14 novembre 2020; 21 novembre 2020]

On va représenter des couleurs différentes en fonction du taux d’évolution découpées en 5 classes

  • taux d’évolution négatif
  • taux compris entre 0 et \(50\%\)
  • taux compris entre \(50\%\) et \(100\%\)
  • taux compris entre \(100\%\) et \(200\%\)
  • taux supérieur à \(200\%\)

2 Données relatives aux résultats des tests virologiques COVID-19

On met à jour les données chaque jour :

On va calculer quelques chiffres clés pour mesurer la situation des régions sur une fenêtre de 7 jours [19 novembre 2020; 25 novembre 2020]. On ne peut pas représenter les 7 derniers jours car les données ne sont pas encore diffusées.

On représente par région:

region total semaine moyenne jour evolution en %
Auvergne-Rhône-Alpes 14677 2097 -51.6
Ile-de-France 11884 1698 -44.8
Hauts-de-France 7824 1118 -48.6
Grand Est 7681 1097 -40.9
Provence-Alpes-Côte d’Azur 6089 870 -48.9
Nouvelle-Aquitaine 5063 723 -45.1
Bourgogne-Franche-Comté 4855 694 -46.3
Occitanie 4825 689 -51.1
Normandie 3192 456 -46
Pays de la Loire 2961 423 -52
Centre-Val de Loire 2674 382 -47.2
Bretagne 1593 228 -51.4
DOM-TOM 961 137 -30.9
Corse 108 15 -63
France entière 74387 10627 -47.7

2.1 Représentation des testés positifs par tranche d’âge en fonction du temps

On représente les testés positifs par tranche d’age:

On représente les testés positifs par région :

2.2 Graphique d’évolution du nombre de détectés positifs par départements groupés par région

2.3 Choix du décallage

Hypothèse: on suppose que le nombre d’admis en hospitalisations à la semaine t0 dépend du nombre de cas testés positifs sur une fenêtre de 7 jours qui aura commencé 10 jours avant la semaine t0. Exemple: la semaine t0 est [22 novembre 2020; 28 novembre 2020], on va l’expliquer par le nombre de personnes testées positive du [12 novembre 2020; 18 novembre 2020].

3 Préparation des données pour la modélisation

On prépare ici les données pour l’étape de modélisation:

3.1 Représentation du lien entre entre le nombre d’hospitalisations et le nombre de testés positifs

Dans un premier temps, on va rerésenter les départements par des cercles de taille proportionnelle aux nombres de testés positifs la semaine du [12 novembre 2020; 18 novembre 2020]. La couleur dépend du nombre d’hospitalisations observés la semaine du [22 novembre 2020; 28 novembre 2020].

On représente le nombre de nouvelles hospitalisations par semaine et par département en fonction du nombre de personnes testées positives quelques jours auparavant et on constate un lien très fort.

4 Prédire le nombre de testés positifs

On rappelle que les données sur le nombre de testés positifs ne sont disponible que jusqu’au 25 novembre 2020. Notre objectif est de prédire le nombre de testés positifs du 26 novembre 2020 au 02 décembre 2020 en utilisant des modèles de séries temporelles. En utilisant un modèle de série temporelle on suppose que ce qu’on observe à la date \(j\) dépend de ce qu’il s’est passé les dates antérieures. On va utiliser 3 modèles différents et en fonction de leur performence (sur les données passées), on va leur donner plus ou moins d’importance.

4.1 Modèle de type Box-Jenkins

Ici, on considère les données journalières, et non hebdomadaires. On va expliquer \(y_{d, t}^a\), le nombre de testés positifs le jour \(t\) dans le département \(d\) et dans la tranche d’âge \(a\). La stratégie utilisée est la suivante :

  • on différencie chaque série pour les rendre stationnaire (on ne vérifiera pas l’hypothèse de stationarité après la différenciation car on modélise énormément de modèle, ici on a \(A\times D\) séries où \(A\) est le nombre de classe d’âge et \(D\) le nombre de département et notre but est d’avoir une procédure automatique)

  • on cherche le meilleur modèle \(ARIMA(p,d,q)\) selon le critère AIC, à l’aide de la fonction auto.arima() (package forecast)

  • on prédit sur les 7 prochains jours à venir et on cummule ces prédictions pour avoir une prédiction du nombre de cas positifs sur la semaine à venir.

4.2 Modèle de type Lissage exponentiel

On va appliquer deux modèles de lissage exponentiels:

  • un modèle journalier qui va permettre de modéliser \(y_{d, t}^a\), le nombre de testés positifs le jour \(t\) dans le département \(d\) et dans la tranche d’âge \(a\) afin de prédire le nombre de testés positifs dans les 7 jours.

  • un modèle hebdomadaire qui va permettre de modéliser \(y_{d, s}^a\), le nombre de testés positifs la semaine \(s\) dans le département \(d\) et dans la tranche d’âge \(a\) afin de prédire le nombre de testés positifs la semaine à venir.

4.3 Combinaison des prédictions

On apprentit les modèles ci-dessus en enlevant la dernière semaine observée dans le but de donner des poids différents aux trois modèles de prédictions utilisés. Ainsi, on donnera davantage de poids aux modèles qui ont mieux prédit la dernière semaine observée.

# prediction par department 
nom_dep <- my_basis[my_basis$semaine == "semaine_t00", "dep"]
pred_cas <- numeric(length(nom_dep))
my_tab <- data.frame(true_P = numeric(0), pred_1 = numeric(0), pred_2 = numeric(0), pred_3 = numeric(0))
      
# apprentissage
for (k in length(nom_dep):1) {
  if (nom_dep[k] %in% c("975", "977", "978")) {
    my_basis <- rbind(data_k, my_basis)
  } else {
    for (age in c(0, 9, 19, 29, 39, 49, 59, 69, 79, 89, 90)) {
      
      # apprentissage
      temp <- test[test$dep == nom_dep[k] & test$cl_age90 == age & test$jour <= max(test$jour) - 7, ]
      my_ts <- zoo(temp$P, temp$jour)
      
      # Methode 1 : ARIMA
      my_ts_diff <- diff(my_ts)
      # tseries::adf.test(my_ts) 
      # tseries::adf.test(my_ts_diff)
      my_mod <- forecast::auto.arima(my_ts_diff)
      forecast_my_mod <- as.numeric(forecast(my_mod)$mean)
      pred_1 <- round(sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:7]), 0)
      # Méthode 2 : lissage exponentiel
      my_mod_exp <- ets(my_ts)
      forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp)$mean)
      forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
      pred_2 <- round(sum(forecast_my_mod_exp[1:7]), 0)
      # Méthode 3 : lissage exponentiel sur données hebdomadaires
      temp <- my_basis[my_basis$dep == nom_dep[k] & !(my_basis$semaine %in% c("semaine_t0-1", "semaine_t0-2")), ]
      my_ts_exp <- zoo(temp[ , paste0("tranche_", age)], temp$jour)
      my_mod_exp_2 <- ets(my_ts_exp)
      forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp_2)$mean)
      forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
      pred_3 <- round(forecast_my_mod_exp[1], 0)
      
      true_P <- sum(test[which(test$dep == nom_dep[k] & test$cl_age90 == age & 
                               test$jour > (max(test$jour) - 7)), "P"])
      my_tab <- rbind(my_tab, data.frame(true_P = true_P, pred_1 = pred_1, pred_2 = pred_2, pred_3 = pred_3))
    }
  }
}

res_lm_cas <- lm(true_P ~ pred_1 + pred_2 + pred_3, data = my_tab)

for (k in length(nom_dep):1) {
  data_k <- data.frame(dep = nom_dep[k], semaine = "semaine_t0-2", hosp = NA, rea = NA, rad = NA, dc = NA,
                         jour = to_day + 13, region = dep_region[match(nom_dep[k], dep_region$dep) , "region"],
                         tranche_9 = NA, tranche_19 = NA,  tranche_29 = NA,  tranche_39 = NA,  tranche_49 = NA, 
                         tranche_59 = NA,  tranche_69 = NA,  tranche_79 = NA,  tranche_89 = NA,  tranche_90 = NA, 
                         tranche_0  = NA)
  data_k_2 <- data.frame(dep = nom_dep[k], semaine = "semaine_t0-2", hosp = NA, rea = NA, rad = NA, dc = NA,
                         jour = to_day + 13, region = dep_region[match(nom_dep[k], dep_region$dep) , "region"],
                         tranche_9 = NA, tranche_19 = NA,  tranche_29 = NA,  tranche_39 = NA,  tranche_49 = NA, 
                         tranche_59 = NA,  tranche_69 = NA,  tranche_79 = NA,  tranche_89 = NA,  tranche_90 = NA, 
                         tranche_0  = NA)
      
  if (nom_dep[k] %in% c("975", "977", "978")) {
    my_basis <- rbind(data_k, my_basis)
  } else {
    for (age in c(0, 9, 19, 29, 39, 49, 59, 69, 79, 89, 90)) {
      # modèle journaliers 
      temp <- test[test$dep == nom_dep[k] & test$cl_age90 == age, ]
      my_ts <- zoo(temp$P, temp$jour)
      my_ts_diff <- diff(my_ts)
      # tseries::adf.test(my_ts) 
      # tseries::adf.test(my_ts_diff)
      my_mod <- forecast::auto.arima(my_ts_diff)
      forecast_my_mod <- as.numeric(forecast(my_mod)$mean)
      pred_1 <- round(sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:7]), 0)
      # modèles exponentiels
      # Méthode 2 : lissage exponentiel
      my_mod_exp <- ets(my_ts)
      forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp)$mean)
      forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
      pred_2 <- round(sum(forecast_my_mod_exp[1:7]), 0)
      # Méthode 3 : lissage exponentiel sur données hebdomadaires
      temp <- my_basis[my_basis$dep == nom_dep[k] & !(my_basis$semaine %in% c("semaine_t0-2")), ]
      my_ts_exp <- zoo(temp[ , paste0("tranche_", age)], temp$jour)
      my_mod_exp_2 <- ets(my_ts_exp)
      forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp_2)$mean)
      forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
      pred_3 <- round(forecast_my_mod_exp[1], 0)
      data_k[ , paste0("tranche_", age)] <- predict(res_lm_cas, newdata = data.frame(pred_1 = pred_1,
                                                                                     pred_2 = pred_2,
                                                                                     pred_3 = pred_3))
    }
    my_basis <- rbind(data_k, my_basis)
  }
}

On représente les testés positifs par région en ajoutant les valeurs de la semaine prédite:

On aggrège les données à la France entière:

5 Prédire le nombre d’hospitalisation de la semaine à venir

5.1 Modèle linéaire 1 (sur les départements) en fonction du nombre de cas détectés positifs : 1 modèle par région

Ici, pour chaque région \(r\), le modèle est de la forme

\[y_{i,t}^r=\beta_0^r+\beta_1^rx_{i,t'}^r+\epsilon_{i,t}^r\] avec:

  • \(y_{i,t}\) le nombre d’entrées à l’hôpital dans le département \(i\in r\) sur la période \(t\), où \(t\) est une fenêtre de 7 jours.
  • \(x_{i,t'}\) est le nombre de testés positifs dans le département \(i\in r\) sur la période \(t'\) où \(t'\) correspond à la fenêtre \(t\), décalé de 10 jours.

En d’autres termes, on fait ici un modèle de régression par région. Cela suppose que le lien entre les tests virologiques et le nombre d’hospitalisation est homogène à l’intérieur d’une région et peut différer d’une région à une autre.

Apprentissage:

On modélise sur les observations des semaines précédentes:

Dependent variable:
hosp
regionAuvergne-Rhône-Alpes 8.093***
(1.629)
regionBourgogne-Franche-Comté 3.550*
(2.017)
regionBretagne 1.707
(3.037)
regionCentre-Val de Loire 4.142*
(2.398)
regionCorse 1.535
(4.346)
regionDOM-TOM 2.091
(3.035)
regionGrand Est 9.509***
(1.801)
regionHauts-de-France 12.043***
(2.526)
regionIle-de-France 22.963***
(2.298)
regionNormandie 2.994
(2.609)
regionNouvelle-Aquitaine 1.361
(1.667)
regionOccitanie 1.288
(1.586)
regionPays de la Loire 7.253***
(2.655)
regionProvence-Alpes-Côte d’Azur 4.404*
(2.371)
regionAuvergne-Rhône-Alpes:tranche_0 0.062***
(0.001)
regionBourgogne-Franche-Comté:tranche_0 0.072***
(0.003)
regionBretagne:tranche_0 0.051***
(0.004)
regionCentre-Val de Loire:tranche_0 0.054***
(0.004)
regionCorse:tranche_0 0.037*
(0.020)
regionDOM-TOM:tranche_0 0.114***
(0.008)
regionGrand Est:tranche_0 0.057***
(0.002)
regionHauts-de-France:tranche_0 0.053***
(0.001)
regionIle-de-France:tranche_0 0.056***
(0.001)
regionNormandie:tranche_0 0.066***
(0.003)
regionNouvelle-Aquitaine:tranche_0 0.057***
(0.002)
regionOccitanie:tranche_0 0.053***
(0.002)
regionPays de la Loire:tranche_0 0.046***
(0.003)
regionProvence-Alpes-Côte d’Azur:tranche_0 0.085***
(0.001)
Observations 2,727
R2 0.943
Adjusted R2 0.943
Residual Std. Error 27.049 (df = 2699)
F Statistic 1,608.024*** (df = 28; 2699)
Note: p<0.1; p<0.05; p<0.01

On représente comme si on on avait fait un modèle par région pour faciliter la lecture des coefficients :

Dependent variable:
hosp
Auvergne-Rhône-Alpes Hauts-de-France Provence-Alpes-Côte d’Azur Grand Est Occitanie Normandie Nouvelle-Aquitaine Centre-Val de Loire Bourgogne-Franche-Comté Bretagne Corse Pays de la Loire Ile-de-France DOM-TOM
(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) (14)
tranche_0 0.062*** 0.053*** 0.085*** 0.057*** 0.053*** 0.066*** 0.057*** 0.054*** 0.072*** 0.051*** 0.037*** 0.046*** 0.056*** 0.114***
(0.001) (0.001) (0.002) (0.001) (0.001) (0.002) (0.001) (0.002) (0.002) (0.002) (0.003) (0.002) (0.001) (0.007)
Constant 8.093*** 12.043*** 4.404 9.509*** 1.288* 2.994* 1.361* 4.142*** 3.550** 1.707 1.535** 7.253*** 22.963*** 2.091
(2.141) (3.852) (4.483) (1.382) (0.763) (1.776) (0.736) (1.286) (1.529) (1.323) (0.661) (1.602) (3.724) (2.598)
Observations 324 135 162 270 351 135 324 162 216 108 54 135 216 135
R2 0.948 0.948 0.916 0.858 0.932 0.914 0.896 0.819 0.865 0.888 0.728 0.873 0.916 0.680
Adjusted R2 0.948 0.948 0.915 0.858 0.932 0.914 0.896 0.818 0.864 0.887 0.723 0.872 0.915 0.677
Residual Std. Error 35.553 (df = 322) 41.253 (df = 133) 51.149 (df = 160) 20.762 (df = 268) 13.018 (df = 349) 18.417 (df = 133) 11.941 (df = 322) 14.506 (df = 160) 20.505 (df = 214) 11.781 (df = 106) 4.117 (df = 52) 16.322 (df = 133) 43.833 (df = 214) 23.156 (df = 133)
F Statistic 5,833.399*** (df = 1; 322) 2,432.405*** (df = 1; 133) 1,739.055*** (df = 1; 160) 1,625.486*** (df = 1; 268) 4,819.329*** (df = 1; 349) 1,416.269*** (df = 1; 133) 2,782.368*** (df = 1; 322) 723.260*** (df = 1; 160) 1,366.579*** (df = 1; 214) 838.188*** (df = 1; 106) 139.393*** (df = 1; 52) 912.887*** (df = 1; 133) 2,323.321*** (df = 1; 214) 282.039*** (df = 1; 133)
Note: p<0.1; p<0.05; p<0.01

Test:

On teste le modèle sur les données de la semaine actuelle:

L’écart quadratique moyen est égal ici à :

## [1] 1122.359

5.2 Modèle linéaire 2 (sur les région) : 1 modèle par classe d’âge

Ici, on va faire un modèle qui prend en compte les classes d’âges. Les données d’hospitalisation par classe d’âge ne sont disponibles que par région. Le modèle est de la forme

\[y_{i,t}^a=\beta_0^a+\beta_1^ax_{i,t'}^a+\epsilon_{i,t}^a\] avec:

  • \(y_{i,t}\) le nombre d’entrées de la classe d’âge \(a\) à l’hôpital dans la région \(i\) sur la période \(t\), où \(t\) est une fenêtre de 7 jours.

  • \(x_{i,t'}\) est le nombre de testés positifs de la classe d’âge \(a\) dans la région \(i\) sur la période \(t'\) où \(t'\) correspond à la fenêtre \(t\), décalé de 10 jours.

En d’autres termes, on fait ici un modèle de régression par classe d’âge, toute région confondue. Cela suppose que le lien entre les tests virologiques et le nombre d’hospitalisation est homogène dans une classe d’âge quelque soit les régions.

On merge avec le nombre de test positifs:

Apprentissage:

On modélise sur les observations des semaines précédentes et on représente les résultats tranche d’âge par tranche d’âge

Dependent variable:
hosp_9 hosp_19 hosp_29 hosp_39 hosp_49 hosp_59 hosp_69 hosp_79 hosp_89 hosp_90
(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)
tranche_9 0.021***
(0.001)
tranche_19 0.004***
(0.0002)
tranche_29 0.008***
(0.0002)
tranche_39 0.014***
(0.0003)
tranche_49 0.022***
(0.0004)
tranche_59 0.047***
(0.001)
tranche_69 0.106***
(0.001)
tranche_79 0.222***
(0.004)
tranche_89 0.349***
(0.006)
tranche_90 0.287***
(0.007)
Constant 0.985*** 1.021*** 2.698*** 3.951*** 4.583*** 5.268*** 8.208*** 10.930*** 14.916*** 9.035***
(0.274) (0.243) (0.450) (0.535) (0.673) (1.183) (1.566) (2.698) (2.988) (1.766)
Observations 378 378 378 378 378 378 378 378 378 378
R2 0.529 0.527 0.774 0.853 0.894 0.914 0.933 0.888 0.894 0.835
Adjusted R2 0.528 0.525 0.773 0.853 0.894 0.913 0.933 0.887 0.894 0.835
Residual Std. Error (df = 376) 4.725 4.244 7.828 9.496 11.976 21.049 27.885 48.119 53.591 31.830
F Statistic (df = 1; 376) 422.792*** 418.274*** 1,288.279*** 2,182.063*** 3,176.875*** 3,977.026*** 5,253.455*** 2,968.046*** 3,181.602*** 1,906.240***
Note: p<0.1; p<0.05; p<0.01

Test:

On teste le modèle sur les données de la semaine actuelle. On revient sur les données départementales, on suppose donc que les modèles estimés pour chaque tranche d’âge sur les régions est valable aussi pour les départements.

On a donc une prédiction par tranche d’âge et pour obtenir la prédiction finale, il faut donc faire la somme sur les différentes prédictions :

L’écart quadratique moyen est égal ici à :

## [1] 7439.268

5.3 Modèle de série temporelle

On utilise la même stratégie que celle présentée pour prédire le nombre de cas positifs.

Etape d’apprentissage : on entraîne l’agorithme sur les données passées en enlevant la dernière semaine observée et on prédit sur cette semaine afin de calculer les écarts quadratiques avec les valeurs observées.

On obtient le graphique suivant de valeurs prédites/valeurs observées :

L’écart quadratique moyen est égal ici à :

## [1] 3965.265 2684.802 4756.416

Les 3 prédictions sont très proches et on va choisir un algorithme de type stepwise sur les prédictions pour choisir la meilleure combinaison des modèles de séries temporelles.

## Start:  AIC=640.13
## my_basis[my_basis$semaine == "semaine_t00", "hosp"] ~ pred_3a + 
##     pred_3b + pred_3c - 1
## 
##           Df Sum of Sq   RSS    AIC
## - pred_3b  1      25.2 53861 638.18
## <none>                 53835 640.13
## - pred_3c  1    4751.7 58587 646.68
## - pred_3a  1    8415.0 62250 652.80
## 
## Step:  AIC=638.18
## my_basis[my_basis$semaine == "semaine_t00", "hosp"] ~ pred_3a + 
##     pred_3c - 1
## 
##           Df Sum of Sq   RSS    AIC
## <none>                 53861 638.18
## - pred_3c  1    6261.4 60122 647.29
## - pred_3a  1   10826.0 64687 654.68
## [1] 533.2727

On a donc 3 prédictions obtenues selon :

  • modèle par région
  • modèle par classe d’âge
  • modèle de séries temporelles (lui-même combinaison de plusieurs méthodes)

5.4 Combinaison des prédictions

Combinaison des prédictions:

  • pour prédire les nouvelles hospitalisations la semaine à venir, on va faire un panaché des trois prédictions en donnant plus de poids à la prédiction qui a le mieux marcher sur la semaine \(t_0\). Autrement dit, on fait un modèle linéaire (avec une procédure stepwise) du nombre d’hospitalisation en fonction des 3 méthodes de prédictions. On calcule l’écart quatratique moyen de la combinaison des prédictions.
## [1] 477.737
  • pour prédire les nouvelles hospitalisations la semaine d’après, on va utiliser une autre pondération en utilisant la même procédure que précédemment, mais dans une optique de prédire à deux semaines.
## Start:  AIC=1502.62
## y_true ~ pred_3a_s2 + pred_3b_s2 + pred_3c_s2 - 1
## 
## 
## Step:  AIC=1502.62
## y_true ~ pred_3a_s2 + pred_3b_s2 - 1
## 
## 
## Step:  AIC=1502.62
## y_true ~ pred_3a_s2 - 1
## 
##              Df Sum of Sq     RSS    AIC
## <none>                     340133 1502.6
## - pred_3a_s2  1   5427596 5767729 2072.4

5.5 Prédiction

On prédit le nombre d’hospitalisations :

  • du [29 novembre 2020; 05 décembre 2020] en utilisant les vrais valeurs du nombre de testé positifs la semaine du [19 novembre 2020; 25 novembre 2020].

  • du [06 décembre 2020; 12 décembre 2020] en utilisant les valeurs prédites du nombre de testé positifs la semaine du [26 novembre 2020; 02 décembre 2020].

Avant de faire cela, on actualise en incluant dans l’étape d’apprentissage les données de la dernière semaine observée:

# modèle 1
res_lm <- lm(hosp ~  tranche_0, data = my_basis[!(my_basis$semaine %in% "semaine_t0-1"), ])
# modèle 2
apprentissage_sample <- my_basis_age_wide[!(my_basis_age_wide$semaine %in% 
                                              c("semaine_t0-1")), ]
res_lm_9 <- lm(hosp_9 ~  tranche_9, data = apprentissage_sample)
res_lm_19 <- lm(hosp_19 ~  tranche_19, data = apprentissage_sample)
res_lm_29 <- lm(hosp_29 ~  tranche_29, data = apprentissage_sample)
res_lm_39 <- lm(hosp_39 ~  tranche_39, data = apprentissage_sample)
res_lm_49 <- lm(hosp_49 ~  tranche_49, data = apprentissage_sample)
res_lm_59 <- lm(hosp_59 ~  tranche_59, data = apprentissage_sample)
res_lm_69 <- lm(hosp_69 ~  tranche_69, data = apprentissage_sample)
res_lm_79 <- lm(hosp_79 ~  tranche_79, data = apprentissage_sample)
res_lm_89 <- lm(hosp_89 ~  tranche_89, data = apprentissage_sample)
res_lm_90 <- lm(hosp_90 ~  tranche_90, data = apprentissage_sample)

# On prédit avec la méthode 1 
new_data <- my_basis[my_basis$semaine %in% c("semaine_t0-1", "semaine_t0-2"), ]
pred_1 <- predict(res_lm, newdata = new_data)
# On prédit avec la méthode 2
test_sample <- my_basis[my_basis$semaine %in% c("semaine_t0-1", "semaine_t0-2"), ]
pred_9 <- predict(res_lm_9, newdata = test_sample)
pred_19 <- predict(res_lm_19, newdata = test_sample)
pred_29 <- predict(res_lm_29, newdata = test_sample)
pred_39 <- predict(res_lm_39, newdata = test_sample)
pred_49 <- predict(res_lm_49, newdata = test_sample)
pred_59 <- predict(res_lm_59, newdata = test_sample)
pred_69 <- predict(res_lm_69, newdata = test_sample)
pred_79 <- predict(res_lm_79, newdata = test_sample)
pred_89 <- predict(res_lm_89, newdata = test_sample)
pred_90 <- predict(res_lm_90, newdata = test_sample)
pred_2 <- pred_9 + pred_19 + pred_29 + pred_39 + pred_49 + pred_59 + pred_69 + 
  pred_79 + pred_89 + pred_90

# on prédit avec le modèle 3, mais on actualise les prédictions semaine par semaine
pred_3 <- matrix(0, length(nom_dep), 2)
for (k in 1:length(nom_dep)) {
 if (nom_dep[k] %in% c("975", "977", "978")) {
    pred_3[k, ] <- NA
  } else {
    temp <- hospital[!(hospital$semaine %in% c("semaine_t0-2", "semaine_t0-1")) & 
                     hospital$dep == nom_dep[k], ]
    my_ts <- zoo(temp$incid_hosp, temp$jour)
    my_ts_diff <- diff(my_ts)
    # tseries::adf.test(my_ts) 
    # tseries::adf.test(my_ts_diff)
    # predictions à 7 jours
    my_mod <- forecast::auto.arima(my_ts_diff)
    forecast_my_mod <- as.numeric(forecast(my_mod, h = 14)$mean)
    pred_3a_s1 <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:7])
    # prediction à 14 jours
    pred_3a_s2 <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:14]) - pred_3a_s1
    # Lissage exponentiel
    my_mod_exp <- ets(my_ts)
    forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp, h = 14)$mean)
    forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
    pred_3b_s1 <- round(sum(forecast_my_mod_exp[1:7]), 0)
    pred_3b_s2 <- round(sum(forecast_my_mod_exp[1:14]), 0) - pred_3b_s1
    # Méthode 3 : lissage exponentiel sur données hebdomadaires
    temp <- my_basis[my_basis$dep == nom_dep[k] & !(my_basis$semaine %in% 
                                          c("semaine_t0-2", "semaine_t0-1")), ]
    my_ts_exp <- zoo(temp$hosp, temp$jour)
    my_mod_exp_2 <- ets(my_ts_exp)
    forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp_2)$mean)
    forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
    pred_3c_s1 <- round(forecast_my_mod_exp[1], 0)
    pred_3c_s2 <- round(forecast_my_mod_exp[2], 0)
    
    # prédictions des time series
    pred_3[k, 1] <- predict(lm_3_ts, newdata = data.frame(pred_3a = pred_3a_s1,
                                                 pred_3b = pred_3b_s1,
                                                 pred_3c = pred_3c_s1))
    
    # prediction à 14 jours
    pred_3[k, 2] <- predict(lm_3b_ts, newdata = data.frame(pred_3a_s2 = pred_3a_s2,
                                                 pred_3b_s2 = pred_3b_s2,
                                                 pred_3c_s2 = pred_3c_s2))
  }
}
pred_3 <- as.vector(pred_3)
# On fait le mélande des deux prédictions
res_pred_a <- predict(lm_3, newdata = data.frame(pred_1 = pred_1[new_data$semaine == "semaine_t0-1"], 
                                               pred_2 = pred_2[test_sample$semaine == "semaine_t0-1"],
                                               pred_3 = pred_3[1:(length(pred_3) / 2)]))
res_pred_b <- predict(lm_3b, newdata = data.frame(
  pred_1_s2 = pred_1[new_data$semaine == "semaine_t0-2"],
  pred_2_s2 = pred_2[test_sample$semaine == "semaine_t0-2"],
  pred_3_s2 = pred_3[((length(pred_3) / 2) + 1):length(pred_3)]))
new_data <- my_basis[my_basis$semaine %in% "semaine_t0-1", ]
my_basis[my_basis$semaine %in% "semaine_t0-1", "hosp"] <- res_pred_a
my_basis[my_basis$semaine %in% "semaine_t0-2", "hosp"] <- res_pred_b
new_data$next_week <- res_pred_a
new_data$next_two_week <- res_pred_b

On va représenter l’évolution du nombre de nouveaux patients hospitalisés dans un intervalle de temps de 4 semaines :

  • la semaine du [15 novembre 2020; 21 novembre 2020]
  • les 7 derniers jours passés : [22 novembre 2020; 28 novembre 2020]
  • la semaine à venir : [29 novembre 2020; 05 décembre 2020]
  • la semaine suivante à venir : [06 décembre 2020; 12 décembre 2020]

On aggrège les données à la France entière:

6 Prédire le nombre de réanimations

L’idée est d’expliquer le nombre de nouvelles réanimations la semaine \(t\) par le nombre de nouvelles hospitalisations la semaine \(t-1\).

Ainsi on sera en mesure de prédire le nombre de nouvelles réanimations d’une part la semaine à venir, mais aussi la semaine d’après si on utilise les prédictions du nombre d’hospitalisation la semaine à venir.

On prépare les données et on représente le nombre de nouvelles réanimations par semaine et par département en fonction du nombre de nouvelles hospitalisations la semaine d’avant et on constate un lien très fort.

On ne va faire que deux modèles :

On n’utilise pas le modèle qui utilise les classes d’âges car c’est difficile d’avoir le nombre de nouvelles réanimations par jour/département par classe d’âge. Il se peut donc que les prédictions soient sous-estimées dans le cas où la distribution des patients hospitalisés agés évolue positivement au cours du temps.

6.1 Modèle 1 : modèle linéaire

Apprentissage:

On modélise sur les observations des semaines précédentes:

Dependent variable:
rea
regionAuvergne-Rhône-Alpes 0.094
(0.388)
regionBourgogne-Franche-Comté 0.598
(0.479)
regionBretagne 0.833
(0.715)
regionCentre-Val de Loire 0.249
(0.578)
regionCorse 0.356
(1.064)
regionDOM-TOM 2.095***
(0.687)
regionGrand Est 0.011
(0.449)
regionHauts-de-France 0.829
(0.607)
regionIle-de-France -0.208
(0.569)
regionNormandie 0.268
(0.617)
regionNouvelle-Aquitaine 0.212
(0.393)
regionOccitanie 0.292
(0.374)
regionPays de la Loire 0.526
(0.656)
regionProvence-Alpes-Côte d’Azur 0.166
(0.558)
regionAuvergne-Rhône-Alpes:hosp 0.153***
(0.002)
regionBourgogne-Franche-Comté:hosp 0.130***
(0.008)
regionBretagne:hosp 0.133***
(0.018)
regionCentre-Val de Loire:hosp 0.177***
(0.016)
regionCorse:hosp 0.150
(0.114)
regionDOM-TOM:hosp 0.093***
(0.013)
regionGrand Est:hosp 0.155***
(0.008)
regionHauts-de-France:hosp 0.187***
(0.003)
regionIle-de-France:hosp 0.190***
(0.003)
regionNormandie:hosp 0.147***
(0.009)
regionNouvelle-Aquitaine:hosp 0.154***
(0.010)
regionOccitanie:hosp 0.209***
(0.007)
regionPays de la Loire:hosp 0.151***
(0.013)
regionProvence-Alpes-Côte d’Azur:hosp 0.155***
(0.003)
Observations 2,626
R2 0.897
Adjusted R2 0.896
Residual Std. Error 6.303 (df = 2598)
F Statistic 807.884*** (df = 28; 2598)
Note: p<0.1; p<0.05; p<0.01

Test:

On teste le modèle sur les données de la semaine actuelle:

L’écart quadratique moyen est égal ici à :

## [1] 184.7122

6.2 Modèle 2 : série temporelle

On utilise la même stratégie que celle présentée pour prédire le nombre de nouveaux cas positifs et de nouvelles réanimations.

Etape d’apprentissage : on entraîne l’agorithme sur les données passées en enlevant la dernière semaine observée et on prédit sur cette semaine afin de calculer les écarts quadratiques avec les valeurs observées.

On observe le graphique des valeurs prédites/valeurs observées :

L’écart quadratique moyen est égal ici à :

## [1] 111.13572  97.66337 172.35644

Les 3 prédictions sont très proches et on va choisir un algorithme de type stepwise sur les prédictions pour choisir la meilleure combinaison et ne conserver qu’une seule prédiction basée sur les séries temporelles:

## Start:  AIC=365.3
## my_basis[my_basis$semaine == "semaine_t00", "rea"] ~ pred_rea_2a + 
##     pred_rea_2b + pred_rea_2c - 1
## 
##               Df Sum of Sq    RSS    AIC
## - pred_rea_2a  1      5.70 3547.9 363.46
## - pred_rea_2b  1     11.08 3553.3 363.61
## <none>                     3542.2 365.30
## - pred_rea_2c  1   1425.92 4968.1 397.46
## 
## Step:  AIC=363.46
## my_basis[my_basis$semaine == "semaine_t00", "rea"] ~ pred_rea_2b + 
##     pred_rea_2c - 1
## 
##               Df Sum of Sq    RSS    AIC
## - pred_rea_2b  1     57.75 3605.6 363.09
## <none>                     3547.9 363.46
## - pred_rea_2c  1   1461.00 5008.9 396.29
## 
## Step:  AIC=363.09
## my_basis[my_basis$semaine == "semaine_t00", "rea"] ~ pred_rea_2c - 
##     1
## 
##               Df Sum of Sq   RSS    AIC
## <none>                      3606 363.09
## - pred_rea_2c  1     48983 52589 631.77
## [1] 35.69947

6.3 Combinaison des prédictions

Combinaison des prédictions: on peut envisager de faire un panaché des deux prédictions. Autrement dit, on fait un modèle linéaire (avec une procédure stepwise) du nombre de réanimations observée la semaine t0 en fonction des 2 méthodes de prédictions. On obtient l’écart-quadratique moyen suivant:

## [1] 33.88473

On adapte le poids des prédictions en fonction de la semaine à prédire

semaine_to_drop <- c("semaine_t0-3", "semaine_t0-2", "semaine_t0-1", "semaine_t00")
pred_rea_1_s2 <- numeric(0)
pred_rea_2a_s2 <- numeric(0)
pred_rea_2b_s2 <- numeric(0)
pred_rea_2c_s2 <- numeric(0)
pred_rea_1_s3 <- numeric(0)
pred_rea_2a_s3 <- numeric(0)
pred_rea_2b_s3 <- numeric(0)
pred_rea_2c_s3 <- numeric(0)

y_true <- numeric(0)

for (j in 0:1) {
  
  semaine_to_estim <-  paste0("semaine_t0", j)
  y_true <- c(y_true, my_basis[my_basis$semaine == semaine_to_estim, "rea"])
  semaine_to_drop <- c(semaine_to_drop, paste0("semaine_t0", j + 1))
  res_lm_rea_1 <- lm(rea ~  region + hosp:region - 1, 
             data = my_basis_rea[!(my_basis_rea$semaine %in% semaine_to_drop), ])
  pred_rea_1_s2 <- c(pred_rea_1_s2, 
     round(predict(res_lm_rea_1, newdata = my_basis_rea[my_basis_rea$semaine == semaine_to_estim, ])))
  pred_rea_2a_temp <- numeric(length(nom_dep))
  pred_rea_2b_temp <- numeric(length(nom_dep))
  pred_rea_2c_temp <- numeric(length(nom_dep))

  for (k in 1:length(nom_dep)) {
    temp <- hospital[!(hospital$semaine %in% semaine_to_drop) & 
                     hospital$dep == nom_dep[k], ]
    my_ts <- zoo(temp$incid_rea, temp$jour)
    my_ts_diff <- diff(my_ts)
    # tseries::adf.test(my_ts) 
    # tseries::adf.test(my_ts_diff)
    if (nom_dep[k] %in% c("975", "977", "978")) {
      pred_rea_2[k] <- NA
    } else {
      my_mod <- forecast::auto.arima(my_ts_diff)
      forecast_my_mod <- as.numeric(forecast(my_mod, h = 14)$mean)
      temp <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:7])
      pred_rea_2a_temp[k] <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:14]) - temp
      # modèles exponentiels
      # Méthode 2 : lissage exponentiel
      my_mod_exp <- ets(my_ts)
      forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp, h = 14)$mean)
      forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
      temp <- round(sum(forecast_my_mod_exp[1:7]), 0)
      pred_rea_2b_temp[k] <- round(sum(forecast_my_mod_exp[1:14]), 0) - temp
      # Méthode 3 : lissage exponentiel sur données hebdomadaires
      temp <- my_basis[my_basis$dep == nom_dep[k] & !(my_basis$semaine %in% semaine_to_drop), ]
      my_ts_exp <- zoo(temp$rea, temp$jour)
      if (all(my_ts_exp == 0)) {
        pred_rea_2c_temp[k] <- 0
      } else {
        my_mod_exp_2 <- ets(my_ts_exp)
        forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp_2)$mean)
        forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
        pred_rea_2c_temp[k] <- round(forecast_my_mod_exp[2], 0)
      }
    }
    
  }
   pred_rea_2a_s2 <- c(pred_rea_2a_s2, pred_rea_2a_temp)
   pred_rea_2b_s2 <- c(pred_rea_2b_s2, pred_rea_2b_temp)
   pred_rea_2c_s2 <- c(pred_rea_2c_s2, pred_rea_2c_temp)
}


lm_2_rea_ts_s1 <- step(lm(y_true ~ pred_rea_2a_s2 + pred_rea_2b_s2 + pred_rea_2c_s2 - 1))
## Start:  AIC=827.63
## y_true ~ pred_rea_2a_s2 + pred_rea_2b_s2 + pred_rea_2c_s2 - 1
## 
##                  Df Sum of Sq   RSS    AIC
## - pred_rea_2a_s2  1      0.32 11799 825.64
## <none>                        11799 827.63
## - pred_rea_2c_s2  1    195.82 11995 828.96
## - pred_rea_2b_s2  1    936.35 12735 841.06
## 
## Step:  AIC=825.64
## y_true ~ pred_rea_2b_s2 + pred_rea_2c_s2 - 1
## 
##                  Df Sum of Sq   RSS    AIC
## <none>                        11799 825.64
## - pred_rea_2c_s2  1     195.6 11995 826.96
## - pred_rea_2b_s2  1    3090.1 14889 870.63
pred_rea_2_s2 <- predict(lm_2_rea_ts_s1)

lm_rea_3b <- lm(y_true ~ pred_rea_1_s2 + pred_rea_2_s2 - 1)



#######

semaine_to_drop <- c("semaine_t0-3", "semaine_t0-2", "semaine_t0-1", "semaine_t00", "semaine_t01")
pred_rea_1_s3 <- numeric(0)
pred_rea_2a_s3 <- numeric(0)
pred_rea_2b_s3 <- numeric(0)
pred_rea_2c_s3 <- numeric(0)

y_true <- numeric(0)

for (j in 0:1) {
  
  semaine_to_estim <-  paste0("semaine_t0", j)
  y_true <- c(y_true, my_basis[my_basis$semaine == semaine_to_estim, "rea"])
  semaine_to_drop <- c(semaine_to_drop, paste0("semaine_t0", j + 2))
  
  res_lm_rea_1 <- lm(rea ~  region + hosp:region - 1, 
             data = my_basis_rea[!(my_basis_rea$semaine %in% semaine_to_drop), ])
  pred_rea_1_s3 <- c(pred_rea_1_s3, round(predict(res_lm_rea_1, newdata = my_basis_rea[my_basis_rea$semaine == semaine_to_estim, ])))
  
  pred_rea_2a_temp <- numeric(length(nom_dep))
  pred_rea_2b_temp <- numeric(length(nom_dep))
  pred_rea_2c_temp <- numeric(length(nom_dep))

  for (k in 1:length(nom_dep)) {
     temp <- hospital[!(hospital$semaine %in% semaine_to_drop) & hospital$dep == nom_dep[k], ]
     my_ts <- zoo(temp$incid_rea, temp$jour)
     my_ts_diff <- diff(my_ts)
     # tseries::adf.test(my_ts) 
     # tseries::adf.test(my_ts_diff)
     if (nom_dep[k] %in% c("975", "977", "978")) {
       pred_rea_2[k] <- NA
     } else {
       my_mod <- forecast::auto.arima(my_ts_diff)
       forecast_my_mod <- as.numeric(forecast(my_mod, h = 21)$mean)
       temp1 <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:7])
       temp2 <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:14]) - temp1
       pred_rea_2a_temp[k] <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:21]) - temp1 - temp2
       # modèles exponentiels
       # Méthode 2 : lissage exponentiel
       my_mod_exp <- ets(my_ts)
       forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp, h = 21)$mean)
       forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
       temp1 <- round(sum(forecast_my_mod_exp[1:7]), 0)
       temp2 <- round(sum(forecast_my_mod_exp[1:14]), 0) - temp1
       pred_rea_2b_temp[k] <- round(sum(forecast_my_mod_exp[1:21]), 0) - temp2 - temp1
       # Méthode 3 : lissage exponentiel sur données hebdomadaires
       temp <- my_basis[my_basis$dep == nom_dep[k] & !(my_basis$semaine %in% semaine_to_drop), ]
       my_ts_exp <- zoo(temp$rea, temp$jour)
       if (all(my_ts_exp == 0)) {
         pred_rea_2c_temp[k] <- 0
       } else {
         my_mod_exp_2 <- ets(my_ts_exp)
         forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp_2)$mean)
         forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
         pred_rea_2c_temp[k] <- round(forecast_my_mod_exp[3], 0)
       }
    }
}

  pred_rea_2a_s3 <- c(pred_rea_2a_s3, pred_rea_2a_temp)
  pred_rea_2b_s3 <- c(pred_rea_2b_s3, pred_rea_2b_temp)
  pred_rea_2c_s3 <- c(pred_rea_2c_s3, pred_rea_2c_temp)
}

lm_2_rea_ts_s2 <- step(lm(y_true ~ pred_rea_2a_s3 + pred_rea_2b_s3 + pred_rea_2c_s3 - 1))
## Start:  AIC=961.49
## y_true ~ pred_rea_2a_s3 + pred_rea_2b_s3 + pred_rea_2c_s3 - 1
## 
##                  Df Sum of Sq   RSS    AIC
## - pred_rea_2c_s3  1     65.80 22954 960.07
## - pred_rea_2a_s3  1    225.91 23114 961.47
## <none>                        22889 961.49
## - pred_rea_2b_s3  1   1622.33 24511 973.32
## 
## Step:  AIC=960.07
## y_true ~ pred_rea_2a_s3 + pred_rea_2b_s3 - 1
## 
##                  Df Sum of Sq   RSS    AIC
## <none>                        22954 960.07
## - pred_rea_2a_s3  1    337.94 23292 961.02
## - pred_rea_2b_s3  1   1872.76 24827 973.91

6.4 Prédiction

On prédit:

  • le nombre de réanimations à venir du [29 novembre 2020; 05 décembre 2020] en utilisant les nouvelles hospitalisations du [22 novembre 2020; 28 novembre 2020]

  • le nombre de réanimations à venir du [06 décembre 2020; 12 décembre 2020] en utilisant la prédiction des hospitalisations à venir du [29 novembre 2020; 05 décembre 2020]

  • le nombre de réanimations à venir du [13 décembre 2020; 19 décembre 2020] en utilisant la prédiction des hospitalisations à venir du [06 décembre 2020; 12 décembre 2020]

Pour cela, on actualise le modèle, c’est-à-dire qu’on inclut la dernière semaine observée:

res_lm <- lm(rea ~  region + hosp:region - 1, 
             data = my_basis_rea[!(my_basis_rea$semaine %in% c("semaine_t0-2", "semaine_t0-1")), ])

# semaine t+1
new_data_rea_1 <- my_basis_rea[my_basis_rea$semaine %in% c("semaine_t0-1", "semaine_t0-2", "semaine_t0-3"),  ]
pred_rea_1 <- predict(res_lm, newdata = new_data_rea_1)

pred_rea_2 <- matrix(0, length(nom_dep), 3)
pred_rea_2a <- matrix(0, length(nom_dep), 3)
pred_rea_2b <- matrix(0, length(nom_dep), 3)
pred_rea_2c <- matrix(0, length(nom_dep), 3)
for (k in 1:length(nom_dep)) {
  temp <- hospital[!(hospital$semaine %in% c("semaine_t0-3", "semaine_t0-2", "semaine_t0-1")) & 
                     hospital$dep == nom_dep[k], ]
  my_ts <- zoo(temp$incid_rea, temp$jour)
  my_ts_diff <- diff(my_ts)
  # tseries::adf.test(my_ts) 
  # tseries::adf.test(my_ts_diff)
  if (nom_dep[k] %in% c("975", "977", "978")) {
    pred_rea_2[k] <- NA
  } else {
    my_mod <- forecast::auto.arima(my_ts_diff)
    forecast_my_mod <- as.numeric(forecast(my_mod, h = 21)$mean)
    pred_3a_s1 <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:7])
    pred_3a_s2 <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:14]) - pred_3a_s1
    pred_3a_s3 <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:21]) - 
      pred_3a_s1 - pred_3a_s2
    # Lissage exponentiel
    my_mod_exp <- ets(my_ts)
    forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp, h = 21)$mean)
    forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
    pred_3b_s1 <- round(sum(forecast_my_mod_exp[1:7]), 0)
    pred_3b_s2 <- round(sum(forecast_my_mod_exp[1:14]), 0) - pred_3b_s1
    pred_3b_s3 <- round(sum(forecast_my_mod_exp[1:21]), 0) - pred_3b_s2 - pred_3b_s1
    # Méthode 3 : lissage exponentiel sur données hebdomadaires
    temp <- my_basis[my_basis$dep == nom_dep[k] & !(my_basis$semaine %in% 
                                          c("semaine_t0-3", "semaine_t0-2", "semaine_t0-1")), ]
    my_ts_exp <- zoo(temp$rea, temp$jour)
    my_mod_exp_2 <- ets(my_ts_exp)
    forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp_2)$mean)
    forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
    pred_3c_s1 <- round(forecast_my_mod_exp[1], 0)
    pred_3c_s2 <- round(forecast_my_mod_exp[2], 0)
    pred_3c_s3 <- round(forecast_my_mod_exp[3], 0)    
    
    pred_rea_2[k, 1] <- predict(lm_2_rea_ts, newdata = data.frame(pred_rea_2a = pred_3a_s1,
                                                 pred_rea_2b = pred_3b_s1,
                                                 pred_rea_2c = pred_3c_s1))
    pred_rea_2[k, 2] <- predict(lm_2_rea_ts_s1, newdata = data.frame(pred_rea_2a_s2 = pred_3a_s2,
                                                 pred_rea_2b_s2 = pred_3b_s2,
                                                 pred_rea_2c_s2 = pred_3c_s2))
    pred_rea_2[k, 3] <- predict(lm_2_rea_ts_s2, newdata = data.frame(pred_rea_2a_s3 = pred_3a_s3,
                                                 pred_rea_2b_s3 = pred_3b_s3,
                                                 pred_rea_2c_s3 = pred_3c_s3))
  }
}

pred_rea_a <- predict(lm_rea_3, newdata = data.frame(pred_rea_1 = 
                        pred_rea_1[new_data_rea_1$semaine == "semaine_t0-1"],
                        pred_rea_2 = as.vector(pred_rea_2)[1:(length(pred_rea_1) / 3)]))

pred_rea_b <- predict(lm_rea_3b, newdata = data.frame(
  pred_rea_1_s2 = pred_rea_1[new_data_rea_1$semaine == "semaine_t0-2"],
  pred_rea_2_s2 = as.vector(pred_rea_2)[((length(pred_rea_1) / 3) + 1):(2 * length(pred_rea_1) / 3)]))

pred_rea_c <- predict(lm_rea_3c, newdata = data.frame(
  pred_rea_1_s3 = pred_rea_1[new_data_rea_1$semaine == "semaine_t0-3"],
  pred_rea_2_s3 = as.vector(pred_rea_2)[(2 * length(pred_rea_1) / 3 + 1):length(pred_rea_1)]))

# on synthétise les résultats
new_data <- my_basis[my_basis$semaine %in% "semaine_t0-1", ]
new_data$this_week <- my_basis_rea[my_basis_rea$semaine == "semaine_t00", "rea"]
new_data$next_week <- pred_rea_a
new_data$next_two_week <- pred_rea_b
new_data$next_three_week <- pred_rea_c
my_basis[my_basis$semaine %in% "semaine_t0-1", "rea"] <- pred_rea_a
my_basis[my_basis$semaine %in% "semaine_t0-2", "rea"] <- pred_rea_b
my_basis[my_basis$semaine %in% "semaine_t0-3", "rea"] <- pred_rea_c

On va représenter l’évolution du nombre de patients en réanimations dans un intervalle de temps de 4 semaines :

  • les 7 derniers jours passés : [22 novembre 2020; 28 novembre 2020]
  • la semaine à venir : [29 novembre 2020; 05 décembre 2020]
  • la 2ème semaine à venir : [06 décembre 2020; 12 décembre 2020]
  • la 3ème semaine à venir : [13 décembre 2020; 19 décembre 2020]

On aggrège les données à la France entière:

7 Prédire le nombre de décès

L’idée est d’expliquer le nombre de nouveaux décès la semaine \(t\) par les nouvelles réanimations la semaine \(t-1\).

Ainsi on sera en mesure de prédire le nombre de nouveaux décès la semaine à venir, mais aussi les trois semaines suivantes en utilisant les prédictions des hospitalisations, des réanimations et cas positifs.

On prépare les données et on représente le nombre de nouveaux décès par semaine et par département en fonction du nombre de nouvelles réanimations la semaine d’avant et on constate un lien très fort.

On va faire deux modèles : un modèle régional où on explique les nouvelles réanimations des départements au sein d’une même région ainsi qu’un modèle de série temporelle département par département.

7.1 Modèle 1 : modèle linéaire

Apprentissage:

On modélise sur les observations des semaines précédentes:

Dependent variable:
dc
regionAuvergne-Rhône-Alpes 1.603***
(0.354)
regionBourgogne-Franche-Comté 0.937**
(0.436)
regionBretagne 0.226
(0.660)
regionCentre-Val de Loire 1.296**
(0.516)
regionCorse 0.267
(0.950)
regionDOM-TOM -0.495
(0.676)
regionGrand Est 1.206***
(0.407)
regionHauts-de-France 1.044*
(0.560)
regionIle-de-France 5.023***
(0.498)
regionNormandie -0.062
(0.565)
regionNouvelle-Aquitaine 0.762**
(0.359)
regionOccitanie 0.606*
(0.342)
regionPays de la Loire -0.229
(0.594)
regionProvence-Alpes-Côte d’Azur 0.853*
(0.510)
regionAuvergne-Rhône-Alpes:rea 0.920***
(0.014)
regionBourgogne-Franche-Comté:rea 0.875***
(0.054)
regionBretagne:rea 0.677***
(0.100)
regionCentre-Val de Loire:rea 0.467***
(0.073)
regionCorse:rea 0.738
(0.486)
regionDOM-TOM:rea 0.586***
(0.093)
regionGrand Est:rea 0.843***
(0.047)
regionHauts-de-France:rea 0.807***
(0.017)
regionIle-de-France:rea 0.561***
(0.013)
regionNormandie:rea 1.145***
(0.057)
regionNouvelle-Aquitaine:rea 0.633***
(0.059)
regionOccitanie:rea 0.566***
(0.029)
regionPays de la Loire:rea 0.884***
(0.071)
regionProvence-Alpes-Côte d’Azur:rea 0.829***
(0.017)
Observations 2,626
R2 0.864
Adjusted R2 0.862
Residual Std. Error 5.850 (df = 2598)
F Statistic 588.602*** (df = 28; 2598)
Note: p<0.1; p<0.05; p<0.01

Test:

On teste le modèle sur les données de la semaine actuelle:

L’écart quadratique moyen est égal ici à :

## [1] 145.9601

7.2 Modèle 2 : série temporelle

On utilise la même stratégie que celle présentée pour prédire le nombre de nouveaux cas, de nouvelles hospitalisations et de nouvelles réanimations.

Etape d’apprentissage : on entraîne l’agorithme sur les données passées en enlevant la dernière semaine observée et on prédit sur cette semaine afin de calculer les écarts quadratiques avec les valeurs observées.

On représente le graphique des valeurs prédites / valeurs observées :

L’écart quadratique moyen est égal ici à :

## [1] 102.4050 107.6931 119.4059

Les 3 prédictions sont très proches et on va choisir un algorithme de type stepwise sur les prédictions pour choisir la meilleure combinaison et ne garder qu’une prédiction de type série temporelle:

## Start:  AIC=434.96
## my_basis[my_basis$semaine == "semaine_t00", "dc"] ~ pred_dc_2a + 
##     pred_dc_2b + pred_dc_2c - 1
## 
##              Df Sum of Sq    RSS    AIC
## - pred_dc_2a  1     83.56 7144.1 434.15
## <none>                    7060.5 434.96
## - pred_dc_2b  1    201.01 7261.6 435.80
## - pred_dc_2c  1   1213.03 8273.6 448.98
## 
## Step:  AIC=434.15
## my_basis[my_basis$semaine == "semaine_t00", "dc"] ~ pred_dc_2b + 
##     pred_dc_2c - 1
## 
##              Df Sum of Sq    RSS    AIC
## - pred_dc_2b  1    125.07 7269.2 433.90
## <none>                    7144.1 434.15
## - pred_dc_2c  1   1178.87 8323.0 447.58
## 
## Step:  AIC=433.9
## my_basis[my_basis$semaine == "semaine_t00", "dc"] ~ pred_dc_2c - 
##     1
## 
##              Df Sum of Sq    RSS    AIC
## <none>                      7269 433.90
## - pred_dc_2c  1    113437 120706 715.69
## [1] 71.97194

7.3 Combinaison des prédictions

Combinaison des prédictions: on peut envisager de faire un panaché des deux prédictions en régressant (avec un algorithme de type stepwise) le nombre de décés observé la semaine t0 en fonction des deux méthodes de régression. On obtient l’écart moyen quadratique suivant :

## [1] 61.31425

On adapte le poids des prédictions en fonction de la semaine à prédire

semaine_to_drop <- c("semaine_t0-4", "semaine_t0-3", "semaine_t0-2", "semaine_t0-1", "semaine_t00")
pred_dc_1_s2 <- numeric(0)
pred_dc_2a_s2 <- numeric(0)
pred_dc_2b_s2 <- numeric(0)
pred_dc_2c_s2 <- numeric(0)

y_true <- numeric(0)

for (j in 0:2) {
  
  semaine_to_estim <-  paste0("semaine_t0", j)
  y_true <- c(y_true, my_basis[my_basis$semaine == semaine_to_estim, "dc"])
  semaine_to_drop <- c(semaine_to_drop, paste0("semaine_t0", j + 1))
  res_lm_dc_1 <- lm(dc ~  region + rea:region - 1, 
             data = my_basis_dc[!(my_basis_dc$semaine %in% semaine_to_drop), ])
  pred_dc_1_s2 <- c(pred_dc_1_s2, 
     round(predict(res_lm_dc_1, newdata = my_basis_dc[my_basis_dc$semaine == semaine_to_estim, ])))
  pred_dc_2a_temp <- numeric(length(nom_dep))
  pred_dc_2b_temp <- numeric(length(nom_dep))
  pred_dc_2c_temp <- numeric(length(nom_dep))

  for (k in 1:length(nom_dep)) {
    temp <- hospital[!(hospital$semaine %in% semaine_to_drop) & 
                     hospital$dep == nom_dep[k], ]
    my_ts <- zoo(temp$incid_dc, temp$jour)
    my_ts_diff <- diff(my_ts)
    # tseries::adf.test(my_ts) 
    # tseries::adf.test(my_ts_diff)
    if (nom_dep[k] %in% c("975", "977", "978")) {
      pred_dc_2[k] <- NA
    } else {
      my_mod <- forecast::auto.arima(my_ts_diff)
      forecast_my_mod <- as.numeric(forecast(my_mod, h = 14)$mean)
      temp <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:7])
      pred_dc_2a_temp[k] <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:14]) - temp
      # modèles exponentiels
      # Méthode 2 : lissage exponentiel
      my_mod_exp <- ets(my_ts)
      forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp, h = 14)$mean)
      forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
      temp <- round(sum(forecast_my_mod_exp[1:7]), 0)
      pred_dc_2b_temp[k] <- round(sum(forecast_my_mod_exp[1:14]), 0) - temp
      # Méthode 3 : lissage exponentiel sur données hebdomadaires
      temp <- my_basis[my_basis$dep == nom_dep[k] & !(my_basis$semaine %in% semaine_to_drop), ]
      my_ts_exp <- zoo(temp$dc, temp$jour)
      if (all(my_ts_exp == 0)) {
        pred_dc_2c_temp[k] <- 0
      } else {
        my_mod_exp_2 <- ets(my_ts_exp)
        forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp_2)$mean)
        forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
        pred_dc_2c_temp[k] <- round(forecast_my_mod_exp[2], 0)
      }
    }
    
  }
   pred_dc_2a_s2 <- c(pred_dc_2a_s2, pred_dc_2a_temp)
   pred_dc_2b_s2 <- c(pred_dc_2b_s2, pred_dc_2b_temp)
   pred_dc_2c_s2 <- c(pred_dc_2c_s2, pred_dc_2c_temp)
}


lm_2_dc_ts_s1 <- step(lm(y_true ~ pred_dc_2a_s2 + pred_dc_2b_s2 + pred_dc_2c_s2 - 1))
## Start:  AIC=1623.86
## y_true ~ pred_dc_2a_s2 + pred_dc_2b_s2 + pred_dc_2c_s2 - 1
## 
##                 Df Sum of Sq   RSS    AIC
## <none>                       63147 1623.9
## - pred_dc_2b_s2  1     739.5 63886 1625.4
## - pred_dc_2c_s2  1    1879.1 65026 1630.8
## - pred_dc_2a_s2  1    3868.7 67015 1639.9
pred_dc_2_s2 <- predict(lm_2_dc_ts_s1)

lm_dc_3b <- lm(y_true ~ pred_dc_1_s2 + pred_dc_2_s2 - 1)

#######
# Semaine + 2
#######

semaine_to_drop <- c("semaine_t0-4", "semaine_t0-3", "semaine_t0-2", "semaine_t0-1", "semaine_t00", "semaine_t01")
pred_dc_1_s3 <- pred_dc_1_s2
pred_dc_2a_s3 <- pred_dc_2a_s2
pred_dc_2b_s3 <- pred_dc_2b_s2
pred_dc_2c_s3 <- pred_dc_2c_s2


for (j in 0:1) {
  
  semaine_to_estim <-  paste0("semaine_t0", j)
  y_true <- c(y_true, my_basis[my_basis$semaine == semaine_to_estim, "dc"])
  semaine_to_drop <- c(semaine_to_drop, paste0("semaine_t0", j + 2))
  
  res_lm_dc_1 <- lm(dc ~  region + rea:region - 1, 
             data = my_basis_dc[!(my_basis_dc$semaine %in% semaine_to_drop), ])
  pred_dc_1_s3 <- c(pred_dc_1_s3, round(predict(res_lm_dc_1, 
                      newdata = my_basis_dc[my_basis_dc$semaine == semaine_to_estim, ])))
  
  pred_dc_2a_temp <- numeric(length(nom_dep))
  pred_dc_2b_temp <- numeric(length(nom_dep))
  pred_dc_2c_temp <- numeric(length(nom_dep))

  for (k in 1:length(nom_dep)) {
     temp <- hospital[!(hospital$semaine %in% semaine_to_drop) & hospital$dep == nom_dep[k], ]
     my_ts <- zoo(temp$incid_dc, temp$jour)
     my_ts_diff <- diff(my_ts)
     # tseries::adf.test(my_ts) 
     # tseries::adf.test(my_ts_diff)
     if (nom_dep[k] %in% c("975", "977", "978")) {
       pred_dc_2[k] <- NA
     } else {
       my_mod <- forecast::auto.arima(my_ts_diff)
       forecast_my_mod <- as.numeric(forecast(my_mod, h = 21)$mean)
       temp1 <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:7])
       temp2 <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:14]) - temp1
       temp3 <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:21]) - temp1 - temp2
       pred_dc_2a_temp[k] <- ifelse(temp3 > 0, round(temp3), 0)
       # modèles exponentiels
       # Méthode 2 : lissage exponentiel
       my_mod_exp <- ets(my_ts)
       forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp, h = 21)$mean)
       forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
       temp1 <- round(sum(forecast_my_mod_exp[1:7]), 0)
       temp2 <- round(sum(forecast_my_mod_exp[1:14]), 0) - temp1
       pred_dc_2b_temp[k] <- round(sum(forecast_my_mod_exp[1:21]), 0) - temp2 - temp1
       # Méthode 3 : lissage exponentiel sur données hebdomadaires
       temp <- my_basis[my_basis$dep == nom_dep[k] & !(my_basis$semaine %in% semaine_to_drop), ]
       my_ts_exp <- zoo(temp$dc, temp$jour)
       if (all(my_ts_exp == 0)) {
         pred_dc_2c_temp[k] <- 0
       } else {
         my_mod_exp_2 <- ets(my_ts_exp)
         forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp_2)$mean)
         forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
         pred_dc_2c_temp[k] <- round(forecast_my_mod_exp[3], 0)
       }
    }
}

  pred_dc_2a_s3 <- c(pred_dc_2a_s3, pred_dc_2a_temp)
  pred_dc_2b_s3 <- c(pred_dc_2b_s3, pred_dc_2b_temp)
  pred_dc_2c_s3 <- c(pred_dc_2c_s3, pred_dc_2c_temp)
}

lm_2_dc_ts_s2 <- step(lm(y_true ~ pred_dc_2a_s3 + pred_dc_2b_s3 + pred_dc_2c_s3 - 1))
## Start:  AIC=2755.51
## y_true ~ pred_dc_2a_s3 + pred_dc_2b_s3 + pred_dc_2c_s3 - 1
## 
##                 Df Sum of Sq    RSS    AIC
## <none>                       116907 2755.5
## - pred_dc_2b_s3  1     770.0 117677 2756.8
## - pred_dc_2c_s3  1    2148.5 119056 2762.7
## - pred_dc_2a_s3  1   14934.6 131842 2814.2
pred_dc_2_s3 <- predict(lm_2_dc_ts_s2)

lm_dc_3c <- lm(y_true ~ pred_dc_1_s3 + pred_dc_2_s3 - 1)

#######
# Semaine T + 3
#######

semaine_to_drop <- c("semaine_t0-4", "semaine_t0-3", "semaine_t0-2", "semaine_t0-1", 
                     "semaine_t00", "semaine_t01", "semaine_t02")
pred_dc_1_s4 <- pred_dc_1_s3
pred_dc_2a_s4 <- pred_dc_2a_s3
pred_dc_2b_s4 <- pred_dc_2b_s3
pred_dc_2c_s4 <- pred_dc_2c_s3


for (j in 0:1) {
  
  semaine_to_estim <-  paste0("semaine_t0", j)
  y_true <- c(y_true, my_basis[my_basis$semaine == semaine_to_estim, "dc"])
  semaine_to_drop <- c(semaine_to_drop, paste0("semaine_t0", j + 3))
  
  res_lm_dc_1 <- lm(dc ~  region + rea:region - 1, 
             data = my_basis_dc[!(my_basis_dc$semaine %in% semaine_to_drop), ])
  pred_dc_1_s4 <- c(pred_dc_1_s4, round(predict(res_lm_dc_1, 
                      newdata = my_basis_dc[my_basis_dc$semaine == semaine_to_estim, ])))
  
  pred_dc_2a_temp <- numeric(length(nom_dep))
  pred_dc_2b_temp <- numeric(length(nom_dep))
  pred_dc_2c_temp <- numeric(length(nom_dep))

  for (k in 1:length(nom_dep)) {
     temp <- hospital[!(hospital$semaine %in% semaine_to_drop) & hospital$dep == nom_dep[k], ]
     my_ts <- zoo(temp$incid_dc, temp$jour)
     my_ts_diff <- diff(my_ts)
     # tseries::adf.test(my_ts) 
     # tseries::adf.test(my_ts_diff)
     if (nom_dep[k] %in% c("975", "977", "978")) {
       pred_dc_2[k] <- NA
     } else {
       my_mod <- forecast::auto.arima(my_ts_diff)
       forecast_my_mod <- as.numeric(forecast(my_mod, h = 28)$mean)
       temp1 <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:7])
       temp2 <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:14]) - temp1
       temp3 <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:21]) - temp1 - temp2 
       temp4 <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:28]) - temp1 - temp2 - temp3
       pred_dc_2a_temp[k] <- ifelse(temp4 > 0, round(temp4), 0)
       # modèles exponentiels
       # Méthode 2 : lissage exponentiel
       my_mod_exp <- ets(my_ts)
       forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp, h = 28)$mean)
       forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
       temp1 <- round(sum(forecast_my_mod_exp[1:7]), 0)
       temp2 <- round(sum(forecast_my_mod_exp[1:14]), 0) - temp1
       temp3 <- round(sum(forecast_my_mod_exp[1:21]), 0) - temp1 - temp2
       pred_dc_2b_temp[k] <- round(sum(forecast_my_mod_exp[1:28]), 0) - temp1 - temp2 - temp3
       # Méthode 3 : lissage exponentiel sur données hebdomadaires
       temp <- my_basis[my_basis$dep == nom_dep[k] & !(my_basis$semaine %in% semaine_to_drop), ]
       my_ts_exp <- zoo(temp$dc, temp$jour)
       if (all(my_ts_exp == 0)) {
         pred_dc_2c_temp[k] <- 0
       } else {
         my_mod_exp_2 <- ets(my_ts_exp)
         forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp_2)$mean)
         forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
         pred_dc_2c_temp[k] <- round(forecast_my_mod_exp[4], 0)
       }
    }
}

  pred_dc_2a_s4 <- c(pred_dc_2a_s4, pred_dc_2a_temp)
  pred_dc_2b_s4 <- c(pred_dc_2b_s4, pred_dc_2b_temp)
  pred_dc_2c_s4 <- c(pred_dc_2c_s4, pred_dc_2c_temp)
}

lm_2_dc_ts_s3 <- step(lm(y_true ~ pred_dc_2a_s4 + pred_dc_2b_s4 + pred_dc_2c_s4 - 1))
## Start:  AIC=4009.95
## y_true ~ pred_dc_2a_s4 + pred_dc_2b_s4 + pred_dc_2c_s4 - 1
## 
##                 Df Sum of Sq    RSS    AIC
## <none>                       203685 4010.0
## - pred_dc_2b_s4  1     958.1 204643 4011.3
## - pred_dc_2c_s4  1    1483.1 205168 4013.1
## - pred_dc_2a_s4  1   21959.8 225644 4080.3

7.4 Prédiction

On prédit:

  • le nombre de décès à venir du [29 novembre 2020; 05 décembre 2020] en utilisant les nouvelles réanimations du [22 novembre 2020; 28 novembre 2020]

  • le nombre de décès à venir du [06 décembre 2020; 12 décembre 2020] en utilisant la prédiction des réanimations à venir du [29 novembre 2020; 05 décembre 2020]

  • le nombre de décès à venir du [13 décembre 2020; 19 décembre 2020] en utilisant la prédiction des réanimations à venir du [06 décembre 2020; 12 décembre 2020]

  • le nombre de décès à venir du [20 décembre 2020; 26 décembre 2020] en utilisant la prédiction des réanimations à venir du [06 décembre 2020; 12 décembre 2020]

Pour cela, on actualise le modèle, c’est-à-dire qu’on inclut la dernière semaine observée:

res_lm <- lm(dc ~  region + rea:region - 1, 
             data = my_basis_dc[!(my_basis_dc$semaine %in% c("semaine_t0-2", "semaine_t0-1")), ])

# semaine t+1
new_data_dc_1 <- my_basis_dc[my_basis_dc$semaine %in% c("semaine_t0-1", "semaine_t0-2", "semaine_t0-3", "semaine_t0-4"),  ]
pred_dc_1 <- predict(res_lm, newdata = new_data_dc_1)

pred_dc_2 <- matrix(0, length(nom_dep), 4)
pred_dc_2a <- matrix(0, length(nom_dep), 4)
pred_dc_2b <- matrix(0, length(nom_dep), 4)
pred_dc_2c <- matrix(0, length(nom_dep), 4)
for (k in 1:length(nom_dep)) {
  temp <- hospital[!(hospital$semaine %in% c("semaine_t0-4", "semaine_t0-3", "semaine_t0-2", "semaine_t0-1")) & 
                     hospital$dep == nom_dep[k], ]
  my_ts <- zoo(temp$incid_dc, temp$jour)
  my_ts_diff <- diff(my_ts)
  # tseries::adf.test(my_ts) 
  # tseries::adf.test(my_ts_diff)
  if (nom_dep[k] %in% c("975", "977", "978")) {
    pred_rea_2[k] <- NA
  } else {
    my_mod <- forecast::auto.arima(my_ts_diff)
    forecast_my_mod <- as.numeric(forecast(my_mod, h = 28)$mean)
    pred_3a_s1 <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:7])
    pred_3a_s2 <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:14]) - pred_3a_s1
    pred_3a_s3 <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:21]) - 
      pred_3a_s1 - pred_3a_s2
    pred_3a_s4 <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:28]) - 
      pred_3a_s1 - pred_3a_s2 - pred_3a_s3
        # Lissage exponentiel
    my_mod_exp <- ets(my_ts)
    forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp, h = 28)$mean)
    forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
    pred_3b_s1 <- round(sum(forecast_my_mod_exp[1:7]), 0)
    pred_3b_s2 <- round(sum(forecast_my_mod_exp[1:14]), 0) - pred_3b_s1
    pred_3b_s3 <- round(sum(forecast_my_mod_exp[1:21]), 0) - pred_3b_s2 - pred_3b_s1
    pred_3b_s4 <- round(sum(forecast_my_mod_exp[1:28]), 0) - pred_3b_s3 - pred_3b_s2 - pred_3b_s1
    # Méthode 3 : lissage exponentiel sur données hebdomadaires
    temp <- my_basis[my_basis$dep == nom_dep[k] & !(my_basis$semaine %in% 
                                          c("semaine_t0-4", "semaine_t0-3", "semaine_t0-2", "semaine_t0-1")), ]
    my_ts_exp <- zoo(temp$dc, temp$jour)
    my_mod_exp_2 <- ets(my_ts_exp)
    forecast_my_mod_exp <- as.numeric(forecast(my_mod_exp_2)$mean)
    forecast_my_mod_exp <- ifelse(forecast_my_mod_exp > 0, forecast_my_mod_exp, 0)
    pred_3c_s1 <- round(forecast_my_mod_exp[1], 0)
    pred_3c_s2 <- round(forecast_my_mod_exp[2], 0)
    pred_3c_s3 <- round(forecast_my_mod_exp[3], 0)    
    pred_3c_s4 <- round(forecast_my_mod_exp[4], 0) 
      
    pred_dc_2[k, 1] <- predict(lm_2_dc_ts, newdata = data.frame(pred_dc_2a = pred_3a_s1,
                                                 pred_dc_2b = pred_3b_s1,
                                                 pred_dc_2c = pred_3c_s1))
    pred_dc_2[k, 2] <- predict(lm_2_dc_ts_s1, newdata = data.frame(pred_dc_2a_s2 = pred_3a_s2,
                                                 pred_dc_2b_s2 = pred_3b_s2,
                                                 pred_dc_2c_s2 = pred_3c_s2))
    pred_dc_2[k, 3] <- predict(lm_2_dc_ts_s2, newdata = data.frame(pred_dc_2a_s3 = pred_3a_s3,
                                                 pred_dc_2b_s3 = pred_3b_s3,
                                                 pred_dc_2c_s3 = pred_3c_s3))
    pred_dc_2[k, 4] <- predict(lm_2_dc_ts_s3, newdata = data.frame(pred_dc_2a_s4 = pred_3a_s4,
                                                 pred_dc_2b_s4 = pred_3b_s4,
                                                 pred_dc_2c_s4 = pred_3c_s4))  
  }
}

n_4 <- 4 * length(nom_dep) 

pred_dc_a <- predict(lm_dc_3, newdata = data.frame(pred_dc_1 = 
                        pred_dc_1[new_data_dc_1$semaine == "semaine_t0-1"],
                        pred_dc_2 = as.vector(pred_dc_2)[1:(n_4 / 4)]))

pred_dc_b <- predict(lm_dc_3b, newdata = data.frame(
  pred_dc_1_s2 = pred_dc_1[new_data_dc_1$semaine == "semaine_t0-2"],
  pred_dc_2_s2 = as.vector(pred_dc_2)[((n_4 / 4) + 1):(2 * n_4 / 4)]))

pred_dc_c <- predict(lm_dc_3c, newdata = data.frame(
  pred_dc_1_s3 = pred_dc_1[new_data_dc_1$semaine == "semaine_t0-3"],
  pred_dc_2_s3 = as.vector(pred_dc_2)[(2 * n_4 / 4 + 1):(3 * n_4 / 4)]))

pred_dc_d <- predict(lm_dc_4c, newdata = data.frame(
  pred_dc_1_s4 = pred_dc_1[new_data_dc_1$semaine == "semaine_t0-4"],
  pred_dc_2_s4 = as.vector(pred_dc_2)[(3 * n_4 / 4 + 1):n_4]))



# on synthétise les résultats
new_data <- my_basis[my_basis$semaine %in% "semaine_t0-1", ]
new_data$this_week <- my_basis_dc[my_basis_dc$semaine == "semaine_t00", "dc"]
new_data$next_week <- pred_dc_a
new_data$next_two_week <- pred_dc_b
new_data$next_three_week <- pred_dc_c
new_data$next_four_week <- pred_dc_d
my_basis[my_basis$semaine %in% "semaine_t0-1", "dc"] <- pred_dc_a
my_basis[my_basis$semaine %in% "semaine_t0-2", "dc"] <- pred_dc_b
my_basis[my_basis$semaine %in% "semaine_t0-3", "dc"] <- pred_dc_c
my_basis[my_basis$semaine %in% "semaine_t0-4", "dc"] <- pred_dc_d

On va représenter l’évolution du nombre de décès dans un intervalle de temps de 5 semaines :

  • les 7 derniers jours passés : [22 novembre 2020; 28 novembre 2020]
  • la semaine à venir : [29 novembre 2020; 05 décembre 2020]
  • la 2ème semaine à venir : [06 décembre 2020; 12 décembre 2020]
  • la 3ème semaine à venir : [13 décembre 2020; 19 décembre 2020]
  • la 4ème semaine à venir : [13 décembre 2020; 19 décembre 2020]

On aggrège les données à la France entière:

Soit un nombre de décès par jour de :

## # A tibble: 5 x 2
##   semaine               dc
##   <fct>              <dbl>
## 1 [22 nov.; 28 nov.]  358 
## 2 [29 nov.; 05 déc.]  284.
## 3 [06 déc.; 12 déc.]  216.
## 4 [13 déc.; 19 déc.]  181.
## 5 [20 déc.; 26 déc.]  154.