Skip to content
Permalink
Branch: master
Find file Copy path
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
116 lines (90 sloc) 2.67 KB
---
title: "Estimating R"
author: "Pierre Nouvellet & Anne Cori"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
word_document:
fig_width: 7
fig_height: 5
fig_caption: true
highlight: "tango"
reference_docx: word_style.docx
---
```{r}
### fix the random seed
set.seed(1)
library(branchr)
library(ggplot2)
library(tidyr)
library(gridExtra)
```
---
# Simulating data
---
```{r sim_data, eval = FALSE}
R=.99
s=1
rho=1/1e3
t_max=1e3
N_rep <- 1e2
n_imp <- c(300,500)
n_rep_inmp <- length(n_imp)
Res_sim <-list()
for (n in 1:n_rep_inmp){
Res_sim[[paste0('imp_',n_imp[n])]] <- list(true_sizes=matrix(NA,n_imp[n],N_rep),
obs_sizes=matrix(NA,n_imp[n],N_rep))
}
for (n in 1:n_rep_inmp){
print(n)
for (k in 1:N_rep){
# print(k)
# Sim_I <- simulate_negbin(R,n_imp[n],s,rho,t_max,overdisp = .5)
Sim_I <- simulate_poisson(R,n_imp[n],s,rho,t_max)
Res_sim[[paste0('imp_',n_imp[n])]]$true_sizes[,k] <- Sim_I$true_size
Res_sim[[paste0('imp_',n_imp[n])]]$obs_sizes[,k] <- Sim_I$reported_size
}
}
save.image(file = paste0('sim_polio_R',R,'.RData'))
```
```{r plot_data}
# layout(matrix(1:2,1,2))
for (n in 1:n_rep_inmp){
x <- as.data.frame( Res_sim[[paste0('imp_',n_imp[n])]]$true_sizes)
y <- gather(x)
p1 <- ggplot(y,aes(value, fill = key))+ geom_density(alpha = 0.2,show.legend=FALSE) + xlim(c(0,200))
x <- as.data.frame( Res_sim[[paste0('imp_',n_imp[n])]]$obs_sizes)
y <- gather(x)
p2<- ggplot(y,aes(value, fill = key)) +
geom_density(alpha = 0.2,show.legend=FALSE) +
xlim(c(0,10)) +
ggtitle(paste0('E[# outbreaks obs]=',mean(colSums(x>0)),' (',n_imp[n],' importations)'))
grid.arrange(p1,p2,ncol=2)
}
```
# Inference
## first infering the reproduction number from observed outbreak final sizes
```{r inferR}
# layout(matrix(1:2,1,2))
for (n in 1:n_rep_inmp){
Y_obs <- Res_sim[[paste0('imp_',n_imp[n])]]$obs_sizes
for (k in 1:N_rep){
y_obs <- Y_obs[,k]
y_obs <- y_obs[which(y_obs>0)]
profile <- profile_likelihood(y_obs = y_obs,
rho = rho,
accuracy = 0.01,
max_R = 20)
R_estimate <- theta_max_likelihood(theta = profile$theta,
likelihood = profile$Likelihood,
threshold_CI = 0.95)
as.numeric(R_estimate[c(1,3,4)])
import <- import(y_obs = y_obs,
rho = rho,
profile = profile,
threshold_z = 1e4,
threshold_import = 5e3,
CI = 0.95)
length(y_obs) + as.numeric(import[c(1,3,4)])
}
}
```
You can’t perform that action at this time.