inst/tree_sampling/tree-sample-test-many-roots.R

In skgallagher/TBornotTB: Analyze Tuberculosis Patients in Maryland

## SKG
## 3/19/2020
## Calculating likelhood  and making trees
## Now with many root nodes
##


devtools::load_all()

## Generate partition list
#part_list <- generate_part_list(n = 25)
#g_weight_list <- get_weight_list(part_list)


K <- 1000
inf_params <- c("beta_0" = -2, "beta_1" = 1.5)
max_size <- 50
smear_pos_prob <- .5
start_at_outsider <- FALSE

forest <- simulate_outbreak(K = K,
                                    inf_params = inf_params,
                                    smear_pos_prob = smear_pos_prob,
                                    max_size = max_size,
                                    start_at_outsider = start_at_outsider)

bark <- summarize_clusters(forest)
sum(bark$freq)
print(bark, n = 200)


n_trials <- 1000
cluster_summary <- bark
t <- proc.time()[3]
sampled_trees <- sample_uniform_trees(n_vec = cluster_summary$n,
                                      x_vec = cluster_summary$n_pos,
                                      B = n_trials,
                                      use_one_index = !start_at_outsider)
print(proc.time()[3] - t)
par <- c("beta_0" = -3, "beta_1" = 2)


best_pars <- optim(par, fn = loglike_fast,
                   data = cluster_summary,
                   sampled_trees_summary = sampled_trees,
                   return_neg = TRUE,
                   hessian = TRUE)
best_pars$par
sqrt(diag(solve(best_pars$hessian)))


##############################

bad <- rep(0, nrow(cluster_summary))
for(ii in 1:nrow(cluster_summary)){
  x <- cluster_summary$n[ii]
  y <- cluster_summary$n_pos[ii]
  ################################
  out <- max_g %>% dplyr::filter(cluster_size == x, cluster_pos == y)

  max_g2 <- sampled_trees %>% dplyr::group_by(cluster_size, cluster_pos, cluster_id) %>%
    dplyr::summarize(max_g = max(gen))

  out2 <- max_g2 %>% dplyr::filter(cluster_size == x, cluster_pos == y)

  #par(mfrow = c(1,2))
  #hist(out2$max_g)
  #hist(out$max_g)
  dog <- ks.test(out$max_g, out2$max_g)
  print(dog)
  if(dog$p.value < .05){
    bad[ii] <- 1
  }
}
## TODO: Need to sample trees of given n, n_pos and compare distributions of i_pos, g, etc.










fast2$par
sqrt(diag(solve(fast2$hessian)))