R/parameters.R
In Geidelberg/cotonou:
Defines functions
generate_parameters
lhs_parameters_parallel
delete_parameter_sets
lhs_parameters
fix_parameters
FormerFSWtoGPF
Documented in
fix_parameters
#' @export
#' @useDynLib cotonou
FormerFSWtoGPF <- function(x) {
if(!is.null(dim(x))) {
x[4,] = x[3,]
x[,4] = x[,3]
}
return(x)
}
#' parameters which depend on others, etc
#'
#' Lily's function to be documented
#' @export
#' @useDynLib cotonou
fix_parameters <- function(y, Ncat, Nage, par_seq, condom_seq, groups_seq, years_seq) {
useless_parameter = 1
# why not working
if(y$ignore_ranges_fc_c == 0) {
# y$PrEP_reinit_OnOff_y = c(0, 0, y$PrEP_reinits_on, y$PrEP_reinits_on)
if(y$prep_dropout - y$rate_leave_pro_FSW - y$muF - 1/45 - 0.008 > 0)
{
y$PrEP_loss_to_follow_up <- y$prep_dropout - y$rate_leave_pro_FSW - y$muF - 1/45 - 0.008
} else {
y$PrEP_loss_to_follow_up = 0
}
y$frac_women_exFSW = y$frac_women_ProFSW
# getting client - low FSW n to be same as for pro fsw
# y$n_y_comm_1985_LowFSW_Client = y$n_y_comm_1985_ProFSW_Client
y$n_y_noncomm_2002_LowFSW_Client = y$n_y_noncomm_2002_ProFSW_Client
y$n_y_noncomm_2015_LowFSW_Client = y$n_y_noncomm_2015_ProFSW_Client
y$n_y_noncomm_1985_GPF_Client = y$n_y_noncomm_1985_GPF_GPM
y$n_y_noncomm_1998_GPF_Client = y$n_y_noncomm_1998_GPF_GPM
y$n_y_noncomm_2011_GPF_Client = y$n_y_noncomm_2011_GPF_GPM
# getting low fsw condom with client for comm and noncomm partnerships same as pro fsw
# y$fc_y_comm_1985_LowFSW_Client = y$fc_y_comm_1985_ProFSW_Client
# y$fc_y_comm_1993_LowFSW_Client = y$fc_y_comm_1993_ProFSW_Client
# y$fc_y_comm_1998_LowFSW_Client = y$fc_y_comm_1998_ProFSW_Client
# y$fc_y_comm_2002_LowFSW_Client = y$fc_y_comm_2002_ProFSW_Client
# y$fc_y_noncomm_1985_LowFSW_Client = y$fc_y_noncomm_1985_ProFSW_Client
# y$fc_y_noncomm_2002_LowFSW_Client = y$fc_y_noncomm_2002_ProFSW_Client
# getting client condom with GPF same as GPF-GPM
y$fc_y_noncomm_1985_GPF_Client = y$fc_y_noncomm_1985_GPF_GPM
y$fc_y_noncomm_1998_GPF_Client = y$fc_y_noncomm_1998_GPF_GPM
y$fc_y_noncomm_2011_GPF_Client = y$fc_y_noncomm_2011_GPF_GPM
# CONDOMS
what_we_got_condom = c()
for(year in 1:length(years_seq))
{for(group in 1:length(groups_seq))
{for(group2 in 1:length(groups_seq))
{for(par in 1:length(condom_seq))
{if(paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group], "_", groups_seq[group2]) %in% names(y)){
# average if the opposite value is there too - note must be same year...
if(paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group2], "_", groups_seq[group]) %in% names(y)) {
av=(y[[which(paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group], "_", groups_seq[group2]) == names(y))]] +
y[[which(paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group2], "_", groups_seq[group]) == names(y))]]) / 2;
y[[which(paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group], "_", groups_seq[group2]) == names(y))]] = av;
y[[which(paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group2], "_", groups_seq[group]) == names(y))]] = av;
}
# replacing the non-group specific parameter to include these precisions
y[paste0(condom_seq[par], "_", years_seq[year])][[paste0(condom_seq[par], "_", years_seq[year])]][group, group2] = y[paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group], "_", groups_seq[group2])][[1]]
y[paste0(condom_seq[par], "_", years_seq[year])][[paste0(condom_seq[par], "_", years_seq[year])]][group2, group] = y[paste0(condom_seq[par], "_", years_seq[year])][[paste0(condom_seq[par], "_", years_seq[year])]][group, group2]
what_we_got_condom = rbind(what_we_got_condom, c(par, year, group, group2))
}}}}}
if (length(what_we_got_condom[,1]) > 0)
{
colnames(what_we_got_condom) = c("par", "year", "group", "group2")
# now we have to fill in the rest of the years... IF 2 OR MORE YEARS FOR SAME GROUP AND PARM
# all years before earliest one is same as earliest estimate,
# and all years after last one is same as last...
# and everything in between is interpolated!
d <- data.frame(what_we_got_condom)
par_counts = plyr::count(d, c('par', 'group', 'group2'))
for(i in 1:length(par_counts[,1]))
{
if(par_counts[i,"freq"] > 0)
{
the_years = subset(d, c(par == par_counts[i, "par"] & group == par_counts[i, "group"] & group2 == par_counts[i, "group2"]), year)
#before
if(min(the_years) > 1)
{
for(year in 1:(min(the_years)-1))
{
if(paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year]) %in% names(y))
{y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group"], par_counts[i, "group2"]] =
y[paste0(condom_seq[par_counts[i,"par"]], "_", years_seq[min(the_years)], "_", groups_seq[par_counts[i, "group"]], "_", groups_seq[par_counts[i, "group2"]])][[1]]
# equalising the complementary condom use
y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group2"], par_counts[i, "group"]] =
y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group"], par_counts[i, "group2"]]
}
}
}
#after
if(max(the_years) < length(years_seq))
{
for(year in (max(the_years)+1):length(years_seq))
{
if(paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year]) %in% names(y))
{y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group"], par_counts[i, "group2"]] =
y[paste0(condom_seq[par_counts[i,"par"]], "_", years_seq[max(the_years)], "_", groups_seq[par_counts[i, "group"]], "_", groups_seq[par_counts[i, "group2"]])][[1]]
# equalising the complementary condom use
y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group2"], par_counts[i, "group"]] =
y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group"], par_counts[i, "group2"]]
}
}
}
# interpolate between
# take the years that have values
# calculate the slopes between them
# apply the slope and time difference to the pars in between
if(par_counts[i,"freq"] > 1) {
for(j in 1:(length(unlist(the_years))-1)) {
slope = (y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j+1])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j+1])]][par_counts[i, "group"], par_counts[i, "group2"]] -
y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])]][par_counts[i, "group"], par_counts[i, "group2"]]) /
(years_seq[unlist(the_years)][j+1] - years_seq[unlist(the_years)][j])
for(k in (years_seq[unlist(the_years)][j]+1):(years_seq[unlist(the_years)][j+1]-1))
{
if(paste0(condom_seq[par_counts[i, "par"]], "_", k) %in% names(y))
{ y[paste0(condom_seq[par_counts[i, "par"]], "_", k)][[paste0(condom_seq[par_counts[i, "par"]], "_", k)]][par_counts[i, "group"], par_counts[i, "group2"]] = slope *
(k - years_seq[unlist(the_years)][j]) +
y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])]][par_counts[i, "group"], par_counts[i, "group2"]]
# equalising complementary condom rate
y[paste0(condom_seq[par_counts[i, "par"]], "_", k)][[paste0(condom_seq[par_counts[i, "par"]], "_", k)]][par_counts[i, "group2"], par_counts[i, "group"]] =
y[paste0(condom_seq[par_counts[i, "par"]], "_", k)][[paste0(condom_seq[par_counts[i, "par"]], "_", k)]][par_counts[i, "group"], par_counts[i, "group2"]]}
}
}
}
}
}
}
# MAKING FORMER FSW IN COTONOU EQUAL CONDOM USE TO GPF _ NOTE DO I NEED TO DO THIS FOR ALL GPF PARAMETERS???
y[names(y)[grep("fc_y_", names(y))]] = lapply(invisible(y[names(y)[grep("fc_y_", names(y))]]), FormerFSWtoGPF)
y[names(y)[grep("n_y_", names(y))]] = lapply(invisible(y[names(y)[grep("n_y_", names(y))]]), FormerFSWtoGPF)
##########################################################################
# c_comm and c_noncomm
# what_we_got = data.frame(par = c(), year = c(), group = c())
what_we_got = c()
for(year in 1:length(years_seq))
{for(group in 1:length(groups_seq))
{for(par in 1:length(par_seq))
# {if(grepl(years_seq[year], names(y)) && grepl(groups_seq[group], names(y)) && grepl(par_seq[par], names(y))){
{if(paste0(paste0(par_seq[par], "_", years_seq[year], "_", groups_seq[group])) %in% names(y)){
y[paste0(par_seq[par], "_", years_seq[year])][[paste0(par_seq[par], "_", years_seq[year])]][group] = y[paste0(par_seq[par], "_", years_seq[year], "_", groups_seq[group])][[1]]
# what_we_got = c(what_we_got, paste0(par_seq[par], "_", years_seq[year], "_", groups_seq[group]))
what_we_got = rbind(what_we_got, c(par, year, group))
# print(c(par_seq[par], groups_seq[group], years_seq[year]))
}}}}
# print(list(what_we_got, years_seq, groups_seq, par_seq, what_we_got_condom))
if(length(what_we_got) > 0)
{colnames(what_we_got) = c("par", "year", "group")
# now we have to fill in the rest of the years... IF 2 OR MORE YEARS FOR SAME GROUP AND PARM
# all years before earliest one is same as earliest estimate,
# and all years after last one is same as last...
# and everything in between is interpolated!
d <- data.frame(what_we_got)
par_counts = plyr::count(d, c('par', 'group'))
for(i in 1:length(par_counts[,1]))
{
if(par_counts[i,"freq"] > 0)
{
the_years = subset(d, c(par == par_counts[i, "par"] & group == par_counts[i, "group"]), year)
#before
if(min(the_years) > 1)
{
for(year in 1:(min(the_years)-1))
{
if(paste0(par_seq[par_counts[i, "par"]], "_", years_seq[year]) %in% names(y))
{y[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group"]] =
y[paste0(par_seq[par_counts[i,"par"]], "_", years_seq[min(the_years)], "_", groups_seq[par_counts[i, "group"]])][[1]]}
}
}
#after
if(max(the_years) < length(years_seq))
{
for(year in (max(the_years)+1):length(years_seq))
{
if(paste0(par_seq[par_counts[i, "par"]], "_", years_seq[year]) %in% names(y))
{y[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group"]] =
y[paste0(par_seq[par_counts[i,"par"]], "_", years_seq[max(the_years)], "_", groups_seq[par_counts[i, "group"]])][[1]]}
}
}
# interoplate between
# take the years that have values
# calculate the slopes between them
# apply the slope and time difference to the pars in between
if(par_counts[i,"freq"] > 1) {
for(j in 1:(length(unlist(the_years))-1)) {
slope = (y[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j+1])][[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j+1])]][par_counts[i, "group"]] -
y[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])][[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])]][par_counts[i, "group"]]) /
(years_seq[unlist(the_years)][j+1] - years_seq[unlist(the_years)][j])
for(k in (years_seq[unlist(the_years)][j]+1):(years_seq[unlist(the_years)][j+1]-1))
{
if(paste0(par_seq[par_counts[i, "par"]], "_", k) %in% names(y))
y[paste0(par_seq[par_counts[i, "par"]], "_", k)][[paste0(par_seq[par_counts[i, "par"]], "_", k)]][par_counts[i, "group"]] = slope *
(k - years_seq[unlist(the_years)][j]) +
y[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])][[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])]][par_counts[i, "group"]]
}
}
}
}
}
}
}
# print("after")
# print(y)
# number of clients at initial conditions shaped by commercial partner change rates and the N of Pro FSW!
# partnerships offered by ProFSW / P.C.R of clients / fraction of men * fraction of men that are clients
if(y$init_clientN_from_PCR == 1)
y$frac_men_client = y$fraction_F * y$frac_women_ProFSW * y$initial_Ntot * y$c_comm_1985[1] / y$c_comm_1985[5] / (y$initial_Ntot * (1-y$fraction_F))
# fractions of groups at initialisation
y$N_init = y$initial_Ntot*c(y$fraction_F*y$frac_women_ProFSW,
y$fraction_F*y$frac_women_ProFSW*y$frac_women_LowFSW,
y$fraction_F - (y$fraction_F*y$frac_women_ProFSW + y$fraction_F*y$frac_women_ProFSW*y$frac_women_LowFSW + y$fraction_F*y$frac_women_virgin + y$fraction_F*y$frac_women_exFSW),
y$fraction_F*y$frac_women_exFSW,
(1-y$fraction_F)*y$frac_men_client,
(1 - y$fraction_F - ((1-y$fraction_F)*y$frac_men_client + (1-y$fraction_F) * y$frac_men_virgin)),
y$fraction_F*y$frac_women_virgin,
(1-y$fraction_F) * y$frac_men_virgin,
0
)
y$betaMtoF_FSW = y$betaMtoF_baseline * (1 + (y$RR_beta_HSV2_a_FSW - 1) * y$prev_HSV2_FSW)
y$betaFtoM_client = y$betaMtoF_baseline * y$RR_beta_FtM * y$RR_beta_circum * (1 + (y$RR_beta_HSV2_a_client - 1) * y$prev_HSV2_Client)
y$betaMtoF_GPF = y$betaMtoF_baseline * (1 + (y$RR_beta_HSV2_a_GPF - 1) * y$prev_HSV2_GPF)
y$betaFtoM_GPM = y$betaMtoF_baseline * y$RR_beta_FtM * y$RR_beta_circum * (1 + (y$RR_beta_HSV2_a_GPM - 1) * y$prev_HSV2_GPM)
# y$betaMtoF_noncomm_FSW = y$betaMtoF_baseline * (1 + (y$RR_beta_HSV2_noncomm_a - 1) * y$prev_HSV2_FSW + (y$RR_beta_HSV2_noncomm_t - 1) * y$prev_HSV2_Client)
# y$betaFtoM_noncomm_client = y$betaMtoF_baseline * y$RR_beta_FtM * y$RR_beta_circum * (1 + (y$RR_beta_HSV2_noncomm_a - 1) * y$prev_HSV2_Client + (y$RR_beta_HSV2_noncomm_t - 1) * y$prev_HSV2_FSW)
#
# y$betaMtoF_noncomm = y$betaMtoF_baseline * (1 + (y$RR_beta_HSV2_noncomm_a - 1) * y$prev_HSV2_GPF + (y$RR_beta_HSV2_noncomm_t - 1) * y$prev_HSV2_GPM)
# y$betaFtoM_noncomm = y$betaMtoF_baseline * y$RR_beta_FtM * y$RR_beta_circum * (1 + (y$RR_beta_HSV2_noncomm_a - 1) * y$prev_HSV2_GPM + (y$RR_beta_HSV2_noncomm_t - 1) * y$prev_HSV2_GPF)
#
# y$betaMtoF_comm = y$betaMtoF_baseline * (1 + (y$RR_beta_HSV2_comm_a - 1) * y$prev_HSV2_FSW + (y$RR_beta_HSV2_comm_t - 1) * y$prev_HSV2_Client)
# y$betaFtoM_comm = y$betaMtoF_baseline * y$RR_beta_FtM * y$RR_beta_circum * (1 + (y$RR_beta_HSV2_comm_a - 1) * y$prev_HSV2_Client + (y$RR_beta_HSV2_comm_t - 1) * y$prev_HSV2_FSW)
# if any beta becomes > 1, then make them all zero and flag it
# if(y$betaMtoF_noncomm_FSW * y$infect_acute >= 1 || y$betaFtoM_noncomm_client * y$infect_acute >= 1 || y$betaMtoF_noncomm * y$infect_acute >= 1 || y$betaMtoF_comm * y$infect_acute >= 1 || y$betaFtoM_noncomm * y$infect_acute >= 1 || y$betaFtoM_comm * y$infect_acute >= 1)
if(y$betaMtoF_FSW * y$infect_acute >= 1 || y$betaFtoM_client * y$infect_acute >= 1 || y$betaMtoF_GPF * y$infect_acute >= 1 || y$betaFtoM_GPM * y$infect_acute >= 1)
{
# y$betaMtoF_noncomm_FSW = 0
# y$betaFtoM_noncomm_client = 0
# y$betaMtoF_noncomm = 0
# y$betaFtoM_noncomm = 0
# y$betaMtoF_comm = 0
# y$betaFtoM_comm = 0
y$betaMtoF_FSW = 0
y$betaFtoM_client = 0
y$betaMtoF_GPF = 0
y$betaFtoM_GPM = 0
y$beta_above_1 = 1
}
y$epsilon_y = c(y$epsilon_1985, y$epsilon_1992, y$epsilon_2002, y$epsilon_2013, y$epsilon_2016)
init_prev = if(Ncat == 9) c(y$prev_init_FSW, rep_len(y$prev_init_rest, 5), 0, 0, 0) else c(y$prev_init_FSW, rep_len(y$prev_init_rest, Ncat-1))
y$S0_init = (1-init_prev) * y$N_init
y$I01_init = init_prev * y$N_init
N = y$S0_init + y$I01_init
# BIOLOGICAL
#dur_primary_phase, dur_200_349 input as a DURATION
#SC_to_death input as a duration
y$ART_RR_mort = y$ART_RR_prog
y$ART_RR = y$ART_RR_prog
y$gamma01 = 1/y$dur_primary_phase
y$gamma02 = 2/(y$SC_to_death - y$dur_primary_phase - y$dur_200_349 - y$dur_below_200)
y$gamma03 = y$gamma02
y$gamma04 = 1/y$dur_200_349
y$gamma01 <- rep_len(y$gamma01, Ncat)
y$gamma02 <- rep_len(y$gamma02, Ncat)
y$gamma03 <- rep_len(y$gamma03, Ncat)
y$gamma04 <- rep_len(y$gamma04, Ncat)
y$gamma11 <- y$gamma01 # loss to follow up same as normal movement?
y$gamma22 <- y$gamma02
y$gamma23 <- y$gamma03
y$gamma24 <- y$gamma04
y$gamma42 <- y$gamma02
y$gamma43 <- y$gamma03
y$gamma44 <- y$gamma04
# progression is slowed by ART_RR
y$gamma32_without_supp <- (y$gamma02)/y$ART_RR_prog
y$gamma33_without_supp <- (y$gamma03)/y$ART_RR_prog
y$gamma34_without_supp <- (y$gamma04)/y$ART_RR_prog
# mortality
y$alpha03 <- rep_len(y$alpha03, Ncat)
y$alpha04 <- rep_len(y$alpha04, Ncat)
y$alpha05 <- 1/y$dur_below_200
y$alpha05 <- rep_len(y$alpha05, Ncat)
y$alpha23 <- y$alpha03
y$alpha24 <- y$alpha04
y$alpha25 <- y$alpha05
y$alpha43 <- y$alpha03
y$alpha44 <- y$alpha04
y$alpha45 <- y$alpha05
y$alpha33_without_supp <- (y$alpha03)/y$ART_RR_mort
y$alpha34_without_supp <- (y$alpha04)/y$ART_RR_mort
y$alpha35_without_supp <- (y$alpha05)/y$ART_RR_mort
# PCR
y$c_t_comm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016)
y$c_t_noncomm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016)
y$c_y_comm = rbind(y$c_comm_1985, y$c_comm_1993, y$c_comm_1995, y$c_comm_1998, y$c_comm_2002,
y$c_comm_2005, y$c_comm_2008, y$c_comm_2012, y$c_comm_2015, y$c_comm_2016)
y$c_y_noncomm = rbind(y$c_noncomm_1985, y$c_noncomm_1993, y$c_noncomm_1995, y$c_noncomm_1998, y$c_noncomm_2002,
y$c_noncomm_2005, y$c_noncomm_2008, y$c_noncomm_2012, y$c_noncomm_2015, y$c_noncomm_2016)
y$fc_y_comm = array(data = c(y$fc_y_comm_1985, y$fc_y_comm_1993,
y$fc_y_comm_1995, y$fc_y_comm_1998,
y$fc_y_comm_2002, y$fc_y_comm_2005,
y$fc_y_comm_2008, y$fc_y_comm_2012,
y$fc_y_comm_2015, y$fc_y_comm_2016), dim=c(Ncat, Ncat, 10))
y$fc_y_comm = aperm(y$fc_y_comm, c(3, 1, 2))
y$fc_y_noncomm = array(data = c(y$fc_y_noncomm_1985, y$fc_y_noncomm_1993, y$fc_y_noncomm_1998, y$fc_y_noncomm_2002,
y$fc_y_noncomm_2008, y$fc_y_noncomm_2011, y$fc_y_noncomm_2015,
y$fc_y_noncomm_2016), dim=c(Ncat, Ncat, 8))
y$fc_y_noncomm = aperm(y$fc_y_noncomm, c(3, 1, 2))
y$n_y_comm = array(data = c(y$n_y_comm_1985, y$n_y_comm_2002,
y$n_y_comm_2015, y$n_y_comm_2016), dim=c(Ncat, Ncat, 4))
y$n_y_noncomm = array(data = c(y$n_y_noncomm_1985, y$n_y_noncomm_1998, y$n_y_noncomm_2002, y$n_y_noncomm_2011,
y$n_y_noncomm_2015, y$n_y_noncomm_2016), dim=c(Ncat, Ncat, 6))
y$n_y_comm = aperm(y$n_y_comm, c(3, 1, 2))
y$n_y_noncomm = aperm(y$n_y_noncomm, c(3, 1, 2))
y$omega = N/sum(N)
# when sampling pars, sometimes instead of length Ncat, I only sampled length one
if(length(y$ec) == 1)
y$ec = rep_len(y$ec, 9)
if(length(y$eP0) == 1)
y$eP0 = c(y$eP0, rep_len(0, 8))
if(length(y$eP1a) == 1)
y$eP1a = c(y$eP1a, rep_len(0, 8))
if(length(y$eP1b) == 1)
y$eP1b = c(y$eP1b, rep_len(0, 8))
if(length(y$eP1c) == 1)
y$eP1c = c(y$eP1c, rep_len(0, 8))
if(length(y$eP1d) == 1)
y$eP1d = c(y$eP1d, rep_len(0, 8))
if(length(y$psia) == 1)
y$psia = c(y$psia, rep_len(0, 8))
if(length(y$psib) == 1)
y$psib = c(y$psib, rep_len(0, 8))
y$phi2 = c(y$dropout_rate_FSW, rep(y$dropout_rate_not_FSW, (Ncat-1)))
y$phi3 = y$phi2
y$phi4 = y$phi2
y$phi5 = y$phi2
if (Ncat == 9) {
y$who_believe_comm = round(y$who_believe_comm)
y$omega = c(
y$fraction_F * y$frac_women_ProFSW * y$fraction_FSW_foreign, # some FSW come from outside Cotonou
y$fraction_F * y$frac_women_LowFSW * y$frac_women_ProFSW * y$fraction_FSW_foreign,
y$fraction_F - y$fraction_F * (1-y$fraction_sexually_active_15_F) - y$fraction_F * y$frac_women_ProFSW * y$fraction_FSW_foreign - y$fraction_F * y$frac_women_LowFSW * y$frac_women_ProFSW * y$fraction_FSW_foreign,
0,
0,
(1 - y$fraction_F) - (1 - y$fraction_F) * (1 - y$fraction_sexually_active_15_M),
y$fraction_F * (1-y$fraction_sexually_active_15_F),
(1 - y$fraction_F) * (1 - y$fraction_sexually_active_15_M),
0)
#
# MIXING
###############################################
y$M_noncomm = matrix(c(0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 1, 0, 0, 0,
0, 0, 0, 0, 1, 1, 0, 0, 0,
1, 1, 1, 1, 0, 0, 0, 0, 0,
0, 0, 1, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0),
nrow = 9, ncol = 9, byrow = T)
y$M_comm = matrix(c(0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
1, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0),
nrow = 9, ncol = 9, byrow = T)
# MOVEMENT
###############################################
y$rate_leave_low_FSW = y$rate_leave_pro_FSW
# virgin movement
y$rate_move_out[7] = - y$rate_enter_sexual_pop_F
y$rate_move_out[8] = - y$rate_enter_sexual_pop_M
y$rate_move_in[3,7] = y$rate_enter_sexual_pop_F
y$rate_move_in[6,8] = y$rate_enter_sexual_pop_M
#this is just to show what happens when you increase movement
# y$rate_leave_pro_FSW = y$rate_leave_pro_FSW * 10
# y$rate_leave_low_FSW = y$rate_leave_low_FSW * 10
# y$rate_leave_client = y$rate_leave_client * 10
y$prop_client_GPM = y$N_init[5] / y$N_init[6]
y$prop_pro_FSW_GPF = y$N_init[1] / y$N_init[3]
y$prop_low_FSW_GPF = y$N_init[2] / y$N_init[3]
# FEMALE MOVEMENT
y$rate_move_out[1] = - y$rate_leave_pro_FSW
y$rate_move_out[2] = - y$rate_leave_low_FSW
y$rate_move_out[3] = - (y$rate_leave_pro_FSW + y$muF + y$nu) * y$prop_pro_FSW_GPF * (1-y$fraction_FSW_foreign) -
(y$rate_leave_low_FSW + y$muF + y$nu) * y$prop_low_FSW_GPF* (1-y$fraction_FSW_foreign)
y$rate_move_in[1,3] = (y$rate_leave_pro_FSW + y$muF + y$nu) * y$prop_pro_FSW_GPF * (1-y$fraction_FSW_foreign) # moving from GPF to pro-FSW
y$rate_move_in[2,3] = (y$rate_leave_low_FSW + y$muF + y$nu) * y$prop_low_FSW_GPF * (1-y$fraction_FSW_foreign) # moving from GPF to low-FSW
y$rate_move_in[4,1] = y$rate_leave_pro_FSW * (1 - y$fraction_FSW_foreign) # moving from pro-FSW to former FSW in Cot
y$rate_move_in[4,2] = y$rate_leave_low_FSW * (1 - y$fraction_FSW_foreign) # moving from low-FSW to former FSW in Cot
y$rate_move_in[9,1] = y$rate_leave_pro_FSW * y$fraction_FSW_foreign # moving from low-FSW to former FSW NOT in Cot
y$rate_move_in[9,2] = y$rate_leave_low_FSW * y$fraction_FSW_foreign # moving from low-FSW to former FSW NOT in Cot
# MALE MOVEMENT
y$rate_move_out[5] = - y$rate_leave_client
y$rate_move_out[6] = - y$prop_client_GPM * (y$rate_leave_client + y$muM + y$nu)
y$rate_move_in[6,5] = y$rate_leave_client # moving from client to GPM
y$rate_move_in[5,6] = y$prop_client_GPM * (y$rate_leave_client + y$muM + y$nu) # moving from GPM to client
y$rate_move_out_PrEP = y$rate_move_out
# y$rate_move_out_PrEP[1] = -y$rate_move_out_PrEP_FSW
y$beta_comm = c(y$betaMtoF_FSW, y$betaMtoF_GPF, y$betaMtoF_GPF, y$betaMtoF_GPF, y$betaFtoM_client, y$betaFtoM_GPM, 0, 0, 0)
y$beta_noncomm = c(y$betaMtoF_FSW, y$betaMtoF_GPF, y$betaMtoF_GPF, y$betaMtoF_GPF, y$betaFtoM_client, y$betaFtoM_GPM, 0, 0, 0)
# c_t_comm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016),
# y$c_y_comm =
y$mu = c(y$muF, y$muF, y$muF, y$muF, y$muM, y$muM, y$muF, y$muM, y$muF)
} else {
# y$omega = y$omega/sum(y$omega)
}
if(y$movement == 0) {
y$rate_move_in = y$rate_move_in * 0
y$rate_move_out = y$rate_move_out * 0}
# extending all interpolating parameters beyond 2016 to the max of time
y$c_y_comm <- rbind(y$c_y_comm, y$c_y_comm[nrow(y$c_y_comm),])
y$c_t_comm <- c(y$c_t_comm, max(y$time))
y$c_y_noncomm <- rbind(y$c_y_noncomm, y$c_y_noncomm[nrow(y$c_y_noncomm),])
y$c_t_noncomm <- c(y$c_t_noncomm, max(y$time))
y$fc_y_comm <- abind::abind(y$fc_y_comm, y$fc_y_comm[length(y$fc_y_comm[,1,1]),,], along = 1)
y$fc_t_comm <- c(y$fc_t_comm, max(y$time))
y$fc_y_noncomm <- abind::abind(y$fc_y_noncomm, y$fc_y_noncomm[length(y$fc_y_noncomm[,1,1]),,], along = 1)
y$fc_t_noncomm <- c(y$fc_t_noncomm, max(y$time))
y$n_y_comm <- abind::abind(y$n_y_comm, y$n_y_comm[length(y$n_y_comm[,1,1]),,], along = 1)
y$n_t_comm <- c(y$n_t_comm, max(y$time))
y$n_y_noncomm <- abind::abind(y$n_y_noncomm, y$n_y_noncomm[length(y$n_y_noncomm[,1,1]),,], along = 1)
y$n_t_noncomm <- c(y$n_t_noncomm, max(y$time))
y$zetaa_t <- c(y$zetaa_t, max(y$time))
y$zetab_t <- c(y$zetab_t, max(y$time))
y$zetac_t <- c(y$zetac_t, max(y$time))
y$zetaa_y <- rbind(y$zetaa_y, y$zetaa_y[nrow(y$zetaa_y),])
y$zetab_y <- rbind(y$zetab_y, y$zetab_y[nrow(y$zetab_y),])
y$zetac_y <- rbind(y$zetac_y, y$zetac_y[nrow(y$zetac_y),])
y$testing_prob_y[which(y$testing_prob_t == 2006), c(2,3,4)] <- y$testing_prob_women_2006
y$testing_prob_y[which(y$testing_prob_t == 2008), c(2,3,4)] <- y$testing_prob_women_2008
y$testing_prob_y[which(y$testing_prob_t == 2012), c(2,3,4)] <- y$testing_prob_women_2012
y$testing_prob_y[which(y$testing_prob_t == 2013), c(2,3,4)] <- y$testing_prob_women_2012
y$testing_prob_y[which(y$testing_prob_t == 2015), c(2,3,4)] <- y$testing_prob_women_2012
y$testing_prob_y[which(y$testing_prob_t == 2016), c(2,3,4)] <- y$testing_prob_women_2012
y$testing_prob_y[which(y$testing_prob_t == 2006), c(5,6)] <- y$testing_prob_men_2006
y$testing_prob_y[which(y$testing_prob_t == 2008), c(5,6)] <- y$testing_prob_men_2008
y$testing_prob_y[which(y$testing_prob_t == 2012), c(5,6)] <- y$testing_prob_men_2012
y$testing_prob_y[which(y$testing_prob_t == 2013), c(5,6)] <- y$testing_prob_men_2012
y$testing_prob_y[which(y$testing_prob_t == 2015), c(5,6)] <- y$testing_prob_men_2012
y$testing_prob_y[which(y$testing_prob_t == 2016), c(5,6)] <- y$testing_prob_men_2012
y$testing_prob_y <- rbind(y$testing_prob_y, y$testing_prob_y[nrow(y$testing_prob_y),])
y$testing_prob_t <- c(y$testing_prob_t, max(y$time))
y$ART_prob_y <- rbind(y$ART_prob_y, y$ART_prob_y[nrow(y$ART_prob_y),])
y$ART_prob_t <- c(y$ART_prob_t, max(y$time))
y$fPb = 1 - y$fPa - y$fPc
# ART and efficacy
y$viral_supp_y[1,] = y$viral_supp_y_1986_rest
y$viral_supp_y[2,] = c(y$viral_supp_y_2015_ProFSW, rep_len(y$viral_supp_y_1986_rest, 8))
y$viral_supp_y[3,] = y$viral_supp_y[2,]
y$infect_ART_y = 1-(y$viral_supp_y * y$ART_eff)
y$infect_ART_t = y$viral_supp_t
y$rho = c(y$ART_recruit_rate_FSW, y$ART_recruit_rate_rest,
y$ART_recruit_rate_rest,
y$ART_recruit_rate_rest,
y$ART_recruit_rate_rest * y$ART_init_ratio_MF,
y$ART_recruit_rate_rest * y$ART_init_ratio_MF, 0, 0, 0)
y$iota = c(y$ART_reinit_rate_FSW, y$ART_reinit_rate_rest,
y$ART_reinit_rate_rest,
y$ART_reinit_rate_rest,
y$ART_reinit_rate_rest * y$ART_init_ratio_MF,
y$ART_reinit_rate_rest * y$ART_init_ratio_MF, 0, 0, 0)
if(y$long_intervention == 0)
{
y$tau_intervention_y = matrix(c(rep(0, 9), y$intervention_testing_increase, c(rep(0, 8)), y$intervention_testing_increase, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T)
y$rho_intervention_y = matrix(c(rep(0, 9), y$intervention_ART_increase, c(rep(0, 8)), y$intervention_ART_increase, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T)
y$prep_intervention_y = matrix(c(rep(0, 9), y$prep_offering_rate , rep(0, 9-1), y$prep_offering_rate , rep(0, 9-1), rep(0, 9)), ncol = 9, byrow = T)
} else {
y$tau_intervention_y = matrix(c(rep(0, 9), y$intervention_testing_increase, c(rep(0, 8)), y$intervention_testing_increase, c(rep(0, 8)), y$intervention_testing_increase, c(rep(0, 8))), ncol = 9, nrow = 4, byrow = T)
y$rho_intervention_y = matrix(c(rep(0, 9), y$intervention_ART_increase, c(rep(0, 8)), y$intervention_ART_increase, c(rep(0, 8)), y$intervention_ART_increase, c(rep(0, 8))), ncol = 9, nrow = 4, byrow = T)
y$prep_intervention_y = matrix(c(rep(0, 9), y$prep_offering_rate , rep(0, 9-1), y$prep_offering_rate , rep(0, 9-1), y$prep_offering_rate , rep(0, 9-1)), ncol = 9, byrow = T)
}
y$kappaa = c(y$PrEP_loss_to_follow_up, rep_len(0,(8)))
# y$kappaa = c(y$prep_dropout, rep_len(0,(8)))
y$kappab = y$kappaa
y$kappac = y$kappaa
y$kappa1 = y$kappaa
y$above_500_by_group = c(y$FSW_eligible, y$GP_eligible, y$GP_eligible,
y$GP_eligible, y$GP_eligible, y$GP_eligible,
y$GP_eligible, y$GP_eligible, y$GP_eligible)
# y$W1 = y$W1
# y$W2 = y$W2
# y$W3 = y$W3
# y$W4 = y$W4
y$pfFSW_y = matrix(c(rep(0, 9), y$prev_non_ben_fsw_1993, c(rep(0, 8)), y$prev_non_ben_fsw_2015, c(rep(0, 8)), y$prev_non_ben_fsw_2015, c(rep(0, 8))), ncol = 9, nrow = 4, byrow = T)
return(y)
}
# N_Ncat7 = function(y) return S0_init[1] + I01_init[1] + I02_init[1] + I03_init[1] + I04_init[1] + I05_init[1]
#
# the parameters below will be sampled from an LHS and will replace their respective defaults
# unless I put something in the args of the function, eg sample = mu
#' @export
#' @useDynLib cotonou
lhs_parameters <- function(n, sample = NULL, Ncat = 9, Nage = 1, ..., set_pars = list(...), forced_pars = list(...), set_null= list(...), ranges = NULL, par_seq, condom_seq, groups_seq, years_seq) {
set_pars <- modifyList(set_pars, forced_pars)
# parameters that I have defined in ranges will be removed from set_pars and fixed pars (fixed pars below)
if(length(which(names(set_pars) %in% rownames(ranges))) > 0)
set_pars <- set_pars[-which(names(set_pars) %in% rownames(ranges))]
#fixed pars list i think for the fix parameters function
fixed_pars = list(
# test_rate_prep = c(4, 0, 0, 0, 0, 0, 0, 0, 0),
fraction_sexually_active_15_F = 0,
fraction_sexually_active_15_M = 0,
initial_Ntot = 286114,
frac_women_ProFSW = 0.0024,
frac_women_LowFSW = 3,
frac_women_exFSW = 0.0024,
frac_men_client = 0.2,
frac_women_virgin = 0.1,
frac_men_virgin = 0.1,
fraction_F = 0.516,
rate_move_in = matrix(0, nrow = Ncat, ncol = Ncat),
rate_move_out = rep_len(0, Ncat),
rate_move_out_PrEP = rep_len(0, Ncat),
epsilon_y = 0,
rate_enter_sexual_pop_F = 0.4,
rate_enter_sexual_pop_M = 0.4,
fraction_FSW_foreign = 0,
movement = 1,
dur_below_200 = rep_len(0.3,Ncat),
fc_y_comm_1985 = matrix(0, Ncat, Ncat),
fc_y_comm_1993 = matrix(0, Ncat, Ncat),
fc_y_comm_1995 = matrix(0, Ncat, Ncat),
fc_y_comm_1998 = matrix(0, Ncat, Ncat),
fc_y_comm_2002 = matrix(0, Ncat, Ncat),
fc_y_comm_2005 = matrix(0, Ncat, Ncat),
fc_y_comm_2008 = matrix(0, Ncat, Ncat),
fc_y_comm_2012 = matrix(0, Ncat, Ncat),
fc_y_comm_2015 = matrix(0, Ncat, Ncat),
fc_y_comm_2016 = matrix(0, Ncat, Ncat),
fc_y_noncomm_1985 = matrix(0, Ncat, Ncat),
fc_y_noncomm_1993 = matrix(0, Ncat, Ncat),
fc_y_noncomm_1998 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2002 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2008 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2011 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2015 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2016 = matrix(0, Ncat, Ncat),
n_y_comm_1985 = matrix(1.02, Ncat, Ncat),
n_y_comm_2002 = matrix(1.02, Ncat, Ncat),
n_y_comm_2015 = matrix(1.02, Ncat, Ncat),
n_y_comm_2016 = matrix(1.02, Ncat, Ncat),
n_y_noncomm_1985 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_1998 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_2011 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_2002 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_2015 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_2016 = matrix(1.03, Ncat, Ncat),
c_comm_1985 = rep_len(2, Ncat),
c_comm_1993 = rep_len(2, Ncat),
c_comm_1995 = rep_len(2, Ncat),
c_comm_1998 = rep_len(2, Ncat),
c_comm_2002 = rep_len(2, Ncat),
c_comm_2005 = rep_len(2, Ncat),
c_comm_2008 = rep_len(2, Ncat),
c_comm_2012 = rep_len(2, Ncat),
c_comm_2015 = rep_len(2, Ncat),
c_comm_2016 = rep_len(2, Ncat),
c_noncomm_1985 = rep_len(1, Ncat),
c_noncomm_1993 = rep_len(1, Ncat),
c_noncomm_1995 = rep_len(1, Ncat),
c_noncomm_1998 = rep_len(1, Ncat),
c_noncomm_2002 = rep_len(1, Ncat),
c_noncomm_2005 = rep_len(1, Ncat),
c_noncomm_2008 = rep_len(1, Ncat),
c_noncomm_2012 = rep_len(1, Ncat),
c_noncomm_2015 = rep_len(1, Ncat),
c_noncomm_2016 = rep_len(1, Ncat),
betaMtoF_noncomm = 0.00193,
betaFtoM_noncomm = 0.00867,
betaMtoF_comm = 0.00193,
betaFtoM_comm = 0.00867,
betaMtoF_baseline = 0.00081,
muF = 0.02597403,
muM = 0.02739726,
RR_beta_FtM = 1,
RR_beta_circum = 0.44,
prev_HSV2_FSW = 1,
prev_HSV2_Client = 1,
prev_HSV2_GPF = 1,
prev_HSV2_GPM = 1,
RR_beta_HSV2_comm_t = 1,
RR_beta_HSV2_noncomm_t = 1,
RR_beta_HSV2_comm_a = 2,
RR_beta_HSV2_noncomm_a = 2,
who_believe_comm = 0,
init_clientN_from_PCR = 0,
beta_above_1 = 0,
ignore_ranges_fc_c = 0,
dropout_rate_not_FSW = 0.025,
dropout_rate_FSW = 0.025,
delete = 0,
nu = 0.02222222,
PrEPOnOff = 0,
fPa = 0.5,
fPc = 0.4,
iota = c(0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0, 0, 0),
viral_supp_t = c(1986, 2015, 2016),
viral_supp_y = matrix(0, nrow = 3, ncol = 9),
viral_supp_y_1986_rest = 0.6,
viral_supp_y_2015_ProFSW = 0.7,
ART_eff = 0.98,
rho = c(0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0, 0, 0),
ART_recruit_rate = 0.2,
ART_reinit_rate = 0.1,
RR_test_CD4200 = 5.4,
tau_intervention_t = c(1986, 2015, 2016, 2017),
rho_intervention_t = c(1986, 2015, 2016, 2017),
tau_intervention_y = matrix(c(rep(0, 9), 1.5, c(rep(0, 8)), 1.5, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T),
rho_intervention_y = matrix(c(rep(0, 9), 6, c(rep(0, 8)), 6, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T),
intervention_testing_increase = 1.4,
intervention_ART_increase = 6.5,
psia = rep_len(0.1,9),
psib = rep_len(0.1,9),
prep_offering_rate = 1,
kappaa = c(0.2, rep_len(0,(9-1))),
kappab = c(0.2, rep_len(0,(9-1))),
kappac = c(0.2, rep_len(0,(9-1))),
kappa1 = c(0.2, rep_len(0,(9-1))),
prep_dropout = 2,
eP1a = 0.9,
dur_FSW = 30,
cost_ART_initiation_study_FSW = 1,
cost_1_year_ART_study_FSW = 1,
cost_1_year_ART_rest = 1,
W0 = 1,
W1 = 1,
W2 = 1,
W3 = 1,
cost_Initiation_of_ART_study_FSW = 1,
cost_Initiation_of_ART_government_FSW = 1,
cost_1_year_of_ART_study_FSW = 1,
cost_1_year_of_ART_government_FSW = 1,
cost_Initiation_ART_rest_of_population = 1,
cost_1_year_of_ART_rest_of_population = 1,
cost_FSW_initiation_ART_Patient_costs = 1,
cost_FSW_1_year_ART_Patient_costs = 1,
cost_Initiation_of_PrEP_study = 1,
cost_1_year_PrEP_perfect_adherence_study = 1,
cost_1_year_PrEP_intermediate_adherence_study = 1,
cost_1_year_PrEP_non_adherence_study = 1,
cost_Initiation_of_PrEP_government = 1,
cost_1_year_PrEP_perfect_adherence_government =1,
cost_1_year_PrEP_intermediate_adherence_government = 1,
cost_1_year_PrEP_non_adherence_government = 1,
cost_PREP_initiation_Patient_costs = 1,
cost_PREP_1_year_ART_Patient_costs = 1,
TasP_testing = 0,
long_intervention = 0,
above_500_by_group = c(1, 1, 1, 1, 1, 1, 1, 1, 1),
FSW_eligible = 0,
GP_eligible = 1,
prev_non_ben_fsw_1993 = 0.55,
prev_non_ben_fsw_2015 = 0.2,
infected_FSW_incoming = 1,
ART_recruit_rate_rest = 0.1,
ART_reinit_rate_rest = 0.1,
ART_init_ratio_MF = 2,
RR_beta_HSV2_a_FSW = 1,
RR_beta_HSV2_a_client = 1,
RR_beta_HSV2_a_GPF = 1,
RR_beta_HSV2_a_GPM = 1,
testing_prob_men_2006 = 0.111,
testing_prob_men_2008 = 0.112,
testing_prob_men_2012 = 0.113,
testing_prob_women_2006 = 0.114,
testing_prob_women_2008 = 0.115,
testing_prob_women_2012 = 0.116,
ART_recruit_rate_FSW = 1,
ART_reinit_rate_FSW = 1,
PrEP_reinit_OnOff_t = c(1985, 2017, 2017.01, 2060),
PrEP_reinit_OnOff_y = c(0, 0, 0, 0),
PrEP_reinits_on = 1,
old_VS_assumption = 1
)
if(length(which(names(fixed_pars) %in% rownames(ranges))) > 0)
fixed_pars <- fixed_pars[-which(names(fixed_pars) %in% rownames(ranges))]
mu <- matrix(rep(c(1/50, 1/42), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("mu", Ncat), NULL))
omega <- if(Ncat == 9) matrix(c(0.0017, 0.0067, 0, 0, 0, 0, 0, 0, 0.1, 0.2, 0, 0, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("omega", Ncat), NULL)) else
matrix(rep(c(0.4, 0.6), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("omega", Ncat), NULL))
# c_y_comm <- if(Ncat == 9) matrix(c(300, 1400, 40, 64, 0, 0, 0, 0, 18.67, 37.5, 0, 0, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_y_comm", Ncat), NULL)) else c(1, 3)
#these parameters need to be here so fix_parameters works? below are fixed...
S0_init = matrix(rep(c(4000, 4000), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("S0_init", Ncat), NULL))
I01_init = matrix(rep(c(1000, 1000), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("I01_init", Ncat), NULL))
N_init = if(Ncat == 9) matrix(c(672, 672, 757, 757, 130895, 130895, 672, 672, 27091, 27091, 100335, 100335, 14544, 14544, 11148, 11148, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("N_init", Ncat), NULL)) else c(300000, 300000)
# c_comm = if(Ncat == 9) matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL)) else
# matrix(rep(c(1,3), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL))
# c_comm = if(Ncat == 9) matrix(c(272, 1439, 40, 64, 0, 0, 0, 0, 18.67, 37.5, 0, 0, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL)) else
# matrix(rep(c(1,3), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL))
#
# c_noncomm = if(Ncat == 9) matrix(c(0.2729358, 0.4682779, 0.2729358, 0.4682779, 0.90, 1.02, 0.90, 1.02, 1.21, 2.5, 1.28, 1.40, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_noncomm", Ncat), NULL)) else
# matrix(rep(c(1,3), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_noncomm", Ncat), NULL))
#
old_ranges <- rbind(
# c_y_comm,
epsilon_1985 = c(0.059346131*1.5, 0.059346131*1.5),
epsilon_1992 = c(0.053594832*1.5, 0.053594832*1.5),
epsilon_2002 = c(0.026936907*1.5, 0.026936907*1.5),
epsilon_2013 = c(0.026936907*1.5, 0.026936907*1.5),
epsilon_2016 = c(0.026936907*1.5, 0.026936907*1.5),
rate_leave_pro_FSW = c(0.2, 0.2),
rate_leave_low_FSW = c(0.1, 0.1),
rate_leave_client = c(0.05, 0.05),
betaMtoF = c(0.00086, 0.00433),
betaFtoM = c(0.00279, 0.02701),
prev_init_FSW = c(0.01318836, 0.06592892),
prev_init_rest = c(0.0003134459, 0.0029420363),
# c_comm,
# c_noncomm,
mu,
eff_ART = c(0.96, 1), # infectiousness RR when on ART
infect_acute = c(4, 18), # RR for acute phase
# dur_primary_phase = c(1/0.5, 1/0.16), # rate
dur_primary_phase = c(0.16, 0.5), # from Mathieu's parameters IN YEARS
# dur_primary_phase = c(2, 6.25), # from Mathieu's parameters
SC_to_death = c(8.7, 12.3),
dur_200_349 = c(3.9, 5), # from Mathieu's parameters IN YEARS
# dur_200_349 = c(3.9, 5), # from Mathieu's parameters
ART_RR_prog = c(8.8, 12.1),
# ART_RR_mort = c(2, 3),
# omega,
#below are fixed...
S0_init,
I01_init,
N_init
)
if (is.null(ranges)) {
ranges <- old_ranges
}
if (!is.null(sample)) {
ranges <- ranges[rownames(ranges) %in% sample, drop=FALSE]
}
samples <- tgp::lhs(n, ranges)
nms <- rownames(ranges)
i <- split(seq_along(nms), nms)
f <- function(x) {
lapply(i, function(j) x[j])
}
samples_list <- apply(samples, 1, f)
# print("ehre:")
# print(samples_list)
samples_list <- lapply(samples_list, function(x) modifyList(x, fixed_pars)) # btw earlier, fixed pars is replaced by ranges where they overlap
# HERE!
# samples_list <- lapply(samples_list, after_LHS, set_pars = set_pars)
samples_list <- lapply(samples_list, function(x) modifyList(x, set_pars)) # set pars before fixed pars in order to get right fixed
# print("ehre2:")
# print(samples_list)
# samples_list_test <<- samples_list
samples_list <- lapply(samples_list, fix_parameters, Ncat = Ncat, par_seq = par_seq, condom_seq = condom_seq, groups_seq = groups_seq, years_seq = years_seq)
# samples_list_test_fixed <<- samples_list
# samples_list <- lapply(samples_list, function(x) modifyList(x, set_pars)) # set pars after fixed pars in order to get right set pars
samples_list <- lapply(samples_list, function(x) modifyList(x, forced_pars)) # set pars after fixed pars in order to get right set pars
samples_list <- lapply(samples_list, delete_parameter_sets)
samples_list_test <<- samples_list
if(n > 1)
samples_list <- samples_list[!unlist(lapply(samples_list, function(x) x$delete))] # removing those runs which do not pass initial criteria
lapply(samples_list, function(x) generate_parameters(parameters = x, Ncat = Ncat, set_null = set_null))
}
#' @export
#' @useDynLib cotonou
delete_parameter_sets <- function(x) {
x$delete <- ifelse(x$beta_above_1 == 1, 1, # if beta is above 1 due to the various RRs
# if believing FSW partner change rates, then remove runs where clients pop size too big or too small ONLY IF client size is fixed at start from PCR balancing
ifelse(((x$c_comm_1985[1] * x$N_init[1] + x$c_comm_1985[2] * x$N_init[2])/ (x$c_comm_1985[5])) > (0.4 * sum(c(x$N_init[5], x$N_init[6], x$N_init[8]))) && (x$init_clientN_from_PCR == 1), 1,
# if believing client partner change rates, we are now also believing pro FSW and altering low FSW. killing runs where frac of low fsw is too low
ifelse(((x$c_comm_1985[5] * x$N_init[5] - x$c_comm_1985[2] * x$N_init[2])/ (x$c_comm_1985[1])) < (0* sum(x$N_init)) && (x$init_clientN_from_PCR == 0), 1, 0)))
return(x)
}
#' @export
#' @useDynLib cotonou
lhs_parameters_parallel <- function(n, sample = NULL, Ncat = 9, Nage = 1, ..., set_pars = list(...), forced_pars = list(...), set_null= list(...), ranges = NULL, par_seq, condom_seq, groups_seq, years_seq) {
set_pars <- modifyList(set_pars, forced_pars)
# parameters that I have defined in ranges will be removed from set_pars and fixed pars (fixed pars below)
if(length(which(names(set_pars) %in% rownames(ranges))) > 0)
set_pars <- set_pars[-which(names(set_pars) %in% rownames(ranges))]
#fixed pars list i think for the fix parameters function
fixed_pars = list(
# test_rate_prep = c(4, 0, 0, 0, 0, 0, 0, 0, 0),
fraction_sexually_active_15_F = 0,
fraction_sexually_active_15_M = 0,
initial_Ntot = 286114,
frac_women_ProFSW = 0.0024,
frac_women_LowFSW = 3,
frac_women_exFSW = 0.0024,
frac_men_client = 0.2,
frac_women_virgin = 0.1,
frac_men_virgin = 0.1,
fraction_F = 0.516,
rate_move_in = matrix(0, nrow = Ncat, ncol = Ncat),
rate_move_out = rep_len(0, Ncat),
rate_move_out_PrEP = rep_len(0, Ncat),
epsilon_y = 0,
rate_enter_sexual_pop_F = 0.4,
rate_enter_sexual_pop_M = 0.4,
fraction_FSW_foreign = 0,
movement = 1,
dur_below_200 = rep_len(0.3,Ncat),
fc_y_comm_1985 = matrix(0, Ncat, Ncat),
fc_y_comm_1993 = matrix(0, Ncat, Ncat),
fc_y_comm_1995 = matrix(0, Ncat, Ncat),
fc_y_comm_1998 = matrix(0, Ncat, Ncat),
fc_y_comm_2002 = matrix(0, Ncat, Ncat),
fc_y_comm_2005 = matrix(0, Ncat, Ncat),
fc_y_comm_2008 = matrix(0, Ncat, Ncat),
fc_y_comm_2012 = matrix(0, Ncat, Ncat),
fc_y_comm_2015 = matrix(0, Ncat, Ncat),
fc_y_comm_2016 = matrix(0, Ncat, Ncat),
fc_y_noncomm_1985 = matrix(0, Ncat, Ncat),
fc_y_noncomm_1993 = matrix(0, Ncat, Ncat),
fc_y_noncomm_1998 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2002 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2008 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2011 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2015 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2016 = matrix(0, Ncat, Ncat),
n_y_comm_1985 = matrix(1.02, Ncat, Ncat),
n_y_comm_2002 = matrix(1.02, Ncat, Ncat),
n_y_comm_2015 = matrix(1.02, Ncat, Ncat),
n_y_comm_2016 = matrix(1.02, Ncat, Ncat),
n_y_noncomm_1985 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_2002 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_2015 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_1998 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_2011 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_2016 = matrix(1.03, Ncat, Ncat),
c_comm_1985 = rep_len(2, Ncat),
c_comm_1993 = rep_len(2, Ncat),
c_comm_1995 = rep_len(2, Ncat),
c_comm_1998 = rep_len(2, Ncat),
c_comm_2002 = rep_len(2, Ncat),
c_comm_2005 = rep_len(2, Ncat),
c_comm_2008 = rep_len(2, Ncat),
c_comm_2012 = rep_len(2, Ncat),
c_comm_2015 = rep_len(2, Ncat),
c_comm_2016 = rep_len(2, Ncat),
c_noncomm_1985 = rep_len(1, Ncat),
c_noncomm_1993 = rep_len(1, Ncat),
c_noncomm_1995 = rep_len(1, Ncat),
c_noncomm_1998 = rep_len(1, Ncat),
c_noncomm_2002 = rep_len(1, Ncat),
c_noncomm_2005 = rep_len(1, Ncat),
c_noncomm_2008 = rep_len(1, Ncat),
c_noncomm_2012 = rep_len(1, Ncat),
c_noncomm_2015 = rep_len(1, Ncat),
c_noncomm_2016 = rep_len(1, Ncat),
betaMtoF_noncomm = 0.00193,
betaFtoM_noncomm = 0.00867,
betaMtoF_comm = 0.00193,
betaFtoM_comm = 0.00867,
betaMtoF_baseline = 0.00081,
muF = 0.02597403,
muM = 0.02739726,
RR_beta_FtM = 1,
RR_beta_circum = 0.44,
prev_HSV2_FSW = 1,
prev_HSV2_Client = 1,
prev_HSV2_GPF = 1,
prev_HSV2_GPM = 1,
RR_beta_HSV2_comm_t = 1,
RR_beta_HSV2_noncomm_t = 1,
RR_beta_HSV2_comm_a = 2,
RR_beta_HSV2_noncomm_a = 2,
who_believe_comm = 0,
init_clientN_from_PCR = 0,
beta_above_1 = 0,
ignore_ranges_fc_c = 0,
dropout_rate_not_FSW = 0.025,
dropout_rate_FSW = 0.025,
delete = 0,
nu = 0.02222222,
PrEPOnOff = 0,
fPa = 0.5,
fPc = 0.4,
iota = c(0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0, 0, 0),
viral_supp_t = c(1986, 2015, 2016),
viral_supp_y = matrix(0, nrow = 3, ncol = 9),
viral_supp_y_1986_rest = 0.6,
viral_supp_y_2015_ProFSW = 0.7,
ART_eff = 0.98,
rho = c(0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0, 0, 0),
ART_recruit_rate = 0.2,
ART_reinit_rate = 0.1,
RR_test_CD4200 = 5.4,
tau_intervention_t = c(1986, 2015, 2016, 2017),
rho_intervention_t = c(1986, 2015, 2016, 2017),
tau_intervention_y = matrix(c(rep(0, 9), 0, c(rep(0, 8)), 0, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T),
rho_intervention_y = matrix(c(rep(0, 9), 0, c(rep(0, 8)), 0, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T),
intervention_testing_increase = 1.4,
intervention_ART_increase = 6.5,
psia = rep_len(0.1,9),
psib = rep_len(0.1,9),
prep_offering_rate = 1,
kappaa = c(0.2, rep_len(0,(9-1))),
kappab = c(0.2, rep_len(0,(9-1))),
kappac = c(0.2, rep_len(0,(9-1))),
kappa1 = c(0.2, rep_len(0,(9-1))),
prep_dropout = 2,
eP1a = 0.9,
dur_FSW = 30,
cost_ART_initiation_study_FSW = 1,
cost_1_year_ART_study_FSW = 1,
cost_1_year_ART_rest = 1,
W0 = 1,
W1 = 1,
W2 = 1,
W3 = 1,
cost_Initiation_of_ART_study_FSW = 1,
cost_Initiation_of_ART_government_FSW = 1,
cost_1_year_of_ART_study_FSW = 1,
cost_1_year_of_ART_government_FSW = 1,
cost_Initiation_ART_rest_of_population = 1,
cost_1_year_of_ART_rest_of_population = 1,
cost_FSW_initiation_ART_Patient_costs = 1,
cost_FSW_1_year_ART_Patient_costs = 1,
cost_Initiation_of_PrEP_study = 1,
cost_1_year_PrEP_perfect_adherence_study = 1,
cost_1_year_PrEP_intermediate_adherence_study = 1,
cost_1_year_PrEP_non_adherence_study = 1,
cost_Initiation_of_PrEP_government = 1,
cost_1_year_PrEP_perfect_adherence_government =1,
cost_1_year_PrEP_intermediate_adherence_government = 1,
cost_1_year_PrEP_non_adherence_government = 1,
cost_PREP_initiation_Patient_costs = 1,
cost_PREP_1_year_ART_Patient_costs = 1,
TasP_testing = 0,
long_intervention = 0,
above_500_by_group = c(1, 1, 1, 1, 1, 1, 1, 1, 1),
FSW_eligible = 0,
GP_eligible = 1,
prev_non_ben_fsw_1993 = 0.55,
prev_non_ben_fsw_2015 = 0.2,
infected_FSW_incoming = 1,
ART_recruit_rate_rest = 0.1,
ART_reinit_rate_rest = 0.1,
ART_init_ratio_MF = 2,
RR_beta_HSV2_a_FSW = 1,
RR_beta_HSV2_a_client = 1,
RR_beta_HSV2_a_GPF = 1,
RR_beta_HSV2_a_GPM = 1,
testing_prob_men_2006 = 0.111,
testing_prob_men_2008 = 0.112,
testing_prob_men_2012 = 0.113,
testing_prob_women_2006 = 0.114,
testing_prob_women_2008 = 0.115,
testing_prob_women_2012 = 0.116,
ART_recruit_rate_FSW = 1,
ART_reinit_rate_FSW = 1,
PrEP_reinit_OnOff_t = c(1985, 2017, 2017.01, 2060),
PrEP_reinit_OnOff_y = c(0, 0, 0, 0),
PrEP_reinits_on = 1,
old_VS_assumption = 1
)
if(length(which(names(fixed_pars) %in% rownames(ranges))) > 0)
fixed_pars <- fixed_pars[-which(names(fixed_pars) %in% rownames(ranges))]
mu <- matrix(rep(c(1/50, 1/42), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("mu", Ncat), NULL))
omega <- if(Ncat == 9) matrix(c(0.0017, 0.0067, 0, 0, 0, 0, 0, 0, 0.1, 0.2, 0, 0, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("omega", Ncat), NULL)) else
matrix(rep(c(0.4, 0.6), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("omega", Ncat), NULL))
# c_y_comm <- if(Ncat == 9) matrix(c(300, 1400, 40, 64, 0, 0, 0, 0, 18.67, 37.5, 0, 0, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_y_comm", Ncat), NULL)) else c(1, 3)
#these parameters need to be here so fix_parameters works? below are fixed...
S0_init = matrix(rep(c(4000, 4000), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("S0_init", Ncat), NULL))
I01_init = matrix(rep(c(1000, 1000), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("I01_init", Ncat), NULL))
N_init = if(Ncat == 9) matrix(c(672, 672, 757, 757, 130895, 130895, 672, 672, 27091, 27091, 100335, 100335, 14544, 14544, 11148, 11148, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("N_init", Ncat), NULL)) else c(300000, 300000)
# c_comm = if(Ncat == 9) matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL)) else
# matrix(rep(c(1,3), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL))
# c_comm = if(Ncat == 9) matrix(c(272, 1439, 40, 64, 0, 0, 0, 0, 18.67, 37.5, 0, 0, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL)) else
# matrix(rep(c(1,3), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL))
#
# c_noncomm = if(Ncat == 9) matrix(c(0.2729358, 0.4682779, 0.2729358, 0.4682779, 0.90, 1.02, 0.90, 1.02, 1.21, 2.5, 1.28, 1.40, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_noncomm", Ncat), NULL)) else
# matrix(rep(c(1,3), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_noncomm", Ncat), NULL))
#
old_ranges <- rbind(
# c_y_comm,
epsilon_1985 = c(0.059346131*1.5, 0.059346131*1.5),
epsilon_1992 = c(0.053594832*1.5, 0.053594832*1.5),
epsilon_2002 = c(0.026936907*1.5, 0.026936907*1.5),
epsilon_2013 = c(0.026936907*1.5, 0.026936907*1.5),
epsilon_2016 = c(0.026936907*1.5, 0.026936907*1.5),
rate_leave_pro_FSW = c(0.2, 0.2),
rate_leave_low_FSW = c(0.1, 0.1),
rate_leave_client = c(0.05, 0.05),
betaMtoF = c(0.00086, 0.00433),
betaFtoM = c(0.00279, 0.02701),
prev_init_FSW = c(0.01318836, 0.06592892),
prev_init_rest = c(0.0003134459, 0.0029420363),
# c_comm,
# c_noncomm,
mu,
eff_ART = c(0.96, 1), # infectiousness RR when on ART
infect_acute = c(4, 18), # RR for acute phase
# dur_primary_phase = c(1/0.5, 1/0.16), # rate
dur_primary_phase = c(0.16, 0.5), # from Mathieu's parameters IN YEARS
# dur_primary_phase = c(2, 6.25), # from Mathieu's parameters
SC_to_death = c(8.7, 12.3),
dur_200_349 = c(3.9, 5), # from Mathieu's parameters IN YEARS
# dur_200_349 = c(3.9, 5), # from Mathieu's parameters
ART_RR_prog = c(8.8, 12.1),
# ART_RR_mort = c(2, 3),
# omega,
#below are fixed...
S0_init,
I01_init,
N_init
)
if (is.null(ranges)) {
ranges <- old_ranges
}
if (!is.null(sample)) {
ranges <- ranges[rownames(ranges) %in% sample, drop=FALSE]
}
samples <- tgp::lhs(n, ranges)
nms <- rownames(ranges)
i <- split(seq_along(nms), nms)
f <- function(x) {
lapply(i, function(j) x[j])
}
samples_list <- apply(samples, 1, f)
# print("ehre:")
# print(samples_list)
samples_list <- parallel::parLapply(NULL, samples_list, function(x) modifyList(x, fixed_pars)) # btw earlier, fixed pars is replaced by ranges where they overlap
# HERE!
# samples_list <- lapply(samples_list, after_LHS, set_pars = set_pars)
samples_list <- parallel::parLapply(NULL, samples_list, function(x) modifyList(x, set_pars)) # set pars before fixed pars in order to get right fixed
# print("ehre2:")
# print(samples_list)
# samples_list_test <<- samples_list
samples_list <- parallel::parLapply(NULL, samples_list, fix_parameters, Ncat = Ncat, par_seq = par_seq, condom_seq = condom_seq, groups_seq = groups_seq, years_seq = years_seq)
# samples_list_test_fixed <<- samples_list
# samples_list <- lapply(samples_list, function(x) modifyList(x, set_pars)) # set pars after fixed pars in order to get right set pars
samples_list <- parallel::parLapply(NULL, samples_list, function(x) modifyList(x, forced_pars)) # set pars after fixed pars in order to get right set pars
samples_list <- parallel::parLapply(NULL, samples_list, delete_parameter_sets)
if(n > 1)
samples_list <- samples_list[!unlist(parallel::parLapply(NULL, samples_list, function(x) x$delete))] # removing those runs which do not pass initial criteria
parallel::parLapply(NULL, samples_list, function(x) generate_parameters(parameters = x, Ncat = Ncat, set_null = set_null))
}
#' @export
#' @useDynLib cotonou
generate_parameters <- function(..., parameters = list(...), set_null = list(...), Ncat = 9, Nage = 1) {
defaults <- list(Ncat = Ncat,
Nage = Nage,
N_init = 300000, #ignore that this is weird - I replace with a vector later
initial_Ntot = 286114,
frac_women_ProFSW = 0.0024,
frac_women_LowFSW = 3,
frac_women_exFSW = 0.0024,
frac_men_client = 0.2,
frac_women_virgin = 0.1,
frac_men_virgin = 0.1,
above_500_by_group = c(1, 1, 1, 1, 1, 1, 1, 1, 1),
ART_eligible_CD4_above_500_t = c(1985, 2002, 2012, 2015, 2016),
ART_eligible_CD4_350_500_t = c(1985, 2002, 2012, 2015, 2016),
ART_eligible_CD4_200_349_t = c(1985, 2002, 2012, 2015, 2016),
ART_eligible_CD4_below_200_t = c(1985, 2002, 2012, 2015, 2016),
ART_eligible_CD4_above_500_y = c(0, 0, 0, 0, 1),
ART_eligible_CD4_350_500_y = c(0, 0, 0, 1, 1),
ART_eligible_CD4_200_349_y = c(0, 0.1, 1, 1, 1),
ART_eligible_CD4_below_200_y = c(0, 1, 1, 1, 1),
# epsilon_t_comm = c(1985, 1991.99, 1992, 2001.99, 2002, 2012.99, 2013, 2016),
# epsilon_y_comm = c(0.059346131, 0.059346131, 0.053594832, 0.053594832, 0.026936907, 0.026936907, 0.026936907, 0.026936907),
epsilon_t = c(1985, 1992, 2002, 2013, 2016),
epsilon_y = c(0.059346131, 0.053594832, 0.026936907, 0.026936907, 0.026936907),
epsilon_1985 = 0.059346131,
epsilon_1992 = 0.053594832,
epsilon_2002 = 0.026936907,
epsilon_2013 = 0.026936907,
epsilon_2016 = 0.026936907,
S0_init = rep_len(2000, Ncat),
S1a_init = rep_len(0, Ncat),
S1b_init = rep_len(0, Ncat),
S1c_init = rep_len(0, Ncat),
S1d_init = rep_len(0, Ncat),
I01_init = rep_len(1000, Ncat),
I11_init = rep_len(0, Ncat),
I02_init = rep_len(0, Ncat),
I03_init = rep_len(0, Ncat),
I04_init = rep_len(0, Ncat),
I05_init = rep_len(0, Ncat),
I22_init = rep_len(0, Ncat),
I23_init = rep_len(0, Ncat),
I24_init = rep_len(0, Ncat),
I25_init = rep_len(0, Ncat),
I32_init = rep_len(0, Ncat),
I33_init = rep_len(0, Ncat),
I34_init = rep_len(0, Ncat),
I35_init = rep_len(0, Ncat),
I42_init = rep_len(0, Ncat),
I43_init = rep_len(0, Ncat),
I44_init = rep_len(0, Ncat),
I45_init = rep_len(0, Ncat),
cumuInf_init = rep_len(0, Ncat),
mu = rep_len(0.02, Ncat),
dur_primary_phase = 0.3,
gamma02 = rep_len(0.2, Ncat),
gamma03 = rep_len(0.2, Ncat),
dur_200_349 = 3,
gamma11 = rep_len(0.2, Ncat),
gamma22 = rep_len(0.2, Ncat),
gamma23 = rep_len(0.2, Ncat),
gamma24 = rep_len(0.2, Ncat),
gamma32 = rep_len(0.2, Ncat),
gamma33 = rep_len(0.2, Ncat),
gamma34 = rep_len(0.2, Ncat),
gamma42 = rep_len(0.2, Ncat),
gamma43 = rep_len(0.2, Ncat),
gamma44 = rep_len(0.2, Ncat),
ART_prob_t = c(1985, 2005, 2012, 2015),
ART_prob_y = matrix(
rep(c(0, 0, 0.2, 0.4), Ncat), ncol = Ncat),
rho = c(0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0, 0, 0),
phi2 = rep_len(0.004,Ncat), # sort out later
phi3 = rep_len(0.004,Ncat),
phi4 = rep_len(0.004,Ncat),
phi5 = rep_len(0.004,Ncat),
psia = rep_len(0.1,Ncat),
psib = rep_len(0.1,Ncat),
testing_prob_t = c(1985, 2001, 2005, 2006, 2008, 2012, 2013, 2015, 2016),
testing_prob_y = matrix(
rep(c(0, 0.1, 0.2, 0.4, 0.5, 0.7, 0.5, 0.8, 0.7), Ncat), ncol = Ncat),
test_rate_prep = c(4, 0, 0, 0, 0, 0, 0, 0, 0),
sigma = c(0.86, 0, 0, 0, 0, 0, 0, 0, 0),
prep_intervention_t = c(1985, 2015, 2016, 2017),
prep_intervention_y = matrix(c(rep(0, Ncat), .1, rep(0, Ncat-1), .1, rep(0, Ncat-1), rep(0, Ncat)), ncol = Ncat, byrow = T), # offering rate
PrEPOnOff = 0,
RR_test_onPrEP = 2,
RR_test_CD4200 = 5.4,
RR_ART_CD4200 = 2,
tau01 = rep_len(1,Ncat),
tau11 = rep_len(1,Ncat),
tau2 = rep_len(1,Ncat),
tau3 = rep_len(1,Ncat),
tau4 = rep_len(1,Ncat),
tau5 = rep_len(1,Ncat),
# PREP
zetaa_t = c(1985, 2013, 2015, 2016),
zetaa_y = matrix(c(rep(0, Ncat), 0.0075, rep(0, Ncat-1), rep(0, Ncat), rep(0, Ncat)), ncol = Ncat, byrow = T),
zetab_t = c(1985, 2013, 2015, 2016),
zetab_y = matrix(c(rep(0, Ncat), 0.0075, rep(0, Ncat-1), rep(0, Ncat), rep(0, Ncat)), ncol = Ncat, byrow = T),
zetac_t = c(1985, 2013, 2015, 2016),
zetac_y = matrix(c(rep(0, Ncat), 0.0075, rep(0, Ncat-1), rep(0, Ncat), rep(0, Ncat)), ncol = Ncat, byrow = T),
eP = rep_len(0.6,Ncat),
eP0 = c(0, rep_len(0, (Ncat-1))),
eP1a = c(0.9, rep_len(0, (Ncat-1))),
eP1b = c(0.45, rep_len(0, (Ncat-1))),
eP1c = c(0, rep_len(0, (Ncat-1))),
eP1d = c(0, rep_len(0, (Ncat-1))),
# partner change rate
c_comm_1985 = rep_len(2, Ncat),
c_comm_1993 = rep_len(2, Ncat),
c_comm_1995 = rep_len(2, Ncat),
c_comm_1998 = rep_len(2, Ncat),
c_comm_2002 = rep_len(2, Ncat),
c_comm_2005 = rep_len(2, Ncat),
c_comm_2008 = rep_len(2, Ncat),
c_comm_2012 = rep_len(2, Ncat),
c_comm_2015 = rep_len(2, Ncat),
c_comm_2016 = rep_len(2, Ncat),
c_noncomm_1985 = rep_len(1, Ncat),
c_noncomm_1993 = rep_len(1, Ncat),
c_noncomm_1995 = rep_len(1, Ncat),
c_noncomm_1998 = rep_len(1, Ncat),
c_noncomm_2002 = rep_len(1, Ncat),
c_noncomm_2005 = rep_len(1, Ncat),
c_noncomm_2008 = rep_len(1, Ncat),
c_noncomm_2012 = rep_len(1, Ncat),
c_noncomm_2015 = rep_len(1, Ncat),
c_noncomm_2016 = rep_len(1, Ncat),
c_t_comm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016),
c_t_noncomm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016),
c_y_comm = matrix(2, nrow = 10, ncol = Ncat),
c_y_noncomm = matrix(1, nrow = 10, ncol = Ncat),
# condoms
# fc_t_noncomm = c(1985, 1990, 1998, 2016),
# fc_y_noncomm = matrix(
# rep(c(0, 0, 0.3, 0.5), Ncat), ncol = Ncat),
fc_y_comm_1985 = matrix(0.2, Ncat, Ncat),
fc_y_comm_1993 = matrix(0.5, Ncat, Ncat),
fc_y_comm_1995 = matrix(0.7, Ncat, Ncat),
fc_y_comm_1998 = matrix(0.3, Ncat, Ncat),
fc_y_comm_2002 = matrix(0.4, Ncat, Ncat),
fc_y_comm_2005 = matrix(0.1, Ncat, Ncat),
fc_y_comm_2008 = matrix(0.1, Ncat, Ncat),
fc_y_comm_2012 = matrix(0.6, Ncat, Ncat),
fc_y_comm_2015 = matrix(0, Ncat, Ncat),
fc_y_comm_2016 = matrix(0, Ncat, Ncat),
fc_y_noncomm_1985 = matrix(0.2, Ncat, Ncat),
fc_y_noncomm_1993 = matrix(0.2, Ncat, Ncat),
fc_y_noncomm_1998 = matrix(0.4, Ncat, Ncat),
fc_y_noncomm_2002 = matrix(0.3, Ncat, Ncat),
fc_y_noncomm_2008 = matrix(0.3, Ncat, Ncat),
# fc_y_noncomm_2011 = matrix(0.3, Ncat, Ncat),
# fc_y_noncomm_2015 = matrix(0.5, Ncat, Ncat),
# fc_y_noncomm_2016 = matrix(0.5, Ncat, Ncat),
fc_y_comm = array(0.5,dim=c(10,Ncat,Ncat)),
fc_y_noncomm = array(0.5,dim=c(7,Ncat,Ncat)),
fc_t_comm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016),
# fc_y_comm = matrix(
# rep(c(0.5, 0.5, 0.9, 0.99), Ncat), ncol = Ncat),
fc_t_noncomm = c(1985, 1993, 1998, 2002, 2008, 2011, 2015, 2016),
# fc_y_comm = matrix(
# rep(c(0.5, 0.5, 0.9, 0.99), Ncat), ncol = Ncat),
kappaa = c(0.2, rep_len(0,(Ncat-1))),
kappab = c(0.2, rep_len(0,(Ncat-1))),
kappac = c(0.2, rep_len(0,(Ncat-1))),
kappa1 = c(0.2, rep_len(0,(Ncat-1))),
alpha01 = rep_len(0,Ncat),
alpha02 = rep_len(0,Ncat),
alpha03 = 0.03,
alpha04 = 0.07,
# dur_below_200 = rep_len(0.3, Ncat),
alpha11 = rep_len(0.01,Ncat),
# alpha21 = rep_len(0.01,Ncat),
alpha22 = rep_len(0.01,Ncat),
alpha23 = rep_len(0.01,Ncat),
alpha24 = rep_len(0.01,Ncat),
alpha25 = rep_len(0.3,Ncat),
alpha32 = rep_len(0.01,Ncat),
# alpha33 = rep_len(0.01,Ncat),
# alpha34 = rep_len(0.01,Ncat),
# alpha35 = rep_len(0.3,Ncat),
alpha42 = rep_len(0.01,Ncat),
alpha43 = rep_len(0.01,Ncat),
alpha44 = rep_len(0.01,Ncat),
alpha45 = rep_len(0.3,Ncat),
# beta = rep_len(0.005,Ncat),
# betaMtoF = 0.00193,
# betaFtoM = 0.00867,
# betaMtoF_noncomm = 0.00193,
# betaFtoM_noncomm = 0.00867,
# betaMtoF_comm = 0.00193,
# betaFtoM_comm = 0.00867,
#beta = 0,
# p_comm = matrix(1, ncol = Ncat, nrow = Ncat),
# p_noncomm = matrix(1, ncol = Ncat, nrow = Ncat),
ec = rep_len(0.8,Ncat),
# ec = rep_len(1,1),
# epsilon = 0.001,
# c_comm = rep_len(2, Ncat),
# c_noncomm = rep_len(2, Ncat),
# c_t_comm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016),
# c_y_comm = matrix(rep(c(1, 2, 3, 1, 3, 4, 2, 1, 2, 4), Ncat), ncol = Ncat),
#
# c_t_noncomm = c(1985, 1990, 1998, 2016),
# c_y_noncomm = matrix(rep(c(0.5, 1, 0.7, 0.9), Ncat), ncol = Ncat),
fP_t_comm = c(1985, 2014, 2015, 2016, 2060),
fP_y_comm = matrix(
rep(c(1, 1, 1, 1, 1), Ncat), ncol = Ncat),
fP_t_noncomm = c(1985, 2014, 2015, 2016, 2060),
fP_y_noncomm = matrix(
rep(c(1, 1, 1, 1, 1), Ncat), ncol = Ncat),
# n_y_comm_1985 = matrix(1.02, Ncat, Ncat),
# n_y_comm_2002 = matrix(1.02, Ncat, Ncat),
# n_y_comm_2015 = matrix(1.02, Ncat, Ncat),
# n_y_comm_2016 = matrix(1.02, Ncat, Ncat),
# n_y_noncomm_1985 = matrix(1.03, Ncat, Ncat),
# n_y_noncomm_2002 = matrix(1.03, Ncat, Ncat),
# n_y_noncomm_2015 = matrix(1.03, Ncat, Ncat),
# n_y_noncomm_1998 = matrix(1.03, Ncat, Ncat),
# n_y_noncomm_2011 = matrix(1.03, Ncat, Ncat),
# n_y_noncomm_2016 = matrix(1.03, Ncat, Ncat),
#n = 0,
R = 1,
omega = rep_len(1/Ncat, Ncat),
theta = matrix(0.5, ncol = Ncat, nrow = Ncat),
M_comm = matrix(1/Ncat, Ncat, Ncat),
M_noncomm = matrix(1/Ncat, Ncat, Ncat),
# A_F = matrix(1/NAge, NAge, NAge),
# A_M = matrix(1/NAge, NAge, NAge),
# p = matrix(1, NAge, NAge),
# ART_RR_prog = 10, # survival extension cofactor
# ART_RR_mort = 2.5, # survival extension cofactor
# infect_ART = c(0, rep_len(0, 8)),
infect_acute = 9, # RR for acute phase
infect_AIDS = 7.27, # RR for AIDS phase
dur_FSW = 30,
# OnPrEP_init = rep_len(0, Ncat),
# have to include the follow in the function for it to work just using generate_parameters(), and not lhs_parameters()
# SC_to_death = 10,
# prev_init_FSW = 0.04,
# prev_init_rest = 0.0008,
rate_leave_pro_FSW = 0.2,
# rate_leave_client = 0.05,
# rate_leave_low_FSW = 0.1,
# prop_client_GPM = 0.2430057, # 27091/111483
# prop_pro_FSW_GPF = 0.004620494, # 672 / 145439
# prop_low_FSW_GPF = 0.005204931, # 757 / 145439
rate_move_in = matrix(0, ncol = Ncat, nrow = Ncat),
rate_move_out = rep_len(0, Ncat),
rate_move_out_PrEP = rep_len(0, Ncat),
# rate_enter_sexual_pop_F = 1,
# rate_enter_sexual_pop_M = 1,
# fraction_F = 0.51,
fraction_FSW_foreign = 0.5,
replaceDeaths = 0,
# movement = 1,
beta_comm = rep_len(0.002, Ncat),
beta_noncomm = rep_len(0.001, Ncat),
# muF = 0.02597403,
# muM = 0.02739726,
# RR_beta_FtM = 1,
# RR_beta_circum = 0.44,
# prev_HSV2_FSW = 1,
# prev_HSV2_Client = 1,
# prev_HSV2_GPF = 1,
# prev_HSV2_GPM = 1,
# RR_beta_HSV2_comm_t = 1,
# RR_beta_HSV2_noncomm_t = 1,
# RR_beta_HSV2_comm_a = 2,
# RR_beta_HSV2_noncomm_a = 2,
# betaMtoF_baseline = 0.00081,
who_believe_comm = 0,
# init_clientN_from_PCR = 0,
# fraction_sexually_active_15_F = 0,
# fraction_sexually_active_15_M = 0,
# beta_above_1 = 0,
# ignore_ranges_fc_c = 0,
# dropout_rate_not_FSW = 0.025,
# dropout_rate_FSW = 0.025,
# delete = 0,
nu = 0.02222222,
fPa = 0.5,
fPc = 0.4,
iota = c(0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0, 0, 0),
viral_supp_t = c(1985, 2015, 2016),
viral_supp_y = matrix(0, nrow = 3, ncol = 9),
# viral_supp_y_1986_rest = 0.6,
# viral_supp_y_2015_ProFSW = 0.7,
# ART_eff = 0.98,
# ART_recruit_rate = 0.2,
# ART_reinit_rate = 0.1,
tau_intervention_t = c(1986, 2015, 2016, 2017),
rho_intervention_t = c(1986, 2015, 2016, 2017),
tau_intervention_y = matrix(c(rep(0, 9), 1.5, c(rep(0, 8)), 1.5, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T),
rho_intervention_y = matrix(c(rep(0, 9), 6, c(rep(0, 8)), 6, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T),
# intervention_testing_increase = 1.4,
# intervention_ART_increase = 6.5,
# prep_offering_rate = 1,
FSW_ONLY = c(1, 1, rep(0, 7)),
# cost_ART_initiation_study_FSW = 1,
# cost_1_year_ART_study_FSW = 1,
# cost_1_year_ART_rest = 1,
W0 = 1,
W1 = 1,
W2 = 1,
W3 = 1,
cost_Initiation_of_ART_study_FSW = 1,
cost_Initiation_of_ART_government_FSW = 1,
cost_1_year_of_ART_study_FSW = 1,
cost_1_year_of_ART_government_FSW = 1,
cost_Initiation_ART_rest_of_population = 1,
cost_1_year_of_ART_rest_of_population = 1,
cost_FSW_initiation_ART_Patient_costs = 1,
cost_FSW_1_year_ART_Patient_costs = 1,
cost_Initiation_of_PrEP_study = 1,
cost_1_year_PrEP_perfect_adherence_study = 1,
cost_1_year_PrEP_intermediate_adherence_study = 1,
cost_1_year_PrEP_non_adherence_study = 1,
cost_Initiation_of_PrEP_government = 1,
cost_1_year_PrEP_perfect_adherence_government =1,
cost_1_year_PrEP_intermediate_adherence_government = 1,
cost_1_year_PrEP_non_adherence_government = 1,
cost_PREP_initiation_Patient_costs = 1,
cost_PREP_1_year_ART_Patient_costs = 1,
TasP_testing = 0,
# long_intervention = 0,
# FSW_eligible = 0,
# GP_eligible = 1,
pfFSW_t = c(1985, 1993, 2015, 2060),
pfFSW_y = matrix(c(rep(0, 9), 0.55, c(rep(0, 8)), 0.2, c(rep(0, 8)), 0.2, c(rep(0, 8))), ncol = 9, nrow = 4, byrow = T),
# prev_non_ben_fsw_1993 = 0.55,
# prev_non_ben_fsw_2015 = 0.2,
infected_FSW_incoming = 1,
prep_efficacious_y = c(0, 1, 1, 0, 0),
prep_efficacious_t = c(1985, 2015, 2017, 2017.01, 2060),
# PrEP_loss_to_follow_up = 0.1,
# ART_recruit_rate_rest = 0.1,
# ART_reinit_rate_rest = 0.1,
# ART_init_ratio_MF = 2,
# RR_beta_HSV2_a_FSW = 1,
# RR_beta_HSV2_a_client = 1,
# RR_beta_HSV2_a_GPF = 1,
# RR_beta_HSV2_a_GPM = 1,
# testing_prob_men_2006 = 0.111,
# testing_prob_men_2008 = 0.112,
# testing_prob_men_2012 = 0.113,
#
# testing_prob_women_2006 = 0.114,
# testing_prob_women_2008 = 0.115,
# testing_prob_women_2012 = 0.116,
# ART_recruit_rate_FSW = 1,
# ART_reinit_rate_FSW = 1,
PrEP_reinit_OnOff_t = c(1985, 2017, 2017.01, 2060),
PrEP_reinit_OnOff_y = c(0, 0, 0, 0),
# PrEP_reinits_on = 1,
count_PrEP_1a = 1,
count_PrEP_1b = 1,
count_PrEP_1c = 1,
re_init_interruption_parm_t = c(1985, 2017, 2017.01, 2060),
re_init_interruption_parm_y = matrix(1, nrow = 4, ncol = 9),
art_dropout_interruption_parm_t = c(1985, 2017, 2017.01, 2060),
art_dropout_interruption_parm_y = matrix(1, nrow = 4, ncol = 9),
art_initiation_interruption_parm_t = c(1985, 2017, 2017.01, 2060),
art_initiation_interruption_parm_y = matrix(1, nrow = 4, ncol = 9),
old_VS_assumption = 1
)
if (length(parameters) == 0L) {
return(defaults)
}
if (is.null(names(parameters)) || !all(nzchar(names(parameters)))) {
stop("All arguments must be named")
}
# extra <- setdiff(names(parameters), names(defaults))
# if (length(extra) > 0L) {
# stop("Unknown arguments: ", extra)
# }
# list of parameters that depend on others
ret <- modifyList(defaults, parameters)
if (length(set_null) > 0L) {
ret = modifyList(ret, lapply(ret[match(set_null, names(ret))], function(x) x*0))
}
# browser()
ret <- ret[which(!names(ret) %in% c("N_init", "initial_Ntot", "frac_women_ProFSW", "frac_women_LowFSW", "frac_women_exFSW", "frac_men_client", "frac_women_virgin", "frac_men_virgin", "epsilon_1985", "epsilon_1992", "epsilon_2002", "epsilon_2013", "epsilon_2016", "dur_primary_phase", "dur_200_349", "gamma32", "gamma33", "gamma34", "ART_prob_t", "ART_prob_y", "RR_test_onPrEP", "RR_ART_CD4200", "tau01", "tau11", "tau2", "tau3", "tau4", "tau5", "zetaa_t", "zetaa_y", "zetab_t", "zetab_y", "zetac_t", "zetac_y", "c_comm_1985", "c_comm_1993", "c_comm_1995", "c_comm_1998", "c_comm_2002", "c_comm_2005", "c_comm_2008", "c_comm_2012", "c_comm_2015", "c_comm_2016", "c_noncomm_1985", "c_noncomm_1993", "c_noncomm_1995", "c_noncomm_1998", "c_noncomm_2002", "c_noncomm_2005", "c_noncomm_2008", "c_noncomm_2012", "c_noncomm_2015", "c_noncomm_2016", "fc_y_comm_1985", "fc_y_comm_1993", "fc_y_comm_1995", "fc_y_comm_1998", "fc_y_comm_2002", "fc_y_comm_2005", "fc_y_comm_2008", "fc_y_comm_2012", "fc_y_comm_2015", "fc_y_comm_2016", "fc_y_noncomm_1985", "fc_y_noncomm_1993", "fc_y_noncomm_1998", "fc_y_noncomm_2002", "fc_y_noncomm_2008", "fc_y_noncomm_20"))]
ret <- ret[which(!names(ret) %in% c("ART_eff", "ART_recruit_rate_FSW", "ART_recruit_rate_rest", "ART_reinit_rate_FSW", "ART_reinit_rate_rest", "ART_RR_prog", "betaMtoF_baseline", "c_comm_1993_ProFSW", "c_comm_1998_Client", "c_comm_2005_ProFSW", "c_comm_2012_Client", "c_comm_2015_Client", "c_non_comm_1985_Client", "c_non_comm_1985_LowFSW", "c_non_comm_1985_ProFSW", "c_noncomm_1998_GPF", "c_noncomm_1998_GPM", "c_noncomm_2008_GPF", "c_noncomm_2008_GPM", "dropout_rate_FSW", "dropout_rate_not_FSW", "dur_below_200", "eff_ART", "fc_y_comm_1985_ProFSW_Client", "fc_y_comm_1993_ProFSW_Client", "fc_y_comm_2002_ProFSW_Client", "fc_y_noncomm_1985_GPF_GPM", "fc_y_noncomm_1985_ProFSW_Client", "fc_y_noncomm_1998_GPF_GPM", "fc_y_noncomm_2002_ProFSW_Client", "fc_y_noncomm_2011_GPF_GPM", "fraction_F", "fraction_sexually_active_15_F", "fraction_sexually_active_15_M", "init_clientN_from_PCR", "muF", "muM", "n_y_comm_1985_Client_LowFSW", "n_y_comm_1985_Client_ProFSW", "n_y_comm_1985_LowFSW_Client", "n_y_comm_1985_ProFSW_Client", "n_y_noncomm_1985_GPF_GPM", "n_y_noncomm_1985_GPM_GPF"))]
ret <- ret[which(!names(ret) %in% c(unlist(strsplit("n_y_noncomm_2002_ProFSW_Client, n_y_noncomm_2015_ProFSW_Client, prev_HSV2_Client, prev_HSV2_FSW, prev_HSV2_GPF, prev_HSV2_GPM, prev_init_FSW, prev_init_rest, rate_enter_sexual_pop_F, rate_enter_sexual_pop_M, rate_leave_client, rate_leave_low_FSW, RR_beta_circum, RR_beta_FtM, RR_beta_HSV2_comm, RR_beta_HSV2_noncomm, SC_to_death, viral_supp_y_1986_rest, viral_supp_y_2015_ProFSW, movement, fc_y_noncomm_2011, fc_y_noncomm_2015, fc_y_noncomm_2016, n_y_comm_1985, n_y_comm_2002, n_y_comm_2015, n_y_comm_2016, n_y_noncomm_1985, n_y_noncomm_1998, n_y_noncomm_2011, n_y_noncomm_2002, n_y_noncomm_2015, n_y_noncomm_2016, betaMtoF_noncomm, betaFtoM_noncomm, betaMtoF_comm, betaFtoM_comm, RR_beta_HSV2_comm_t, RR_beta_HSV2_noncomm_t, RR_beta_HSV2_comm_a, RR_beta_HSV2_noncomm_a, beta_above_1, ignore_ranges_fc_c, delete, ART_recruit_rate, ART_reinit_rate, prep_dropout", ", "))))]
ret <- ret[which(!names(ret) %in% c(unlist(strsplit("cost_ART_initiation_study_FSW, cost_1_year_ART_study_FSW, cost_1_year_ART_rest, long_intervention, FSW_eligible, GP_eligible, prev_non_ben_fsw_1993, prev_non_ben_fsw_2015, ART_init_ratio_MF, RR_beta_HSV2_a_FSW, RR_beta_HSV2_a_client, RR_beta_HSV2_a_GPF, RR_beta_HSV2_a_GPM, testing_prob_men_2006, testing_prob_men_2008, testing_prob_men_2012, testing_prob_women_2006, testing_prob_women_2008, testing_prob_women_2012, PrEP_reinits_on, infect_ART, alpha33, alpha34, alpha35, FSW_leave_Cotonou_fraction, time, PrEP_loss_to_follow_up, n_y_noncomm_2002_LowFSW_Client, n_y_noncomm_2015_LowFSW_Client, n_y_noncomm_1985_GPF_Client, fc_y_noncomm_1985_GPF_Client, fc_y_noncomm_1998_GPF_Client, fc_y_noncomm_2011_GPF_Client, betaMtoF_FSW, betaFtoM_client, betaMtoF_GPF, betaFtoM_GPM, ART_RR_mort, prop_client_GPM, prop_pro_FSW_GPF, prop_low_FSW_GPF", ", "))))]
ret
}
Geidelberg/cotonou documentation built on Aug. 29, 2020, 4:22 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions generate_parameters lhs_parameters_parallel delete_parameter_sets lhs_parameters fix_parameters FormerFSWtoGPF
Documented in fix_parameters
#' @export
#' @useDynLib cotonou
FormerFSWtoGPF <- function(x) {
if(!is.null(dim(x))) {
x[4,] = x[3,]
x[,4] = x[,3]
}
return(x)
}
#' parameters which depend on others, etc
#'
#' Lily's function to be documented
#' @export
#' @useDynLib cotonou
fix_parameters <- function(y, Ncat, Nage, par_seq, condom_seq, groups_seq, years_seq) {
useless_parameter = 1
# why not working
if(y$ignore_ranges_fc_c == 0) {
# y$PrEP_reinit_OnOff_y = c(0, 0, y$PrEP_reinits_on, y$PrEP_reinits_on)
if(y$prep_dropout - y$rate_leave_pro_FSW - y$muF - 1/45 - 0.008 > 0)
{
y$PrEP_loss_to_follow_up <- y$prep_dropout - y$rate_leave_pro_FSW - y$muF - 1/45 - 0.008
} else {
y$PrEP_loss_to_follow_up = 0
}
y$frac_women_exFSW = y$frac_women_ProFSW
# getting client - low FSW n to be same as for pro fsw
# y$n_y_comm_1985_LowFSW_Client = y$n_y_comm_1985_ProFSW_Client
y$n_y_noncomm_2002_LowFSW_Client = y$n_y_noncomm_2002_ProFSW_Client
y$n_y_noncomm_2015_LowFSW_Client = y$n_y_noncomm_2015_ProFSW_Client
y$n_y_noncomm_1985_GPF_Client = y$n_y_noncomm_1985_GPF_GPM
y$n_y_noncomm_1998_GPF_Client = y$n_y_noncomm_1998_GPF_GPM
y$n_y_noncomm_2011_GPF_Client = y$n_y_noncomm_2011_GPF_GPM
# getting low fsw condom with client for comm and noncomm partnerships same as pro fsw
# y$fc_y_comm_1985_LowFSW_Client = y$fc_y_comm_1985_ProFSW_Client
# y$fc_y_comm_1993_LowFSW_Client = y$fc_y_comm_1993_ProFSW_Client
# y$fc_y_comm_1998_LowFSW_Client = y$fc_y_comm_1998_ProFSW_Client
# y$fc_y_comm_2002_LowFSW_Client = y$fc_y_comm_2002_ProFSW_Client
# y$fc_y_noncomm_1985_LowFSW_Client = y$fc_y_noncomm_1985_ProFSW_Client
# y$fc_y_noncomm_2002_LowFSW_Client = y$fc_y_noncomm_2002_ProFSW_Client
# getting client condom with GPF same as GPF-GPM
y$fc_y_noncomm_1985_GPF_Client = y$fc_y_noncomm_1985_GPF_GPM
y$fc_y_noncomm_1998_GPF_Client = y$fc_y_noncomm_1998_GPF_GPM
y$fc_y_noncomm_2011_GPF_Client = y$fc_y_noncomm_2011_GPF_GPM
# CONDOMS
what_we_got_condom = c()
for(year in 1:length(years_seq))
{for(group in 1:length(groups_seq))
{for(group2 in 1:length(groups_seq))
{for(par in 1:length(condom_seq))
{if(paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group], "_", groups_seq[group2]) %in% names(y)){
# average if the opposite value is there too - note must be same year...
if(paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group2], "_", groups_seq[group]) %in% names(y)) {
av=(y[[which(paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group], "_", groups_seq[group2]) == names(y))]] +
y[[which(paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group2], "_", groups_seq[group]) == names(y))]]) / 2;
y[[which(paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group], "_", groups_seq[group2]) == names(y))]] = av;
y[[which(paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group2], "_", groups_seq[group]) == names(y))]] = av;
}
# replacing the non-group specific parameter to include these precisions
y[paste0(condom_seq[par], "_", years_seq[year])][[paste0(condom_seq[par], "_", years_seq[year])]][group, group2] = y[paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group], "_", groups_seq[group2])][[1]]
y[paste0(condom_seq[par], "_", years_seq[year])][[paste0(condom_seq[par], "_", years_seq[year])]][group2, group] = y[paste0(condom_seq[par], "_", years_seq[year])][[paste0(condom_seq[par], "_", years_seq[year])]][group, group2]
what_we_got_condom = rbind(what_we_got_condom, c(par, year, group, group2))
}}}}}
if (length(what_we_got_condom[,1]) > 0)
{
colnames(what_we_got_condom) = c("par", "year", "group", "group2")
# now we have to fill in the rest of the years... IF 2 OR MORE YEARS FOR SAME GROUP AND PARM
# all years before earliest one is same as earliest estimate,
# and all years after last one is same as last...
# and everything in between is interpolated!
d <- data.frame(what_we_got_condom)
par_counts = plyr::count(d, c('par', 'group', 'group2'))
for(i in 1:length(par_counts[,1]))
{
if(par_counts[i,"freq"] > 0)
{
the_years = subset(d, c(par == par_counts[i, "par"] & group == par_counts[i, "group"] & group2 == par_counts[i, "group2"]), year)
#before
if(min(the_years) > 1)
{
for(year in 1:(min(the_years)-1))
{
if(paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year]) %in% names(y))
{y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group"], par_counts[i, "group2"]] =
y[paste0(condom_seq[par_counts[i,"par"]], "_", years_seq[min(the_years)], "_", groups_seq[par_counts[i, "group"]], "_", groups_seq[par_counts[i, "group2"]])][[1]]
# equalising the complementary condom use
y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group2"], par_counts[i, "group"]] =
y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group"], par_counts[i, "group2"]]
}
}
}
#after
if(max(the_years) < length(years_seq))
{
for(year in (max(the_years)+1):length(years_seq))
{
if(paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year]) %in% names(y))
{y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group"], par_counts[i, "group2"]] =
y[paste0(condom_seq[par_counts[i,"par"]], "_", years_seq[max(the_years)], "_", groups_seq[par_counts[i, "group"]], "_", groups_seq[par_counts[i, "group2"]])][[1]]
# equalising the complementary condom use
y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group2"], par_counts[i, "group"]] =
y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group"], par_counts[i, "group2"]]
}
}
}
# interpolate between
# take the years that have values
# calculate the slopes between them
# apply the slope and time difference to the pars in between
if(par_counts[i,"freq"] > 1) {
for(j in 1:(length(unlist(the_years))-1)) {
slope = (y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j+1])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j+1])]][par_counts[i, "group"], par_counts[i, "group2"]] -
y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])]][par_counts[i, "group"], par_counts[i, "group2"]]) /
(years_seq[unlist(the_years)][j+1] - years_seq[unlist(the_years)][j])
for(k in (years_seq[unlist(the_years)][j]+1):(years_seq[unlist(the_years)][j+1]-1))
{
if(paste0(condom_seq[par_counts[i, "par"]], "_", k) %in% names(y))
{ y[paste0(condom_seq[par_counts[i, "par"]], "_", k)][[paste0(condom_seq[par_counts[i, "par"]], "_", k)]][par_counts[i, "group"], par_counts[i, "group2"]] = slope *
(k - years_seq[unlist(the_years)][j]) +
y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])]][par_counts[i, "group"], par_counts[i, "group2"]]
# equalising complementary condom rate
y[paste0(condom_seq[par_counts[i, "par"]], "_", k)][[paste0(condom_seq[par_counts[i, "par"]], "_", k)]][par_counts[i, "group2"], par_counts[i, "group"]] =
y[paste0(condom_seq[par_counts[i, "par"]], "_", k)][[paste0(condom_seq[par_counts[i, "par"]], "_", k)]][par_counts[i, "group"], par_counts[i, "group2"]]}
}
}
}
}
}
}
# MAKING FORMER FSW IN COTONOU EQUAL CONDOM USE TO GPF _ NOTE DO I NEED TO DO THIS FOR ALL GPF PARAMETERS???
y[names(y)[grep("fc_y_", names(y))]] = lapply(invisible(y[names(y)[grep("fc_y_", names(y))]]), FormerFSWtoGPF)
y[names(y)[grep("n_y_", names(y))]] = lapply(invisible(y[names(y)[grep("n_y_", names(y))]]), FormerFSWtoGPF)
##########################################################################
# c_comm and c_noncomm
# what_we_got = data.frame(par = c(), year = c(), group = c())
what_we_got = c()
for(year in 1:length(years_seq))
{for(group in 1:length(groups_seq))
{for(par in 1:length(par_seq))
# {if(grepl(years_seq[year], names(y)) && grepl(groups_seq[group], names(y)) && grepl(par_seq[par], names(y))){
{if(paste0(paste0(par_seq[par], "_", years_seq[year], "_", groups_seq[group])) %in% names(y)){
y[paste0(par_seq[par], "_", years_seq[year])][[paste0(par_seq[par], "_", years_seq[year])]][group] = y[paste0(par_seq[par], "_", years_seq[year], "_", groups_seq[group])][[1]]
# what_we_got = c(what_we_got, paste0(par_seq[par], "_", years_seq[year], "_", groups_seq[group]))
what_we_got = rbind(what_we_got, c(par, year, group))
# print(c(par_seq[par], groups_seq[group], years_seq[year]))
}}}}
# print(list(what_we_got, years_seq, groups_seq, par_seq, what_we_got_condom))
if(length(what_we_got) > 0)
{colnames(what_we_got) = c("par", "year", "group")
# now we have to fill in the rest of the years... IF 2 OR MORE YEARS FOR SAME GROUP AND PARM
# all years before earliest one is same as earliest estimate,
# and all years after last one is same as last...
# and everything in between is interpolated!
d <- data.frame(what_we_got)
par_counts = plyr::count(d, c('par', 'group'))
for(i in 1:length(par_counts[,1]))
{
if(par_counts[i,"freq"] > 0)
{
the_years = subset(d, c(par == par_counts[i, "par"] & group == par_counts[i, "group"]), year)
#before
if(min(the_years) > 1)
{
for(year in 1:(min(the_years)-1))
{
if(paste0(par_seq[par_counts[i, "par"]], "_", years_seq[year]) %in% names(y))
{y[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group"]] =
y[paste0(par_seq[par_counts[i,"par"]], "_", years_seq[min(the_years)], "_", groups_seq[par_counts[i, "group"]])][[1]]}
}
}
#after
if(max(the_years) < length(years_seq))
{
for(year in (max(the_years)+1):length(years_seq))
{
if(paste0(par_seq[par_counts[i, "par"]], "_", years_seq[year]) %in% names(y))
{y[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group"]] =
y[paste0(par_seq[par_counts[i,"par"]], "_", years_seq[max(the_years)], "_", groups_seq[par_counts[i, "group"]])][[1]]}
}
}
# interoplate between
# take the years that have values
# calculate the slopes between them
# apply the slope and time difference to the pars in between
if(par_counts[i,"freq"] > 1) {
for(j in 1:(length(unlist(the_years))-1)) {
slope = (y[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j+1])][[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j+1])]][par_counts[i, "group"]] -
y[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])][[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])]][par_counts[i, "group"]]) /
(years_seq[unlist(the_years)][j+1] - years_seq[unlist(the_years)][j])
for(k in (years_seq[unlist(the_years)][j]+1):(years_seq[unlist(the_years)][j+1]-1))
{
if(paste0(par_seq[par_counts[i, "par"]], "_", k) %in% names(y))
y[paste0(par_seq[par_counts[i, "par"]], "_", k)][[paste0(par_seq[par_counts[i, "par"]], "_", k)]][par_counts[i, "group"]] = slope *
(k - years_seq[unlist(the_years)][j]) +
y[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])][[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])]][par_counts[i, "group"]]
}
}
}
}
}
}
}
# print("after")
# print(y)
# number of clients at initial conditions shaped by commercial partner change rates and the N of Pro FSW!
# partnerships offered by ProFSW / P.C.R of clients / fraction of men * fraction of men that are clients
if(y$init_clientN_from_PCR == 1)
y$frac_men_client = y$fraction_F * y$frac_women_ProFSW * y$initial_Ntot * y$c_comm_1985[1] / y$c_comm_1985[5] / (y$initial_Ntot * (1-y$fraction_F))
# fractions of groups at initialisation
y$N_init = y$initial_Ntot*c(y$fraction_F*y$frac_women_ProFSW,
y$fraction_F*y$frac_women_ProFSW*y$frac_women_LowFSW,
y$fraction_F - (y$fraction_F*y$frac_women_ProFSW + y$fraction_F*y$frac_women_ProFSW*y$frac_women_LowFSW + y$fraction_F*y$frac_women_virgin + y$fraction_F*y$frac_women_exFSW),
y$fraction_F*y$frac_women_exFSW,
(1-y$fraction_F)*y$frac_men_client,
(1 - y$fraction_F - ((1-y$fraction_F)*y$frac_men_client + (1-y$fraction_F) * y$frac_men_virgin)),
y$fraction_F*y$frac_women_virgin,
(1-y$fraction_F) * y$frac_men_virgin,
0
)
y$betaMtoF_FSW = y$betaMtoF_baseline * (1 + (y$RR_beta_HSV2_a_FSW - 1) * y$prev_HSV2_FSW)
y$betaFtoM_client = y$betaMtoF_baseline * y$RR_beta_FtM * y$RR_beta_circum * (1 + (y$RR_beta_HSV2_a_client - 1) * y$prev_HSV2_Client)
y$betaMtoF_GPF = y$betaMtoF_baseline * (1 + (y$RR_beta_HSV2_a_GPF - 1) * y$prev_HSV2_GPF)
y$betaFtoM_GPM = y$betaMtoF_baseline * y$RR_beta_FtM * y$RR_beta_circum * (1 + (y$RR_beta_HSV2_a_GPM - 1) * y$prev_HSV2_GPM)
# y$betaMtoF_noncomm_FSW = y$betaMtoF_baseline * (1 + (y$RR_beta_HSV2_noncomm_a - 1) * y$prev_HSV2_FSW + (y$RR_beta_HSV2_noncomm_t - 1) * y$prev_HSV2_Client)
# y$betaFtoM_noncomm_client = y$betaMtoF_baseline * y$RR_beta_FtM * y$RR_beta_circum * (1 + (y$RR_beta_HSV2_noncomm_a - 1) * y$prev_HSV2_Client + (y$RR_beta_HSV2_noncomm_t - 1) * y$prev_HSV2_FSW)
#
# y$betaMtoF_noncomm = y$betaMtoF_baseline * (1 + (y$RR_beta_HSV2_noncomm_a - 1) * y$prev_HSV2_GPF + (y$RR_beta_HSV2_noncomm_t - 1) * y$prev_HSV2_GPM)
# y$betaFtoM_noncomm = y$betaMtoF_baseline * y$RR_beta_FtM * y$RR_beta_circum * (1 + (y$RR_beta_HSV2_noncomm_a - 1) * y$prev_HSV2_GPM + (y$RR_beta_HSV2_noncomm_t - 1) * y$prev_HSV2_GPF)
#
# y$betaMtoF_comm = y$betaMtoF_baseline * (1 + (y$RR_beta_HSV2_comm_a - 1) * y$prev_HSV2_FSW + (y$RR_beta_HSV2_comm_t - 1) * y$prev_HSV2_Client)
# y$betaFtoM_comm = y$betaMtoF_baseline * y$RR_beta_FtM * y$RR_beta_circum * (1 + (y$RR_beta_HSV2_comm_a - 1) * y$prev_HSV2_Client + (y$RR_beta_HSV2_comm_t - 1) * y$prev_HSV2_FSW)
# if any beta becomes > 1, then make them all zero and flag it
# if(y$betaMtoF_noncomm_FSW * y$infect_acute >= 1 || y$betaFtoM_noncomm_client * y$infect_acute >= 1 || y$betaMtoF_noncomm * y$infect_acute >= 1 || y$betaMtoF_comm * y$infect_acute >= 1 || y$betaFtoM_noncomm * y$infect_acute >= 1 || y$betaFtoM_comm * y$infect_acute >= 1)
if(y$betaMtoF_FSW * y$infect_acute >= 1 || y$betaFtoM_client * y$infect_acute >= 1 || y$betaMtoF_GPF * y$infect_acute >= 1 || y$betaFtoM_GPM * y$infect_acute >= 1)
{
# y$betaMtoF_noncomm_FSW = 0
# y$betaFtoM_noncomm_client = 0
# y$betaMtoF_noncomm = 0
# y$betaFtoM_noncomm = 0
# y$betaMtoF_comm = 0
# y$betaFtoM_comm = 0
y$betaMtoF_FSW = 0
y$betaFtoM_client = 0
y$betaMtoF_GPF = 0
y$betaFtoM_GPM = 0
y$beta_above_1 = 1
}
y$epsilon_y = c(y$epsilon_1985, y$epsilon_1992, y$epsilon_2002, y$epsilon_2013, y$epsilon_2016)
init_prev = if(Ncat == 9) c(y$prev_init_FSW, rep_len(y$prev_init_rest, 5), 0, 0, 0) else c(y$prev_init_FSW, rep_len(y$prev_init_rest, Ncat-1))
y$S0_init = (1-init_prev) * y$N_init
y$I01_init = init_prev * y$N_init
N = y$S0_init + y$I01_init
# BIOLOGICAL
#dur_primary_phase, dur_200_349 input as a DURATION
#SC_to_death input as a duration
y$ART_RR_mort = y$ART_RR_prog
y$ART_RR = y$ART_RR_prog
y$gamma01 = 1/y$dur_primary_phase
y$gamma02 = 2/(y$SC_to_death - y$dur_primary_phase - y$dur_200_349 - y$dur_below_200)
y$gamma03 = y$gamma02
y$gamma04 = 1/y$dur_200_349
y$gamma01 <- rep_len(y$gamma01, Ncat)
y$gamma02 <- rep_len(y$gamma02, Ncat)
y$gamma03 <- rep_len(y$gamma03, Ncat)
y$gamma04 <- rep_len(y$gamma04, Ncat)
y$gamma11 <- y$gamma01 # loss to follow up same as normal movement?
y$gamma22 <- y$gamma02
y$gamma23 <- y$gamma03
y$gamma24 <- y$gamma04
y$gamma42 <- y$gamma02
y$gamma43 <- y$gamma03
y$gamma44 <- y$gamma04
# progression is slowed by ART_RR
y$gamma32_without_supp <- (y$gamma02)/y$ART_RR_prog
y$gamma33_without_supp <- (y$gamma03)/y$ART_RR_prog
y$gamma34_without_supp <- (y$gamma04)/y$ART_RR_prog
# mortality
y$alpha03 <- rep_len(y$alpha03, Ncat)
y$alpha04 <- rep_len(y$alpha04, Ncat)
y$alpha05 <- 1/y$dur_below_200
y$alpha05 <- rep_len(y$alpha05, Ncat)
y$alpha23 <- y$alpha03
y$alpha24 <- y$alpha04
y$alpha25 <- y$alpha05
y$alpha43 <- y$alpha03
y$alpha44 <- y$alpha04
y$alpha45 <- y$alpha05
y$alpha33_without_supp <- (y$alpha03)/y$ART_RR_mort
y$alpha34_without_supp <- (y$alpha04)/y$ART_RR_mort
y$alpha35_without_supp <- (y$alpha05)/y$ART_RR_mort
# PCR
y$c_t_comm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016)
y$c_t_noncomm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016)
y$c_y_comm = rbind(y$c_comm_1985, y$c_comm_1993, y$c_comm_1995, y$c_comm_1998, y$c_comm_2002,
y$c_comm_2005, y$c_comm_2008, y$c_comm_2012, y$c_comm_2015, y$c_comm_2016)
y$c_y_noncomm = rbind(y$c_noncomm_1985, y$c_noncomm_1993, y$c_noncomm_1995, y$c_noncomm_1998, y$c_noncomm_2002,
y$c_noncomm_2005, y$c_noncomm_2008, y$c_noncomm_2012, y$c_noncomm_2015, y$c_noncomm_2016)
y$fc_y_comm = array(data = c(y$fc_y_comm_1985, y$fc_y_comm_1993,
y$fc_y_comm_1995, y$fc_y_comm_1998,
y$fc_y_comm_2002, y$fc_y_comm_2005,
y$fc_y_comm_2008, y$fc_y_comm_2012,
y$fc_y_comm_2015, y$fc_y_comm_2016), dim=c(Ncat, Ncat, 10))
y$fc_y_comm = aperm(y$fc_y_comm, c(3, 1, 2))
y$fc_y_noncomm = array(data = c(y$fc_y_noncomm_1985, y$fc_y_noncomm_1993, y$fc_y_noncomm_1998, y$fc_y_noncomm_2002,
y$fc_y_noncomm_2008, y$fc_y_noncomm_2011, y$fc_y_noncomm_2015,
y$fc_y_noncomm_2016), dim=c(Ncat, Ncat, 8))
y$fc_y_noncomm = aperm(y$fc_y_noncomm, c(3, 1, 2))
y$n_y_comm = array(data = c(y$n_y_comm_1985, y$n_y_comm_2002,
y$n_y_comm_2015, y$n_y_comm_2016), dim=c(Ncat, Ncat, 4))
y$n_y_noncomm = array(data = c(y$n_y_noncomm_1985, y$n_y_noncomm_1998, y$n_y_noncomm_2002, y$n_y_noncomm_2011,
y$n_y_noncomm_2015, y$n_y_noncomm_2016), dim=c(Ncat, Ncat, 6))
y$n_y_comm = aperm(y$n_y_comm, c(3, 1, 2))
y$n_y_noncomm = aperm(y$n_y_noncomm, c(3, 1, 2))
y$omega = N/sum(N)
# when sampling pars, sometimes instead of length Ncat, I only sampled length one
if(length(y$ec) == 1)
y$ec = rep_len(y$ec, 9)
if(length(y$eP0) == 1)
y$eP0 = c(y$eP0, rep_len(0, 8))
if(length(y$eP1a) == 1)
y$eP1a = c(y$eP1a, rep_len(0, 8))
if(length(y$eP1b) == 1)
y$eP1b = c(y$eP1b, rep_len(0, 8))
if(length(y$eP1c) == 1)
y$eP1c = c(y$eP1c, rep_len(0, 8))
if(length(y$eP1d) == 1)
y$eP1d = c(y$eP1d, rep_len(0, 8))
if(length(y$psia) == 1)
y$psia = c(y$psia, rep_len(0, 8))
if(length(y$psib) == 1)
y$psib = c(y$psib, rep_len(0, 8))
y$phi2 = c(y$dropout_rate_FSW, rep(y$dropout_rate_not_FSW, (Ncat-1)))
y$phi3 = y$phi2
y$phi4 = y$phi2
y$phi5 = y$phi2
if (Ncat == 9) {
y$who_believe_comm = round(y$who_believe_comm)
y$omega = c(
y$fraction_F * y$frac_women_ProFSW * y$fraction_FSW_foreign, # some FSW come from outside Cotonou
y$fraction_F * y$frac_women_LowFSW * y$frac_women_ProFSW * y$fraction_FSW_foreign,
y$fraction_F - y$fraction_F * (1-y$fraction_sexually_active_15_F) - y$fraction_F * y$frac_women_ProFSW * y$fraction_FSW_foreign - y$fraction_F * y$frac_women_LowFSW * y$frac_women_ProFSW * y$fraction_FSW_foreign,
0,
0,
(1 - y$fraction_F) - (1 - y$fraction_F) * (1 - y$fraction_sexually_active_15_M),
y$fraction_F * (1-y$fraction_sexually_active_15_F),
(1 - y$fraction_F) * (1 - y$fraction_sexually_active_15_M),
0)
#
# MIXING
###############################################
y$M_noncomm = matrix(c(0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 1, 0, 0, 0,
0, 0, 0, 0, 1, 1, 0, 0, 0,
1, 1, 1, 1, 0, 0, 0, 0, 0,
0, 0, 1, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0),
nrow = 9, ncol = 9, byrow = T)
y$M_comm = matrix(c(0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
1, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0),
nrow = 9, ncol = 9, byrow = T)
# MOVEMENT
###############################################
y$rate_leave_low_FSW = y$rate_leave_pro_FSW
# virgin movement
y$rate_move_out[7] = - y$rate_enter_sexual_pop_F
y$rate_move_out[8] = - y$rate_enter_sexual_pop_M
y$rate_move_in[3,7] = y$rate_enter_sexual_pop_F
y$rate_move_in[6,8] = y$rate_enter_sexual_pop_M
#this is just to show what happens when you increase movement
# y$rate_leave_pro_FSW = y$rate_leave_pro_FSW * 10
# y$rate_leave_low_FSW = y$rate_leave_low_FSW * 10
# y$rate_leave_client = y$rate_leave_client * 10
y$prop_client_GPM = y$N_init[5] / y$N_init[6]
y$prop_pro_FSW_GPF = y$N_init[1] / y$N_init[3]
y$prop_low_FSW_GPF = y$N_init[2] / y$N_init[3]
# FEMALE MOVEMENT
y$rate_move_out[1] = - y$rate_leave_pro_FSW
y$rate_move_out[2] = - y$rate_leave_low_FSW
y$rate_move_out[3] = - (y$rate_leave_pro_FSW + y$muF + y$nu) * y$prop_pro_FSW_GPF * (1-y$fraction_FSW_foreign) -
(y$rate_leave_low_FSW + y$muF + y$nu) * y$prop_low_FSW_GPF* (1-y$fraction_FSW_foreign)
y$rate_move_in[1,3] = (y$rate_leave_pro_FSW + y$muF + y$nu) * y$prop_pro_FSW_GPF * (1-y$fraction_FSW_foreign) # moving from GPF to pro-FSW
y$rate_move_in[2,3] = (y$rate_leave_low_FSW + y$muF + y$nu) * y$prop_low_FSW_GPF * (1-y$fraction_FSW_foreign) # moving from GPF to low-FSW
y$rate_move_in[4,1] = y$rate_leave_pro_FSW * (1 - y$fraction_FSW_foreign) # moving from pro-FSW to former FSW in Cot
y$rate_move_in[4,2] = y$rate_leave_low_FSW * (1 - y$fraction_FSW_foreign) # moving from low-FSW to former FSW in Cot
y$rate_move_in[9,1] = y$rate_leave_pro_FSW * y$fraction_FSW_foreign # moving from low-FSW to former FSW NOT in Cot
y$rate_move_in[9,2] = y$rate_leave_low_FSW * y$fraction_FSW_foreign # moving from low-FSW to former FSW NOT in Cot
# MALE MOVEMENT
y$rate_move_out[5] = - y$rate_leave_client
y$rate_move_out[6] = - y$prop_client_GPM * (y$rate_leave_client + y$muM + y$nu)
y$rate_move_in[6,5] = y$rate_leave_client # moving from client to GPM
y$rate_move_in[5,6] = y$prop_client_GPM * (y$rate_leave_client + y$muM + y$nu) # moving from GPM to client
y$rate_move_out_PrEP = y$rate_move_out
# y$rate_move_out_PrEP[1] = -y$rate_move_out_PrEP_FSW
y$beta_comm = c(y$betaMtoF_FSW, y$betaMtoF_GPF, y$betaMtoF_GPF, y$betaMtoF_GPF, y$betaFtoM_client, y$betaFtoM_GPM, 0, 0, 0)
y$beta_noncomm = c(y$betaMtoF_FSW, y$betaMtoF_GPF, y$betaMtoF_GPF, y$betaMtoF_GPF, y$betaFtoM_client, y$betaFtoM_GPM, 0, 0, 0)
# c_t_comm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016),
# y$c_y_comm =
y$mu = c(y$muF, y$muF, y$muF, y$muF, y$muM, y$muM, y$muF, y$muM, y$muF)
} else {
# y$omega = y$omega/sum(y$omega)
}
if(y$movement == 0) {
y$rate_move_in = y$rate_move_in * 0
y$rate_move_out = y$rate_move_out * 0}
# extending all interpolating parameters beyond 2016 to the max of time
y$c_y_comm <- rbind(y$c_y_comm, y$c_y_comm[nrow(y$c_y_comm),])
y$c_t_comm <- c(y$c_t_comm, max(y$time))
y$c_y_noncomm <- rbind(y$c_y_noncomm, y$c_y_noncomm[nrow(y$c_y_noncomm),])
y$c_t_noncomm <- c(y$c_t_noncomm, max(y$time))
y$fc_y_comm <- abind::abind(y$fc_y_comm, y$fc_y_comm[length(y$fc_y_comm[,1,1]),,], along = 1)
y$fc_t_comm <- c(y$fc_t_comm, max(y$time))
y$fc_y_noncomm <- abind::abind(y$fc_y_noncomm, y$fc_y_noncomm[length(y$fc_y_noncomm[,1,1]),,], along = 1)
y$fc_t_noncomm <- c(y$fc_t_noncomm, max(y$time))
y$n_y_comm <- abind::abind(y$n_y_comm, y$n_y_comm[length(y$n_y_comm[,1,1]),,], along = 1)
y$n_t_comm <- c(y$n_t_comm, max(y$time))
y$n_y_noncomm <- abind::abind(y$n_y_noncomm, y$n_y_noncomm[length(y$n_y_noncomm[,1,1]),,], along = 1)
y$n_t_noncomm <- c(y$n_t_noncomm, max(y$time))
y$zetaa_t <- c(y$zetaa_t, max(y$time))
y$zetab_t <- c(y$zetab_t, max(y$time))
y$zetac_t <- c(y$zetac_t, max(y$time))
y$zetaa_y <- rbind(y$zetaa_y, y$zetaa_y[nrow(y$zetaa_y),])
y$zetab_y <- rbind(y$zetab_y, y$zetab_y[nrow(y$zetab_y),])
y$zetac_y <- rbind(y$zetac_y, y$zetac_y[nrow(y$zetac_y),])
y$testing_prob_y[which(y$testing_prob_t == 2006), c(2,3,4)] <- y$testing_prob_women_2006
y$testing_prob_y[which(y$testing_prob_t == 2008), c(2,3,4)] <- y$testing_prob_women_2008
y$testing_prob_y[which(y$testing_prob_t == 2012), c(2,3,4)] <- y$testing_prob_women_2012
y$testing_prob_y[which(y$testing_prob_t == 2013), c(2,3,4)] <- y$testing_prob_women_2012
y$testing_prob_y[which(y$testing_prob_t == 2015), c(2,3,4)] <- y$testing_prob_women_2012
y$testing_prob_y[which(y$testing_prob_t == 2016), c(2,3,4)] <- y$testing_prob_women_2012
y$testing_prob_y[which(y$testing_prob_t == 2006), c(5,6)] <- y$testing_prob_men_2006
y$testing_prob_y[which(y$testing_prob_t == 2008), c(5,6)] <- y$testing_prob_men_2008
y$testing_prob_y[which(y$testing_prob_t == 2012), c(5,6)] <- y$testing_prob_men_2012
y$testing_prob_y[which(y$testing_prob_t == 2013), c(5,6)] <- y$testing_prob_men_2012
y$testing_prob_y[which(y$testing_prob_t == 2015), c(5,6)] <- y$testing_prob_men_2012
y$testing_prob_y[which(y$testing_prob_t == 2016), c(5,6)] <- y$testing_prob_men_2012
y$testing_prob_y <- rbind(y$testing_prob_y, y$testing_prob_y[nrow(y$testing_prob_y),])
y$testing_prob_t <- c(y$testing_prob_t, max(y$time))
y$ART_prob_y <- rbind(y$ART_prob_y, y$ART_prob_y[nrow(y$ART_prob_y),])
y$ART_prob_t <- c(y$ART_prob_t, max(y$time))
y$fPb = 1 - y$fPa - y$fPc
# ART and efficacy
y$viral_supp_y[1,] = y$viral_supp_y_1986_rest
y$viral_supp_y[2,] = c(y$viral_supp_y_2015_ProFSW, rep_len(y$viral_supp_y_1986_rest, 8))
y$viral_supp_y[3,] = y$viral_supp_y[2,]
y$infect_ART_y = 1-(y$viral_supp_y * y$ART_eff)
y$infect_ART_t = y$viral_supp_t
y$rho = c(y$ART_recruit_rate_FSW, y$ART_recruit_rate_rest,
y$ART_recruit_rate_rest,
y$ART_recruit_rate_rest,
y$ART_recruit_rate_rest * y$ART_init_ratio_MF,
y$ART_recruit_rate_rest * y$ART_init_ratio_MF, 0, 0, 0)
y$iota = c(y$ART_reinit_rate_FSW, y$ART_reinit_rate_rest,
y$ART_reinit_rate_rest,
y$ART_reinit_rate_rest,
y$ART_reinit_rate_rest * y$ART_init_ratio_MF,
y$ART_reinit_rate_rest * y$ART_init_ratio_MF, 0, 0, 0)
if(y$long_intervention == 0)
{
y$tau_intervention_y = matrix(c(rep(0, 9), y$intervention_testing_increase, c(rep(0, 8)), y$intervention_testing_increase, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T)
y$rho_intervention_y = matrix(c(rep(0, 9), y$intervention_ART_increase, c(rep(0, 8)), y$intervention_ART_increase, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T)
y$prep_intervention_y = matrix(c(rep(0, 9), y$prep_offering_rate , rep(0, 9-1), y$prep_offering_rate , rep(0, 9-1), rep(0, 9)), ncol = 9, byrow = T)
} else {
y$tau_intervention_y = matrix(c(rep(0, 9), y$intervention_testing_increase, c(rep(0, 8)), y$intervention_testing_increase, c(rep(0, 8)), y$intervention_testing_increase, c(rep(0, 8))), ncol = 9, nrow = 4, byrow = T)
y$rho_intervention_y = matrix(c(rep(0, 9), y$intervention_ART_increase, c(rep(0, 8)), y$intervention_ART_increase, c(rep(0, 8)), y$intervention_ART_increase, c(rep(0, 8))), ncol = 9, nrow = 4, byrow = T)
y$prep_intervention_y = matrix(c(rep(0, 9), y$prep_offering_rate , rep(0, 9-1), y$prep_offering_rate , rep(0, 9-1), y$prep_offering_rate , rep(0, 9-1)), ncol = 9, byrow = T)
}
y$kappaa = c(y$PrEP_loss_to_follow_up, rep_len(0,(8)))
# y$kappaa = c(y$prep_dropout, rep_len(0,(8)))
y$kappab = y$kappaa
y$kappac = y$kappaa
y$kappa1 = y$kappaa
y$above_500_by_group = c(y$FSW_eligible, y$GP_eligible, y$GP_eligible,
y$GP_eligible, y$GP_eligible, y$GP_eligible,
y$GP_eligible, y$GP_eligible, y$GP_eligible)
# y$W1 = y$W1
# y$W2 = y$W2
# y$W3 = y$W3
# y$W4 = y$W4
y$pfFSW_y = matrix(c(rep(0, 9), y$prev_non_ben_fsw_1993, c(rep(0, 8)), y$prev_non_ben_fsw_2015, c(rep(0, 8)), y$prev_non_ben_fsw_2015, c(rep(0, 8))), ncol = 9, nrow = 4, byrow = T)
return(y)
}
# N_Ncat7 = function(y) return S0_init[1] + I01_init[1] + I02_init[1] + I03_init[1] + I04_init[1] + I05_init[1]
#
# the parameters below will be sampled from an LHS and will replace their respective defaults
# unless I put something in the args of the function, eg sample = mu
#' @export
#' @useDynLib cotonou
lhs_parameters <- function(n, sample = NULL, Ncat = 9, Nage = 1, ..., set_pars = list(...), forced_pars = list(...), set_null= list(...), ranges = NULL, par_seq, condom_seq, groups_seq, years_seq) {
set_pars <- modifyList(set_pars, forced_pars)
# parameters that I have defined in ranges will be removed from set_pars and fixed pars (fixed pars below)
if(length(which(names(set_pars) %in% rownames(ranges))) > 0)
set_pars <- set_pars[-which(names(set_pars) %in% rownames(ranges))]
#fixed pars list i think for the fix parameters function
fixed_pars = list(
# test_rate_prep = c(4, 0, 0, 0, 0, 0, 0, 0, 0),
fraction_sexually_active_15_F = 0,
fraction_sexually_active_15_M = 0,
initial_Ntot = 286114,
frac_women_ProFSW = 0.0024,
frac_women_LowFSW = 3,
frac_women_exFSW = 0.0024,
frac_men_client = 0.2,
frac_women_virgin = 0.1,
frac_men_virgin = 0.1,
fraction_F = 0.516,
rate_move_in = matrix(0, nrow = Ncat, ncol = Ncat),
rate_move_out = rep_len(0, Ncat),
rate_move_out_PrEP = rep_len(0, Ncat),
epsilon_y = 0,
rate_enter_sexual_pop_F = 0.4,
rate_enter_sexual_pop_M = 0.4,
fraction_FSW_foreign = 0,
movement = 1,
dur_below_200 = rep_len(0.3,Ncat),
fc_y_comm_1985 = matrix(0, Ncat, Ncat),
fc_y_comm_1993 = matrix(0, Ncat, Ncat),
fc_y_comm_1995 = matrix(0, Ncat, Ncat),
fc_y_comm_1998 = matrix(0, Ncat, Ncat),
fc_y_comm_2002 = matrix(0, Ncat, Ncat),
fc_y_comm_2005 = matrix(0, Ncat, Ncat),
fc_y_comm_2008 = matrix(0, Ncat, Ncat),
fc_y_comm_2012 = matrix(0, Ncat, Ncat),
fc_y_comm_2015 = matrix(0, Ncat, Ncat),
fc_y_comm_2016 = matrix(0, Ncat, Ncat),
fc_y_noncomm_1985 = matrix(0, Ncat, Ncat),
fc_y_noncomm_1993 = matrix(0, Ncat, Ncat),
fc_y_noncomm_1998 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2002 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2008 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2011 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2015 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2016 = matrix(0, Ncat, Ncat),
n_y_comm_1985 = matrix(1.02, Ncat, Ncat),
n_y_comm_2002 = matrix(1.02, Ncat, Ncat),
n_y_comm_2015 = matrix(1.02, Ncat, Ncat),
n_y_comm_2016 = matrix(1.02, Ncat, Ncat),
n_y_noncomm_1985 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_1998 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_2011 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_2002 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_2015 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_2016 = matrix(1.03, Ncat, Ncat),
c_comm_1985 = rep_len(2, Ncat),
c_comm_1993 = rep_len(2, Ncat),
c_comm_1995 = rep_len(2, Ncat),
c_comm_1998 = rep_len(2, Ncat),
c_comm_2002 = rep_len(2, Ncat),
c_comm_2005 = rep_len(2, Ncat),
c_comm_2008 = rep_len(2, Ncat),
c_comm_2012 = rep_len(2, Ncat),
c_comm_2015 = rep_len(2, Ncat),
c_comm_2016 = rep_len(2, Ncat),
c_noncomm_1985 = rep_len(1, Ncat),
c_noncomm_1993 = rep_len(1, Ncat),
c_noncomm_1995 = rep_len(1, Ncat),
c_noncomm_1998 = rep_len(1, Ncat),
c_noncomm_2002 = rep_len(1, Ncat),
c_noncomm_2005 = rep_len(1, Ncat),
c_noncomm_2008 = rep_len(1, Ncat),
c_noncomm_2012 = rep_len(1, Ncat),
c_noncomm_2015 = rep_len(1, Ncat),
c_noncomm_2016 = rep_len(1, Ncat),
betaMtoF_noncomm = 0.00193,
betaFtoM_noncomm = 0.00867,
betaMtoF_comm = 0.00193,
betaFtoM_comm = 0.00867,
betaMtoF_baseline = 0.00081,
muF = 0.02597403,
muM = 0.02739726,
RR_beta_FtM = 1,
RR_beta_circum = 0.44,
prev_HSV2_FSW = 1,
prev_HSV2_Client = 1,
prev_HSV2_GPF = 1,
prev_HSV2_GPM = 1,
RR_beta_HSV2_comm_t = 1,
RR_beta_HSV2_noncomm_t = 1,
RR_beta_HSV2_comm_a = 2,
RR_beta_HSV2_noncomm_a = 2,
who_believe_comm = 0,
init_clientN_from_PCR = 0,
beta_above_1 = 0,
ignore_ranges_fc_c = 0,
dropout_rate_not_FSW = 0.025,
dropout_rate_FSW = 0.025,
delete = 0,
nu = 0.02222222,
PrEPOnOff = 0,
fPa = 0.5,
fPc = 0.4,
iota = c(0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0, 0, 0),
viral_supp_t = c(1986, 2015, 2016),
viral_supp_y = matrix(0, nrow = 3, ncol = 9),
viral_supp_y_1986_rest = 0.6,
viral_supp_y_2015_ProFSW = 0.7,
ART_eff = 0.98,
rho = c(0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0, 0, 0),
ART_recruit_rate = 0.2,
ART_reinit_rate = 0.1,
RR_test_CD4200 = 5.4,
tau_intervention_t = c(1986, 2015, 2016, 2017),
rho_intervention_t = c(1986, 2015, 2016, 2017),
tau_intervention_y = matrix(c(rep(0, 9), 1.5, c(rep(0, 8)), 1.5, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T),
rho_intervention_y = matrix(c(rep(0, 9), 6, c(rep(0, 8)), 6, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T),
intervention_testing_increase = 1.4,
intervention_ART_increase = 6.5,
psia = rep_len(0.1,9),
psib = rep_len(0.1,9),
prep_offering_rate = 1,
kappaa = c(0.2, rep_len(0,(9-1))),
kappab = c(0.2, rep_len(0,(9-1))),
kappac = c(0.2, rep_len(0,(9-1))),
kappa1 = c(0.2, rep_len(0,(9-1))),
prep_dropout = 2,
eP1a = 0.9,
dur_FSW = 30,
cost_ART_initiation_study_FSW = 1,
cost_1_year_ART_study_FSW = 1,
cost_1_year_ART_rest = 1,
W0 = 1,
W1 = 1,
W2 = 1,
W3 = 1,
cost_Initiation_of_ART_study_FSW = 1,
cost_Initiation_of_ART_government_FSW = 1,
cost_1_year_of_ART_study_FSW = 1,
cost_1_year_of_ART_government_FSW = 1,
cost_Initiation_ART_rest_of_population = 1,
cost_1_year_of_ART_rest_of_population = 1,
cost_FSW_initiation_ART_Patient_costs = 1,
cost_FSW_1_year_ART_Patient_costs = 1,
cost_Initiation_of_PrEP_study = 1,
cost_1_year_PrEP_perfect_adherence_study = 1,
cost_1_year_PrEP_intermediate_adherence_study = 1,
cost_1_year_PrEP_non_adherence_study = 1,
cost_Initiation_of_PrEP_government = 1,
cost_1_year_PrEP_perfect_adherence_government =1,
cost_1_year_PrEP_intermediate_adherence_government = 1,
cost_1_year_PrEP_non_adherence_government = 1,
cost_PREP_initiation_Patient_costs = 1,
cost_PREP_1_year_ART_Patient_costs = 1,
TasP_testing = 0,
long_intervention = 0,
above_500_by_group = c(1, 1, 1, 1, 1, 1, 1, 1, 1),
FSW_eligible = 0,
GP_eligible = 1,
prev_non_ben_fsw_1993 = 0.55,
prev_non_ben_fsw_2015 = 0.2,
infected_FSW_incoming = 1,
ART_recruit_rate_rest = 0.1,
ART_reinit_rate_rest = 0.1,
ART_init_ratio_MF = 2,
RR_beta_HSV2_a_FSW = 1,
RR_beta_HSV2_a_client = 1,
RR_beta_HSV2_a_GPF = 1,
RR_beta_HSV2_a_GPM = 1,
testing_prob_men_2006 = 0.111,
testing_prob_men_2008 = 0.112,
testing_prob_men_2012 = 0.113,
testing_prob_women_2006 = 0.114,
testing_prob_women_2008 = 0.115,
testing_prob_women_2012 = 0.116,
ART_recruit_rate_FSW = 1,
ART_reinit_rate_FSW = 1,
PrEP_reinit_OnOff_t = c(1985, 2017, 2017.01, 2060),
PrEP_reinit_OnOff_y = c(0, 0, 0, 0),
PrEP_reinits_on = 1,
old_VS_assumption = 1
)
if(length(which(names(fixed_pars) %in% rownames(ranges))) > 0)
fixed_pars <- fixed_pars[-which(names(fixed_pars) %in% rownames(ranges))]
mu <- matrix(rep(c(1/50, 1/42), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("mu", Ncat), NULL))
omega <- if(Ncat == 9) matrix(c(0.0017, 0.0067, 0, 0, 0, 0, 0, 0, 0.1, 0.2, 0, 0, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("omega", Ncat), NULL)) else
matrix(rep(c(0.4, 0.6), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("omega", Ncat), NULL))
# c_y_comm <- if(Ncat == 9) matrix(c(300, 1400, 40, 64, 0, 0, 0, 0, 18.67, 37.5, 0, 0, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_y_comm", Ncat), NULL)) else c(1, 3)
#these parameters need to be here so fix_parameters works? below are fixed...
S0_init = matrix(rep(c(4000, 4000), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("S0_init", Ncat), NULL))
I01_init = matrix(rep(c(1000, 1000), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("I01_init", Ncat), NULL))
N_init = if(Ncat == 9) matrix(c(672, 672, 757, 757, 130895, 130895, 672, 672, 27091, 27091, 100335, 100335, 14544, 14544, 11148, 11148, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("N_init", Ncat), NULL)) else c(300000, 300000)
# c_comm = if(Ncat == 9) matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL)) else
# matrix(rep(c(1,3), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL))
# c_comm = if(Ncat == 9) matrix(c(272, 1439, 40, 64, 0, 0, 0, 0, 18.67, 37.5, 0, 0, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL)) else
# matrix(rep(c(1,3), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL))
#
# c_noncomm = if(Ncat == 9) matrix(c(0.2729358, 0.4682779, 0.2729358, 0.4682779, 0.90, 1.02, 0.90, 1.02, 1.21, 2.5, 1.28, 1.40, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_noncomm", Ncat), NULL)) else
# matrix(rep(c(1,3), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_noncomm", Ncat), NULL))
#
old_ranges <- rbind(
# c_y_comm,
epsilon_1985 = c(0.059346131*1.5, 0.059346131*1.5),
epsilon_1992 = c(0.053594832*1.5, 0.053594832*1.5),
epsilon_2002 = c(0.026936907*1.5, 0.026936907*1.5),
epsilon_2013 = c(0.026936907*1.5, 0.026936907*1.5),
epsilon_2016 = c(0.026936907*1.5, 0.026936907*1.5),
rate_leave_pro_FSW = c(0.2, 0.2),
rate_leave_low_FSW = c(0.1, 0.1),
rate_leave_client = c(0.05, 0.05),
betaMtoF = c(0.00086, 0.00433),
betaFtoM = c(0.00279, 0.02701),
prev_init_FSW = c(0.01318836, 0.06592892),
prev_init_rest = c(0.0003134459, 0.0029420363),
# c_comm,
# c_noncomm,
mu,
eff_ART = c(0.96, 1), # infectiousness RR when on ART
infect_acute = c(4, 18), # RR for acute phase
# dur_primary_phase = c(1/0.5, 1/0.16), # rate
dur_primary_phase = c(0.16, 0.5), # from Mathieu's parameters IN YEARS
# dur_primary_phase = c(2, 6.25), # from Mathieu's parameters
SC_to_death = c(8.7, 12.3),
dur_200_349 = c(3.9, 5), # from Mathieu's parameters IN YEARS
# dur_200_349 = c(3.9, 5), # from Mathieu's parameters
ART_RR_prog = c(8.8, 12.1),
# ART_RR_mort = c(2, 3),
# omega,
#below are fixed...
S0_init,
I01_init,
N_init
)
if (is.null(ranges)) {
ranges <- old_ranges
}
if (!is.null(sample)) {
ranges <- ranges[rownames(ranges) %in% sample, drop=FALSE]
}
samples <- tgp::lhs(n, ranges)
nms <- rownames(ranges)
i <- split(seq_along(nms), nms)
f <- function(x) {
lapply(i, function(j) x[j])
}
samples_list <- apply(samples, 1, f)
# print("ehre:")
# print(samples_list)
samples_list <- lapply(samples_list, function(x) modifyList(x, fixed_pars)) # btw earlier, fixed pars is replaced by ranges where they overlap
# HERE!
# samples_list <- lapply(samples_list, after_LHS, set_pars = set_pars)
samples_list <- lapply(samples_list, function(x) modifyList(x, set_pars)) # set pars before fixed pars in order to get right fixed
# print("ehre2:")
# print(samples_list)
# samples_list_test <<- samples_list
samples_list <- lapply(samples_list, fix_parameters, Ncat = Ncat, par_seq = par_seq, condom_seq = condom_seq, groups_seq = groups_seq, years_seq = years_seq)
# samples_list_test_fixed <<- samples_list
# samples_list <- lapply(samples_list, function(x) modifyList(x, set_pars)) # set pars after fixed pars in order to get right set pars
samples_list <- lapply(samples_list, function(x) modifyList(x, forced_pars)) # set pars after fixed pars in order to get right set pars
samples_list <- lapply(samples_list, delete_parameter_sets)
samples_list_test <<- samples_list
if(n > 1)
samples_list <- samples_list[!unlist(lapply(samples_list, function(x) x$delete))] # removing those runs which do not pass initial criteria
lapply(samples_list, function(x) generate_parameters(parameters = x, Ncat = Ncat, set_null = set_null))
}
#' @export
#' @useDynLib cotonou
delete_parameter_sets <- function(x) {
x$delete <- ifelse(x$beta_above_1 == 1, 1, # if beta is above 1 due to the various RRs
# if believing FSW partner change rates, then remove runs where clients pop size too big or too small ONLY IF client size is fixed at start from PCR balancing
ifelse(((x$c_comm_1985[1] * x$N_init[1] + x$c_comm_1985[2] * x$N_init[2])/ (x$c_comm_1985[5])) > (0.4 * sum(c(x$N_init[5], x$N_init[6], x$N_init[8]))) && (x$init_clientN_from_PCR == 1), 1,
# if believing client partner change rates, we are now also believing pro FSW and altering low FSW. killing runs where frac of low fsw is too low
ifelse(((x$c_comm_1985[5] * x$N_init[5] - x$c_comm_1985[2] * x$N_init[2])/ (x$c_comm_1985[1])) < (0* sum(x$N_init)) && (x$init_clientN_from_PCR == 0), 1, 0)))
return(x)
}
#' @export
#' @useDynLib cotonou
lhs_parameters_parallel <- function(n, sample = NULL, Ncat = 9, Nage = 1, ..., set_pars = list(...), forced_pars = list(...), set_null= list(...), ranges = NULL, par_seq, condom_seq, groups_seq, years_seq) {
set_pars <- modifyList(set_pars, forced_pars)
# parameters that I have defined in ranges will be removed from set_pars and fixed pars (fixed pars below)
if(length(which(names(set_pars) %in% rownames(ranges))) > 0)
set_pars <- set_pars[-which(names(set_pars) %in% rownames(ranges))]
#fixed pars list i think for the fix parameters function
fixed_pars = list(
# test_rate_prep = c(4, 0, 0, 0, 0, 0, 0, 0, 0),
fraction_sexually_active_15_F = 0,
fraction_sexually_active_15_M = 0,
initial_Ntot = 286114,
frac_women_ProFSW = 0.0024,
frac_women_LowFSW = 3,
frac_women_exFSW = 0.0024,
frac_men_client = 0.2,
frac_women_virgin = 0.1,
frac_men_virgin = 0.1,
fraction_F = 0.516,
rate_move_in = matrix(0, nrow = Ncat, ncol = Ncat),
rate_move_out = rep_len(0, Ncat),
rate_move_out_PrEP = rep_len(0, Ncat),
epsilon_y = 0,
rate_enter_sexual_pop_F = 0.4,
rate_enter_sexual_pop_M = 0.4,
fraction_FSW_foreign = 0,
movement = 1,
dur_below_200 = rep_len(0.3,Ncat),
fc_y_comm_1985 = matrix(0, Ncat, Ncat),
fc_y_comm_1993 = matrix(0, Ncat, Ncat),
fc_y_comm_1995 = matrix(0, Ncat, Ncat),
fc_y_comm_1998 = matrix(0, Ncat, Ncat),
fc_y_comm_2002 = matrix(0, Ncat, Ncat),
fc_y_comm_2005 = matrix(0, Ncat, Ncat),
fc_y_comm_2008 = matrix(0, Ncat, Ncat),
fc_y_comm_2012 = matrix(0, Ncat, Ncat),
fc_y_comm_2015 = matrix(0, Ncat, Ncat),
fc_y_comm_2016 = matrix(0, Ncat, Ncat),
fc_y_noncomm_1985 = matrix(0, Ncat, Ncat),
fc_y_noncomm_1993 = matrix(0, Ncat, Ncat),
fc_y_noncomm_1998 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2002 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2008 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2011 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2015 = matrix(0, Ncat, Ncat),
fc_y_noncomm_2016 = matrix(0, Ncat, Ncat),
n_y_comm_1985 = matrix(1.02, Ncat, Ncat),
n_y_comm_2002 = matrix(1.02, Ncat, Ncat),
n_y_comm_2015 = matrix(1.02, Ncat, Ncat),
n_y_comm_2016 = matrix(1.02, Ncat, Ncat),
n_y_noncomm_1985 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_2002 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_2015 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_1998 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_2011 = matrix(1.03, Ncat, Ncat),
n_y_noncomm_2016 = matrix(1.03, Ncat, Ncat),
c_comm_1985 = rep_len(2, Ncat),
c_comm_1993 = rep_len(2, Ncat),
c_comm_1995 = rep_len(2, Ncat),
c_comm_1998 = rep_len(2, Ncat),
c_comm_2002 = rep_len(2, Ncat),
c_comm_2005 = rep_len(2, Ncat),
c_comm_2008 = rep_len(2, Ncat),
c_comm_2012 = rep_len(2, Ncat),
c_comm_2015 = rep_len(2, Ncat),
c_comm_2016 = rep_len(2, Ncat),
c_noncomm_1985 = rep_len(1, Ncat),
c_noncomm_1993 = rep_len(1, Ncat),
c_noncomm_1995 = rep_len(1, Ncat),
c_noncomm_1998 = rep_len(1, Ncat),
c_noncomm_2002 = rep_len(1, Ncat),
c_noncomm_2005 = rep_len(1, Ncat),
c_noncomm_2008 = rep_len(1, Ncat),
c_noncomm_2012 = rep_len(1, Ncat),
c_noncomm_2015 = rep_len(1, Ncat),
c_noncomm_2016 = rep_len(1, Ncat),
betaMtoF_noncomm = 0.00193,
betaFtoM_noncomm = 0.00867,
betaMtoF_comm = 0.00193,
betaFtoM_comm = 0.00867,
betaMtoF_baseline = 0.00081,
muF = 0.02597403,
muM = 0.02739726,
RR_beta_FtM = 1,
RR_beta_circum = 0.44,
prev_HSV2_FSW = 1,
prev_HSV2_Client = 1,
prev_HSV2_GPF = 1,
prev_HSV2_GPM = 1,
RR_beta_HSV2_comm_t = 1,
RR_beta_HSV2_noncomm_t = 1,
RR_beta_HSV2_comm_a = 2,
RR_beta_HSV2_noncomm_a = 2,
who_believe_comm = 0,
init_clientN_from_PCR = 0,
beta_above_1 = 0,
ignore_ranges_fc_c = 0,
dropout_rate_not_FSW = 0.025,
dropout_rate_FSW = 0.025,
delete = 0,
nu = 0.02222222,
PrEPOnOff = 0,
fPa = 0.5,
fPc = 0.4,
iota = c(0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0, 0, 0),
viral_supp_t = c(1986, 2015, 2016),
viral_supp_y = matrix(0, nrow = 3, ncol = 9),
viral_supp_y_1986_rest = 0.6,
viral_supp_y_2015_ProFSW = 0.7,
ART_eff = 0.98,
rho = c(0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0, 0, 0),
ART_recruit_rate = 0.2,
ART_reinit_rate = 0.1,
RR_test_CD4200 = 5.4,
tau_intervention_t = c(1986, 2015, 2016, 2017),
rho_intervention_t = c(1986, 2015, 2016, 2017),
tau_intervention_y = matrix(c(rep(0, 9), 0, c(rep(0, 8)), 0, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T),
rho_intervention_y = matrix(c(rep(0, 9), 0, c(rep(0, 8)), 0, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T),
intervention_testing_increase = 1.4,
intervention_ART_increase = 6.5,
psia = rep_len(0.1,9),
psib = rep_len(0.1,9),
prep_offering_rate = 1,
kappaa = c(0.2, rep_len(0,(9-1))),
kappab = c(0.2, rep_len(0,(9-1))),
kappac = c(0.2, rep_len(0,(9-1))),
kappa1 = c(0.2, rep_len(0,(9-1))),
prep_dropout = 2,
eP1a = 0.9,
dur_FSW = 30,
cost_ART_initiation_study_FSW = 1,
cost_1_year_ART_study_FSW = 1,
cost_1_year_ART_rest = 1,
W0 = 1,
W1 = 1,
W2 = 1,
W3 = 1,
cost_Initiation_of_ART_study_FSW = 1,
cost_Initiation_of_ART_government_FSW = 1,
cost_1_year_of_ART_study_FSW = 1,
cost_1_year_of_ART_government_FSW = 1,
cost_Initiation_ART_rest_of_population = 1,
cost_1_year_of_ART_rest_of_population = 1,
cost_FSW_initiation_ART_Patient_costs = 1,
cost_FSW_1_year_ART_Patient_costs = 1,
cost_Initiation_of_PrEP_study = 1,
cost_1_year_PrEP_perfect_adherence_study = 1,
cost_1_year_PrEP_intermediate_adherence_study = 1,
cost_1_year_PrEP_non_adherence_study = 1,
cost_Initiation_of_PrEP_government = 1,
cost_1_year_PrEP_perfect_adherence_government =1,
cost_1_year_PrEP_intermediate_adherence_government = 1,
cost_1_year_PrEP_non_adherence_government = 1,
cost_PREP_initiation_Patient_costs = 1,
cost_PREP_1_year_ART_Patient_costs = 1,
TasP_testing = 0,
long_intervention = 0,
above_500_by_group = c(1, 1, 1, 1, 1, 1, 1, 1, 1),
FSW_eligible = 0,
GP_eligible = 1,
prev_non_ben_fsw_1993 = 0.55,
prev_non_ben_fsw_2015 = 0.2,
infected_FSW_incoming = 1,
ART_recruit_rate_rest = 0.1,
ART_reinit_rate_rest = 0.1,
ART_init_ratio_MF = 2,
RR_beta_HSV2_a_FSW = 1,
RR_beta_HSV2_a_client = 1,
RR_beta_HSV2_a_GPF = 1,
RR_beta_HSV2_a_GPM = 1,
testing_prob_men_2006 = 0.111,
testing_prob_men_2008 = 0.112,
testing_prob_men_2012 = 0.113,
testing_prob_women_2006 = 0.114,
testing_prob_women_2008 = 0.115,
testing_prob_women_2012 = 0.116,
ART_recruit_rate_FSW = 1,
ART_reinit_rate_FSW = 1,
PrEP_reinit_OnOff_t = c(1985, 2017, 2017.01, 2060),
PrEP_reinit_OnOff_y = c(0, 0, 0, 0),
PrEP_reinits_on = 1,
old_VS_assumption = 1
)
if(length(which(names(fixed_pars) %in% rownames(ranges))) > 0)
fixed_pars <- fixed_pars[-which(names(fixed_pars) %in% rownames(ranges))]
mu <- matrix(rep(c(1/50, 1/42), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("mu", Ncat), NULL))
omega <- if(Ncat == 9) matrix(c(0.0017, 0.0067, 0, 0, 0, 0, 0, 0, 0.1, 0.2, 0, 0, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("omega", Ncat), NULL)) else
matrix(rep(c(0.4, 0.6), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("omega", Ncat), NULL))
# c_y_comm <- if(Ncat == 9) matrix(c(300, 1400, 40, 64, 0, 0, 0, 0, 18.67, 37.5, 0, 0, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_y_comm", Ncat), NULL)) else c(1, 3)
#these parameters need to be here so fix_parameters works? below are fixed...
S0_init = matrix(rep(c(4000, 4000), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("S0_init", Ncat), NULL))
I01_init = matrix(rep(c(1000, 1000), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("I01_init", Ncat), NULL))
N_init = if(Ncat == 9) matrix(c(672, 672, 757, 757, 130895, 130895, 672, 672, 27091, 27091, 100335, 100335, 14544, 14544, 11148, 11148, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("N_init", Ncat), NULL)) else c(300000, 300000)
# c_comm = if(Ncat == 9) matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL)) else
# matrix(rep(c(1,3), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL))
# c_comm = if(Ncat == 9) matrix(c(272, 1439, 40, 64, 0, 0, 0, 0, 18.67, 37.5, 0, 0, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL)) else
# matrix(rep(c(1,3), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL))
#
# c_noncomm = if(Ncat == 9) matrix(c(0.2729358, 0.4682779, 0.2729358, 0.4682779, 0.90, 1.02, 0.90, 1.02, 1.21, 2.5, 1.28, 1.40, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_noncomm", Ncat), NULL)) else
# matrix(rep(c(1,3), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_noncomm", Ncat), NULL))
#
old_ranges <- rbind(
# c_y_comm,
epsilon_1985 = c(0.059346131*1.5, 0.059346131*1.5),
epsilon_1992 = c(0.053594832*1.5, 0.053594832*1.5),
epsilon_2002 = c(0.026936907*1.5, 0.026936907*1.5),
epsilon_2013 = c(0.026936907*1.5, 0.026936907*1.5),
epsilon_2016 = c(0.026936907*1.5, 0.026936907*1.5),
rate_leave_pro_FSW = c(0.2, 0.2),
rate_leave_low_FSW = c(0.1, 0.1),
rate_leave_client = c(0.05, 0.05),
betaMtoF = c(0.00086, 0.00433),
betaFtoM = c(0.00279, 0.02701),
prev_init_FSW = c(0.01318836, 0.06592892),
prev_init_rest = c(0.0003134459, 0.0029420363),
# c_comm,
# c_noncomm,
mu,
eff_ART = c(0.96, 1), # infectiousness RR when on ART
infect_acute = c(4, 18), # RR for acute phase
# dur_primary_phase = c(1/0.5, 1/0.16), # rate
dur_primary_phase = c(0.16, 0.5), # from Mathieu's parameters IN YEARS
# dur_primary_phase = c(2, 6.25), # from Mathieu's parameters
SC_to_death = c(8.7, 12.3),
dur_200_349 = c(3.9, 5), # from Mathieu's parameters IN YEARS
# dur_200_349 = c(3.9, 5), # from Mathieu's parameters
ART_RR_prog = c(8.8, 12.1),
# ART_RR_mort = c(2, 3),
# omega,
#below are fixed...
S0_init,
I01_init,
N_init
)
if (is.null(ranges)) {
ranges <- old_ranges
}
if (!is.null(sample)) {
ranges <- ranges[rownames(ranges) %in% sample, drop=FALSE]
}
samples <- tgp::lhs(n, ranges)
nms <- rownames(ranges)
i <- split(seq_along(nms), nms)
f <- function(x) {
lapply(i, function(j) x[j])
}
samples_list <- apply(samples, 1, f)
# print("ehre:")
# print(samples_list)
samples_list <- parallel::parLapply(NULL, samples_list, function(x) modifyList(x, fixed_pars)) # btw earlier, fixed pars is replaced by ranges where they overlap
# HERE!
# samples_list <- lapply(samples_list, after_LHS, set_pars = set_pars)
samples_list <- parallel::parLapply(NULL, samples_list, function(x) modifyList(x, set_pars)) # set pars before fixed pars in order to get right fixed
# print("ehre2:")
# print(samples_list)
# samples_list_test <<- samples_list
samples_list <- parallel::parLapply(NULL, samples_list, fix_parameters, Ncat = Ncat, par_seq = par_seq, condom_seq = condom_seq, groups_seq = groups_seq, years_seq = years_seq)
# samples_list_test_fixed <<- samples_list
# samples_list <- lapply(samples_list, function(x) modifyList(x, set_pars)) # set pars after fixed pars in order to get right set pars
samples_list <- parallel::parLapply(NULL, samples_list, function(x) modifyList(x, forced_pars)) # set pars after fixed pars in order to get right set pars
samples_list <- parallel::parLapply(NULL, samples_list, delete_parameter_sets)
if(n > 1)
samples_list <- samples_list[!unlist(parallel::parLapply(NULL, samples_list, function(x) x$delete))] # removing those runs which do not pass initial criteria
parallel::parLapply(NULL, samples_list, function(x) generate_parameters(parameters = x, Ncat = Ncat, set_null = set_null))
}
#' @export
#' @useDynLib cotonou
generate_parameters <- function(..., parameters = list(...), set_null = list(...), Ncat = 9, Nage = 1) {
defaults <- list(Ncat = Ncat,
Nage = Nage,
N_init = 300000, #ignore that this is weird - I replace with a vector later
initial_Ntot = 286114,
frac_women_ProFSW = 0.0024,
frac_women_LowFSW = 3,
frac_women_exFSW = 0.0024,
frac_men_client = 0.2,
frac_women_virgin = 0.1,
frac_men_virgin = 0.1,
above_500_by_group = c(1, 1, 1, 1, 1, 1, 1, 1, 1),
ART_eligible_CD4_above_500_t = c(1985, 2002, 2012, 2015, 2016),
ART_eligible_CD4_350_500_t = c(1985, 2002, 2012, 2015, 2016),
ART_eligible_CD4_200_349_t = c(1985, 2002, 2012, 2015, 2016),
ART_eligible_CD4_below_200_t = c(1985, 2002, 2012, 2015, 2016),
ART_eligible_CD4_above_500_y = c(0, 0, 0, 0, 1),
ART_eligible_CD4_350_500_y = c(0, 0, 0, 1, 1),
ART_eligible_CD4_200_349_y = c(0, 0.1, 1, 1, 1),
ART_eligible_CD4_below_200_y = c(0, 1, 1, 1, 1),
# epsilon_t_comm = c(1985, 1991.99, 1992, 2001.99, 2002, 2012.99, 2013, 2016),
# epsilon_y_comm = c(0.059346131, 0.059346131, 0.053594832, 0.053594832, 0.026936907, 0.026936907, 0.026936907, 0.026936907),
epsilon_t = c(1985, 1992, 2002, 2013, 2016),
epsilon_y = c(0.059346131, 0.053594832, 0.026936907, 0.026936907, 0.026936907),
epsilon_1985 = 0.059346131,
epsilon_1992 = 0.053594832,
epsilon_2002 = 0.026936907,
epsilon_2013 = 0.026936907,
epsilon_2016 = 0.026936907,
S0_init = rep_len(2000, Ncat),
S1a_init = rep_len(0, Ncat),
S1b_init = rep_len(0, Ncat),
S1c_init = rep_len(0, Ncat),
S1d_init = rep_len(0, Ncat),
I01_init = rep_len(1000, Ncat),
I11_init = rep_len(0, Ncat),
I02_init = rep_len(0, Ncat),
I03_init = rep_len(0, Ncat),
I04_init = rep_len(0, Ncat),
I05_init = rep_len(0, Ncat),
I22_init = rep_len(0, Ncat),
I23_init = rep_len(0, Ncat),
I24_init = rep_len(0, Ncat),
I25_init = rep_len(0, Ncat),
I32_init = rep_len(0, Ncat),
I33_init = rep_len(0, Ncat),
I34_init = rep_len(0, Ncat),
I35_init = rep_len(0, Ncat),
I42_init = rep_len(0, Ncat),
I43_init = rep_len(0, Ncat),
I44_init = rep_len(0, Ncat),
I45_init = rep_len(0, Ncat),
cumuInf_init = rep_len(0, Ncat),
mu = rep_len(0.02, Ncat),
dur_primary_phase = 0.3,
gamma02 = rep_len(0.2, Ncat),
gamma03 = rep_len(0.2, Ncat),
dur_200_349 = 3,
gamma11 = rep_len(0.2, Ncat),
gamma22 = rep_len(0.2, Ncat),
gamma23 = rep_len(0.2, Ncat),
gamma24 = rep_len(0.2, Ncat),
gamma32 = rep_len(0.2, Ncat),
gamma33 = rep_len(0.2, Ncat),
gamma34 = rep_len(0.2, Ncat),
gamma42 = rep_len(0.2, Ncat),
gamma43 = rep_len(0.2, Ncat),
gamma44 = rep_len(0.2, Ncat),
ART_prob_t = c(1985, 2005, 2012, 2015),
ART_prob_y = matrix(
rep(c(0, 0, 0.2, 0.4), Ncat), ncol = Ncat),
rho = c(0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0, 0, 0),
phi2 = rep_len(0.004,Ncat), # sort out later
phi3 = rep_len(0.004,Ncat),
phi4 = rep_len(0.004,Ncat),
phi5 = rep_len(0.004,Ncat),
psia = rep_len(0.1,Ncat),
psib = rep_len(0.1,Ncat),
testing_prob_t = c(1985, 2001, 2005, 2006, 2008, 2012, 2013, 2015, 2016),
testing_prob_y = matrix(
rep(c(0, 0.1, 0.2, 0.4, 0.5, 0.7, 0.5, 0.8, 0.7), Ncat), ncol = Ncat),
test_rate_prep = c(4, 0, 0, 0, 0, 0, 0, 0, 0),
sigma = c(0.86, 0, 0, 0, 0, 0, 0, 0, 0),
prep_intervention_t = c(1985, 2015, 2016, 2017),
prep_intervention_y = matrix(c(rep(0, Ncat), .1, rep(0, Ncat-1), .1, rep(0, Ncat-1), rep(0, Ncat)), ncol = Ncat, byrow = T), # offering rate
PrEPOnOff = 0,
RR_test_onPrEP = 2,
RR_test_CD4200 = 5.4,
RR_ART_CD4200 = 2,
tau01 = rep_len(1,Ncat),
tau11 = rep_len(1,Ncat),
tau2 = rep_len(1,Ncat),
tau3 = rep_len(1,Ncat),
tau4 = rep_len(1,Ncat),
tau5 = rep_len(1,Ncat),
# PREP
zetaa_t = c(1985, 2013, 2015, 2016),
zetaa_y = matrix(c(rep(0, Ncat), 0.0075, rep(0, Ncat-1), rep(0, Ncat), rep(0, Ncat)), ncol = Ncat, byrow = T),
zetab_t = c(1985, 2013, 2015, 2016),
zetab_y = matrix(c(rep(0, Ncat), 0.0075, rep(0, Ncat-1), rep(0, Ncat), rep(0, Ncat)), ncol = Ncat, byrow = T),
zetac_t = c(1985, 2013, 2015, 2016),
zetac_y = matrix(c(rep(0, Ncat), 0.0075, rep(0, Ncat-1), rep(0, Ncat), rep(0, Ncat)), ncol = Ncat, byrow = T),
eP = rep_len(0.6,Ncat),
eP0 = c(0, rep_len(0, (Ncat-1))),
eP1a = c(0.9, rep_len(0, (Ncat-1))),
eP1b = c(0.45, rep_len(0, (Ncat-1))),
eP1c = c(0, rep_len(0, (Ncat-1))),
eP1d = c(0, rep_len(0, (Ncat-1))),
# partner change rate
c_comm_1985 = rep_len(2, Ncat),
c_comm_1993 = rep_len(2, Ncat),
c_comm_1995 = rep_len(2, Ncat),
c_comm_1998 = rep_len(2, Ncat),
c_comm_2002 = rep_len(2, Ncat),
c_comm_2005 = rep_len(2, Ncat),
c_comm_2008 = rep_len(2, Ncat),
c_comm_2012 = rep_len(2, Ncat),
c_comm_2015 = rep_len(2, Ncat),
c_comm_2016 = rep_len(2, Ncat),
c_noncomm_1985 = rep_len(1, Ncat),
c_noncomm_1993 = rep_len(1, Ncat),
c_noncomm_1995 = rep_len(1, Ncat),
c_noncomm_1998 = rep_len(1, Ncat),
c_noncomm_2002 = rep_len(1, Ncat),
c_noncomm_2005 = rep_len(1, Ncat),
c_noncomm_2008 = rep_len(1, Ncat),
c_noncomm_2012 = rep_len(1, Ncat),
c_noncomm_2015 = rep_len(1, Ncat),
c_noncomm_2016 = rep_len(1, Ncat),
c_t_comm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016),
c_t_noncomm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016),
c_y_comm = matrix(2, nrow = 10, ncol = Ncat),
c_y_noncomm = matrix(1, nrow = 10, ncol = Ncat),
# condoms
# fc_t_noncomm = c(1985, 1990, 1998, 2016),
# fc_y_noncomm = matrix(
# rep(c(0, 0, 0.3, 0.5), Ncat), ncol = Ncat),
fc_y_comm_1985 = matrix(0.2, Ncat, Ncat),
fc_y_comm_1993 = matrix(0.5, Ncat, Ncat),
fc_y_comm_1995 = matrix(0.7, Ncat, Ncat),
fc_y_comm_1998 = matrix(0.3, Ncat, Ncat),
fc_y_comm_2002 = matrix(0.4, Ncat, Ncat),
fc_y_comm_2005 = matrix(0.1, Ncat, Ncat),
fc_y_comm_2008 = matrix(0.1, Ncat, Ncat),
fc_y_comm_2012 = matrix(0.6, Ncat, Ncat),
fc_y_comm_2015 = matrix(0, Ncat, Ncat),
fc_y_comm_2016 = matrix(0, Ncat, Ncat),
fc_y_noncomm_1985 = matrix(0.2, Ncat, Ncat),
fc_y_noncomm_1993 = matrix(0.2, Ncat, Ncat),
fc_y_noncomm_1998 = matrix(0.4, Ncat, Ncat),
fc_y_noncomm_2002 = matrix(0.3, Ncat, Ncat),
fc_y_noncomm_2008 = matrix(0.3, Ncat, Ncat),
# fc_y_noncomm_2011 = matrix(0.3, Ncat, Ncat),
# fc_y_noncomm_2015 = matrix(0.5, Ncat, Ncat),
# fc_y_noncomm_2016 = matrix(0.5, Ncat, Ncat),
fc_y_comm = array(0.5,dim=c(10,Ncat,Ncat)),
fc_y_noncomm = array(0.5,dim=c(7,Ncat,Ncat)),
fc_t_comm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016),
# fc_y_comm = matrix(
# rep(c(0.5, 0.5, 0.9, 0.99), Ncat), ncol = Ncat),
fc_t_noncomm = c(1985, 1993, 1998, 2002, 2008, 2011, 2015, 2016),
# fc_y_comm = matrix(
# rep(c(0.5, 0.5, 0.9, 0.99), Ncat), ncol = Ncat),
kappaa = c(0.2, rep_len(0,(Ncat-1))),
kappab = c(0.2, rep_len(0,(Ncat-1))),
kappac = c(0.2, rep_len(0,(Ncat-1))),
kappa1 = c(0.2, rep_len(0,(Ncat-1))),
alpha01 = rep_len(0,Ncat),
alpha02 = rep_len(0,Ncat),
alpha03 = 0.03,
alpha04 = 0.07,
# dur_below_200 = rep_len(0.3, Ncat),
alpha11 = rep_len(0.01,Ncat),
# alpha21 = rep_len(0.01,Ncat),
alpha22 = rep_len(0.01,Ncat),
alpha23 = rep_len(0.01,Ncat),
alpha24 = rep_len(0.01,Ncat),
alpha25 = rep_len(0.3,Ncat),
alpha32 = rep_len(0.01,Ncat),
# alpha33 = rep_len(0.01,Ncat),
# alpha34 = rep_len(0.01,Ncat),
# alpha35 = rep_len(0.3,Ncat),
alpha42 = rep_len(0.01,Ncat),
alpha43 = rep_len(0.01,Ncat),
alpha44 = rep_len(0.01,Ncat),
alpha45 = rep_len(0.3,Ncat),
# beta = rep_len(0.005,Ncat),
# betaMtoF = 0.00193,
# betaFtoM = 0.00867,
# betaMtoF_noncomm = 0.00193,
# betaFtoM_noncomm = 0.00867,
# betaMtoF_comm = 0.00193,
# betaFtoM_comm = 0.00867,
#beta = 0,
# p_comm = matrix(1, ncol = Ncat, nrow = Ncat),
# p_noncomm = matrix(1, ncol = Ncat, nrow = Ncat),
ec = rep_len(0.8,Ncat),
# ec = rep_len(1,1),
# epsilon = 0.001,
# c_comm = rep_len(2, Ncat),
# c_noncomm = rep_len(2, Ncat),
# c_t_comm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016),
# c_y_comm = matrix(rep(c(1, 2, 3, 1, 3, 4, 2, 1, 2, 4), Ncat), ncol = Ncat),
#
# c_t_noncomm = c(1985, 1990, 1998, 2016),
# c_y_noncomm = matrix(rep(c(0.5, 1, 0.7, 0.9), Ncat), ncol = Ncat),
fP_t_comm = c(1985, 2014, 2015, 2016, 2060),
fP_y_comm = matrix(
rep(c(1, 1, 1, 1, 1), Ncat), ncol = Ncat),
fP_t_noncomm = c(1985, 2014, 2015, 2016, 2060),
fP_y_noncomm = matrix(
rep(c(1, 1, 1, 1, 1), Ncat), ncol = Ncat),
# n_y_comm_1985 = matrix(1.02, Ncat, Ncat),
# n_y_comm_2002 = matrix(1.02, Ncat, Ncat),
# n_y_comm_2015 = matrix(1.02, Ncat, Ncat),
# n_y_comm_2016 = matrix(1.02, Ncat, Ncat),
# n_y_noncomm_1985 = matrix(1.03, Ncat, Ncat),
# n_y_noncomm_2002 = matrix(1.03, Ncat, Ncat),
# n_y_noncomm_2015 = matrix(1.03, Ncat, Ncat),
# n_y_noncomm_1998 = matrix(1.03, Ncat, Ncat),
# n_y_noncomm_2011 = matrix(1.03, Ncat, Ncat),
# n_y_noncomm_2016 = matrix(1.03, Ncat, Ncat),
#n = 0,
R = 1,
omega = rep_len(1/Ncat, Ncat),
theta = matrix(0.5, ncol = Ncat, nrow = Ncat),
M_comm = matrix(1/Ncat, Ncat, Ncat),
M_noncomm = matrix(1/Ncat, Ncat, Ncat),
# A_F = matrix(1/NAge, NAge, NAge),
# A_M = matrix(1/NAge, NAge, NAge),
# p = matrix(1, NAge, NAge),
# ART_RR_prog = 10, # survival extension cofactor
# ART_RR_mort = 2.5, # survival extension cofactor
# infect_ART = c(0, rep_len(0, 8)),
infect_acute = 9, # RR for acute phase
infect_AIDS = 7.27, # RR for AIDS phase
dur_FSW = 30,
# OnPrEP_init = rep_len(0, Ncat),
# have to include the follow in the function for it to work just using generate_parameters(), and not lhs_parameters()
# SC_to_death = 10,
# prev_init_FSW = 0.04,
# prev_init_rest = 0.0008,
rate_leave_pro_FSW = 0.2,
# rate_leave_client = 0.05,
# rate_leave_low_FSW = 0.1,
# prop_client_GPM = 0.2430057, # 27091/111483
# prop_pro_FSW_GPF = 0.004620494, # 672 / 145439
# prop_low_FSW_GPF = 0.005204931, # 757 / 145439
rate_move_in = matrix(0, ncol = Ncat, nrow = Ncat),
rate_move_out = rep_len(0, Ncat),
rate_move_out_PrEP = rep_len(0, Ncat),
# rate_enter_sexual_pop_F = 1,
# rate_enter_sexual_pop_M = 1,
# fraction_F = 0.51,
fraction_FSW_foreign = 0.5,
replaceDeaths = 0,
# movement = 1,
beta_comm = rep_len(0.002, Ncat),
beta_noncomm = rep_len(0.001, Ncat),
# muF = 0.02597403,
# muM = 0.02739726,
# RR_beta_FtM = 1,
# RR_beta_circum = 0.44,
# prev_HSV2_FSW = 1,
# prev_HSV2_Client = 1,
# prev_HSV2_GPF = 1,
# prev_HSV2_GPM = 1,
# RR_beta_HSV2_comm_t = 1,
# RR_beta_HSV2_noncomm_t = 1,
# RR_beta_HSV2_comm_a = 2,
# RR_beta_HSV2_noncomm_a = 2,
# betaMtoF_baseline = 0.00081,
who_believe_comm = 0,
# init_clientN_from_PCR = 0,
# fraction_sexually_active_15_F = 0,
# fraction_sexually_active_15_M = 0,
# beta_above_1 = 0,
# ignore_ranges_fc_c = 0,
# dropout_rate_not_FSW = 0.025,
# dropout_rate_FSW = 0.025,
# delete = 0,
nu = 0.02222222,
fPa = 0.5,
fPc = 0.4,
iota = c(0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0, 0, 0),
viral_supp_t = c(1985, 2015, 2016),
viral_supp_y = matrix(0, nrow = 3, ncol = 9),
# viral_supp_y_1986_rest = 0.6,
# viral_supp_y_2015_ProFSW = 0.7,
# ART_eff = 0.98,
# ART_recruit_rate = 0.2,
# ART_reinit_rate = 0.1,
tau_intervention_t = c(1986, 2015, 2016, 2017),
rho_intervention_t = c(1986, 2015, 2016, 2017),
tau_intervention_y = matrix(c(rep(0, 9), 1.5, c(rep(0, 8)), 1.5, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T),
rho_intervention_y = matrix(c(rep(0, 9), 6, c(rep(0, 8)), 6, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T),
# intervention_testing_increase = 1.4,
# intervention_ART_increase = 6.5,
# prep_offering_rate = 1,
FSW_ONLY = c(1, 1, rep(0, 7)),
# cost_ART_initiation_study_FSW = 1,
# cost_1_year_ART_study_FSW = 1,
# cost_1_year_ART_rest = 1,
W0 = 1,
W1 = 1,
W2 = 1,
W3 = 1,
cost_Initiation_of_ART_study_FSW = 1,
cost_Initiation_of_ART_government_FSW = 1,
cost_1_year_of_ART_study_FSW = 1,
cost_1_year_of_ART_government_FSW = 1,
cost_Initiation_ART_rest_of_population = 1,
cost_1_year_of_ART_rest_of_population = 1,
cost_FSW_initiation_ART_Patient_costs = 1,
cost_FSW_1_year_ART_Patient_costs = 1,
cost_Initiation_of_PrEP_study = 1,
cost_1_year_PrEP_perfect_adherence_study = 1,
cost_1_year_PrEP_intermediate_adherence_study = 1,
cost_1_year_PrEP_non_adherence_study = 1,
cost_Initiation_of_PrEP_government = 1,
cost_1_year_PrEP_perfect_adherence_government =1,
cost_1_year_PrEP_intermediate_adherence_government = 1,
cost_1_year_PrEP_non_adherence_government = 1,
cost_PREP_initiation_Patient_costs = 1,
cost_PREP_1_year_ART_Patient_costs = 1,
TasP_testing = 0,
# long_intervention = 0,
# FSW_eligible = 0,
# GP_eligible = 1,
pfFSW_t = c(1985, 1993, 2015, 2060),
pfFSW_y = matrix(c(rep(0, 9), 0.55, c(rep(0, 8)), 0.2, c(rep(0, 8)), 0.2, c(rep(0, 8))), ncol = 9, nrow = 4, byrow = T),
# prev_non_ben_fsw_1993 = 0.55,
# prev_non_ben_fsw_2015 = 0.2,
infected_FSW_incoming = 1,
prep_efficacious_y = c(0, 1, 1, 0, 0),
prep_efficacious_t = c(1985, 2015, 2017, 2017.01, 2060),
# PrEP_loss_to_follow_up = 0.1,
# ART_recruit_rate_rest = 0.1,
# ART_reinit_rate_rest = 0.1,
# ART_init_ratio_MF = 2,
# RR_beta_HSV2_a_FSW = 1,
# RR_beta_HSV2_a_client = 1,
# RR_beta_HSV2_a_GPF = 1,
# RR_beta_HSV2_a_GPM = 1,
# testing_prob_men_2006 = 0.111,
# testing_prob_men_2008 = 0.112,
# testing_prob_men_2012 = 0.113,
#
# testing_prob_women_2006 = 0.114,
# testing_prob_women_2008 = 0.115,
# testing_prob_women_2012 = 0.116,
# ART_recruit_rate_FSW = 1,
# ART_reinit_rate_FSW = 1,
PrEP_reinit_OnOff_t = c(1985, 2017, 2017.01, 2060),
PrEP_reinit_OnOff_y = c(0, 0, 0, 0),
# PrEP_reinits_on = 1,
count_PrEP_1a = 1,
count_PrEP_1b = 1,
count_PrEP_1c = 1,
re_init_interruption_parm_t = c(1985, 2017, 2017.01, 2060),
re_init_interruption_parm_y = matrix(1, nrow = 4, ncol = 9),
art_dropout_interruption_parm_t = c(1985, 2017, 2017.01, 2060),
art_dropout_interruption_parm_y = matrix(1, nrow = 4, ncol = 9),
art_initiation_interruption_parm_t = c(1985, 2017, 2017.01, 2060),
art_initiation_interruption_parm_y = matrix(1, nrow = 4, ncol = 9),
old_VS_assumption = 1
)
if (length(parameters) == 0L) {
return(defaults)
}
if (is.null(names(parameters)) || !all(nzchar(names(parameters)))) {
stop("All arguments must be named")
}
# extra <- setdiff(names(parameters), names(defaults))
# if (length(extra) > 0L) {
# stop("Unknown arguments: ", extra)
# }
# list of parameters that depend on others
ret <- modifyList(defaults, parameters)
if (length(set_null) > 0L) {
ret = modifyList(ret, lapply(ret[match(set_null, names(ret))], function(x) x*0))
}
# browser()
ret <- ret[which(!names(ret) %in% c("N_init", "initial_Ntot", "frac_women_ProFSW", "frac_women_LowFSW", "frac_women_exFSW", "frac_men_client", "frac_women_virgin", "frac_men_virgin", "epsilon_1985", "epsilon_1992", "epsilon_2002", "epsilon_2013", "epsilon_2016", "dur_primary_phase", "dur_200_349", "gamma32", "gamma33", "gamma34", "ART_prob_t", "ART_prob_y", "RR_test_onPrEP", "RR_ART_CD4200", "tau01", "tau11", "tau2", "tau3", "tau4", "tau5", "zetaa_t", "zetaa_y", "zetab_t", "zetab_y", "zetac_t", "zetac_y", "c_comm_1985", "c_comm_1993", "c_comm_1995", "c_comm_1998", "c_comm_2002", "c_comm_2005", "c_comm_2008", "c_comm_2012", "c_comm_2015", "c_comm_2016", "c_noncomm_1985", "c_noncomm_1993", "c_noncomm_1995", "c_noncomm_1998", "c_noncomm_2002", "c_noncomm_2005", "c_noncomm_2008", "c_noncomm_2012", "c_noncomm_2015", "c_noncomm_2016", "fc_y_comm_1985", "fc_y_comm_1993", "fc_y_comm_1995", "fc_y_comm_1998", "fc_y_comm_2002", "fc_y_comm_2005", "fc_y_comm_2008", "fc_y_comm_2012", "fc_y_comm_2015", "fc_y_comm_2016", "fc_y_noncomm_1985", "fc_y_noncomm_1993", "fc_y_noncomm_1998", "fc_y_noncomm_2002", "fc_y_noncomm_2008", "fc_y_noncomm_20"))]
ret <- ret[which(!names(ret) %in% c("ART_eff", "ART_recruit_rate_FSW", "ART_recruit_rate_rest", "ART_reinit_rate_FSW", "ART_reinit_rate_rest", "ART_RR_prog", "betaMtoF_baseline", "c_comm_1993_ProFSW", "c_comm_1998_Client", "c_comm_2005_ProFSW", "c_comm_2012_Client", "c_comm_2015_Client", "c_non_comm_1985_Client", "c_non_comm_1985_LowFSW", "c_non_comm_1985_ProFSW", "c_noncomm_1998_GPF", "c_noncomm_1998_GPM", "c_noncomm_2008_GPF", "c_noncomm_2008_GPM", "dropout_rate_FSW", "dropout_rate_not_FSW", "dur_below_200", "eff_ART", "fc_y_comm_1985_ProFSW_Client", "fc_y_comm_1993_ProFSW_Client", "fc_y_comm_2002_ProFSW_Client", "fc_y_noncomm_1985_GPF_GPM", "fc_y_noncomm_1985_ProFSW_Client", "fc_y_noncomm_1998_GPF_GPM", "fc_y_noncomm_2002_ProFSW_Client", "fc_y_noncomm_2011_GPF_GPM", "fraction_F", "fraction_sexually_active_15_F", "fraction_sexually_active_15_M", "init_clientN_from_PCR", "muF", "muM", "n_y_comm_1985_Client_LowFSW", "n_y_comm_1985_Client_ProFSW", "n_y_comm_1985_LowFSW_Client", "n_y_comm_1985_ProFSW_Client", "n_y_noncomm_1985_GPF_GPM", "n_y_noncomm_1985_GPM_GPF"))]
ret <- ret[which(!names(ret) %in% c(unlist(strsplit("n_y_noncomm_2002_ProFSW_Client, n_y_noncomm_2015_ProFSW_Client, prev_HSV2_Client, prev_HSV2_FSW, prev_HSV2_GPF, prev_HSV2_GPM, prev_init_FSW, prev_init_rest, rate_enter_sexual_pop_F, rate_enter_sexual_pop_M, rate_leave_client, rate_leave_low_FSW, RR_beta_circum, RR_beta_FtM, RR_beta_HSV2_comm, RR_beta_HSV2_noncomm, SC_to_death, viral_supp_y_1986_rest, viral_supp_y_2015_ProFSW, movement, fc_y_noncomm_2011, fc_y_noncomm_2015, fc_y_noncomm_2016, n_y_comm_1985, n_y_comm_2002, n_y_comm_2015, n_y_comm_2016, n_y_noncomm_1985, n_y_noncomm_1998, n_y_noncomm_2011, n_y_noncomm_2002, n_y_noncomm_2015, n_y_noncomm_2016, betaMtoF_noncomm, betaFtoM_noncomm, betaMtoF_comm, betaFtoM_comm, RR_beta_HSV2_comm_t, RR_beta_HSV2_noncomm_t, RR_beta_HSV2_comm_a, RR_beta_HSV2_noncomm_a, beta_above_1, ignore_ranges_fc_c, delete, ART_recruit_rate, ART_reinit_rate, prep_dropout", ", "))))]
ret <- ret[which(!names(ret) %in% c(unlist(strsplit("cost_ART_initiation_study_FSW, cost_1_year_ART_study_FSW, cost_1_year_ART_rest, long_intervention, FSW_eligible, GP_eligible, prev_non_ben_fsw_1993, prev_non_ben_fsw_2015, ART_init_ratio_MF, RR_beta_HSV2_a_FSW, RR_beta_HSV2_a_client, RR_beta_HSV2_a_GPF, RR_beta_HSV2_a_GPM, testing_prob_men_2006, testing_prob_men_2008, testing_prob_men_2012, testing_prob_women_2006, testing_prob_women_2008, testing_prob_women_2012, PrEP_reinits_on, infect_ART, alpha33, alpha34, alpha35, FSW_leave_Cotonou_fraction, time, PrEP_loss_to_follow_up, n_y_noncomm_2002_LowFSW_Client, n_y_noncomm_2015_LowFSW_Client, n_y_noncomm_1985_GPF_Client, fc_y_noncomm_1985_GPF_Client, fc_y_noncomm_1998_GPF_Client, fc_y_noncomm_2011_GPF_Client, betaMtoF_FSW, betaFtoM_client, betaMtoF_GPF, betaFtoM_GPM, ART_RR_mort, prop_client_GPM, prop_pro_FSW_GPF, prop_low_FSW_GPF", ", "))))]
ret
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.