24 Modélisation des épidémies
24.1 Overview
Il existe un nombre croissant d’outils pour la modélisation des épidémies qui nous permettent de mener des analyses assez complexes avec un effort minimal.Cette section fournira une aperçu sur la façon d’utiliser ces outils pour :
- estimer le nombre de reproduction effectif Rt et les statistiques connexes. telles que le temps de doublement
- produire des projections à court terme de l’incidence future.
Il ne s’agit pas d’un aperçu des méthodologies et des méthodes statistiques qui sous-tendent ces outils. Veuillez donc vous référer à l’onglet Ressources pour des liens vers des documents traitant de ce sujet. Assurez-vous d’avoir une bonne compréhension des les méthodes avant d’utiliser ces outils ; cela vous permettra d’interpréter correctement leurs résultats.
Voici un exemple de l’un des résultats que nous produirons dans cette section.
24.2 Préparation
Nous allons utiliser deux méthodes et packages différents pour l’estimation Rt, à savoir EpiNow et EpiEstim, ainsi que le package projections pour la prévision de l’incidence des cas.
Ce fragment de code montre le chargement des paquets nécessaires aux analyses. Dans ce manuel, nous mettons l’accent sur p_load()
de pacman, qui installe le paquet si nécessaire et le charge pour l’utiliser. Vous pouvez également charger les paquets installés avec library()
de base R. Voir la page sur R basics pour plus d’informations sur les paquets R.
::p_load(
pacman# Importation de fichiers
rio, # Localisation de fichiers
here, # Gestion des données + graphiques ggplot2
tidyverse, # Analyse des réseaux de transmission
epicontacts, # Estimation de Rt
EpiNow2, # Estimation Rt
EpiEstim, # Projections d'incidence
projections, # Traitement des données d'incidence
incidence2, # Fonctions epi utiles
epitrix, # Distributions discrètes des délais
distcrete )
Nous utiliserons la linelist de cas nettoyée pour toutes les analyses de cette section. Si vous voulez suivre, cliquez pour télécharger la linelist “propre” (en tant que fichier .rds). Consultez la page Télécharger le manuel et les données pour télécharger tous les exemples de données utilisés dans ce manuel.
# Importez la liste de cas nettoyée
<- import("linelist_cleaned.rds") linelist
24.3 Estimation de Rt
EpiNow2 vs. EpiEstim
Le taux de reproduction R est une mesure de la transmissibilité d’une maladie, et est défini comme le nombre attendu de cas secondaires par cas infecté. Dans une population totalement sensible, cette valeur représente le nombre de reproduction de base, R0. Cependant, comme le nombre d’individus sensibles dans une population évolue au cours d’une épidémie ou d’une pandémie, et que diverses mesures de réponse sont mises en œuvre, la mesure la plus couramment utilisée de la transmissibilité est le taux de reproduction effectif, Rt ; il est défini défini comme le nombre attendu de cas secondaires par cas infecté à un moment, t.
Le paquet EpiNow2 fournit le cadre le plus sophistiqué pour l’estimation de Rt. Il présente deux avantages essentiels par rapport à l’autre paquet couramment utilisé, EpiEstim :
- Il tient compte des délais de déclaration et peut donc estimer Rt même lorsque les données récentes sont incomplètes.
- Il estime Rt sur les dates d’infection plutôt que sur les dates de début de déclaration, ce qui signifie que l’effet d’une intervention sera immédiatement reflété dans un changement de Rt, plutôt qu’avec un delai.
Cependant, elle présente également deux inconvénients majeurs :
- Elle nécessite la connaissance de la distribution des temps de génération (c’est-à-dire la distribution des délais entre l’infection d’un cas primaire et d’un cas secondaire), la distribution de la période d’incubation (c’est-à-dire la distribution des délais entre l’infection et l’apparition des symptômes) et toute autre distribution de délai pertinente pour vos données (par exemple, si vous avez des dates de déclaration, vous avez besoin de la distribution des délais entre l’apparition des symptômes et la déclaration, ou la période d’incubation). Bien que cela permette une estimation plus précise de Rt, EpiEstim ne requiert que la distribution de l’intervalle sériel (c’est-à-dire la distribution des délais entre l’apparition des symptômes d’un cas primaire et d’un cas secondaire), qui peut être la seule distribution disponible pour vous.
- EpiNow2 est significativement plus lent que EpiEstim, de manière anecdotique par un facteur de 100 à 1000 ! Par exemple, l’estimation de Rt pour l’échantillon de foyers considéré dans cette section prend environ quatre heures (ceci a été exécuté pour un grand d’itérations pour garantir une grande précision et pourrait probablement être réduite si nécessaire) mais il n’en reste pas moins que l’algorithme est lent en général. Cela peut être irréalisable si vous mettez régulièrement à jour votre base de données pour Rt.
Le paquet que vous choisirez d’utiliser dépendra donc des données, du temps et des ressources informatiques dont vous disposez.
EpiNow2
24.3.0.1 Estimation des distributions de retard
Les distributions de retard requises pour exécuter EpiNow2 dépendent des données dont vous disposez. Essentiellement, vous devez être en mesure de décrire le délai entre la date d’infection à la date de l’événement que vous voulez utiliser pour estimer Rt. Si vous utilisez les dates d’apparition, il s’agit simplement de la distribution de la période d’incubation. Si vous utilisez les dates de déclaration, vous avez besoin du délai entre l’infection et la déclaration. Comme il est peu probable que cette distribution soit connue directement, EpiNow2 vous permet d’enchaîner plusieurs distributions de délai ; dans ce cas, le délai entre l’infection et la déclaration est le même.
Comme nous disposons des dates d’apparition des symptômes pour tous nos cas dans la liste d’exemples, nous n’aurons besoin que de la distribution de la période d’incubation pour déterminer le délai d’apparition des symptômes.Nous pouvons soit estimer cette distribution à partir des données ou utiliser les valeurs de la littérature.
Une estimation de la période d’incubation d’Ebola dans la littérature (tirée de cet article) avec une moyenne de 9,1, un écart-type de 7,3 et une valeur maximale de 30, serait spécifiée comme suit :
<- list(
incubation_period_lit mean = log(9.1),
mean_sd = log(0.1),
sd = log(7.3),
sd_sd = log(0.1),
max = 30
)
Notez que EpiNow2 exige que ces distributions de délais soient fournies sur une échelle log d’où l’appel log
autour de chaque valeur (sauf le paramètre max
qui doit être fourni sur une échelle naturelle). Les paramètres mean_sd
et sd_sd
définissent l’écart type des estimations de la moyenne. Comme ceux-ci ne sont pas connus dans ce cas, nous choisissons la valeur assez arbitraire de 0.1.
Dans cette analyse, nous estimons plutôt la distribution de la période d’incubation à partir de la linelist elle-même en utilisant la fonction bootstrapped_dist_fit
, ce qui va une distribution lognormale aux délais observés entre l’infection et l’apparition de la maladie.
## Estimation de la période d'incubation
<- bootstrapped_dist_fit(
incubation_period $date_onset - linelist$date_infection,
linelistdist = "lognormal",
max_value = 100,
bootstraps = 1
)
L’autre distribution dont nous avons besoin est le temps de génération. Comme nous avons des données sur les temps d’infection et les liens de transmission, nous pouvons estimer cette distribution à partir de la liste de liens en calculant le délai entre les temps d’infection des paires infecteur-infecte. Pour ce faire, nous utilisons la fonction pratique get_pairwise
du paquet epicontacts, qui nous permet de calculer les différences par paire des propriétés de la linelist entre les paires de transmission. Nous créons d’abord un objet epicontacts (voir la page Chaînes de transmission pour plus de détails) :
## générer des contacts
<- linelist %>%
contacts transmute(
from = infector,
to = case_id
%>%
) drop_na()
## générer un objet epicontacts
<- make_epicontacts(
epic linelist = linelist,
contacts = contacts,
directed = TRUE
)
Nous ajustons ensuite la différence de temps d’infection entre les paires de transmission, calculée en utilisant get_pairwise
, à une distribution gamma :
## estimation du temps de génération gamma
<- bootstrapped_dist_fit(
generation_time get_pairwise(epic, "date_infection"),
dist = "gamma",
max_value = 20,
bootstraps = 1
)
Exécution de EpiNow2
Maintenant, il ne nous reste plus qu’à calculer l’incidence journalière à partir de la liste linéaire, ce que nous pouvons faire facilement avec les fonctions dplyr group_by()
et n()
. Notez que EpiNow2 exige que les noms des colonnes soient date
et confirm
.
## Obtenir l'incidence à partir des dates d'apparition
<- linelist %>%
cases group_by(date = date_onset) %>%
summarise(confirm = n())
Nous pouvons ensuite estimer Rt en utilisant la fonction epinow
. Quelques remarques sur les entrées :
- Nous pouvons fournir n’importe quel nombre de distributions de délais “enchaînés” à l’argument
delays
. Nous les insérons simplement à côté de l’objetincubation_period
dans la fonctiondelay_opts
. -
return_output
permet de s’assurer que la sortie est retournée dans R et pas seulement un fichier. -
verbose
spécifie que nous voulons une lecture de la progression. -
horizon
indique pour combien de jours nous voulons projeter l’incidence future. - Nous passons des options supplémentaires à l’argument
stan
pour spécifier combien de temps nous voulons exécuter l’inférence pour. L’augmentation desamples
et dechains
vous donnera une estimation plus précise qui caractérisera mieux l’incertitude.
Cependant, l’exécution sera plus longue.
## exécuter epinow
<- epinow(
epinow_res reported_cases = cases,
generation_time = generation_time,
delays = delay_opts(incubation_period),
return_output = TRUE,
verbose = TRUE,
horizon = 21,
stan = stan_opts(samples = 750, chains = 4)
)
Analyser les sorties
Une fois l’exécution du code terminée, nous pouvons tracer un résumé très facilement comme suit. Faites défiler l’image pour voir l’étendue complète.
## Tracer la figure récapitulative
plot(epinow_res)
Nous pouvons également examiner diverses statistiques sommaires :
## tableau récapitulatif
$summary epinow_res
measure estimate
<char> <char>
1: New confirmed cases by infection date 4 (2 -- 6)
2: Expected change in daily cases Unsure
3: Effective reproduction no. 0.88 (0.73 -- 1.1)
4: Rate of growth -0.012 (-0.028 -- 0.0052)
5: Doubling/halving time (days) -60 (130 -- -25)
numeric_estimate
<list>
1: <data.table[1x9]>
2: 0.56
3: <data.table[1x9]>
4: <data.table[1x9]>
5: <data.table[1x9]>
Pour des analyses plus approfondies et des tracés personnalisés, vous pouvez accéder aux estimations quotidiennes résumées via $estimates$summarised
. Nous allons convertir le tableau par défaut data.table
en un tibble
pour faciliter l’utilisation avec dplyr.
## extraire le résumé et le convertir en tibble
<- as_tibble(epinow_res$estimates$summarised)
estimates estimates
A titre d’exemple, faisons un graphique du temps de doublement et de Rt. Nous n’examinerons que les premiers mois de l’épidémie, lorsque Rt est largement supérieur à un, pour éviter de tracer des temps de doublement extrêmement élevés.
Nous utilisons la formule log(2)/taux de croissance
pour calculer le temps de doublement à partir du taux de croissance estimé.
## faire des df larges pour le tracé de la médiane
<- estimates %>%
df_wide filter(
%in% c("growth_rate", "R"),
variable < as.Date("2014-09-01")
date %>%
) ## convertir les taux de croissance en temps de doublement
mutate(
across(
c(median, lower_90:upper_90),
~ case_when(
== "growth_rate" ~ log(2)/.x,
variable TRUE ~ .x
)
),## renommer la variable pour refléter la transformation
variable = replace(variable, variable == "growth_rate", "doubling_time")
)
## créer des df longs pour le tracé des quantiles
<- df_wide %>%
df_long ## ici, nous faisons correspondre les quantiles (par exemple, lower_90 à upper_90)
pivot_longer(
:upper_90,
lower_90names_to = c(".value", "quantile"),
names_pattern = "(.+)_(.+)"
)
## créer un graphique
ggplot() +
geom_ribbon(
data = df_long,
aes(x = date, ymin = lower, ymax = upper, alpha = quantile),
color = NA
+
) geom_line(
data = df_wide,
aes(x = date, y = median)
+
) ## utiliser label_parsed pour permettre l'utilisation d'une étiquette en indice
facet_wrap(
~ variable,
ncol = 1,
scales = "free_y",
labeller = as_labeller(c(R = "R[t]", doubling_time = "Doubling~time"), label_parsed),
strip.position = 'left'
+
) ## définir manuellement la transparence des quantiles
scale_alpha_manual(
values = c(`20` = 0.7, `50` = 0.4, `90` = 0.2),
labels = function(x) paste0(x, "%")
+
) labs(
x = NULL,
y = NULL,
alpha = "Credible\ninterval"
+
) scale_x_date(
date_breaks = "1 month",
date_labels = "%b %d\n%Y"
+
) theme_minimal(base_size = 14) +
theme(
strip.background = element_blank(),
strip.placement = 'outside'
)
EpiEstim
Pour exécuter EpiEstim, nous devons fournir des données sur l’incidence journalière et spécifier l’intervalle sériel (c’est-à-dire la distribution des délais entre l’apparition des symptômes des cas primaires et secondaires).
Les données d’incidence peuvent être fournies à EpiEstim sous la forme d’un vecteur, d’un cadre de données ou d’un objet incidence
provenant du paquetage original incidence. Vous pouvez même faire la distinction entre les importations et les infections acquises localement ; voir la documentation de ?estimate_R
pour plus de détails.
Nous allons créer l’entrée en utilisant incidence2. Voir la page sur Epidemic curves pour plus d’exemples avec le paquet incidence2. Comme il y a eu des mises à jour du paquet incidence2 qui ne correspondent pas complètement à l’entrée attendue de estimateR()
, quelques étapes supplémentaires mineures sont nécessaires. L’objet incidence consiste en un tibble avec des dates et leurs nombres de cas respectifs. Nous utilisons complete()
de tidyr pour nous assurer que toutes les dates sont incluses (même celles sans cas), puis nous rename()
les colonnes pour les aligner avec ce qui est attendu par estimate_R()
dans une étape ultérieure.
## Obtenir l'incidence à partir de la date d'apparition
<- incidence2::incidence(linelist, date_index = "date_onset") %>% # obtient le nombre de cas par jour
cases ::complete(date_index = seq.Date( # s'assurer que toutes les dates sont représentées
tidyrfrom = min(date_index, na.rm = T),
to = max(date_index, na.rm=T),
by = "day"),
fill = list(count = 0)) %>% # convertit les comptes NA en 0
rename(I = count, # renomme aux noms attendus par estimateR
dates = date_index)
Le paquetage fournit plusieurs options pour spécifier l’intervalle sériel, dont les détails sont fournis dans la documentation de ?estimate_R
.
Utiliser des estimations d’intervalles sériels issues de la littérature
En utilisant l’option method = "parametric_si"
, nous pouvons spécifier manuellement la moyenne et l’écart type de l’intervalle sériel dans la littérature ou dans un objet config
créé à l’aide de la fonction make_config
. Nous utilisons une moyenne et un écart-type de 12.0 et 5.2, respectivement, définis dans cet article :
## créer config
<- make_config(
config_lit mean_si = 12.0,
std_si = 5.2
)
Nous pouvons ensuite estimer Rt avec la fonction estimate_R
:
<- cases %>%
cases filter(!is.na(date))
#créer un cadre de données pour la fonction estimate_R()
<- data.frame(dates = seq.Date(from = min(cases$dates),
cases_incidence to = max(cases$dates),
by = 1))
<- left_join(cases_incidence, cases) %>%
cases_incidence select(dates, I) %>%
mutate(I = ifelse(is.na(I), 0, I))
Joining with `by = join_by(dates)`
<- estimate_R(
epiestim_res_lit incid = cases_incidence,
method = "parametric_si",
config = config_lit
)
Default config will estimate R on weekly sliding windows.
To change this change the t_start and t_end arguments.
et tracer un résumé des résultats :
plot(epiestim_res_lit)
Utilisation d’estimations d’intervalles en série à partir des données
Comme nous avons des données sur les dates d’apparition des symptômes et les liens de transmission, nous pouvons également estimer l’intervalle sériel à partir de la liste de liens en calculant le délai entre les dates d’apparition des symptômes des paires infecteur-infecté. Comme nous l’avons fait dans la section EpiNow2 nous allons utiliser la fonction get_pairwise
du paquet epicontacts qui nous permet de calculer les différences par paires des propriétés de la liste de liens entre les paires de transmission. Nous créons d’abord un objet epicontacts (voir la page Chaînes de transmission pour plus de détails) :
## générer des contacts
<- linelist %>%
contacts transmute(
from = infector,
to = case_id
%>%
) drop_na()
## générer un objet epicontacts
<- make_epicontacts(
epic linelist = linelist,
contacts = contacts,
directed = TRUE
)
Nous ajustons ensuite la différence de dates d’apparition entre les paires de transmissions, calculée en utilisant get_pairwise
, à une distribution gamma. Nous utilisons l’outil pratique fit_disc_gamma
du paquet epitrix pour cette procédure d’ajustement, car nous avons besoin d’une distribution discrète.
## Estimation de l'intervalle sériel gamma
<- fit_disc_gamma(get_pairwise(epic, "date_onset")) serial_interval
Nous passons ensuite ces informations à l’objet config
, exécutons EpiEstim et traçons les résultats :
## faire le config
<- make_config(
config_emp mean_si = serial_interval$mu,
std_si = serial_interval$sd
)
## Exécuter epiestim
<- estimate_R(
epiestim_res_emp incid = cases_incidence,
method = "parametric_si",
config = config_emp
)
Default config will estimate R on weekly sliding windows.
To change this change the t_start and t_end arguments.
## tracer les résultats
plot(epiestim_res_emp)
Spécification des fenêtres temporelles d’estimation
Ces options par défaut fournissent une estimation hebdomadaire glissante et peuvent servir d’avertissement si vous estimez Rt trop tôt dans l’épidémie pour une estimation précise.Vous pouvez changer cela en fixant une date de début ultérieure pour l’estimation de Rt, comme indiqué ci-dessous.
Malheureusement, EpiEstim n’offre qu’une façon très maladroite de spécifier ces temps d’estimation, en ce sens que vous devez fournir un vecteur d’entiers __ se référant aux dates de début et de fin de chaque fenêtre temporelle.
## définir un vecteur de dates commençant le 1er juin
<- seq.Date(
start_dates as.Date("2014-06-01"),
max(cases$dates) - 7,
by = 1
%>%
) ## soustraire la date de départ pour la convertir en numérique
`-`(min(cases$dates)) %>%
## convertir en entier
as.integer()
## ajouter six jours pour une fenêtre glissante d'une semaine
<- start_dates + 6
end_dates
## faire la configuration
<- make_config(
config_partial mean_si = 12.0,
std_si = 5.2,
t_start = start_dates,
t_end = end_dates
)
Maintenant, nous réexécutons EpiEstim et nous pouvons voir que les estimations ne commencent qu’à partir de juin :
## exécuter epiestim
<- estimate_R(
epiestim_res_partial incid = cases_incidence,
method = "parametric_si",
config = config_partial
)
## tracer les résultats
plot(epiestim_res_partial)
Analyser les sorties
Les principales sorties sont accessibles via $R
. A titre d’exemple, nous allons créer un graphe de Rt et une mesure de “potentiel de transmission” donnée par le produit de Rt et du nombre de cas signalés ce jour-là ; cela représente le nombre attendu de cas dans la prochaine génération d’infection.
## créer un cadre de données large pour la médiane
<- epiestim_res_lit$R %>%
df_wide rename_all(clean_labels) %>%
rename(
lower_95_r = quantile_0_025_r,
lower_90_r = quantile_0_05_r,
lower_50_r = quantile_0_25_r,
upper_50_r = quantile_0_75_r,
upper_90_r = quantile_0_95_r,
upper_95_r = quantile_0_975_r,
%>%
) mutate(
## extraire la date médiane de t_start et t_end
dates = epiestim_res_emp$dates[round(map2_dbl(t_start, t_end, median))],
var = "R[t]"
%>%
) ## fusionner les données d'incidence quotidienne
left_join(cases, "dates") %>%
## calculer le risque pour toutes les estimations r
mutate(
across(
:upper_95_r,
lower_95_r~ .x*I,
.names = "{str_replace(.col, '_r', '_risk')}"
)%>%
) ## séparer les estimations de r et les estimations de risque
pivot_longer(
contains("median"),
names_to = c(".value", "variable"),
names_pattern = "(.+)_(.+)"
%>%
) ## Assigner des niveaux de facteurs
mutate(variable = factor(variable, c("risk", "r")))
## créer un cadre de données long à partir des quantiles
<- df_wide %>%
df_long select(-variable, -median) %>%
## séparer les estimations de r/risque et les niveaux de quantile
pivot_longer(
contains(c("lower", "upper")),
names_to = c(".value", "quantile", "variable"),
names_pattern = "(.+)_(.+)_(.+)"
%>%
) mutate(variable = factor(variable, c("risk", "r")))
## créer un graphique
ggplot() +
geom_ribbon(
data = df_long,
aes(x = dates, ymin = lower, ymax = upper, alpha = quantile),
color = NA
+
) geom_line(
data = df_wide,
aes(x = dates, y = median),
alpha = 0.2
+
) ## utiliser label_parsed pour permettre l'utilisation d'une étiquette en indice
facet_wrap(
~ variable,
ncol = 1,
scales = "free_y",
labeller = as_labeller(c(r = "R[t]", risk = "Transmission~potential"), label_parsed),
strip.position = 'left'
+
) ## définir manuellement la transparence des quantiles
scale_alpha_manual(
values = c(`50` = 0.7, `90` = 0.4, `95` = 0.2),
labels = function(x) paste0(x, "%")
+
) labs(
x = NULL,
y = NULL,
alpha = "Credible\ninterval"
+
) scale_x_date(
date_breaks = "1 month",
date_labels = "%b %d\n%Y"
+
) theme_minimal(base_size = 14) +
theme(
strip.background = element_blank(),
strip.placement = 'outside'
)
24.4 Projection de l’incidence
EpiNow2
En plus de l’estimation de Rt, EpiNow2 permet également la prévision de Rt et les projections du nombre de cas par l’intégration avec le paquet EpiSoon sous le capot. Tout ce que vous avez à faire est de spécifier l’argument horizon
dans votre appel de fonction epinow
, indiquant le nombre de jours que vous voulez projeter dans le futur ; voir la section EpiNow2 sous la rubrique “Estimation Rt” pour plus de détails sur la façon de mettre en place EpiNow2. Dans cette section, nous allons simplement tracer les sorties de cette analyse, stockées dans le fichier l’objet epinow_res
.
## définir la date minimale pour le tracé
<- as.Date("2015-03-01")
min_date
## extraire les estimations résumées
<- as_tibble(epinow_res$estimates$summarised)
estimates
## extraire les données brutes sur l'incidence des cas
<- as_tibble(epinow_res$estimates$observations) %>%
observations filter(date > min_date)
## extraire les estimations prévisionnelles du nombre de cas
<- estimates %>%
df_wide filter(
== "reported_cases",
variable == "forecast",
type > min_date
date
)
## convertir en un format encore plus long pour le tracé des quantiles
<- df_wide %>%
df_long ## ici nous faisons correspondre les quantiles (par exemple, lower_90 à upper_90)
pivot_longer(
:upper_90,
lower_90names_to = c(".value", "quantile"),
names_pattern = "(.+)_(.+)"
)
## créer un graphique
ggplot() +
geom_histogram(
data = observations,
aes(x = date, y = confirm),
stat = 'identity',
binwidth = 1
+
) geom_ribbon(
data = df_long,
aes(x = date, ymin = lower, ymax = upper, alpha = quantile),
color = NA
+
) geom_line(
data = df_wide,
aes(x = date, y = median)
+
) geom_vline(xintercept = min(df_long$date), linetype = 2) +
## Définir manuellement la transparence des quantiles
scale_alpha_manual(
values = c(`20` = 0.7, `50` = 0.4, `90` = 0.2),
labels = function(x) paste0(x, "%")
+
) labs(
x = NULL,
y = "Daily reported cases",
alpha = "Credible\ninterval"
+
) scale_x_date(
date_breaks = "1 month",
date_labels = "%b %d\n%Y"
+
) theme_minimal(base_size = 14)
projections
Le paquet projections développé par RECON permet de faire très facilement des prévisions d’incidence à court terme, ne nécessitant que la connaissance du nombre de reproduction effectif de reproduction Rt et de l’intervalle de série. Nous verrons ici comment utiliser des estimations d’intervalle sériel de la littérature et comment utiliser nos propres estimations de la liste de diffusion.
Utiliser les estimations d’intervalles sériels de la littérature
projections nécessite une distribution d’intervalle série discrétisée de la classe distcrete
du paquet distcrete. Nous utiliserons une distribution gamma avec une moyenne de 12,0 et un écart-type de 5,2 définie dans cet article. Pour convertir ces valeurs en paramètres de forme et d’échelle requis pour une distribution gamma. nous utiliserons la fonction gamma_mucv2shapescale
du paquet epitrix.
## obtenir les paramètres de forme et d'échelle à partir du mu moyen et du coefficient de
## variation (par exemple, le rapport entre l'écart type et la moyenne).
<- epitrix::gamma_mucv2shapescale(mu = 12.0, cv = 5.2/12)
shapescale
## fabriquer un objet distcrete
<- distcrete::distcrete(
serial_interval_lit name = "gamma",
interval = 1,
shape = shapescale$shape,
scale = shapescale$scale
)
Voici une vérification rapide pour s’assurer que l’intervalle de série est correct. Nous accédons à la densité de la distribution gamma que nous venons de définir par $d
, ce qui revient à appeler dgamma
:
## vérifiez que l'intervalle série est correct
qplot(
x = 0:50, y = serial_interval_lit$d(0:50), geom = "area",
xlab = "Serial interval", ylab = "Density"
)
Warning: `qplot()` was deprecated in ggplot2 3.4.0.
Utilisation des estimations d’intervalles sériels à partir des données
Comme nous avons des données sur les dates d’apparition des symptômes et les liens de transmission, nous pouvons également estimer l’intervalle sériel à partir de la liste de liens en calculant le délai entre les dates d’apparition des symptômes des paires infecteur-infecté. Comme nous l’avons fait dans la section EpiNow2, nous allons utiliser la fonction get_pairwise
du paquet epicontacts qui nous permet de calculer les différences par paires des propriétés de la liste de liens entre les paires de transmission. Nous créons d’abord un objet epicontacts (voir la page Chaînes de transmission pour plus de détails) :
## générer des contacts
<- linelist %>%
contacts transmute(
from = infector,
to = case_id
%>%
) drop_na()
## générer un objet epicontacts
<- make_epicontacts(
epic linelist = linelist,
contacts = contacts,
directed = TRUE
)
Nous ajustons ensuite la différence de dates d’apparition entre les paires de transmissions, calculée avec get_pairwise
, à une distribution gamma. Nous utilisons l’outil pratique fit_disc_gamma
du paquet epitrix pour cette procédure d’ajustement, car nous avons besoin d’une distribution discrète.
## Estimation de l'intervalle sériel gamma
<- fit_disc_gamma(get_pairwise(epic, "date_onset"))
serial_interval
## inspecter l'estimation
c("mu", "sd")] serial_interval[
$mu
[1] 11.51047
$sd
[1] 7.696056
Projection de l’incidence
Pour projeter l’incidence future, nous devons fournir l’incidence historique sous la forme d’un objet incidence
, ainsi qu’un échantillon de valeurs Rt plausibles. Nous générerons ces valeurs en utilisant les estimations Rt générées par EpiEstim dans la section précédente (sous “Estimation de la valeur de Rt”) et stockées dans l’objet epiestim_res_emp
. Dans le code ci-dessous nous extrayons les estimations de la moyenne et de l’écart type de Rt pour la dernière fenêtre de temps de l’épidémie (en utilisant la fonction tail
pour accéder au dernier élément d’un vecteur), et nous simulons 1000 valeurs à partir d’une distribution gamma en utilisant rgamma
. Vous pouvez également fournir votre propre vecteur de valeurs Rt que vous souhaitez utiliser pour les projections.
## créer un objet d'incidence à partir des dates d'apparition des symptômes
<- incidence::incidence(linelist$date_onset) inc
256 missing observations were removed.
## extraire les valeurs r plausibles de l'estimation la plus récente
<- tail(epiestim_res_emp$R$`Mean(R)`, 1)
mean_r <- tail(epiestim_res_emp$R$`Std(R)`, 1)
sd_r <- gamma_mucv2shapescale(mu = mean_r, cv = sd_r/mean_r)
shapescale <- rgamma(1000, shape = shapescale$shape, scale = shapescale$scale)
plausible_r
## vérifier la distribution
qplot(x = plausible_r, geom = "histogram", xlab = expression(R[t]), ylab = "Counts")
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Nous utilisons ensuite la fonction project()
pour effectuer la prévision réelle. Nous spécifions le nombre de jours pour lesquels nous voulons faire une projection via les arguments n_days
, et nous spécifions le nombre de simulations en utilisant les arguments n_sim
.
## faire une projection
<- project(
proj x = inc,
R = plausible_r,
si = serial_interval$distribution,
n_days = 21,
n_sim = 1000
)
Nous pouvons alors facilement tracer l’incidence et les projections en utilisant les fonctions plot()
et add_projections()
. On peut facilement sous-évaluer l’objet incidence
pour ne montrer que les cas les plus récents en utilisant l’opérateur de crochets.
## Tracer l'incidence et les projections
plot(inc[inc$dates > as.Date("2015-03-01")]) %>%
add_projections(proj)
Vous pouvez également extraire facilement les estimations brutes du nombre de cas quotidiens en convertissant la sortie en un cadre de données.
## convertir en cadre de données pour les données brutes
<- as.data.frame(proj)
proj_df proj_df
24.5 Ressources
Voici un article qui décrit la méthodologie mise en œuvre dans EpiEstim. Voici un article décrivant la méthodologie mise en œuvre dans EpiNow. Voici un article décrivant diverses considérations méthodologiques et pratiques pour l’estimation de Rt.