In reconhub/branchr: R Estimation from Cluster Sizes

### fix the random seed
set.seed(1)
library(branchr)
library(ggplot2)
library(tidyr)
library(gridExtra)

Simulating data

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

# 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

# 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)])
  }
}