Analyse des données de nouvelles hospitalisations/réanimations/décès de la COVID-19

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:

• Données hospitalières relatives à l’épidémie de COVID-19 : https://www.data.gouv.fr/fr/datasets/donnees-hospitalieres-relatives-a-lepidemie-de-covid-19/

• Données relatives aux résultats des tests virologiques COVID-19 : https://www.data.gouv.fr/fr/datasets/donnees-relatives-aux-resultats-des-tests-virologiques-covid-19/

Questions de recherche:

• Peut-on expliquer le nombre de nouvelles entrées par semaine en fonction du nombre de personnes testées positives quelques jours auparavant ?

• Peut-on expliquer les entrées en réanimation par les entrées en hospitalisations ?

• Peut-on expliquer les nouveaux décès par les entrées réanimations ?

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

• Peut-on prédire à court terme les nouveaux cas positifs, les nouvelles hospitalisations, les nouvelles réanimations et les nouveaux décès ?

Packages et fonctions locales à charger:

library(cartogram)
library(cartography)

## Warning in fun(libname, pkgname): rgeos: versions of GEOS runtime 3.9.0-CAPI-1.16.2
## and GEOS at installation 3.8.0-CAPI-1.13.1differ

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.4                         
## year           2021                        
## month          02                          
## day            15                          
## svn rev        80002                       
## language       R                           
## version.string R version 4.0.4 (2021-02-15)
## nickname       Lost Library Book

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.

# Date du jour pour actualiser les données:
to_day <- Sys.Date()
my_url <- "https://www.data.gouv.fr/fr/datasets/r/6fadff46-9efd-4c53-942a-54aca783c30c"
if (!file.exists(paste0(getwd(), "/data/", to_day, ".csv"))) {
    download.file(my_url, destfile = paste0(getwd(), "/data/", to_day, ".csv"))
}
hospital <- read.csv(paste0(getwd(), "/data/", to_day, ".csv"), sep = ";")
# On ajoute le nom des régions:
hospital <- merge(hospital, dep_region, by = "dep")
# On utilise le format date pour coder le jour:
hospital$jour <- as.Date(hospital$jour)

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

temp <- hospital[hospital$dep == "69" & hospital$jour > to_day - 5, ]
kableExtra::kbl(temp)

	dep	jour	incid_hosp	incid_rea	incid_dc	incid_rad	nom_dep	region
25372	69	2021-03-16	60	23	18	59	Rhône	Auvergne-Rhône-Alpes
25433	69	2021-03-17	49	17	9	70	Rhône	Auvergne-Rhône-Alpes
25481	69	2021-03-18	35	4	7	47	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 19 mars 2021 la semaine qui correspond à la semaine t0 correspond à la fenêtre [13 mars 2021; 19 mars 2021]. Dans chaque fenêtre, on calcule le nombre de nouvelles hospitalisations, réanimations et décès par département.

hospital$semaine <- num_semaine(hospital$jour)
# On aggrège les données en fonction de cette fenêtre et on garde tous les départements:
my_basis <- hospital %>%
  group_by(dep, semaine) %>%
  dplyr::summarize(hosp = sum(incid_hosp),
            rea = sum(incid_rea),
            rad = sum(incid_rad),
            dc = sum(incid_dc),
            jour = max(jour),
            region = unique(region))

timeline_data <- data.frame(event = c("Semaine t0", "Semaine t1", "Semaine t2", "Semaine t3"), 
                            start = c(to_day - 8, to_day - 15, to_day - 22, to_day - 29), 
                            end = c(to_day - 1, to_day - 8, to_day - 15, to_day - 22), group = "temps")
vistime(timeline_data)

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 : [13 mars 2021; 19 mars 2021].

# On aggrège les données par région sur la semaine `r paste0("[", format(to_day - 7, '%d %B %Y'), "; ",  format(to_day - 1, '%d %B %Y'), "]")` :
vs_my_basis_t0 <- my_basis %>%
  filter(semaine %in% "semaine_t00") %>%
  group_by(region) %>%
  summarise(hosp = sum(hosp),
            rea = sum(rea),
            dc = sum(dc)) 
# On aggrège les données par région sur la semaine `r paste0("[", format(to_day - 13, '%d %B %Y'), "; ",  format(to_day - 7, '%d %B %Y'), "]")` :
vs_my_basis_t1 <- my_basis %>%
  filter(semaine %in% "semaine_t01") %>%
  group_by(region) %>%
  summarise(hosp = sum(hosp),
            rea = sum(rea),
            dc = sum(dc)) 

1.2.1 Résumé des hospitalisations

On représente par région:

• le nombre total de nouvelles hospitalisations (semaine [13 mars 2021; 19 mars 2021]).

• le nombre moyen journalier de nouvelles hospitalisations (semaine [13 mars 2021; 19 mars 2021]).

• l’évolution (en pourcentage) entre la semaine [06 mars 2021; 12 mars 2021] et la semaine [13 mars 2021; 19 mars 2021].

hosp_region <- vs_my_basis_t0 %>%
  select(region, hosp) %>%
  mutate(`moyenne jour` = hosp / 7,
         `evolution en %` = (vs_my_basis_t0$hosp - vs_my_basis_t1$hosp) / vs_my_basis_t1$hosp * 100) %>%
    rename(`total semaine` = hosp) %>%
  arrange(-`total semaine`)
# On représente les données :
hosp_region[, 3] <- round(hosp_region[, 3])
hosp_region[, 4] <- round(hosp_region[, 4], 1)
hosp_region[2:4] <- lapply(hosp_region[2:4], function(x) {
    cell_spec(x, bold = T, color = spec_color(x, end = 0.9),
              font_size = spec_font_size(x))
})
hosp_region <- rbind(hosp_region, tibble(region = "France entière", 
          `total semaine` = sum(vs_my_basis_t0$hosp), 
          `moyenne jour` = round(sum(vs_my_basis_t0$hosp) / 7, 0), 
          `evolution en %` = round((sum(vs_my_basis_t0$hosp) - sum(vs_my_basis_t1$hosp)) / 
                                     sum(vs_my_basis_t1$hosp) * 100, 1)))
kbl(hosp_region, escape = F, align = "c") %>% kable_classic("striped", full_width = F)

region	total semaine	moyenne jour	evolution en %
Ile-de-France	2879	411	4.4
Hauts-de-France	1453	208	23
Provence-Alpes-Côte d’Azur	1235	176	-3.4
Auvergne-Rhône-Alpes	1022	146	-2.1
Grand Est	919	131	15.5
Occitanie	489	70	12.9
Normandie	460	66	2.2
Nouvelle-Aquitaine	435	62	-1.8
Bourgogne-Franche-Comté	379	54	1.9
Pays de la Loire	372	53	18.5
Centre-Val de Loire	324	46	-3.9
Bretagne	261	37	28.6
DOM-TOM	183	26	16.6
Corse	32	5	23.1
France entière	10443	1492	6.7

On représente par département la carte des nouvelles hospitalisations sur la dernière semaine observée ([13 mars 2021; 19 mars 2021])

# On importe les contours des départements
dep.2015 <- readOGR(dsn="./departements 2015/DEPARTEMENT", layer="DEPARTEMENT")

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

dep.2015@data$CODE_DEPT <- as.character(dep.2015@data$CODE_DEPT) 

# data 
dep.2015_00 <- merge(dep.2015, filter(my_basis, semaine == "semaine_t00"), 
                     by.x = "CODE_DEPT", by.y = "dep")
# On représente les nouvelles hospitalisations sur la dernière semaine observée (`r paste0("[", format(to_day - 7, '%d %B %Y'), "; ",  format(to_day - 1, '%d %B %Y'), "]")`) (Source pour le code : https://rgeomatic.hypotheses.org/1361#more-1361)
# quantization breaks of the rate
bks <- c(0, getBreaks(v = dep.2015_00$hosp, method = "kmeans", nclass = 5))
# correct the breaks to use the global rate as limit of class 
# get a color palette
cols <- carto.pal(pal1 = "green.pal", n1 = 3, pal2 = "wine.pal", n2 = 3)
## Choropleth layer
# set figure margins and background color
par(mar = c(0, 0, 1.2, 0), bg = "lemonchiffon")
# Hospitalisations
choroLayer(spdf = dep.2015_00, var = "hosp", breaks = bks, col = cols,
           border = "khaki", lwd = 0.5, 
           legend.title.txt = "Hospitalisations", 
           legend.pos = 'topleft', legend.values.rnd = 0)
# add a title and layout
layoutLayer(title = paste0("Nouvelles hospitalisations ", 
  paste0("[", format(to_day - 7, '%d %B %Y'), "; ",  format(to_day - 1, '%d %B %Y'), "]")), 
            sources = "", north = TRUE, scale = 50, tabtitle = TRUE,
            theme = "sand.pal", frame = FALSE,  
            author = "")

1.2.2 Résumé des réanimations

On représente par région:

• le nombre total de nouvelles réanimations (semaine [13 mars 2021; 19 mars 2021]).

• le nombre moyen journalier de nouvelles réanimations (semaine [13 mars 2021; 19 mars 2021]).

• l’évolution (en pourcentage) entre la semaine [06 mars 2021; 12 mars 2021] et la semaine [13 mars 2021; 19 mars 2021].

rea_region <- vs_my_basis_t0 %>%
  select(region, rea) %>%
  mutate(`moyenne jour` = rea / 7,
         `evolution en %` = (vs_my_basis_t0$rea - vs_my_basis_t1$rea) / vs_my_basis_t1$rea * 100) %>%
    rename(`total semaine` = rea) %>%
  arrange(-`total semaine`)
# On représente les données :
rea_region[, 3] <- round(rea_region[, 3])
rea_region[, 4] <- round(rea_region[, 4], 1)
rea_region[2:4] <- lapply(rea_region[2:4], function(x) {
    cell_spec(x, bold = T, color = spec_color(x, end = 0.9),
              font_size = spec_font_size(x))
})

rea_region <- rbind(rea_region, tibble(region = "France entière", 
          `total semaine` = sum(vs_my_basis_t0$rea), 
          `moyenne jour` = round(sum(vs_my_basis_t0$rea) / 7, 0), 
          `evolution en %` = round((sum(vs_my_basis_t0$rea) - sum(vs_my_basis_t1$rea)) / 
                                     sum(vs_my_basis_t1$rea) * 100, 1)))

kbl(rea_region, escape = F, align = "c") %>% kable_classic("striped", full_width = F)

region	total semaine	moyenne jour	evolution en %
Ile-de-France	733	105	-3.6
Hauts-de-France	306	44	7
Provence-Alpes-Côte d’Azur	267	38	3.9
Auvergne-Rhône-Alpes	224	32	-7.1
Grand Est	203	29	38.1
Occitanie	111	16	-19.6
Nouvelle-Aquitaine	91	13	-7.1
Normandie	85	12	16.4
Bourgogne-Franche-Comté	72	10	4.3
Centre-Val de Loire	67	10	-6.9
Pays de la Loire	60	9	1.7
DOM-TOM	55	8	41
Bretagne	48	7	50
Corse	10	1	233.3
France entière	2332	333	2.6

On représente par département la carte des nouvelles réanimations sur la dernière semaine observée ([13 mars 2021; 19 mars 2021])

# quantization breaks of the rate
bks <- c(0, getBreaks(v = dep.2015_00$rea, method = "kmeans", nclass = 5))
# correct the breaks to use the global rate as limit of class 
# get a color palette
cols <- carto.pal(pal1 = "green.pal", n1 = 3, pal2 = "wine.pal", n2 = 3)
## Choropleth layer
# set figure margins and background color
par(mar = c(0, 0, 1.2, 0), bg = "lemonchiffon")
# Hospitalisations
choroLayer(spdf = dep.2015_00, var = "rea", breaks = bks, col = cols,
           border = "khaki", lwd = 0.5, 
           legend.title.txt = "Réanimations", 
           legend.pos = 'topleft', legend.values.rnd = 0)
# add a title and layout
layoutLayer(title = paste0("Nouvelles Réanimations ", 
  paste0("[", format(to_day - 7, '%d %B %Y'), "; ",  format(to_day - 1, '%d %B %Y'), "]")), 
            sources = "", north = TRUE, scale = 50, tabtitle = TRUE,
            theme = "sand.pal", frame = FALSE,  
            author = "")

1.2.3 Résumé des décès

On représente par région:

• le nombre total de nouveaux décès (semaine [13 mars 2021; 19 mars 2021]).

• le nombre moyen journalier de nouveaux décès (semaine [13 mars 2021; 19 mars 2021]).

• l’évolution (en pourcentage) entre la semaine [06 mars 2021; 12 mars 2021] et la semaine [13 mars 2021; 19 mars 2021].

dc_region <- vs_my_basis_t0 %>%
  select(region, dc) %>%
  mutate(`moyenne jour` = dc / 7,
         `evolution en %` = (vs_my_basis_t0$dc - vs_my_basis_t1$dc) / vs_my_basis_t1$dc * 100) %>%
    rename(`total semaine` = dc) %>%
  arrange(-`total semaine`)
# On représente les données :
dc_region[, 3] <- round(dc_region[, 3])
dc_region[, 4] <- round(dc_region[, 4], 1)
dc_region[2:4] <- lapply(dc_region[2:4], function(x) {
    cell_spec(x, bold = T, color = spec_color(x, end = 0.9),
              font_size = spec_font_size(x))
})

dc_region <- rbind(dc_region, tibble(region = "France entière", 
          `total semaine` = sum(vs_my_basis_t0$dc), 
          `moyenne jour` = round(sum(vs_my_basis_t0$dc) / 7, 0), 
          `evolution en %` = round((sum(vs_my_basis_t0$dc) - sum(vs_my_basis_t1$dc)) / 
                                     sum(vs_my_basis_t1$dc) * 100, 1)))

kbl(dc_region, escape = F, align = "c") %>% kable_classic("striped", full_width = F)

region	total semaine	moyenne jour	evolution en %
Ile-de-France	388	55	-6.5
Hauts-de-France	242	35	1.7
Provence-Alpes-Côte d’Azur	212	30	-15.2
Auvergne-Rhône-Alpes	203	29	-2.9
Grand Est	171	24	0
Nouvelle-Aquitaine	89	13	6
Normandie	85	12	16.4
Bourgogne-Franche-Comté	81	12	26.6
Occitanie	72	10	-2.7
Centre-Val de Loire	57	8	26.7
Pays de la Loire	51	7	4.1
Bretagne	31	4	-24.4
DOM-TOM	27	4	92.9
Corse	5	1	-16.7
France entière	1714	245	-1.1

On représente par département la carte des nouveaux décès sur la dernière semaine observée ([13 mars 2021; 19 mars 2021])

# quantization breaks of the rate
bks <- c(0, getBreaks(v = dep.2015_00$dc, method = "kmeans", nclass = 5))
# correct the breaks to use the global rate as limit of class 
# get a color palette
cols <- carto.pal(pal1 = "green.pal", n1 = 3, pal2 = "wine.pal", n2 = 3)
## Choropleth layer
# set figure margins and background color
par(mar = c(0, 0, 1.2, 0), bg = "lemonchiffon")
# Hospitalisations
choroLayer(spdf = dep.2015_00, var = "dc", breaks = bks, col = cols,
           border = "khaki", lwd = 0.5, 
           legend.title.txt = "Décès", 
           legend.pos = 'topleft', legend.values.rnd = 0)
# add a title and layout
layoutLayer(title = paste0("Nouveaux décès ", 
  paste0("[", format(to_day - 7, '%d %B %Y'), "; ",  format(to_day - 1, '%d %B %Y'), "]")), 
            sources = "", north = TRUE, scale = 50, tabtitle = TRUE,
            theme = "sand.pal", frame = FALSE,  
            author = "")

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.

my_basis_fr <- my_basis %>%
  group_by(semaine, jour) %>%
  summarise(hosp = sum(hosp),
            rea = sum(rea),
            dc = sum(dc))
p <- ggplot(data = my_basis_fr,
            aes(x = jour, y = hosp)) +
  geom_line() +
  labs(title = "Nouvelles hospitalisations par semaine",
       x = "semaine",
       y = "hospitalisations",
       fill = "Age") 
plotly::ggplotly(p)

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

my_basis$region <- factor(my_basis$region, levels = hosp_region$region)
# On aggrège les données par région et semaine:
my_basis_by_region <- my_basis %>%
  group_by(semaine, region) %>%
  summarize(hosp = sum(hosp),
            rea = sum(rea),
            dc = sum(dc),
            jour = max(jour))
# On représente les shares par tranche d'age 
p <- ggplot(my_basis_by_region) + 
  aes(x = jour, y = hosp, fill = region) +
  geom_area(color = "black") +
  labs(title = "Nouvelles hospitalisations par région",
       x = "Semaine",
       y = "Total",
       fill = "Hospitalisation") +
  scale_fill_brewer(palette = "Set3") +
  theme_minimal() 
plotly::ggplotly(p)

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

# On aggrège les données par age et semaine:
my_basis_by_age_1 <- my_basis_age %>%
  group_by(semaine, cl_age90) %>%
  summarize(hosp = sum(hosp),
            rea = sum(rea),
            dc = sum(dc),
            jour = max(jour))
# On représente les shares par tranche d'age 
p <- ggplot(filter(my_basis_by_age_1, cl_age90 != 0)) + 
  aes(x = jour, y = hosp, fill = factor(cl_age90)) +
  geom_area(color = "black") +
  labs(title = "Nouvelles hospitalisations par classe d'âges",
       x = "Semaine",
       y = "Total",
       fill = "Hospitalisation") +
  scale_fill_brewer(palette = "Set3") +
  theme_minimal() 
plotly::ggplotly(p)

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

# On aggrège les données par age et semaine:
my_basis_by_age_2 <- my_basis_age %>%
  filter(cl_age90 != 0) %>%
  group_by(semaine, cl_age90) %>%
  summarize(hosp = sum(hosp),
            rea = sum(rea),
            dc = sum(dc),
            jour = max(jour)) %>%
  group_by(cl_age90)  %>%
  group_by(semaine) %>%
  mutate(percent_hosp = hosp / sum(hosp),
         percent_rea = rea / sum(rea), 
         percent_dc = dc / sum(dc))
# On représente les shares par tranche d'age 
p <- ggplot(my_basis_by_age_2) + 
  aes(x = jour, y = percent_hosp, fill = factor(cl_age90)) +
  geom_area(color = "black") +
  labs(title = "Répartition des nouvelles hospitalisations par classe d'âges",
       x = "Semaine",
       y = "Total",
       fill = "Hospitalisation") +
  scale_fill_brewer(palette = "Set3") +
  theme_minimal() 
plotly::ggplotly(p)

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.

p <- ggplot(data = my_basis_fr,
            aes(x = jour, y = rea)) +
  geom_line() +
  labs(title = "Nouvelles réanimations par semaine",
       x = "semaine",
       y = "réanimations",
       fill = "Age") 
plotly::ggplotly(p)

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

# On représente les shares par régions
p <- ggplot(my_basis_by_region) + 
  aes(x = jour, y = rea, fill = region) +
  geom_area(color = "black") +
  labs(title = "Nouvelles réanimations par région",
       x = "Semaine",
       y = "Total",
       fill = "Réanimations") +
  scale_fill_brewer(palette = "Set3") +
  theme_minimal() 
plotly::ggplotly(p)

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

p <- ggplot(filter(my_basis_by_age_1, cl_age90 != 0)) + 
  aes(x = jour, y = rea, fill = factor(cl_age90)) +
  geom_area(color = "black") +
  labs(title = "Nouvelles réanimations par classe d'âges",
       x = "Semaine",
       y = "Total",
       fill = "Réanimations") +
  scale_fill_brewer(palette = "Set3") +
  theme_minimal() 
plotly::ggplotly(p)

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

# On représente les shares par tranche d'age 
p <- ggplot(my_basis_by_age_2) + 
  aes(x = jour, y = percent_rea, fill = factor(cl_age90)) +
  geom_area(color = "black") +
  labs(title = "Répartition des nouvelles réanimations par classe d'âges",
       x = "Semaine",
       y = "Total",
       fill = "Réanimations") +
  scale_fill_brewer(palette = "Set3") +
  theme_minimal() 
plotly::ggplotly(p)

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

p <- ggplot(data = my_basis_by_age_2,
            aes(x = jour, y = rea / hosp, col = factor(cl_age90))) +
  geom_line() +
  labs(title = "Ratio réanimations/hospitalisations par semaine",
       x = "semaine",
       y = "Ratio réa / hospi",
       color = "Age") 
plotly::ggplotly(p)

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.

p <- ggplot(data = my_basis_fr,
            aes(x = jour, y = dc)) +
  geom_line() +
  labs(title = "Nouveux décès par semaine",
       x = "semaine",
       y = "décès",
       fill = "Age") 
plotly::ggplotly(p)

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

# On représente les shares par régions
p <- ggplot(my_basis_by_region) + 
  aes(x = jour, y = dc, fill = region) +
  geom_area(color = "black") +
  labs(title = "Nouveux décès par région",
       x = "Semaine",
       y = "Total",
       fill = "Décès") +
  scale_fill_brewer(palette = "Set3") +
  theme_minimal() 
plotly::ggplotly(p)

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

p <- ggplot(filter(my_basis_by_age_1, cl_age90 != 0)) + 
  aes(x = jour, y = dc, fill = factor(cl_age90)) +
  geom_area(color = "black") +
  labs(title = "Nouveux décès par classe d'âges",
       x = "Semaine",
       y = "Total",
       fill = "Décès") +
  scale_fill_brewer(palette = "Set3") +
  theme_minimal() 
plotly::ggplotly(p)

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

# On représente les shares par tranche d'age 
p <- ggplot(my_basis_by_age_2) + 
  aes(x = jour, y = percent_dc, fill = factor(cl_age90)) +
  geom_area(color = "black") +
  labs(title = "Répartition des nouveux décès par classe d'âges",
       x = "Semaine",
       y = "Total",
       fill = "Décès") +
  scale_fill_brewer(palette = "Set3") +
  theme_minimal() 
plotly::ggplotly(p)

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

p <- ggplot(data = my_basis_by_age_2,
            aes(x = jour, y = dc / rea, col = factor(cl_age90))) +
  geom_line() +
  labs(title = "Ratio décès/réanimations par semaine",
       x = "semaine",
       y = "Ratio décès / réanimations",
       color = "Age") 
plotly::ggplotly(p)

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.

p <- ggplot(data = filter(my_basis, region %in% c("Provence-Alpes-Côte d'Azur",
                                                  "Auvergne-Rhône-Alpes", 
                                                  "Hauts-de-France",
                                                   "Ile-de-France")),
            aes(x = jour, y = hosp, color = region, group = dep)) +
  geom_line() + 
  coord_cartesian(ylim = c(0, 1200)) +
  facet_wrap(~ region, nrow = 4)
plotly::ggplotly(p)

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

p <- ggplot(data = filter(my_basis, region %in% c("Occitanie",
                                                  "Grand Est", 
                                                  "Bourgogne-Franche-Comté",
                                                  "Nouvelle-Aquitaine", 
                                                  "Normandie",
                                                  "Pays de la Loire", 
                                                  "Centre-Val de Loire",
                                                  "Bretagne")),
            aes(x = jour, y = hosp, color = region, group = dep)) +
  geom_line()  + 
  coord_cartesian(ylim = c(0, 400)) +
  facet_wrap(~ region, nrow = 8)
plotly::ggplotly(p)

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

p <- ggplot(data = filter(my_basis, region %in% c("DOM-TOM", "Corse")),
            aes(x = jour, y = hosp, color = region, group = dep)) +
  geom_line()  + 
  coord_cartesian(ylim = c(0, 200)) +
  facet_wrap(~ region, nrow = 2)
plotly::ggplotly(p)

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 [13 mars 2021; 19 mars 2021] et sur la période [05 mars 2021; 12 mars 2021]

my_basis_evol <- merge(my_basis %>% 
  filter(semaine == "semaine_t00") %>%
  rename(hosp_t0 = hosp) %>%
  select(dep, hosp_t0, region),
    my_basis %>% 
  filter(semaine == "semaine_t01") %>%
  select(dep, hosp) %>%
  rename(hosp_t1 = hosp),
by = "dep") %>%
  mutate(diff_abs = hosp_t0 - hosp_t1,
         diff_rel = (hosp_t0 - hosp_t1) / hosp_t1) %>%
  mutate(evol = factor(ifelse(diff_rel < 0, "<0", 
                       ifelse(diff_rel >= 0 & diff_rel < 0.5, "[0,50%[",
                           ifelse(diff_rel >= 0.5 & diff_rel < 1, "[50,100%[", 
                                  ifelse(diff_rel >= 1 & diff_rel < 2, "[100,200%[",
                                         ">200%")))), 
                       levels = c("<0", "[0,50%[", "[50,100%[", "[100,200%[", ">200%")))
my_basis_evol_long <- tidyr::pivot_longer(data = my_basis_evol,
                                   col = c(2, 4),
                                   names_to = "semaine",
                                   values_to = "hospitalisations")
my_basis_evol_long$semaine <- factor(my_basis_evol_long$semaine,
      levels = c("hosp_t1", "hosp_t0"),
      labels = c(paste0("[", format(to_day - 14, '%d %b'), "; ",  
                        format(to_day - 8, '%d %b'), "]"),
                 paste0("[", format(to_day - 7, '%d %b'), "; ",  
                        format(to_day - 1, '%d %b'), "]")))

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

p <- ggplot(my_basis_evol_long, aes(x = semaine, y = hospitalisations, colour = evol, group = dep))+
    geom_line() +
  scale_colour_manual(values = c("blue", "#FC9272", "#FB6A4A", "#DE2D26", "#A50F15")) + 
  facet_wrap(~region)
plotly::ggplotly(p)

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

On met à jour les données chaque jour :

# On récupére directement l'url depuis le site du ministère:
url <- "https://www.data.gouv.fr/fr/datasets/r/406c6a23-e283-4300-9484-54e78c8ae675"
if (!file.exists(paste0(getwd(), "/data/test", to_day, ".csv"))) {
  download.file(url, destfile = paste0(getwd(), "/data/test", to_day, ".csv"))
    }
test <- read.csv(paste0(getwd(), "/data/test", to_day, ".csv"), sep = ";")
# on enlève les départements qui ne sont pas présents dans la table hopital
test <- test %>%
  filter(!(dep %in% c("975", "977", "978")))
# On utilise le format date pour coder le jour:
test$jour <- as.Date(test$jour)
test$semaine <- num_semaine(test$jour, begin = max(hospital$jour) - 3, decallage = TRUE)
test_region <- merge(test, dep_region, by = "dep")

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

# On aggrège les données par région sur la semaine `r paste0("[", format(to_day - 7, '%d %B %Y'), "; ",  format(to_day - 1, '%d %B %Y'), "]")` :
vs_my_basis_t0 <- test_region %>%
  filter(semaine %in% "semaine_t0-1", cl_age90 == 0) %>%
  group_by(region) %>%
  summarise(P = sum(P)) 
# On aggrège les données par région sur la semaine `r paste0("[", format(to_day - 13, '%d %B %Y'), "; ",  format(to_day - 7, '%d %B %Y'), "]")` :
vs_my_basis_t1 <- test_region %>%
  filter(semaine %in% "semaine_t00", cl_age90 == 0) %>%
  group_by(region) %>%
  summarise(P = sum(P)) 

On représente par région:

• le nombre total de testés positifs (semaine [10 mars 2021; 16 mars 2021]).

• le nombre moyen journalier de testés positifs (semaine [10 mars 2021; 16 mars 2021]).

• l’évolution (en pourcentage) entre la semaine [03 mars 2021; 09 mars 2021] et la semaine [10 mars 2021; 16 mars 2021].

P_region <- vs_my_basis_t0 %>%
  select(region, P) %>%
  mutate(`moyenne jour` = P / 7,
         `evolution en %` = (vs_my_basis_t0$P - vs_my_basis_t1$P) / vs_my_basis_t1$P * 100) %>%
    rename(`total semaine` = P) %>%
  arrange(-`total semaine`)
# On représente les données :
P_region[, 3] <- round(P_region[, 3])
P_region[, 4] <- round(P_region[, 4], 1)
P_region[2:4] <- lapply(P_region[2:4], function(x) {
    cell_spec(x, bold = T, color = spec_color(x, end = 0.9),
              font_size = spec_font_size(x))
})
P_region <- rbind(P_region, tibble(region = "France entière", 
          `total semaine` = sum(vs_my_basis_t0$P), 
          `moyenne jour` = round(sum(vs_my_basis_t0$P) / 7, 0), 
          `evolution en %` = round((sum(vs_my_basis_t0$P) - sum(vs_my_basis_t1$P)) / 
                                     sum(vs_my_basis_t1$P) * 100, 1)))
kbl(P_region, escape = F, align = "c") %>% kable_classic("striped", full_width = F)

region	total semaine	moyenne jour	evolution en %
Ile-de-France	54778	7825	24
Hauts-de-France	22766	3252	12.4
Auvergne-Rhône-Alpes	18580	2654	24.3
Provence-Alpes-Côte d’Azur	16584	2369	-2.1
Grand Est	12104	1729	14.8
Occitanie	10613	1516	16
Nouvelle-Aquitaine	8640	1234	32.2
Normandie	8563	1223	35.6
Pays de la Loire	6740	963	17.9
Bourgogne-Franche-Comté	6004	858	39.3
Centre-Val de Loire	5165	738	17.1
Bretagne	4722	675	20.5
DOM-TOM	2381	340	9.8
Corse	514	73	-5.7
France entière	178154	25451	18.8

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 aggrège les données par tranche d'age et semaine:
test_by_age <- test %>%
  filter(cl_age90 != 0) %>%
  group_by(semaine, cl_age90) %>%
  summarize(P = sum(P),
            jour = max(jour)) %>%
  group_by(semaine) %>%
  mutate(percent_P = P / sum(P))
p <- ggplot(test_by_age) + 
  aes(x = jour, y = P, fill = factor(cl_age90)) +
  geom_area(color = "black") +
  labs(title = "Nombre de téstés positifs par classe d'âge ",
       x = "Week",
       y = "Effectifs",
       fill = "Age") +
  scale_fill_brewer(palette = "Set3") +
  theme_minimal() 
plotly::ggplotly(p)

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

test_region$region <- factor(test_region$region, levels = hosp_region$region)
test_by_region <- test_region %>%
  filter(cl_age90 == 0) %>%
  group_by(semaine, region) %>%
  summarize(P = sum(P),
            jour = max(jour)) %>%
  group_by(semaine) %>%
  mutate(percent_P = P / sum(P))
p <- ggplot(test_by_region) + 
  aes(x = jour, y = P, fill = factor(region)) +
  geom_area(color = "black") +
  labs(title = "Nombre de testés positifs par région",
       x = "Week",
       y = "Effectifs",
       fill = "Region") +
  scale_fill_brewer(palette = "Set3") +
  theme_minimal() 
plotly::ggplotly(p)

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

test_by_dep <- test_region %>%
  filter(cl_age90 == 0) %>%
  group_by(semaine, region, dep) %>%
  summarize(P = sum(P),
            jour = max(jour)) %>%
  group_by(semaine) %>%
  mutate(percent_P = P / sum(P))

p <- ggplot(data = test_by_dep, 
            aes(x = jour, y = P, color = region, group = dep)) +
  geom_line() +
  coord_cartesian(ylim = c(0, 20000)) +
  facet_wrap(~ region)
plotly::ggplotly(p)

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 [13 mars 2021; 19 mars 2021], on va l’expliquer par le nombre de personnes testées positive du [03 mars 2021; 09 mars 2021].

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

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

# On aggrège les données par département et semaine:
test_by <- test %>%
  filter(cl_age90 != 0) %>%
  group_by(dep, semaine, cl_age90) %>%
  summarize(P = sum(P),
            jour = max(jour)) %>%
  group_by(dep, semaine) %>%
  mutate(percent_P = P / sum(P))
# On passe d'un format long à un format wide les tranches d'âge pour les avoir comme des variables explicatives
test_by$cl_age90 <- paste0("tranche_", test_by$cl_age90)
test_long <- tidyr::pivot_wider(test_by,
                           id_cols = c("semaine", "dep", "jour", "P", "cl_age90"),
                           names_from = "cl_age90",
                           values_from = "P")
# On aggrège les données par département et semaine:
test_by <- test_long %>%
  group_by(dep, semaine) %>%
  summarize(tranche_9 = sum(tranche_9),
            tranche_19 = sum(tranche_19),
            tranche_29 = sum(tranche_29),
            tranche_39 = sum(tranche_39),
            tranche_49 = sum(tranche_49),
            tranche_59 = sum(tranche_59),
            tranche_69 = sum(tranche_69),
            tranche_79 = sum(tranche_79),
            tranche_89 = sum(tranche_89),
            tranche_90 = sum(tranche_90)) %>%
  mutate(tranche_0 = tranche_9 + tranche_19 + tranche_29 + tranche_39 + 
           tranche_49 + tranche_59 + tranche_69 + tranche_79 + tranche_89 + tranche_90)

# On merge les jeux de données :
my_basis <- merge(my_basis, test_by, by = c("dep", "semaine"), all.y = T)
my_basis[which(is.na(my_basis$jour)), "jour"] <- to_day + 6
my_basis[which(is.na(my_basis$region)), "region"] <- dep_region$region[match(my_basis[which(is.na(my_basis$region)), "dep"], dep_region$dep)]

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 [03 mars 2021; 09 mars 2021]. La couleur dépend du nombre d’hospitalisations observés la semaine du [13 mars 2021; 19 mars 2021].

dep.2015_00 <- merge(dep.2015, filter(my_basis, semaine == "semaine_t00"), 
                     by.x = "CODE_DEPT", by.y = "dep")
w <- 1 - (dep.2015_00$tranche_0 / max(dep.2015_00$tranche_0)) 
dep.2015_dorling <- cartogram_dorling(dep.2015_00, "tranche_0", m_weight = w, k = 5)
# dep.2015$tranche_0 <- as.numeric(dep.2015$tranche_0)
# dep.2015_ncont <- cartogram_ncont(dep.2015, "tranche_0")
# set figure margins and background color
par(mar = c(0, 0, 1.2, 0), bg = "lemonchiffon")
# Hospitalisations
bks <- c(0, getBreaks(v = dep.2015_00$hosp, method = "kmeans", nclass = 5))
# correct the breaks to use the global rate as limit of class 
# get a color palette
cols <- carto.pal(pal1 = "green.pal", n1 = 3, pal2 = "wine.pal", n2 = 3)
choroLayer(spdf = dep.2015_dorling, var = "hosp", breaks = bks, col = cols, lwd = 0.5,
           legend.title.txt = "Hospitalisations", 
           legend.pos = 'topleft', legend.values.rnd = 0)
# plot(dep.2015, add = T, border = "khaki")
# add a title and layout
layoutLayer(title = paste0("Nouvelles hospitalisations ", 
  paste0("[", format(to_day - 7, '%d %B %Y'), "; ",  format(to_day - 1, '%d %B %Y'), "]")), 
            sources = "", north = TRUE, scale = 50, tabtitle = TRUE,
            theme = "sand.pal", frame = FALSE,  
            author = "")

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.

p <- ggplot(my_basis) +                
  aes(x = tranche_0, y = hosp) +     
  geom_point() +                  
  geom_smooth(method = "loess") + 
  geom_smooth(method = "lm",     
              col = "red") +
  facet_wrap(~ region)
plotly::ggplotly(p)

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 16 mars 2021. Notre objectif est de prédire le nombre de testés positifs du 17 mars 2021 au 23 mars 2021 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)
  }
}

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

## Warning in ets(my_ts_exp): Missing values encountered. Using longest contiguous
## portion of time series

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

test_by_region <- my_basis %>%
  group_by(semaine, region) %>%
  summarize(P = sum(tranche_0),
            jour = max(jour)) 
test_by_region$region <- factor(test_by_region$region, levels = hosp_region$region)
p <- ggplot(test_by_region) + 
  aes(x = jour - 10, y = P, fill = factor(region)) +
  geom_area(color = "black") +
  labs(title = "Prédictions du nombre de testés positifs par région",
       x = "Week",
       y = "Effectifs",
       fill = "Region") +
  scale_fill_brewer(palette = "Set3") +
  theme_minimal() 
plotly::ggplotly(p)

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

  my_basis_fr <- my_basis %>%
    group_by(semaine, jour) %>%
    summarise(P = sum(tranche_0))
  p <- ggplot(data = filter(my_basis_fr, semaine %in% c("semaine_t0-2", "semaine_t0-1", "semaine_t00")),
              aes(x = jour - 10, y = P)) +
    geom_line(col = "red") +
    geom_line(data = filter(my_basis_fr, !(semaine %in% c("semaine_t0-2"))), 
              aes(x = jour - 10, y = P)) +
  labs(title = "Prédiction des nouveaux cas positifs dans les 7 jours",
       x = "semaine",
       y = "Cas positifs",
       fill = "Age") 
plotly::ggplotly(p)

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:

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

stargazer::stargazer(res_lm, type = "html")

	Dependent variable:
	hosp
regionAuvergne-Rhône-Alpes	14.475***
	(1.719)
regionBourgogne-Franche-Comté	10.344***
	(2.283)
regionBretagne	4.124
	(3.373)
regionCentre-Val de Loire	9.072***
	(2.753)
regionCorse	2.086
	(4.672)
regionDOM-TOM	3.740
	(2.923)
regionGrand Est	12.535***
	(2.032)
regionHauts-de-France	28.005***
	(2.717)
regionIle-de-France	29.832***
	(2.651)
regionNormandie	6.344**
	(2.896)
regionNouvelle-Aquitaine	2.975
	(1.820)
regionOccitanie	3.622**
	(1.677)
regionPays de la Loire	14.125***
	(2.994)
regionProvence-Alpes-Côte d’Azur	9.750***
	(2.567)
regionAuvergne-Rhône-Alpes:tranche_0	0.063***
	(0.001)
regionBourgogne-Franche-Comté:tranche_0	0.073***
	(0.003)
regionBretagne:tranche_0	0.058***
	(0.004)
regionCentre-Val de Loire:tranche_0	0.058***
	(0.004)
regionCorse:tranche_0	0.043*
	(0.024)
regionDOM-TOM:tranche_0	0.084***
	(0.006)
regionGrand Est:tranche_0	0.069***
	(0.002)
regionHauts-de-France:tranche_0	0.056***
	(0.001)
regionIle-de-France:tranche_0	0.058***
	(0.001)
regionNormandie:tranche_0	0.070***
	(0.003)
regionNouvelle-Aquitaine:tranche_0	0.061***
	(0.003)
regionOccitanie:tranche_0	0.053***
	(0.002)
regionPays de la Loire:tranche_0	0.048***
	(0.003)
regionProvence-Alpes-Côte d’Azur:tranche_0	0.081***
	(0.001)
Observations	4,242
R2	0.925
Adjusted R2	0.925
Residual Std. Error	33.644 (df = 4214)
F Statistic	1,867.852*** (df = 28; 4214)
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 :

stargazer::stargazer(res_lm_r1, res_lm_r2, res_lm_r3, res_lm_r4, res_lm_r5, res_lm_r6, res_lm_r7, 
                     res_lm_r8, res_lm_r9, res_lm_r10, res_lm_r11, res_lm_r12, res_lm_r13, res_lm_r14, 
                     type = "html", column.labels = nom_region)

	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.063***	0.056***	0.081***	0.069***	0.053***	0.070***	0.061***	0.058***	0.073***	0.058***	0.043***	0.048***	0.058***	0.084***
	(0.001)	(0.001)	(0.002)	(0.002)	(0.001)	(0.002)	(0.001)	(0.002)	(0.003)	(0.002)	(0.005)	(0.002)	(0.001)	(0.004)
Constant	14.475***	28.005***	9.750**	12.535***	3.622***	6.344***	2.975***	9.072***	10.344***	4.124**	2.086**	14.125***	29.832***	3.740*
	(2.045)	(4.608)	(4.713)	(1.817)	(0.758)	(2.184)	(0.907)	(1.498)	(2.100)	(1.881)	(0.983)	(2.133)	(3.998)	(1.996)
Observations	504	210	252	420	546	210	504	252	336	168	84	210	336	210
R2	0.915	0.893	0.891	0.816	0.899	0.852	0.826	0.740	0.714	0.778	0.478	0.734	0.867	0.672
Adjusted R2	0.914	0.892	0.891	0.816	0.899	0.851	0.825	0.739	0.713	0.776	0.472	0.733	0.867	0.670
Residual Std. Error	40.027 (df = 502)	57.055 (df = 208)	61.775 (df = 250)	30.085 (df = 418)	15.209 (df = 544)	25.373 (df = 208)	16.760 (df = 502)	18.308 (df = 250)	30.949 (df = 334)	18.760 (df = 166)	7.078 (df = 82)	23.976 (df = 208)	50.740 (df = 334)	22.976 (df = 208)
F Statistic	5,370.797*** (df = 1; 502)	1,734.497*** (df = 1; 208)	2,050.474*** (df = 1; 250)	1,856.033*** (df = 1; 418)	4,861.813*** (df = 1; 544)	1,194.232*** (df = 1; 208)	2,376.495*** (df = 1; 502)	711.943*** (df = 1; 250)	834.458*** (df = 1; 334)	580.322*** (df = 1; 166)	75.130*** (df = 1; 82)	573.771*** (df = 1; 208)	2,185.952*** (df = 1; 334)	425.609*** (df = 1; 208)
Note:	p<0.1; p<0.05; p<0.01

Test:

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

pred_1 <- predict(res_lm, newdata = my_basis[my_basis$semaine == "semaine_t00", ])

plot(pred_1, my_basis[my_basis$semaine == "semaine_t00", "hosp"],
     xlab = "valeurs prédites", ylab = "valeurs observées")
abline(a = 0, b = 1)
text(pred_1, my_basis[my_basis$semaine == "semaine_t00", "hosp"], 
     my_basis[my_basis$semaine == "semaine_t00", "dep"], pos = 2)

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

mean((pred_1 - my_basis[my_basis$semaine == "semaine_t00", "hosp"]) ^ 2, na.rm = T)

## [1] 1029.62

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:

# débord on aggrège le nombre de testés positifs par région
test_by <- merge(test_by, dep_region, by = "dep")
test_reg <- aggregate(test_by[, c("tranche_9", "tranche_19", "tranche_29", "tranche_39", 
                                   "tranche_49", "tranche_59", "tranche_69", "tranche_79", 
                                   "tranche_89", "tranche_90", "tranche_0")], 
                      by = list(
  semaine = test_by$semaine,
  region = test_by$region), FUN = sum)
# ensuite, on fait le merge
my_basis_age_wide <- merge(my_basis_age_wide, test_reg, by = c("semaine", "region"))

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

apprentissage_sample <- my_basis_age_wide[!(my_basis_age_wide$semaine %in% 
                                              c("semaine_t0-2", "semaine_t0-1", "semaine_t00")), ]
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)

stargazer::stargazer(res_lm_9, res_lm_19, res_lm_29, res_lm_39, res_lm_49, res_lm_59,
                     res_lm_69, res_lm_79, res_lm_89, res_lm_90, type = "html")

	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.013***
	(0.001)
tranche_19		0.004***
		(0.0001)
tranche_29			0.008***
			(0.0002)
tranche_39				0.015***
				(0.0003)
tranche_49					0.023***
					(0.0004)
tranche_59						0.048***
						(0.001)
tranche_69							0.107***
							(0.001)
tranche_79								0.223***
								(0.004)
tranche_89									0.344***
									(0.005)
tranche_90										0.288***
										(0.006)
Constant	0.793***	1.011***	2.653***	3.335***	4.153***	4.739***	8.716***	12.415***	15.427***	10.174***
	(0.227)	(0.194)	(0.355)	(0.488)	(0.625)	(1.102)	(1.579)	(2.442)	(2.850)	(1.778)
Observations	588	588	588	588	588	588	588	588	588	588
R2	0.389	0.532	0.779	0.835	0.871	0.892	0.900	0.867	0.877	0.814
Adjusted R2	0.388	0.531	0.778	0.835	0.871	0.892	0.900	0.867	0.876	0.814
Residual Std. Error (df = 586)	4.503	3.881	7.263	10.087	12.884	22.635	32.286	49.489	56.962	35.610
F Statistic (df = 1; 586)	372.615***	665.835***	2,060.524***	2,962.246***	3,966.228***	4,847.972***	5,296.178***	3,817.554***	4,161.806***	2,568.788***
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.

test_sample <- my_basis[my_basis$semaine == "semaine_t00", ]
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)

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 :

pred_2 <- pred_9 + pred_19 + pred_29 + pred_39 + pred_49 + pred_59 + pred_69 + 
  pred_79 + pred_89 + pred_90

plot(pred_2, my_basis[my_basis$semaine == "semaine_t00", "hosp"],
     xlab = "valeurs prédites", ylab = "valeurs observées")
abline(a = 0, b = 1)
text(pred_2, my_basis[my_basis$semaine == "semaine_t00", "hosp"], 
     my_basis[my_basis$semaine == "semaine_t00", "dep"], pos = 2)

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

mean((pred_2 - my_basis[my_basis$semaine == "semaine_t00", "hosp"]) ^ 2, na.rm = T)

## [1] 3543.97

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.

nom_dep <- my_basis[my_basis$semaine == "semaine_t00", "dep"]
pred_3a <- numeric(length(nom_dep))
pred_3b <- numeric(length(nom_dep))
pred_3c <- numeric(length(nom_dep))
for (k in 1:length(nom_dep)) {
  temp <- hospital[!(hospital$semaine %in% c("semaine_t0-2", "semaine_t0-1", "semaine_t00")) & 
                     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)
  if (nom_dep[k] %in% c("975", "977", "978")) {
    pred_3[k] <- NA
  } else {
    my_mod <- forecast::auto.arima(my_ts_diff)
    forecast_my_mod <- as.numeric(forecast(my_mod)$mean)
    pred_3a[k] <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:7])
    # 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_3b[k] <- 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", "semaine_t0-1", "semaine_t00")), ]
    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[k] <- round(forecast_my_mod_exp[1], 0)
  }
}

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

op <- par(mfrow = c(1, 3), oma = c(0, 0, 0, 0))
plot(pred_3a, my_basis[my_basis$semaine == "semaine_t00", "hosp"],
     xlab = "valeurs prédites", ylab = "valeurs observées", main = "Box-Jenkins")
abline(a = 0, b = 1)
text(pred_3a, my_basis[my_basis$semaine == "semaine_t00", "hosp"], 
     my_basis[my_basis$semaine == "semaine_t00", "dep"], pos = 2)
plot(pred_3b, my_basis[my_basis$semaine == "semaine_t00", "hosp"],
     xlab = "valeurs prédites", ylab = "valeurs observées", main = "Lissage exponentiel journalier")
abline(a = 0, b = 1)
text(pred_3b, my_basis[my_basis$semaine == "semaine_t00", "hosp"], 
     my_basis[my_basis$semaine == "semaine_t00", "dep"], pos = 2)
plot(pred_3c, my_basis[my_basis$semaine == "semaine_t00", "hosp"],
     xlab = "valeurs prédites", ylab = "valeurs observées", main = "Lissage exponentiel hebdomadaire")
abline(a = 0, b = 1)
text(pred_3c, my_basis[my_basis$semaine == "semaine_t00", "hosp"], 
     my_basis[my_basis$semaine == "semaine_t00", "dep"], pos = 2)

par(op)

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

c(mean((pred_3a - my_basis[my_basis$semaine == "semaine_t00", "hosp"]) ^ 2, na.rm = T),
  mean((pred_3b - my_basis[my_basis$semaine == "semaine_t00", "hosp"]) ^ 2, na.rm = T),
  mean((pred_3c - my_basis[my_basis$semaine == "semaine_t00", "hosp"]) ^ 2, na.rm = T)
)

## [1] 1535.1888  784.6733  612.9604

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.

lm_3_ts <- step(lm(my_basis[my_basis$semaine == "semaine_t00", "hosp"] ~ pred_3a + pred_3b + pred_3c - 1))

## Start:  AIC=649.21
## my_basis[my_basis$semaine == "semaine_t00", "hosp"] ~ pred_3a + 
##     pred_3b + pred_3c - 1
## 
##           Df Sum of Sq   RSS    AIC
## - pred_3b  1       8.4 58908 647.23
## - pred_3a  1     707.3 59607 648.42
## <none>                 58899 649.21
## - pred_3c  1   13014.9 71914 667.38
## 
## Step:  AIC=647.23
## my_basis[my_basis$semaine == "semaine_t00", "hosp"] ~ pred_3a + 
##     pred_3c - 1
## 
##           Df Sum of Sq    RSS    AIC
## - pred_3a  1       707  59615 646.43
## <none>                  58908 647.23
## - pred_3c  1     83901 142809 734.67
## 
## Step:  AIC=646.43
## my_basis[my_basis$semaine == "semaine_t00", "hosp"] ~ pred_3c - 
##     1
## 
##           Df Sum of Sq     RSS     AIC
## <none>                   59615  646.43
## - pred_3c  1   2636118 2695733 1029.40

pred_3 <- predict(lm_3_ts)
mean((pred_3 - na.omit(my_basis[my_basis$semaine == "semaine_t00", "hosp"])) ^ 2)

## [1] 590.2445

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.

lm_3 <- lm(my_basis[my_basis$semaine == "semaine_t00", "hosp"] ~ pred_1 + pred_2 + pred_3 - 1)
mean((predict(lm_3) - na.omit(my_basis[my_basis$semaine == "semaine_t00", "hosp"])) ^ 2)

## [1] 479.0464

• 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=1358.18
## y_true ~ pred_3a_s2 + pred_3b_s2 + pred_3c_s2 - 1
## 
##              Df Sum of Sq    RSS    AIC
## <none>                    163121 1358.2
## - pred_3a_s2  1    4371.3 167493 1361.5
## - pred_3b_s2  1    5167.9 168289 1362.5
## - pred_3c_s2  1   17316.5 180438 1376.6

5.5 Prédiction

On prédit le nombre d’hospitalisations :

• du [20 mars 2021; 26 mars 2021] en utilisant les vrais valeurs du nombre de testé positifs la semaine du [10 mars 2021; 16 mars 2021].

• du [27 mars 2021; 02 avril 2021] en utilisant les valeurs prédites du nombre de testé positifs la semaine du [17 mars 2021; 23 mars 2021].

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 [06 mars 2021; 12 mars 2021]
• les 7 derniers jours passés : [13 mars 2021; 19 mars 2021]
• la semaine à venir : [20 mars 2021; 26 mars 2021]
• la semaine suivante à venir : [27 mars 2021; 02 avril 2021]

new_data[, "last_week"] <- my_basis[my_basis$semaine == "semaine_t01", "hosp"]
new_data[, "this_week"] <- my_basis[my_basis$semaine == "semaine_t00", "hosp"]
new_data_long <- tidyr::pivot_longer(data = select(new_data, dep, region, last_week, this_week, next_week, next_two_week),
                                   col = 3:6,
                                   names_to = "semaine",
                                   values_to = "hospitalisations")
new_data_long$semaine <- factor(new_data_long$semaine,
      levels = c("last_week", "this_week", "next_week", "next_two_week"),
      labels = c(paste0("[", format(to_day - 14, '%d %b'), "; ",  
                        format(to_day - 8, '%d %b'), "]"),
                 paste0("[", format(to_day - 7, '%d %b'), "; ",  
                        format(to_day - 1, '%d %b'), "]"),
                 paste0("[", format(to_day, '%d %b'), "; ",  
                        format(to_day + 6, '%d %b'), "]"),
                    paste0("[", format(to_day + 7, '%d %b'), "; ",  
                        format(to_day + 13, '%d %b'), "]"))
      )
# On ajoute aux données de la semaine dernière :
new_data_long$region <- factor(new_data_long$region, levels = hosp_region$region)
p <- ggplot(new_data_long, aes(x = semaine, y = hospitalisations, group = dep))+
    geom_line() +
  facet_wrap(~region)
plotly::ggplotly(p)

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

  my_basis_fr <- my_basis %>%
    group_by(semaine, jour) %>%
    summarise(hosp = sum(hosp),
              rea = sum(rea),
              dc = sum(dc))
  p <- ggplot(data = filter(my_basis_fr, semaine %in% c("semaine_t0-2", "semaine_t0-1", "semaine_t00")),
              aes(x = jour, y = hosp)) +
    geom_line(col = "red") +
    geom_line(data = filter(my_basis_fr, !(semaine %in% c("semaine_t0-2", "semaine_t0-1"))), 
              aes(x = jour, y = hosp)) +
  labs(title = "Prédiction des nouvelles hospitalisations dans les 14 jours",
       x = "semaine",
       y = "hospitalisations",
       fill = "Age") 
plotly::ggplotly(p)

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.

my_basis_temp <- my_basis %>% 
  filter(semaine == "semaine_t0-2")
my_basis_temp$hosp <- NA
my_basis_temp$semaine <- "semaine_t0-3"
my_basis_temp$jour <- my_basis_temp$jour + 7
my_basis_temp[, paste0("tranche_", c(9, 19, 29, 39, 49, 59, 69, 79, 89, 90, 0))] <- NA
my_basis <- rbind(my_basis_temp, my_basis)
my_basis_rea <- my_basis %>%
  select(dep, semaine, jour, region, rea)
temp_hosp <- my_basis %>%
  select(dep, jour, hosp) %>%
  mutate(jour = jour + 7)
my_basis_rea <- merge(my_basis_rea, temp_hosp, by = c("jour", "dep"), all.x = T)
my_basis_rea <- my_basis_rea %>%
  as.data.frame
my_basis_rea <- my_basis_rea[order(my_basis_rea$jour, decreasing = T), ]

p <- ggplot(my_basis_rea) +                
  aes(x = hosp, y = rea) +     
  geom_point() +                  
  geom_smooth(method = "loess") + 
  geom_smooth(method = "lm",     
              col = "red") +
  facet_wrap(~ region)
plotly::ggplotly(p)

On ne va faire que 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.

• un modèle de série temporelles où on explique dans un département les réanimations à l’instant \(t\) par les valeurs du passé.

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:

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

stargazer::stargazer(res_lm_rea_1, type = "html")

	Dependent variable:
	rea
regionAuvergne-Rhône-Alpes	-0.561
	(0.415)
regionBourgogne-Franche-Comté	1.006*
	(0.548)
regionBretagne	0.845
	(0.785)
regionCentre-Val de Loire	-0.039
	(0.683)
regionCorse	0.647
	(1.068)
regionDOM-TOM	1.636**
	(0.667)
regionGrand Est	-0.265
	(0.499)
regionHauts-de-France	-0.389
	(0.670)
regionIle-de-France	-1.299**
	(0.662)
regionNormandie	0.510
	(0.685)
regionNouvelle-Aquitaine	0.046
	(0.426)
regionOccitanie	-0.168
	(0.399)
regionPays de la Loire	0.647
	(0.763)
regionProvence-Alpes-Côte d’Azur	-1.115*
	(0.601)
regionAuvergne-Rhône-Alpes:hosp	0.160***
	(0.003)
regionBourgogne-Franche-Comté:hosp	0.116***
	(0.007)
regionBretagne:hosp	0.108***
	(0.015)
regionCentre-Val de Loire:hosp	0.171***
	(0.014)
regionCorse:hosp	0.154*
	(0.088)
regionDOM-TOM:hosp	0.121***
	(0.014)
regionGrand Est:hosp	0.151***
	(0.005)
regionHauts-de-France:hosp	0.183***
	(0.003)
regionIle-de-France:hosp	0.206***
	(0.003)
regionNormandie:hosp	0.128***
	(0.008)
regionNouvelle-Aquitaine:hosp	0.153***
	(0.009)
regionOccitanie:hosp	0.209***
	(0.007)
regionPays de la Loire:hosp	0.126***
	(0.012)
regionProvence-Alpes-Côte d’Azur:hosp	0.167***
	(0.003)
Observations	4,141
R2	0.881
Adjusted R2	0.881
Residual Std. Error	7.766 (df = 4113)
F Statistic	1,091.925*** (df = 28; 4113)
Note:	p<0.1; p<0.05; p<0.01

Test:

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

pred_rea_1 <- predict(res_lm_rea_1, newdata = my_basis_rea[my_basis_rea$semaine == "semaine_t00", ])

plot(pred_rea_1, my_basis_rea[my_basis_rea$semaine == "semaine_t00", "rea"],
     xlab = "valeurs prédites", ylab = "valeurs observées")
abline(a = 0, b = 1)
text(pred_rea_1, my_basis_rea[my_basis_rea$semaine == "semaine_t00", "rea"], 
     my_basis_rea[my_basis_rea$semaine == "semaine_t00", "dep"], pos = 2)

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

mean((pred_rea_1 - my_basis_rea[my_basis_rea$semaine == "semaine_t00", "rea"]) ^ 2, na.rm = T)

## [1] 224.1271

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.

nom_dep <- my_basis_rea[my_basis_rea$semaine == "semaine_t00", "dep"]
pred_rea_2a <- numeric(length(nom_dep))
pred_rea_2b <- numeric(length(nom_dep))
pred_rea_2c <- numeric(length(nom_dep))
for (k in 1:length(nom_dep)) {
  temp <- hospital[!(hospital$semaine %in% c("semaine_t0-2", "semaine_t0-1", "semaine_t00")) & 
                     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 {
    # box jenkins
    my_mod <- forecast::auto.arima(my_ts_diff)
    forecast_my_mod <- as.numeric(forecast(my_mod)$mean)
    pred_rea_2a[k] <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:7])
    # 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_rea_2b[k] <- 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-3", "semaine_t0-2", "semaine_t0-1", "semaine_t00")), ]
    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_rea_2c[k] <- round(forecast_my_mod_exp[1], 0)
  }
}

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

op <- par(mfrow = c(1, 3), oma = c(0, 0, 0, 0))
plot(pred_rea_2a, my_basis_rea[my_basis_rea$semaine == "semaine_t00", "rea"],
     xlab = "valeurs prédites", ylab = "valeurs observées", main = "Box-Jenkins")
abline(a = 0, b = 1)
text(pred_rea_2a, my_basis_rea[my_basis_rea$semaine == "semaine_t00", "rea"], 
     my_basis_rea[my_basis_rea$semaine == "semaine_t00", "dep"], pos = 2)
plot(pred_rea_2b, my_basis_rea[my_basis_rea$semaine == "semaine_t00", "rea"],
     xlab = "valeurs prédites", ylab = "valeurs observées", main = "Lissage exponentiel journalier")
abline(a = 0, b = 1)
text(pred_rea_2b, my_basis_rea[my_basis_rea$semaine == "semaine_t00", "rea"], 
     my_basis_rea[my_basis_rea$semaine == "semaine_t00", "dep"], pos = 2)
plot(pred_rea_2c, my_basis_rea[my_basis_rea$semaine == "semaine_t00", "rea"],
     xlab = "valeurs prédites", ylab = "valeurs observées", main = "Lissage exponentiel hebdomadaire")
abline(a = 0, b = 1)
text(pred_rea_2c, my_basis_rea[my_basis_rea$semaine == "semaine_t00", "rea"], 
     my_basis_rea[my_basis_rea$semaine == "semaine_t00", "dep"], pos = 2)

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

c(mean((pred_rea_2a - my_basis_rea[my_basis_rea$semaine == "semaine_t00", "rea"]) ^ 2, na.rm = T), 
  mean((pred_rea_2b - my_basis_rea[my_basis_rea$semaine == "semaine_t00", "rea"]) ^ 2, na.rm = T),
  mean((pred_rea_2c - my_basis_rea[my_basis_rea$semaine == "semaine_t00", "rea"]) ^ 2, na.rm = T)
)

## [1]  97.28744 105.29703  66.52475

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:

lm_2_rea_ts <- step(lm(my_basis[my_basis$semaine == "semaine_t00", "rea"] ~ 
                         pred_rea_2a + pred_rea_2b + pred_rea_2c - 1))

## Start:  AIC=427.51
## 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      2.83 6561.1 425.55
## - pred_rea_2b  1     42.22 6600.5 426.16
## <none>                     6558.3 427.51
## - pred_rea_2c  1   1020.61 7578.9 440.12
## 
## Step:  AIC=425.55
## 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      59.3 6620.4 424.46
## <none>                     6561.1 425.55
## - pred_rea_2c  1    1319.4 7880.5 442.06
## 
## Step:  AIC=424.46
## my_basis[my_basis$semaine == "semaine_t00", "rea"] ~ pred_rea_2c - 
##     1
## 
##               Df Sum of Sq    RSS    AIC
## <none>                       6620 424.46
## - pred_rea_2c  1    165686 172306 751.63

pred_rea_2 <- predict(lm_2_rea_ts)
mean((pred_rea_2 - na.omit(my_basis[my_basis$semaine == "semaine_t00", "rea"])) ^ 2)

## [1] 65.54872

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:

lm_rea_3 <- lm(my_basis_rea[my_basis_rea$semaine == "semaine_t00", "rea"] ~ 
                 pred_rea_1 + pred_rea_2 - 1)
mean((predict(lm_rea_3) - na.omit(my_basis_rea[my_basis_rea$semaine == "semaine_t00", "rea"])) ^ 2)

## [1] 50.02158

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 <- pred_rea_1
pred_rea_2a_s2 <- pred_rea_2a
pred_rea_2b_s2 <- pred_rea_2b
pred_rea_2c_s2 <- pred_rea_2c

y_true <- my_basis_rea[my_basis_rea$semaine == "semaine_t00", "rea"]

for (j in 0:0) {
  
  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=913.81
## 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     53.76 18130 912.41
## - pred_rea_2b_s2  1    131.33 18208 913.27
## <none>                        18076 913.81
## - pred_rea_2c_s2  1   1578.28 19655 928.72
## 
## Step:  AIC=912.41
## y_true ~ pred_rea_2b_s2 + pred_rea_2c_s2 - 1
## 
##                  Df Sum of Sq   RSS    AIC
## - pred_rea_2b_s2  1     78.48 18209 911.28
## <none>                        18130 912.41
## - pred_rea_2c_s2  1   1827.51 19958 929.81
## 
## Step:  AIC=911.28
## y_true ~ pred_rea_2c_s2 - 1
## 
##                  Df Sum of Sq    RSS     AIC
## <none>                         18209  911.28
## - pred_rea_2c_s2  1    326403 344612 1503.27

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 <- pred_rea_1
pred_rea_2a_s3 <- pred_rea_2a
pred_rea_2b_s3 <- pred_rea_2b
pred_rea_2c_s3 <- pred_rea_2c

y_true <- my_basis_rea[my_basis_rea$semaine == "semaine_t00", "rea"]

for (j in 0:0) {
  
  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=935.89
## y_true ~ pred_rea_2a_s3 + pred_rea_2b_s3 + pred_rea_2c_s3 - 1
## 
##                  Df Sum of Sq   RSS    AIC
## - pred_rea_2b_s3  1     184.4 20349 935.73
## <none>                        20165 935.89
## - pred_rea_2a_s3  1    1178.3 21343 945.36
## - pred_rea_2c_s3  1    3605.5 23770 967.12
## 
## Step:  AIC=935.73
## y_true ~ pred_rea_2a_s3 + pred_rea_2c_s3 - 1
## 
##                  Df Sum of Sq   RSS    AIC
## <none>                        20349 935.73
## - pred_rea_2a_s3  1    1590.3 21939 948.93
## - pred_rea_2c_s3  1    3669.7 24019 967.22

pred_rea_2_s3 <- predict(lm_2_rea_ts_s2)

lm_rea_3c <- lm(y_true ~ pred_rea_1_s3 + pred_rea_2_s3 - 1)

6.4 Prédiction

On prédit:

• le nombre de réanimations à venir du [20 mars 2021; 26 mars 2021] en utilisant les nouvelles hospitalisations du [13 mars 2021; 19 mars 2021]

• le nombre de réanimations à venir du [27 mars 2021; 02 avril 2021] en utilisant la prédiction des hospitalisations à venir du [20 mars 2021; 26 mars 2021]

• le nombre de réanimations à venir du [03 avril 2021; 09 avril 2021] en utilisant la prédiction des hospitalisations à venir du [27 mars 2021; 02 avril 2021]

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 : [13 mars 2021; 19 mars 2021]
• la semaine à venir : [20 mars 2021; 26 mars 2021]
• la 2ème semaine à venir : [27 mars 2021; 02 avril 2021]
• la 3ème semaine à venir : [03 avril 2021; 09 avril 2021]

new_data_long <- tidyr::pivot_longer(data = select(new_data, dep, region, 
                                   this_week, next_week,  next_two_week, next_three_week),
                                   col = 3:6,
                                   names_to = "semaine",
                                   values_to = "rea")
new_data_long$semaine <- factor(new_data_long$semaine,
      levels = c("this_week", "next_week", "next_two_week", "next_three_week"),
      labels = c(paste0("[", format(to_day - 7, '%d %b'), "; ",  
                        format(to_day - 1, '%d %b'), "]"),
                 paste0("[", format(to_day, '%d %b'), "; ",  
                        format(to_day + 6, '%d %b'), "]"),
                 paste0("[", format(to_day+7, '%d %b'), "; ",  
                        format(to_day + 13, '%d %b'), "]"), 
                 paste0("[", format(to_day + 14, '%d %b'), "; ",  
                        format(to_day + 20, '%d %b'), "]")
                 )
      )

new_data_long$region <- factor(new_data_long$region, levels = hosp_region$region)
p <- ggplot(new_data_long, aes(x = semaine, y = rea,  group = dep))+
    geom_line() +
  facet_wrap(~region)
plotly::ggplotly(p)

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

my_basis_fr <- my_basis %>%
  group_by(semaine, jour) %>%
  summarise(rea = sum(rea))
p <- ggplot(data = filter(my_basis_fr, semaine %in% 
                      c("semaine_t0-3", "semaine_t0-2", "semaine_t0-1", "semaine_t00")),
            aes(x = jour, y = rea)) +
  geom_line(col = "red") +
  geom_line(data = filter(my_basis_fr, !(semaine %in% 
                        c("semaine_t0-3", "semaine_t0-2", "semaine_t0-1"))), 
            aes(x = jour, y = rea)) +
  labs(title = "Prédiction des nouvelles réanimations dans les 21 jours",
       x = "semaine",
       y = "réanimations",
       fill = "Age") 
plotly::ggplotly(p)

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.

my_basis_temp <- my_basis %>% 
  filter(semaine == "semaine_t0-3")
my_basis_temp$rea <- NA
my_basis_temp$semaine <- "semaine_t0-4"
my_basis_temp$jour <- my_basis_temp$jour + 7
my_basis <- rbind(my_basis_temp, my_basis)
my_basis_dc <- my_basis %>%
  select(dep, semaine, jour, region, dc)
temp_rea <- my_basis %>%
  select(dep, jour, rea) %>%
  mutate(jour = jour + 7)
my_basis_dc <- merge(my_basis_dc, temp_rea, by = c("jour", "dep"), all.x = T)
my_basis_dc <- my_basis_dc %>%
  as.data.frame
my_basis_dc <- my_basis_dc[order(my_basis_dc$jour, decreasing = T), ]

p <- ggplot(my_basis_dc) +                
  aes(x = rea, y = dc) +     
  geom_point() +                  
  geom_smooth(method = "loess") + 
  geom_smooth(method = "lm",     
              col = "red") +
  facet_wrap(~ region)
plotly::ggplotly(p)

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:

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

stargazer::stargazer(res_lm_dc_1, type = "html")

	Dependent variable:
	dc
regionAuvergne-Rhône-Alpes	3.926***
	(0.412)
regionBourgogne-Franche-Comté	2.979***
	(0.537)
regionBretagne	1.436*
	(0.795)
regionCentre-Val de Loire	3.494***
	(0.618)
regionCorse	0.351
	(1.071)
regionDOM-TOM	-0.073
	(0.698)
regionGrand Est	2.515***
	(0.483)
regionHauts-de-France	5.501***
	(0.663)
regionIle-de-France	8.419***
	(0.614)
regionNormandie	1.132*
	(0.686)
regionNouvelle-Aquitaine	2.087***
	(0.410)
regionOccitanie	1.650***
	(0.392)
regionPays de la Loire	2.190***
	(0.724)
regionProvence-Alpes-Côte d’Azur	4.233***
	(0.592)
regionAuvergne-Rhône-Alpes:rea	0.933***
	(0.016)
regionBourgogne-Franche-Comté:rea	0.950***
	(0.050)
regionBretagne:rea	1.070***
	(0.118)
regionCentre-Val de Loire:rea	0.537***
	(0.064)
regionCorse:rea	0.474
	(0.354)
regionDOM-TOM:rea	0.454***
	(0.085)
regionGrand Est:rea	1.170***
	(0.034)
regionHauts-de-France:rea	0.868***
	(0.017)
regionIle-de-France:rea	0.598***
	(0.014)
regionNormandie:rea	1.269***
	(0.061)
regionNouvelle-Aquitaine:rea	0.770***
	(0.048)
regionOccitanie:rea	0.623***
	(0.031)
regionPays de la Loire:rea	0.973***
	(0.075)
regionProvence-Alpes-Côte d’Azur:rea	0.886***
	(0.016)
Observations	4,141
R2	0.863
Adjusted R2	0.862
Residual Std. Error	7.952 (df = 4113)
F Statistic	921.948*** (df = 28; 4113)
Note:	p<0.1; p<0.05; p<0.01

Test:

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

pred_dc_1 <- predict(res_lm_dc_1, newdata = my_basis_dc[my_basis_dc$semaine == "semaine_t00", ])

plot(pred_dc_1, my_basis_dc[my_basis_dc$semaine == "semaine_t00", "dc"],
     xlab = "valeurs prédites", ylab = "valeurs observées")
abline(a = 0, b = 1)
text(pred_dc_1, my_basis_dc[my_basis_dc$semaine == "semaine_t00", "dc"], 
     my_basis_dc[my_basis_dc$semaine == "semaine_t00", "dep"], pos = 2)

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

mean((pred_dc_1 - my_basis_dc[my_basis_dc$semaine == "semaine_t00", "dc"]) ^ 2, na.rm = T)

## [1] 159.3435

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.

pred_dc_2a <- numeric(length(nom_dep))
pred_dc_2b <- numeric(length(nom_dep))
pred_dc_2c <- numeric(length(nom_dep))
for (k in 1:length(nom_dep)) {
  temp <- hospital[!(hospital$semaine %in% c("semaine_t0-1", "semaine_t00")) & 
                     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)$mean)
    pred_dc_2a[k] <- sum((as.numeric(my_ts[length(my_ts)]) + cumsum(forecast_my_mod))[1:7])
     # 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_dc_2b[k] <- 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-4", "semaine_t0-3", 
                                            "semaine_t0-2", "semaine_t0-1", "semaine_t00")), ]
    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_dc_2c[k] <- round(forecast_my_mod_exp[1], 0)
    
    
  }
}

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

par(mfrow = c(1, 3))
plot(pred_dc_2a, my_basis_dc[my_basis_dc$semaine == "semaine_t00", "dc"],
     xlab = "valeurs prédites", ylab = "valeurs observées", main = "Box-Jenkins")
abline(a = 0, b = 1)
text(pred_dc_2a, my_basis_dc[my_basis_dc$semaine == "semaine_t00", "dc"], 
     my_basis_dc[my_basis_dc$semaine == "semaine_t00", "dep"], pos = 2)
plot(pred_dc_2b, my_basis_dc[my_basis_dc$semaine == "semaine_t00", "dc"],
     xlab = "valeurs prédites", ylab = "valeurs observées", main = "Lissage exponentiel journalier")
abline(a = 0, b = 1)
text(pred_dc_2b, my_basis_dc[my_basis_dc$semaine == "semaine_t00", "dc"], 
     my_basis_dc[my_basis_dc$semaine == "semaine_t00", "dep"], pos = 2)
plot(pred_dc_2c, my_basis_dc[my_basis_dc$semaine == "semaine_t00", "dc"],
     xlab = "valeurs prédites", ylab = "valeurs observées", main = "Lissage exponentiel hebdomadaire")
abline(a = 0, b = 1)
text(pred_dc_2c, my_basis_dc[my_basis_dc$semaine == "semaine_t00", "dc"], 
     my_basis_dc[my_basis_dc$semaine == "semaine_t00", "dep"], pos = 2)

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

c(
  mean((pred_dc_2a - my_basis_dc[my_basis_dc$semaine == "semaine_t00", "dc"]) ^ 2, na.rm = T),
  mean((pred_dc_2b - my_basis_dc[my_basis_dc$semaine == "semaine_t00", "dc"]) ^ 2, na.rm = T),
  mean((pred_dc_2c - my_basis_dc[my_basis_dc$semaine == "semaine_t00", "dc"]) ^ 2, na.rm = T)
)

## [1] 124.15947  71.25743  46.19802

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:

lm_2_dc_ts <- step(lm(my_basis[my_basis$semaine == "semaine_t00", "dc"] ~ 
                         pred_dc_2a + pred_dc_2b + pred_dc_2c - 1))

## Start:  AIC=343.53
## my_basis[my_basis$semaine == "semaine_t00", "dc"] ~ pred_dc_2a + 
##     pred_dc_2b + pred_dc_2c - 1
## 
##              Df Sum of Sq    RSS    AIC
## <none>                    2855.4 343.53
## - pred_dc_2b  1    131.37 2986.8 346.07
## - pred_dc_2c  1    637.08 3492.5 361.87
## - pred_dc_2a  1    886.35 3741.8 368.83

pred_dc_2 <- predict(lm_2_dc_ts)
mean((pred_dc_2 - na.omit(my_basis[my_basis$semaine == "semaine_t00", "dc"])) ^ 2)

## [1] 28.2715

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 :

lm_dc_3 <- lm(my_basis_dc[my_basis_dc$semaine == "semaine_t00", "dc"] ~ 
                 pred_dc_1 + pred_dc_2 - 1)
mean((predict(lm_dc_3) - na.omit(my_basis_dc[my_basis_dc$semaine == "semaine_t00", "dc"])) ^ 2)

## [1] 27.8706

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 <- pred_dc_1
pred_dc_2a_s2 <- pred_dc_2a
pred_dc_2b_s2 <- pred_dc_2b
pred_dc_2c_s2 <- pred_dc_2c

y_true <- my_basis[my_basis$semaine == "semaine_t00", "dc"]

for (j in 0:0) {
  
  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=726.83
## y_true ~ pred_dc_2a_s2 + pred_dc_2b_s2 + pred_dc_2c_s2 - 1
## 
##                 Df Sum of Sq    RSS    AIC
## <none>                       7163.2 726.83
## - pred_dc_2b_s2  1    551.07 7714.3 739.80
## - pred_dc_2a_s2  1   1259.22 8422.5 757.54
## - pred_dc_2c_s2  1   1889.23 9052.5 772.11

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:0) {
  
  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=1151.96
## y_true ~ pred_dc_2a_s3 + pred_dc_2b_s3 + pred_dc_2c_s3 - 1
## 
##                 Df Sum of Sq   RSS    AIC
## <none>                       13304 1152.0
## - pred_dc_2b_s3  1     478.2 13782 1160.7
## - pred_dc_2a_s3  1    1190.7 14494 1175.9
## - pred_dc_2c_s3  1    3178.9 16483 1214.9

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:0) {
  
  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=1655.72
## y_true ~ pred_dc_2a_s4 + pred_dc_2b_s4 + pred_dc_2c_s4 - 1
## 
##                 Df Sum of Sq   RSS    AIC
## <none>                       23978 1655.7
## - pred_dc_2b_s4  1     177.4 24155 1656.7
## - pred_dc_2a_s4  1     543.4 24521 1662.8
## - pred_dc_2c_s4  1    4890.7 28868 1728.7

pred_dc_2_s4 <- predict(lm_2_dc_ts_s3)

lm_dc_4c <- lm(y_true ~ pred_dc_1_s4 + pred_dc_2_s4 - 1)

7.4 Prédiction

On prédit:

• le nombre de décès à venir du [20 mars 2021; 26 mars 2021] en utilisant les nouvelles réanimations du [13 mars 2021; 19 mars 2021]

• le nombre de décès à venir du [27 mars 2021; 02 avril 2021] en utilisant la prédiction des réanimations à venir du [20 mars 2021; 26 mars 2021]

• le nombre de décès à venir du [03 avril 2021; 09 avril 2021] en utilisant la prédiction des réanimations à venir du [27 mars 2021; 02 avril 2021]

• le nombre de décès à venir du [10 avril 2021; 16 avril 2021] en utilisant la prédiction des réanimations à venir du [27 mars 2021; 02 avril 2021]

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 : [13 mars 2021; 19 mars 2021]
• la semaine à venir : [20 mars 2021; 26 mars 2021]
• la 2ème semaine à venir : [27 mars 2021; 02 avril 2021]
• la 3ème semaine à venir : [03 avril 2021; 09 avril 2021]
• la 4ème semaine à venir : [03 avril 2021; 09 avril 2021]

new_data_long <- tidyr::pivot_longer(data = select(new_data, dep, region, 
                                    this_week, next_week, next_two_week, next_three_week, next_four_week),
                                   col = 3:7,
                                   names_to = "semaine",
                                   values_to = "dc")
new_data_long$semaine <- factor(new_data_long$semaine,
      levels = c("this_week", "next_week", "next_two_week", "next_three_week", "next_four_week"),
      labels = c(paste0("[", format(to_day - 7, '%d %b'), "; ",  
                        format(to_day - 1, '%d %b'), "]"),
                 paste0("[", format(to_day, '%d %b'), "; ",  
                        format(to_day + 6, '%d %b'), "]"),
                 paste0("[", format(to_day + 7, '%d %b'), "; ",  
                        format(to_day + 13, '%d %b'), "]"),
                 paste0("[", format(to_day + 14, '%d %b'), "; ",  
                        format(to_day + 20, '%d %b'), "]"), 
                 paste0("[", format(to_day + 21, '%d %b'), "; ",  
                        format(to_day + 27, '%d %b'), "]")
                 )
      )

new_data_long$region <- factor(new_data_long$region, levels = hosp_region$region)
p <- ggplot(new_data_long, aes(x = semaine, y = dc,  group = dep))+
    geom_line() +
  facet_wrap(~region)
plotly::ggplotly(p)

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

my_basis_fr <- my_basis %>%
  group_by(semaine, jour) %>%
  summarise(dc = sum(dc))
p <- ggplot(data = filter(my_basis_fr, semaine %in% 
                      c("semaine_t0-4", "semaine_t0-3", "semaine_t0-2", "semaine_t0-1", "semaine_t00")),
            aes(x = jour, y = dc)) +
  geom_line(col = "red") +
  geom_line(data = filter(my_basis_fr, !(semaine %in% 
                        c("semaine_t0-4", "semaine_t0-3", "semaine_t0-2", "semaine_t0-1"))), 
            aes(x = jour, y = dc)) +
  labs(title = "Prédiction des nouveaux décès dans les 28 jours",
       x = "semaine",
       y = "décès",
       fill = "Age") 
plotly::ggplotly(p)

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

new_data_long %>%
  group_by(semaine) %>%
  summarize(dc = sum(dc, na.rm = T) / 7)

## # A tibble: 5 x 2
##   semaine                 dc
##   <fct>                <dbl>
## 1 [13 mars; 19 mars]    245.
## 2 [20 mars; 26 mars]      0 
## 3 [27 mars; 02 avril]     0 
## 4 [03 avril; 09 avril]  235.
## 5 [10 avril; 16 avril]  249.