# CREATE parTab
# ranges ------------------------------------------------------------------
# yup
ranges = rbind(
testing_prob_men_2006 = c(0.0975, 0.21),
testing_prob_men_2008 = c(0.1, 0.26),
testing_prob_men_2012 = c(0.058, 0.26), # NOTE 0.26 is from MICHEL 2008
testing_prob_women_2006 = c(0.142, 0.4),
testing_prob_women_2008 = c(0.21, 0.54),
testing_prob_women_2012 = c(0.331, 0.513),
infected_FSW_incoming = c(1,1),
n_y_noncomm_1998_GPF_GPM = c(34, 44),
n_y_noncomm_2011_GPF_GPM = c(29, 38),
prev_non_ben_fsw_1993 = c(0.027, 0.163),
# prev_non_ben_fsw_1993 = c(0.12, 0.20),
# prev_non_ben_fsw_2015 = c(0.03, 0.157),
prev_non_ben_fsw_2015 = c(0, 0.046),
# prev_non_ben_fsw_2015 = c(0.10, 0.157), # high
# prev_non_ben_fsw_2015 = c(0.03, 0.046),
# MISC
# init_clientN_from_PCR = c(0,0),
who_believe_comm = c(0, 0),
# # growth rates
# epsilon_1985 = c(0.08, 0.08),
# epsilon_1992 = c(0.08, 0.08),
# epsilon_2002 = c(0.06, 0.07),
# epsilon_2013 = c(0.04, 0.06),
# epsilon_2016 = c(0.04, 0.06),
epsilon_1985 = c(0.059, 0.059),
epsilon_1992 = c(0.048, 0.058),
epsilon_2002 = c(0.027, 0.027),
epsilon_2013 = c(0.027, 0.027),
epsilon_2016 = c(0.027, 0.027),
# DEMOGRAPHIC
fraction_F = c(0.512, 0.52), # fraction of population born female
# frac_women_ProFSW = c(0.0024, 0.0036), # fraction of women that are professional FSW
frac_women_ProFSW = c(0.0024, 0.00715), # fraction of women that are professional FSW
frac_women_LowFSW = c(1, 2), # relative abundance of low FSW relative to pro FSW
frac_men_client = c(0.066, 0.3), # fraction of men that are clients
frac_women_virgin = c(0.079, 0.2), # fraction of women that are virgins
frac_men_virgin = c(0.070, 0.17), # fraction of men that are virgins
# prev_init_FSW = c(0.0132, 0.1), # initial prevalence of FSW
prev_init_FSW = c(0.0132, 0.0659), # initial prevalence of FSW
prev_init_rest = c(0.000313, 0.00294), # initial prevalence of the other groups
muF = c(0.0187, 0.02), # female mortality
muM = c(0.0194, 0.022), # male mortality
# rate_leave_pro_FSW = c(0, 0.125),
rate_leave_pro_FSW = c(0, 0.33), # rate of exit of professional sex work
# rate_leave_low_FSW = c(0, 1), # rate of exit of low level sex work
fraction_FSW_foreign = c(0.5, 0.9),
# fraction_FSW_foreign = c(0.1, 0.5),
rate_leave_client = c(0, 0.295), # rate of exit of clients
# rate_leave_client = c(0, 0.2), # rate of exit of clients
rate_enter_sexual_pop_F = c(0.2, 0.5), # rate of entering sexual population women
rate_enter_sexual_pop_M = c(0.2, 0.5), # rate of entering sexual population men
fraction_sexually_active_15_F = c(0.119, 0.17), # fraction of 15 year old women sexually active
fraction_sexually_active_15_M = c(0.18, 0.35), # fraction of 15 year old men sexually active
# BEHAVIOURAL
# commercial partnerships
c_comm_1993_ProFSW = c(192, 1277),
c_comm_1995_ProFSW = c(192, 1277),
c_comm_2005_ProFSW = c(81, 562),
# c_comm_2015_ProFSW = c(71, 501),
c_comm_1993_LowFSW = c(26, 78),
c_comm_1998_Client = c(8.4, 32),
c_comm_2002_Client = c(11.1, 19.8),
# c_comm_2012_Client = c(11.8, 15),
# c_comm_2015_Client = c(14.5, 19.8),
#non commercial partnerships
c_noncomm_1985_ProFSW = c(0.31, 0.86),
c_noncomm_1985_LowFSW = c(0.41, 1.04),
c_noncomm_1985_Client = c(1.6, 3.3),
c_noncomm_1998_GPF = c(0.93, 0.99),
c_noncomm_2008_GPF = c(0.77, 0.82),
c_noncomm_1998_GPM = c(1.25, 1.43),
c_noncomm_2008_GPM = c(0.73, 0.84),
# sex acts per partnership comm
n_y_comm_1985_ProFSW_Client = c(1, 3.3),
# n_y_comm_1985_Client_ProFSW = c(1.45, 11.45),
# n_y_comm_1985_ProFSW_Client = c(1, 4),
# n_y_comm_2002_ProFSW_Client = c(1, 3),
n_y_comm_1985_LowFSW_Client = c(1, 1),
n_y_comm_1985_Client_LowFSW = c(1, 1),
# sex acts per partnership noncomm
n_y_noncomm_2002_ProFSW_Client = c(13, 20),
n_y_noncomm_2015_ProFSW_Client = c(38.2, 60),
# n_y_noncomm_1985_GPF_GPM = c(39, 100),
# n_y_noncomm_1985_GPM_GPF = c(19.4, 46.7),
#BETA
betaMtoF_baseline = c(0.0006, 0.00109), # baseline male to female transmission rate
RR_beta_FtM = c(0.53, 2), # RR for transmission female to male
# RR_beta_HSV2_comm_a = c(1.4, 2.1), # RR for commercial sex acts where the susceptible individual is infected HSV2
# RR_beta_HSV2_noncomm_a = c(2.2, 3.4), # RR for non commercial sex acts where the susceptible individual is infected HSV2
RR_beta_HSV2_a_FSW = c(0.9, 2.3),
RR_beta_HSV2_a_client = c(1.5, 2.2),
RR_beta_HSV2_a_GPF = c(1.8, 3.4),
RR_beta_HSV2_a_GPM = c(2.2, 4.3),
prev_HSV2_FSW = c(0.87, 0.94), # prevalence HSV2 in FSW
prev_HSV2_Client = c(0.18, 0.28), # prevalence HSV2 in clients
prev_HSV2_GPF = c(0.27, 0.32), # prevalence of HSV2 in GPF
prev_HSV2_GPM = c(0.098, 0.14), # prevalence of HSV2 in GPM
RR_beta_circum = c(0.34, 0.72), # RR for transmission if susceptible individual is circumcised
# Progression parameters
infect_acute = c(4.5, 18.8), # RR for transmission rate if infected is acute stage
infect_AIDS = c(4.5, 11.9), # RR for transmission rate if infected is in AIDS stage
ART_eff = c(0.96, 0.99), # infectiousness RR when on ART (efficacy ART assuimed 90% * % undetectable which is 52.3%)
viral_supp_y_1986_rest = c(0.424, 0.85),
viral_supp_y_2015_ProFSW = c(0.75, 0.85),
ec = c(0.58, 0.95), # condom efficacy
# eP1a = c(0.9, 0.9), # prep efficacy perfect adherence
# eP1b = c(0, 0.9), # prep efficacy intermediate adherence
# eP1c = c(0, 0), # prep efficacy poor adherence
SC_to_death = c(8.7, 12.3),
dur_primary_phase = c(0.25, 0.42),
dur_200_349 = c(2.3, 4.4),
dur_below_200 = c(0.58, 3.17),
alpha03 = c(0.01, 0.05),
alpha04 = c(0.03, 0.1),
ART_RR_prog = c(4.82, 10.23),
# intervention_testing_increase = c(1, 2),
# intervention_testing_increase = c(0.5, 2), # keep
intervention_testing_increase = c(0, 0),
RR_test_CD4200 = c(1, 5.4),
# ART_recruit_rate_FSW = c(0.5, 6),
# ART_recruit_rate_FSW = c(0.5, 1.5),
# ART_recruit_rate_FSW = c(0.5, 3),
ART_recruit_rate_FSW = c(0, 5),
# ART_recruit_rate_rest = c(0.5, 1.5),
# ART_recruit_rate_rest = c(0.5, 6),
ART_recruit_rate_rest = c(3, 12), # NEW
ART_init_ratio_MF = c(1.76, 3.8), # NEW
# intervention_ART_increase = c(0, 12),
# intervention_ART_increase = c(0, 24),
# intervention_ART_increase = c(0.5, 5), # keep
intervention_ART_increase = c(0, 0),
dropout_rate_not_FSW = c(0.0233, 0.11),
dropout_rate_FSW = c(0.0233, 0.11),
ART_reinit_rate_FSW = c(0.25, 1.5),
ART_reinit_rate_rest = c(0.25, 1.5),
# condoms
# fc_y_comm_1985_ProFSW_Client = c(0, 0),
# fc_y_comm_1985_ProFSW_Client = c(0.54, 0.69),
# fc_y_comm_1998_ProFSW_Client = c(0.54, 0.99),
fc_y_comm_1985_ProFSW_Client = c(0, 0.18),
fc_y_comm_1993_ProFSW_Client = c(0.18, 0.33),
fc_y_comm_1998_ProFSW_Client = c(0.4, 0.73),
fc_y_comm_2002_ProFSW_Client = c(0.61, 0.99),
fc_y_comm_2008_ProFSW_Client = c(0.86, 0.99),
fc_y_comm_1985_LowFSW_Client = 0,
fc_y_comm_2015_LowFSW_Client = c(0.25, 0.52),
fc_y_noncomm_1985_ProFSW_Client = 0,
fc_y_noncomm_2002_ProFSW_Client = c(0.19, 0.62),
fc_y_noncomm_1985_LowFSW_Client = 0,
fc_y_noncomm_2015_LowFSW_Client = c(0.138, 0.383),
# fc_y_noncomm_1985_GPF_GPM = 0,
fc_y_noncomm_1998_GPF_GPM = c(0.033, 0.05),
fc_y_noncomm_2011_GPF_GPM = c(0.16, 0.26)
)
# ranges end ------------------------------------------------------------------
# who_believe_comm = c(0, 0), !!!!
parTab_create = data.frame(
names = rownames(ranges),
values = ranges[,1],
fixed = 0,
steps = 0.1,
lower_bound = ranges[,1],
upper_bound = ranges[,2]
)
good_previous_fit = get(load("best_pars_combined_0904_top.Rdata"))
good_previous_fit_variedpars = good_previous_fit[[1]][rownames(ranges)]
j=0
for(i in 1:length(rownames(ranges)))
{
if(parTab_create[i,"names"] %in% names(good_previous_fit_variedpars))
{
parTab_create[i, "values"] = as.numeric(good_previous_fit_variedpars[as.character(parTab_create[i,"names"])][[1]][1]);
j=j+1
}
}
sum(names(good_previous_fit_variedpars) == rownames(ranges), na.rm = T)
parTab_create[parTab_create$lower_bound ==parTab_create$upper_bound, "fixed"] = 1
rownames(parTab_create) = NULL
write.csv(parTab_create, file = "parTab.csv")
###
create_lik <- function(parTab, data, PRIOR_FUNC,...){
epi_start = 1986
# epi_end = 2030
epi_end = 2017
# setup -------------------------------------------------------------------
par_seq = c("c_comm", "c_noncomm")
condom_seq = c("fc_y_comm", "fc_y_noncomm", "n_y_comm", "n_y_noncomm")
groups_seq = c("ProFSW", "LowFSW", "GPF", "FormerFSW", "Client", "GPM", "VirginF", "VirginM", "FormerFSWoutside")
years_seq = seq(1985, 2016)
time <- seq(epi_start, epi_end, length.out = epi_end - epi_start + 1)
time_with_mid <- seq(epi_start, epi_end, length.out = (epi_end - epi_start + 0.5)*2)
#####################################################
# this is the best set of parameters (the fixed ones)
# best_set ----------------------------------------------------------------
best_set = list(
init_clientN_from_PCR = 0,
initial_Ntot = 286114,
frac_women_ProFSW = 0.0024,
frac_women_LowFSW = 0.0027,
frac_women_exFSW = 0.0024,
frac_men_client = 0.2,
frac_women_virgin = 0.1,
frac_men_virgin = 0.1,
prev_init_FSW = 0.0326,
prev_init_rest = 0.0012,
nu = 0.022,
# N_init = c(672, 757, 130895, 672, 27124, 100305, 14544, 11145, 0),
# fraction_F = 0.5,
fraction_F = 0.515666224,
epsilon_1985 = 0.08,
epsilon_1992 = 0.08,
epsilon_2002 = 0.026936907 * 1.5,
epsilon_2013 = 0.026936907 * 1.5,
epsilon_2016 = 0.026936907 * 1.5,
# mu = c(0.02597403, 0.02597403, 0.02597403, 0.02597403, 0.02739726, 0.02739726, 0.02597403, 0.02739726, 0.02597403), # women 1/((27 + 50)/2) # men 1/((25 + 48)/2)
# c_comm = c(750, 52, 0, 0, 13.5, 0, 0, 0, 0),
# c_noncomm = c(0.38, 0.38, 0.88, 0.88, 4, 1.065, 0, 0, 0), # partner change rate lowlevel FSW same as pro, others are approximations from various surveys
#
muF = 0.02597403,
muM = 0.02739726,
# PARTNER CHANGE RATE
c_comm_1985 = c(1229.5, 52, 0, 0, 10.15873, 0, 0, 0, 0), # (1020 + 1439)/2
c_comm_1993 = c(1229.5, 52, 0, 0, 10.15873, 0, 0, 0, 0), # (1020 + 1439)/2
c_comm_1995 = c(1280, 52, 0, 0, 10.15873, 0, 0, 0, 0), # (1135 + 1425)/2
c_comm_1998 = c(881, 52, 0, 0, 10.15873, 0, 0, 0, 0), # (757 + 1005)/2
c_comm_2002 = c(598.5, 52, 0, 0, 11.08109, 0, 0, 0, 0), # (498 + 699)/2, (13.387-10.15873)/14 * 4 + 10.15873
c_comm_2005 = c(424, 52, 0, 0, 11.77286, 0, 0, 0, 0), # (366 + 482)/2, (13.387-10.15873)/14 * 7 + 10.15873
c_comm_2008 = c(371.5, 52, 0, 0, 12.46464, 0, 0, 0, 0), # (272 + 471)/2, (13.387-10.15873)/14 * 10 + 10.15873
c_comm_2012 = c(541, 52, 0, 0, 13.387, 0, 0, 0, 0), # (459 + 623)/2
c_comm_2015 = c(400, 52, 0, 0, 17.15294, 0, 0, 0, 0), # (309 + 491)/2
c_comm_2016 = c(400, 52, 0, 0, 17.15294, 0, 0, 0, 0), # (309 + 491)/2
c_noncomm_1985 = c(0.3766285, 0.3766285, 0.9610526, 0.9610526, 2.028986, 1.337444, 0, 0, 0), # (0.4682779 + 0.3886719 + 0.2729358)/3
c_noncomm_1993 = c(0.3766285, 0.3766285, 0.9610526, 0.9610526, 2.028986, 1.337444, 0, 0, 0),
c_noncomm_1995 = c(0.3766285, 0.3766285, 0.9610526, 0.9610526, 2.028986, 1.337444, 0, 0, 0),
c_noncomm_1998 = c(0.3766285, 0.3766285, 0.9610526, 0.9610526, 2.028986, 1.337444, 0, 0, 0),
c_noncomm_2002 = c(0.3766285, 0.3766285, 0.9610526, 0.9610526, 2.028986, 1.337444, 0, 0, 0),
c_noncomm_2005 = c(0.3766285, 0.3766285, 0.9610526, 0.9610526, 2.028986, 1.337444, 0, 0, 0),
c_noncomm_2008 = c(0.3766285, 0.3766285, 0.7943578, 0.7943578, 2.028986, 0.7878543, 0, 0, 0),
c_noncomm_2012 = c(0.3766285, 0.3766285, 0.7943578, 0.7943578, 8.086957, 0.7878543, 0, 0, 0),
c_noncomm_2015 = c(0.3766285, 0.3766285, 0.7943578, 0.7943578, 6.258258, 0.7878543, 0, 0, 0),
c_noncomm_2016 = c(0.3766285, 0.3766285, 0.7943578, 0.7943578, 6.258258, 0.7878543, 0, 0, 0),
#think about transforming to matrix
betaMtoF_comm = 0.00051, # RR circumcision = 0.44
betaFtoM_comm = 0.02442*0.44,
betaMtoF_noncomm = 0.003,
betaFtoM_noncomm = 0.0038*0.44,
infect_acute = 9, # RR for acute phase
infect_AIDS = 2, #7.27, # RR for AIDS phase
infect_ART = 0.9 * 0.523, # infectiousness RR when on ART (efficacy ART assuimed 90% * % undetectable which is 52.3%)
ec = rep_len(0.8, 9), # from kate's paper on nigeria SD couples
eP0 = c(0, rep_len(0, 8)), # assumptions!
eP1a = c(0.9, rep_len(0, 8)),
eP1b = c(0.45, rep_len(0, 8)),
eP1c = c(0, rep_len(0, 8)),
eP1d = c(0, rep_len(0, 8)),
# gamma01 = 0.4166667, #years
# gamma04 = 4.45, #years
#
alpha01 = rep_len(0, 9),
alpha11 = rep_len(0, 9),
alpha02 = rep_len(0, 9),
alpha03 = 0.03,
alpha04 = 0.07,
alpha05 = 2,
alpha11 = rep_len(0, 9),
alpha22 = rep_len(0, 9),
# alpha23 = rep_len(0.05, 9),
# alpha24 = rep_len(0.08, 9),
# alpha25 = rep_len(0.27, 9),
alpha32 = rep_len(0, 9),
# alpha33 = rep_len(0.05, 9),
# alpha34 = rep_len(0.08, 9),
# alpha35 = rep_len(0.27, 9),
alpha42 = rep_len(0, 9),
# alpha43 = rep_len(0.05, 9),
# alpha44 = rep_len(0.08, 9),
# alpha45 = rep_len(0.27, 9),
test_rate_prep = c(4, 0, 0, 0, 0, 0, 0, 0, 0),
sigma = c(0.82, 0, 0, 0, 0, 0, 0, 0, 0),
prep_intervention_t = c(1985, 2015, 2016, 2016.001),
prep_intervention_y = matrix(c(rep(0, 9), 1, rep(0, 9-1), 1, rep(0, 9-1), rep(0, 9)), ncol = 9, byrow = T),
PrEPOnOff = 0,
#PREP
zetaa_t = c(1985, 2013, 2015, 2016),
zetaa_y = matrix(c(rep(0, 9), 0, rep(0, 9-1), rep(0, 9), rep(0, 9)), ncol = 9, byrow = T),
zetab_t = c(1985, 2013, 2015, 2016),
zetab_y = matrix(c(rep(0, 9), 0, rep(0, 9-1), rep(0, 9), rep(0, 9)), ncol = 9, byrow = T),
zetac_t = c(1985, 2013, 2015, 2016),
zetac_y = matrix(c(rep(0, 9), 0, rep(0, 9-1), rep(0, 9), rep(0, 9)), ncol = 9, byrow = T),
# zetac_y = matrix(c(rep(0, 9), 0.0075, rep(0, 9-1), rep(0, 9), rep(0, 9)), ncol = 9, byrow = T),
psia = rep_len(0.1,9),
psib = rep_len(0.1,9),
#TESTING
testing_prob_t = c(1985, 2001, 2005, 2006, 2008, 2012, 2013, 2015, 2016),
# testing_prob_y = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, # 1985 columns are the risk groups
# 0, 0, 0, 0, 0, 0, 0, 0, 0, # 2001
# 0, 0, 0, 0, 0, 0, 0, 0, 0, # 2005
# 0.142, 0.142, 0.142, 0.142, 0.142, 0.142, 0, 0, 0, # 2006 0.653/8 slope
# 0.21, 0.21, 0.21, 0.21, 0.21, 0.21, 0, 0, 0, # 2008 3*0.653/8
# 0.331, 0.331, 0.331, 0.331, 0.331, 0.331, 0, 0, 0, # 2012 7*0.653/8
# 0.331, 0.331, 0.331, 0.331, 0.331, 0.331, 0, 0, 0, # 2013
# 0.331, 0.331, 0.331, 0.331, 0.331, 0.331, 0, 0, 0, # 2015
# 0.331, 0.331, 0.331, 0.331, 0.331, 0.331, 0, 0, 0), # 2016
# nrow = 9, ncol = 9, byrow = T),
testing_prob_y = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, # 1985 columns are the risk groups
0, 0, 0, 0, 0, 0, 0, 0, 0, # 2001
0, 0.118, 0.118, 0.118, 0.08125, 0.08125, 0, 0, 0, # 2005 0.142*5/6 0.0975*5/6
0.081625, 0.142, 0.142, 0.142, 0.0975, 0.0975, 0, 0, 0, # 2006 0.653/8 slope
0.244875, 0.21, 0.21, 0.21, 0.1, 0.1, 0, 0, 0, # 2008 3*0.653/8
0.571375, 0.331, 0.331, 0.331, 0.0582, 0.0582, 0, 0, 0, # 2012 7*0.653/8
0.653, 0.331, 0.331, 0.331, 0.0582, 0.0582, 0, 0, 0, # 2013
0.68, 0.331, 0.331, 0.331, 0.0582, 0.0582, 0, 0, 0, # 2015
0.68, 0.331, 0.331, 0.331, 0.0582, 0.0582, 0, 0, 0), # 2016
nrow = 9, ncol = 9, byrow = T),
#ART
ART_prob_t = c(1985, 2002, 2005, 2016),
# ART_prob_y = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, # 1985
# 0, 0, 0, 0, 0, 0, 0, 0, 0, # 2002
# 0.1448571, 0.1448571, 0.1448571, 0.1448571, 0.1448571, 0.1448571, 0, 0, 0, # 2005 0.676/14 * 3
# 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0, 0, 0),
# nrow = 4, ncol = 9, byrow = T), # 2016 GP: (0.8+0.552)/2
ART_prob_y = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, # 1985
0, 0, 0, 0, 0, 0, 0, 0, 0, # 2002
0, 0.1448571, 0.1448571, 0.1448571, 0.1448571, 0.1448571, 0, 0, 0, # 2005 0.676/14 * 3
0.6739, 0.676, 0.676, 0.676, 0.676, 0.676, 0, 0, 0),
nrow = 4, ncol = 9, byrow = T), # 2016 GP: (0.8+0.552)/2
RR_ART_CD4200 = 5.39,
# phi2 = c(0.105360516, rep_len(0.025,8)), # former sex workers drop out rate??!
# phi3 = c(0.105360516, rep_len(0.025,8)),
# phi4 = c(0.105360516, rep_len(0.025,8)),
# phi5 = c(0.105360516, rep_len(0.025,8)),
#CONDOM
fc_y_comm_1985 = matrix(
c(0, 0, 0, 0, 0.145524, 0, 0, 0, 0, # 0.145524 is using John's FSW condom 1989 as prop of 1993, * our measure of 1993
0, 0, 0, 0, 0.145524, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0.145524, 0.145524, 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),
fc_y_comm_1993 = matrix(
c(0, 0, 0, 0, 0.536, 0, 0, 0, 0,
0, 0, 0, 0, 0.536, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0.536, 0.536, 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),
fc_y_comm_1995 = matrix(
c(0, 0, 0, 0, 0.536, 0, 0, 0, 0,
0, 0, 0, 0, 0.536, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0.536, 0.536, 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),
fc_y_comm_1998 = matrix(
c(0, 0, 0, 0, 0.536, 0, 0, 0, 0,
0, 0, 0, 0, 0.536, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0.536, 0.536, 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),
fc_y_comm_2002 = matrix(
c(0, 0, 0, 0, 0.8, 0, 0, 0, 0,
0, 0, 0, 0, 0.8, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0.8, 0.8, 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),
fc_y_comm_2005 = matrix(
c(0, 0, 0, 0, 0.8, 0, 0, 0, 0,
0, 0, 0, 0, 0.8, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0.8, 0.8, 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),
fc_y_comm_2008 = matrix(
c(0, 0, 0, 0, 0.8, 0, 0, 0, 0,
0, 0, 0, 0, 0.8, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0.8, 0.8, 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),
fc_y_comm_2012 = matrix(
c(0, 0, 0, 0, 0.8, 0, 0, 0, 0,
0, 0, 0, 0, 0.8, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0.8, 0.8, 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),
fc_y_comm_2015 = matrix(
c(0, 0, 0, 0, 0.8, 0, 0, 0, 0,
0, 0, 0, 0, 0.8, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0.8, 0.8, 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),
fc_y_comm_2015 = matrix(
c(0, 0, 0, 0, 0.8, 0, 0, 0, 0,
0, 0, 0, 0, 0.8, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0,
0.8, 0.8, 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),
fc_y_noncomm_1985 = matrix(
c(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, 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),
fc_y_noncomm_1993 = matrix(
c(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, 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),
# 1998
# (0.33 + 0.2705314)/ 2 # average FSW client
# (0.0326087 + 0.2705314)/ 2 # average client GPF
# (0.0326087 + 0.04989035) / 2 # average gpm gpf
fc_y_noncomm_1998 = matrix(
c(0, 0, 0, 0, 0.3002657, 0, 0, 0, 0,
0, 0, 0, 0, 0.3002657, 0, 0, 0, 0,
0, 0, 0, 0, 0.15157, 0.04124952, 0, 0, 0,
0, 0, 0, 0, 0.15157, 0.04124952, 0, 0, 0,
0.3002657, 0.3002657, 0.15157, 0.15157, 0, 0, 0, 0, 0,
0, 0, 0.04124952, 0.04124952, 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),
# 2008
# (0.33 + 0.4)/ 2 # average FSW client (both approx)
# ((0.05042017+0.241404781)/2 + 0.4)/ 2 # average client GPF (gpf averaged from 2 estimtes)
# ((0.05042017+0.241404781)/2 + (0.07103825+0.34838295)/2) / 2 # average gpm gpf
fc_y_noncomm_2002 = matrix(
c(0, 0, 0, 0, 0.365, 0, 0, 0, 0,
0, 0, 0, 0, 0.365, 0, 0, 0, 0,
0, 0, 0, 0, 0.2729562, 0.1778115, 0, 0, 0,
0, 0, 0, 0, 0.2729562, 0.1778115, 0, 0, 0,
0.365, 0.365, 0.2729562, 0.2729562, 0, 0, 0, 0, 0,
0, 0, 0.1778115, 0.1778115, 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),
fc_y_noncomm_2008 = matrix(
c(0, 0, 0, 0, 0.365, 0, 0, 0, 0,
0, 0, 0, 0, 0.365, 0, 0, 0, 0,
0, 0, 0, 0, 0.2729562, 0.1778115, 0, 0, 0,
0, 0, 0, 0, 0.2729562, 0.1778115, 0, 0, 0,
0.365, 0.365, 0.2729562, 0.2729562, 0, 0, 0, 0, 0,
0, 0, 0.1778115, 0.1778115, 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),
fc_y_noncomm_2011 = matrix(
c(0, 0, 0, 0, 0.365, 0, 0, 0, 0,
0, 0, 0, 0, 0.365, 0, 0, 0, 0,
0, 0, 0, 0, 0.2729562, 0.1778115, 0, 0, 0,
0, 0, 0, 0, 0.2729562, 0.1778115, 0, 0, 0,
0.365, 0.365, 0.2729562, 0.2729562, 0, 0, 0, 0, 0,
0, 0, 0.1778115, 0.1778115, 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),
fc_y_noncomm_2015 = matrix(
c(0, 0, 0, 0, 0.365, 0, 0, 0, 0,
0, 0, 0, 0, 0.365, 0, 0, 0, 0,
0, 0, 0, 0, 0.2729562, 0.1778115, 0, 0, 0,
0, 0, 0, 0, 0.2729562, 0.1778115, 0, 0, 0,
0.365, 0.365, 0.2729562, 0.2729562, 0, 0, 0, 0, 0,
0, 0, 0.1778115, 0.1778115, 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),
fc_y_noncomm_2016 = matrix(
c(0, 0, 0, 0, 0.365, 0, 0, 0, 0,
0, 0, 0, 0, 0.365, 0, 0, 0, 0,
0, 0, 0, 0, 0.2729562, 0.1778115, 0, 0, 0,
0, 0, 0, 0, 0.2729562, 0.1778115, 0, 0, 0,
0.365, 0.365, 0.2729562, 0.2729562, 0, 0, 0, 0, 0,
0, 0, 0.1778115, 0.1778115, 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),
fc_t_comm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016),
fc_t_noncomm = c(1985, 1993, 1998, 2002, 2008, 2011, 2015, 2016),
n_y_comm_1985 = matrix(
c(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, 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),
n_y_comm_2002 = matrix(
c(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, 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),
n_y_comm_2015 = matrix(
c(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, 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),
n_y_comm_2016 = matrix(
c(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, 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),
n_t_comm = c(1985, 2002, 2015, 2016),
n_y_noncomm_1985 = matrix(
c(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, 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),
n_y_noncomm_2002 = matrix(
c(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, 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),
n_y_noncomm_1998 = matrix(
c(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, 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),
n_y_noncomm_2011 = matrix(
c(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, 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),
n_y_noncomm_2015 = matrix(
c(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, 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),
n_y_noncomm_2016 = matrix(
c(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, 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),
n_t_noncomm = c(1985, 1998, 2002, 2011, 2015, 2016),
rate_leave_pro_FSW = 0.2,
FSW_leave_Cotonou_fraction = 0.1,
rate_leave_low_FSW = 0.1,
rate_leave_client = 0.05,
dropout_rate_not_FSW = 0.025,
replaceDeaths = 0,
movement = 1,
ART_recruit_rate_rest = 0.25,
ART_recruit_rate_FSW = 0.25,
ART_reinit_rate_FSW = 0.2,
ART_reinit_rate_rest = 0.2
)
#
# ranges ------------------------------------------------------------------
# yup
ranges = rbind(
testing_prob_men_2006 = c(0.0975, 0.21),
testing_prob_men_2008 = c(0.1, 0.26),
testing_prob_men_2012 = c(0.058, 0.26), # NOTE 0.26 is from MICHEL 2008
testing_prob_women_2006 = c(0.142, 0.4),
testing_prob_women_2008 = c(0.21, 0.54),
testing_prob_women_2012 = c(0.331, 0.513),
infected_FSW_incoming = c(1,1),
n_y_noncomm_1998_GPF_GPM = c(34, 44),
n_y_noncomm_2011_GPF_GPM = c(29, 38),
prev_non_ben_fsw_1993 = c(0.027, 0.163),
# prev_non_ben_fsw_1993 = c(0.12, 0.20),
# prev_non_ben_fsw_2015 = c(0.03, 0.157),
prev_non_ben_fsw_2015 = c(0, 0.046),
# prev_non_ben_fsw_2015 = c(0.10, 0.157), # high
# prev_non_ben_fsw_2015 = c(0.03, 0.046),
# MISC
# init_clientN_from_PCR = c(0,0),
who_believe_comm = c(0, 1),
# # growth rates
# epsilon_1985 = c(0.08, 0.08),
# epsilon_1992 = c(0.08, 0.08),
# epsilon_2002 = c(0.06, 0.07),
# epsilon_2013 = c(0.04, 0.06),
# epsilon_2016 = c(0.04, 0.06),
epsilon_1985 = c(0.059, 0.059),
epsilon_1992 = c(0.048, 0.058),
epsilon_2002 = c(0.027, 0.027),
epsilon_2013 = c(0.027, 0.027),
epsilon_2016 = c(0.027, 0.027),
# DEMOGRAPHIC
fraction_F = c(0.512, 0.52), # fraction of population born female
# frac_women_ProFSW = c(0.0024, 0.0036), # fraction of women that are professional FSW
frac_women_ProFSW = c(0.0024, 0.00715), # fraction of women that are professional FSW
frac_women_LowFSW = c(1, 2), # relative abundance of low FSW relative to pro FSW
frac_men_client = c(0.066, 0.3), # fraction of men that are clients
frac_women_virgin = c(0.079, 0.2), # fraction of women that are virgins
frac_men_virgin = c(0.070, 0.17), # fraction of men that are virgins
# prev_init_FSW = c(0.0132, 0.1), # initial prevalence of FSW
prev_init_FSW = c(0.0132, 0.0659), # initial prevalence of FSW
prev_init_rest = c(0.000313, 0.00294), # initial prevalence of the other groups
muF = c(0.0187, 0.02), # female mortality
muM = c(0.0194, 0.022), # male mortality
# rate_leave_pro_FSW = c(0, 0.125),
rate_leave_pro_FSW = c(0, 0.33), # rate of exit of professional sex work
# rate_leave_low_FSW = c(0, 1), # rate of exit of low level sex work
fraction_FSW_foreign = c(0.5, 0.9),
# fraction_FSW_foreign = c(0.1, 0.5),
rate_leave_client = c(0, 0.295), # rate of exit of clients
# rate_leave_client = c(0, 0.2), # rate of exit of clients
rate_enter_sexual_pop_F = c(0.2, 0.5), # rate of entering sexual population women
rate_enter_sexual_pop_M = c(0.2, 0.5), # rate of entering sexual population men
fraction_sexually_active_15_F = c(0.119, 0.17), # fraction of 15 year old women sexually active
fraction_sexually_active_15_M = c(0.18, 0.35), # fraction of 15 year old men sexually active
# BEHAVIOURAL
# commercial partnerships
c_comm_1993_ProFSW = c(192, 1277),
c_comm_1995_ProFSW = c(192, 1277),
c_comm_2005_ProFSW = c(81, 562),
# c_comm_2015_ProFSW = c(71, 501),
c_comm_1993_LowFSW = c(26, 78),
c_comm_1998_Client = c(8.4, 32),
c_comm_2002_Client = c(11.1, 19.8),
# c_comm_2012_Client = c(11.8, 15),
# c_comm_2015_Client = c(14.5, 19.8),
#non commercial partnerships
c_noncomm_1985_ProFSW = c(0.31, 0.86),
c_noncomm_1985_LowFSW = c(0.41, 1.04),
c_noncomm_1985_Client = c(1.6, 3.3),
c_noncomm_1998_GPF = c(0.93, 0.99),
c_noncomm_2008_GPF = c(0.77, 0.82),
c_noncomm_1998_GPM = c(1.25, 1.43),
c_noncomm_2008_GPM = c(0.73, 0.84),
# sex acts per partnership comm
n_y_comm_1985_ProFSW_Client = c(1, 3.3),
# n_y_comm_1985_Client_ProFSW = c(1.45, 11.45),
# n_y_comm_1985_ProFSW_Client = c(1, 4),
# n_y_comm_2002_ProFSW_Client = c(1, 3),
n_y_comm_1985_LowFSW_Client = c(1, 1),
n_y_comm_1985_Client_LowFSW = c(1, 1),
# sex acts per partnership noncomm
n_y_noncomm_2002_ProFSW_Client = c(13, 20),
n_y_noncomm_2015_ProFSW_Client = c(38.2, 60),
# n_y_noncomm_1985_GPF_GPM = c(39, 100),
# n_y_noncomm_1985_GPM_GPF = c(19.4, 46.7),
#BETA
betaMtoF_baseline = c(0.0006, 0.00109), # baseline male to female transmission rate
RR_beta_FtM = c(0.53, 2), # RR for transmission female to male
# RR_beta_HSV2_comm_a = c(1.4, 2.1), # RR for commercial sex acts where the susceptible individual is infected HSV2
# RR_beta_HSV2_noncomm_a = c(2.2, 3.4), # RR for non commercial sex acts where the susceptible individual is infected HSV2
RR_beta_HSV2_a_FSW = c(0.9, 2.3),
RR_beta_HSV2_a_client = c(1.5, 2.2),
RR_beta_HSV2_a_GPF = c(1.8, 3.4),
RR_beta_HSV2_a_GPM = c(2.2, 4.3),
prev_HSV2_FSW = c(0.87, 0.94), # prevalence HSV2 in FSW
prev_HSV2_Client = c(0.18, 0.28), # prevalence HSV2 in clients
prev_HSV2_GPF = c(0.27, 0.32), # prevalence of HSV2 in GPF
prev_HSV2_GPM = c(0.098, 0.14), # prevalence of HSV2 in GPM
RR_beta_circum = c(0.34, 0.72), # RR for transmission if susceptible individual is circumcised
# Progression parameters
infect_acute = c(4.5, 18.8), # RR for transmission rate if infected is acute stage
infect_AIDS = c(4.5, 11.9), # RR for transmission rate if infected is in AIDS stage
ART_eff = c(0.96, 0.99), # infectiousness RR when on ART (efficacy ART assuimed 90% * % undetectable which is 52.3%)
viral_supp_y_1986_rest = c(0.424, 0.85),
viral_supp_y_2015_ProFSW = c(0.75, 0.85),
ec = c(0.58, 0.95), # condom efficacy
# eP1a = c(0.9, 0.9), # prep efficacy perfect adherence
# eP1b = c(0, 0.9), # prep efficacy intermediate adherence
# eP1c = c(0, 0), # prep efficacy poor adherence
SC_to_death = c(8.7, 12.3),
dur_primary_phase = c(0.25, 0.42),
dur_200_349 = c(2.3, 4.4),
dur_below_200 = c(0.58, 3.17),
alpha03 = c(0.01, 0.05),
alpha04 = c(0.03, 0.1),
ART_RR_prog = c(4.82, 10.23),
# intervention_testing_increase = c(1, 2),
# intervention_testing_increase = c(0.5, 2), # keep
intervention_testing_increase = c(0, 0),
RR_test_CD4200 = c(1, 5.4),
# ART_recruit_rate_FSW = c(0.5, 6),
# ART_recruit_rate_FSW = c(0.5, 1.5),
# ART_recruit_rate_FSW = c(0.5, 3),
ART_recruit_rate_FSW = c(0, 5),
# ART_recruit_rate_rest = c(0.5, 1.5),
# ART_recruit_rate_rest = c(0.5, 6),
ART_recruit_rate_rest = c(3, 12), # NEW
ART_init_ratio_MF = c(1.76, 3.8), # NEW
# intervention_ART_increase = c(0, 12),
# intervention_ART_increase = c(0, 24),
# intervention_ART_increase = c(0.5, 5), # keep
intervention_ART_increase = c(0, 0),
dropout_rate_not_FSW = c(0.0233, 0.11),
dropout_rate_FSW = c(0.0233, 0.11),
ART_reinit_rate_FSW = c(0.25, 1.5),
ART_reinit_rate_rest = c(0.25, 1.5),
# condoms
# fc_y_comm_1985_ProFSW_Client = c(0, 0),
# fc_y_comm_1985_ProFSW_Client = c(0.54, 0.69),
# fc_y_comm_1998_ProFSW_Client = c(0.54, 0.99),
fc_y_comm_1985_ProFSW_Client = c(0, 0.18),
fc_y_comm_1993_ProFSW_Client = c(0.18, 0.33),
fc_y_comm_1998_ProFSW_Client = c(0.4, 0.73),
fc_y_comm_2002_ProFSW_Client = c(0.61, 0.99),
fc_y_comm_2008_ProFSW_Client = c(0.86, 0.99),
fc_y_comm_1985_LowFSW_Client = 0,
fc_y_comm_2015_LowFSW_Client = c(0.25, 0.52),
fc_y_noncomm_1985_ProFSW_Client = 0,
fc_y_noncomm_2002_ProFSW_Client = c(0.19, 0.62),
fc_y_noncomm_1985_LowFSW_Client = 0,
fc_y_noncomm_2015_LowFSW_Client = c(0.138, 0.383),
# fc_y_noncomm_1985_GPF_GPM = 0,
fc_y_noncomm_1998_GPF_GPM = c(0.033, 0.05),
fc_y_noncomm_2011_GPF_GPM = c(0.16, 0.26)
)
# outputs -----------------------------------------------------------------
outputs = c("HIV_positive_On_ART","lambda_sum_0","HIV_positive","pc_of_FOI_on_clients_from_pro_FSW", "S0", "S1a", "S1b", "S1c", "S1d", "prev", "frac_N", "Ntot", "epsilon", "rate_leave_client", "alphaItot", "prev_FSW", "prev_LowFSW", "prev_client", "prev_men", "prev_women", "c_comm_balanced", "c_noncomm_balanced", "who_believe_comm", "ART_coverage_FSW", "ART_coverage_men", "ART_coverage_women", "ART_coverage_all", "rho", "n_comm", "n_noncomm", "fc_comm", "fc_noncomm", "N", "cumuHIVDeaths", "lambda_0", "lambda_1a", "lambda_1b", "lambda_1c", "lambda_1d")
# prev_points -------------------------------------------------------------
prev_points = data.frame(time = c(1986, 1987, 1988, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015,
1998, 2002, 2005, 2008, 2012, 2015,
1998, 2008, 2011,
1998, 2008, 2011,
2012, 2015),
variable = c(rep("Pro FSW", 11),
rep("Clients", 6),
rep("Women", 3),
rep("Men", 3),
rep("Low-level FSW", 2)),
value = c(3.3, 8.2, 19.2, 53.3, 48.7, 40.6, 38.9, 34.8, 29.3, 27.4, 18.7,
100*0.084, 9, 6.9, 5.8, 100*0.028, 100*0.016,
100*0.035, 100*0.04, 2.2,
100*0.033, 100*0.02, 1.6,
100*0.084, 100*0.043),
lower = c(3.3, 8.2, 19.2, 48.02, 43.02, 36.58, 31.97, 30.42, 24.93, 23.01, 15.71,
100*0.05898524, 100*0.068218538, 100*0.04293149, 100*0.034772131, 100*0.012660836, 100*0.006039259,
100*0.024181624, 100*0.030073668, 100*0.012980254,
100*0.022857312, 100*0.012427931, 100*0.007517563,
# 100*0.091838441, 100*0.026704897),
100*0.055700329, 100*0.024043597),
upper = c(3.3, 8.2, 19.2, 58.48, 54.42, 44.67, 46.27, 39.38, 33.88, 32.23, 22.01,
100*0.11561791, 100*0.115608811, 100*0.105215792, 100*0.090216628, 100*0.051602442, 100*0.035338436,
100*0.047726245, 100*0.052817187, 100*0.035296286,
100*0.047183668, 100*0.029774338, 100*0.028546718,
100*0.120857355, 100*0.069311506))
prev_points_all = prev_points
prev_points = prev_points[-c(1,2,3),]
prev_points_low_fsw_2015 = prev_points
prev_points_low_fsw_2015[prev_points_low_fsw_2015$time == 2015 & prev_points_low_fsw_2015$variable == "Pro FSW", "lower"] = 13.79
prev_points_extended_low = data.frame(time = c(1986, 1987, 1988, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015,
1998, 2002, 2005, 2008, 2012, 2015,
1998, 2008, 2011,
1998, 2008, 2011,
2012, 2015),
variable = c(rep("Pro FSW", 11),
rep("Clients", 6),
rep("Women", 3),
rep("Men", 3),
rep("Low-level FSW", 2)),
value = c(3.3, 8.2, 19.2, 53.3, 48.7, 40.6, 38.9, 34.8, 29.3, 27.4, 18.7,
100*0.084, 9, 6.9, 5.8, 100*0.028, 100*0.016,
100*0.035, 100*0.04, 2.2,
100*0.033, 100*0.02, 1.6,
100*0.084, 100*0.043),
lower = c(3.3, 8.2, 19.2, 48.02, 43.02, 36.58, 31.97, 30.42, 24.93, 23.01, 15.71,
100*0.05898524, 100*0.068218538, 100*0.04293149, 100*0.034772131, 100*0.012660836, 100*0.006039259,
100*0.024181624, 100*0.030073668, 100*0.012980254,
100*0.022857312, 100*0.012427931, 100*0.007517563,
100*0, 100*0),
upper = c(3.3, 8.2, 19.2, 58.48, 54.42, 44.67, 46.27, 39.38, 33.88, 32.23, 22.01,
100*0.11561791, 100*0.115608811, 100*0.105215792, 100*0.090216628, 100*0.051602442, 100*0.035338436,
100*0.047726245, 100*0.052817187, 100*0.035296286,
100*0.047183668, 100*0.029774338, 100*0.028546718,
100*0.120857355, 100*0.069311506))
prev_points_FSW_only = data.frame(time = c(1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015
),
variable = c(rep("Pro FSW", 8)
),
value = c(53.3, 48.7, 40.6, 38.9, 34.8, 29.3, 27.4, 18.7
),
lower = c(48.02, 43.02, 36.58, 31.97, 30.42, 24.93, 23.01, 15.71
),
upper = c(58.48, 54.42, 44.67, 46.27, 39.38, 33.88, 32.23, 22.01))
prev_points_FSW_only_all_8_LB = prev_points_FSW_only
prev_points_FSW_only_all_8_LB[prev_points_FSW_only_all_8_LB$time == 2015,"lower"] = 13.79
prev_points_FSW_only_no_2012 <- prev_points_FSW_only_all_8_LB[-7,]
prev_points_FSW_only_even_less_2 = prev_points_FSW_only[c(1, 4, 6, 8),]
prev_points_FSW_Cotonou_centrale_lower_bound = prev_points_FSW_only_even_less_2
prev_points_FSW_Cotonou_centrale_lower_bound[prev_points_FSW_Cotonou_centrale_lower_bound$time == 2015,"lower"] = 13.79
prev_points_FSW_Cotonou_centrale_lower_bound_mid_year = prev_points_FSW_Cotonou_centrale_lower_bound
prev_points_FSW_Cotonou_centrale_lower_bound_mid_year$time = prev_points_FSW_Cotonou_centrale_lower_bound_mid_year$time + 0.5
prev_points_FSW_Cotonou_centrale_lower_bound$x = c(195, 74, 122, 116)
prev_points_FSW_Cotonou_centrale_lower_bound$N = c(366, 190, 417, 620)
# frac N data points ------------------------------------------------------
frac_N_data_points = data.frame(time = c(1998, 2014,
1998, 1998,
1998, 2008, 2011,
1998, 2008, 2011),
point = c(1.43*0.515666224, 0.24*0.515666224,
100*0.195738802*(1-0.515666224), 40*(1-0.515666224),
100*0.1292392*0.515666224, 100*0.0972973*0.515666224, 100*0.18*0.515666224,
100*0.124632*(1-0.515666224), 100*0.08840413*(1-0.515666224), 100*0.1175*(1-0.515666224)),
variable = c("Pro FSW", "Pro FSW",
"Clients", "Clients",
"Virgin female", "Virgin female", "Virgin female",
"Virgin male", "Virgin male", "Virgin male"))
# frac_N_discard_points = data.frame(variable = c("Pro FSW", "Clients", "Virgin female", "Virgin male", "Low-level FSW"),
# min = c(0.001237599, 0.1509*(1-0.515666224), 0.07896475*0.515666224, 0.07039551*(1-0.515666224), 2*0.001237599),
# max = c(0.007374027, 0.40 * (1-0.515666224), 0.2*0.515666224, 0.17*(1-0.515666224), 5*0.007374027))
frac_N_discard_points = data.frame(variable = c("Pro FSW", "Clients", "Virgin female", "Virgin male", "Active FSW", "Low Pro Ratio"),
min = c(0.001237599, 0.074*(1-0.515666224), 0.07896475*0.515666224, 0.07039551*(1-0.515666224), 0.0048*0.516, 1),
max = c(0.0143*0.515666224/2, 0.3 * (1-0.515666224), 0.2*0.515666224, 0.17*(1-0.515666224), 0.0143*0.516, 5))
frac_N_discard_points_no_FSW_LB = data.frame(variable = c("Pro FSW", "Clients", "Virgin female", "Virgin male", "Active FSW", "Low Pro Ratio"),
min = c(0, 0.074*(1-0.515666224), 0.07896475*0.515666224, 0.07039551*(1-0.515666224), 0.0048*0.516, 1),
max = c(0.0143*0.515666224/2, 0.3 * (1-0.515666224), 0.2*0.515666224, 0.17*(1-0.515666224), 0.0143*0.516, 5))
frac_N_discard_points_no_FSW_LB
# Ntot data points ------------------------------------------------------
Ntot_data_points = data.frame(time = c(1992, 2002, 2013, 2020, 2030),
point = c(404359.0418, 681559.032, 913029.606, 1128727.062, 1423887.65),
lower = c(343705.15, 579325.15, 776075.5, 959417.95, 1210304.8),
upper = c(465012.85, 783792.85, 1049984.5, 1298036.05, 1637471.2),
colour = c("data", "data", "data", "predicted", "predicted"))
# ART coverage data points ------------------------------------------------
ART_data_points = data.frame(time = c(2011, 2012, 2013, 2014, 2015, 2016, 2017,
2012, 2014, 2015, 2016, 2017
),
Lower = c(0.33, 0.42, 0.44, 0.39, 0.44, 0.51, 0.57,
0.09, 0.14, 0.2, 0.22, 0.3
),
Upper = c(0.52, 0.66, 0.69, 0.61, 0.69, 0.8, 0.8,
0.13, 0.2, 0.28, 0.3, 0.42
),
variable = c("All", "All", "All", "All", "All", "All", "All",
"Pro FSW", "Pro FSW", "Pro FSW", "Pro FSW", "Pro FSW"))
ART_data_points_extra = data.frame(time = c(2011, 2012, 2013, 2014, 2015, 2016, 2017,
2012, 2014, 2015, 2016, 2017, 2017.5
),
Lower = c(0.33, 0.42, 0.44, 0.39, 0.44, 0.51, 0.57,
0.08, 0.14, 0.2, 0.50, 0.56, 0.8
),
Upper = c(0.52, 0.66, 0.69, 0.61, 0.69, 0.8, 0.8,
0.11, 0.2, 0.28, 0.7, 0.79, 0.9
),
variable = c("All", "All", "All", "All", "All", "All", "All",
"Pro FSW", "Pro FSW", "Pro FSW", "Pro FSW", "Pro FSW", "Pro FSW"))
ART_data_points_FSW = ART_data_points[c(8, 9, 10, 11, 12),]
ART_data_points_FSW_last_2 = ART_data_points[c(11, 12),]
ART_data_points_first_and_last_FSW = ART_data_points[c(8, 12),]
# first and last GP, first FSW and 2016, 2017 FSW
ART_data_points_1611 = ART_data_points[c(1, 7, 8, 11, 12),]
# first and last GP, all fsw before intervention
ART_data_points_1911 = ART_data_points[c(1, 7, 8, 9, 10),]
# NUMBERS OF FSW ON ART
ART_data_points_with_numbers = data.frame(time = c(2011, 2017,
2012, 2014, 2015#, 2016, 2017
),
Lower = c(0.33, 0.57,
27, 34, 42#, 45, 62
),
Upper = c(0.52, 0.8,
37, 46, 56#, 83, 107
),
variable = c("All", "All",
"Numbers FSW", "Numbers FSW", "Numbers FSW"#, "Numbers FSW", "Numbers FSW"
))
ART_data_points_with_numbers_reduced = ART_data_points_with_numbers[c(1,2,5),]
ART_data_points_NULL = data.frame(time = c(2011
),
Lower = c(-1
),
Upper = c(1
),
variable = c("All"))
# PrEP_fitting ------------------------------------------------
# PrEP_fitting = data.frame(time = c(2016, 2017, 2016, 2017, 2016, 2017),
# group = c("S1a", "S1a", "S1b", "S1b", "S1c", "S1c"),
# lower = c(50, 40, 50, 40, 50, 40),
# point = c(55, 45, 50, 40, 50, 40),
# upper = c(61, 66, 61, 66, 61, 66)
#
#
# )
#
# PrEP_fitting = data.frame(time = c(2016, 2017),
# group = c("S1a", "S1a"),
# lower = c(50, 40),
# point = c(55, 45),
# upper = c(61, 66)
#
#
# )
PrEP_fitting = NULL
#
#
# ART_data_points_FSW = data.frame(time = c(2010, 2013, 2014, 2015, 2016, 2017
# ),
# Lower = c(0.07, 0.13, 0.17, 0.50, 0.56, 0.80
#
# ),
# Upper = c(0.10, 0.18, 0.24, 0.7, 0.79, 0.9
#
# ),
# variable = c("Pro FSW", "Pro FSW", "Pro FSW", "Pro FSW", "Pro FSW", "Pro FSW"))
# ART_data_points_FSW_last_3 = data.frame(time = c(2015, 2016, 2017
# ),
# Lower = c(0.50, 0.56, 0.80
#
# ),
# Upper = c(0.7, 0.79, 0.9
#
# ),
# variable = c("Pro FSW", "Pro FSW", "Pro FSW"))
# ART_data_points_first_and_last = data.frame(time = c(2010, 2017
# ),
# Lower = c(0.07, 0.80
#
# ),
# Upper = c(0.10, 0.9
#
# ),
# variable = c("Pro FSW", "Pro FSW"))
# ART_data_points_all = data.frame(time = c(2010, 2011, 2012, 2013, 2014, 2015, 2016,
# 2010, 2011, 2012, 2013, 2014, 2015, 2016,
# 2010, 2011, 2012, 2013, 2014, 2015, 2016,
# 2010, 2013, 2014, 2015, 2016
# ),
# Lower = c(0.32, 0.4, 0.43, 0.38, 0.43, 0.49, 0.552,
# 0.32, 0.4, 0.43, 0.38, 0.43, 0.49, 0.552,
# 0.32, 0.4, 0.43, 0.38, 0.43, 0.49, 0.552,
# 7.5, 12.5, 17.4, 50.0, 56.1
#
# ),
# Upper = c(0.522449, 0.653061, 0.702041, 0.620408, 0.702041, 0.8, 0.8,
# 0.522449, 0.653061, 0.702041, 0.620408, 0.702041, 0.8, 0.8,
# 0.522449, 0.653061, 0.702041, 0.620408, 0.702041, 0.8, 0.8,
# 0.522449, 0.620408, 0.702041, 0.8, 0.8
#
# ),
# variable = c("Women", "Women", "Women", "Women", "Women", "Women", "Women",
# "Men", "Men", "Men", "Men", "Men", "Men", "Men",
# "All", "All", "All", "All", "All", "All", "All",
# "Pro FSW", "Pro FSW", "Pro FSW", "Pro FSW", "Pro FSW"))
#####################################################
ART_data_points_lazymcmc = data.frame(time = c(2017, 2017, 2015),
variable = c("Men", "Women","Pro FSW"),
x = c(3908, 9784, 49))
number_simulations=1
par_names <- rownames(ranges)
# pars = parTab_create$values
## using the `...` bit.
likelihood_func <- function(pars){
# put pars into way in which the model can read it
# names(pars) = par_names
ranges_new = cbind(pars, pars)
rownames(ranges_new) = par_names
parameters = cotonou::lhs_parameters(number_simulations, set_pars = best_set, Ncat = 9, time = time,
ranges = ranges_new, par_seq = par_seq, condom_seq = condom_seq, groups_seq = groups_seq, years_seq = years_seq)
#ranges = NULL needs to be changed to incorporate those that vary!
# run model
# res = lapply(parameters, cotonou::return_outputs, cotonou::main_model, time = time, outputs = outputs)
res = cotonou::return_outputs(parameters[[1]], cotonou::main_model, time = time, outputs = outputs)
# calculate likelihood
# likelihood_list = lapply(res, cotonou::likelihood_lazymcmc, time = time,
# prev_points = prev_points_FSW_Cotonou_centrale_lower_bound,
# frac_N_discard_points = frac_N_discard_points_no_FSW_LB,
# Ntot_data_points = Ntot_data_points,
# ART_data_points = ART_data_points_lazymcmc,
# PrEP_fitting = PrEP_fitting)
lik = cotonou::likelihood_lazymcmc(res, time = time,
prev_points = prev_points_FSW_Cotonou_centrale_lower_bound,
frac_N_discard_points = frac_N_discard_points_no_FSW_LB,
Ntot_data_points = Ntot_data_points,
ART_data_points = ART_data_points_lazymcmc,
PrEP_fitting = PrEP_fitting)
# lik = likelihood_list[[1]]
return(lik)
# return(par_names)
}
return(likelihood_func)
}
parTab <- read.csv("parTab.csv")
# f <- create_lik(parTab, NULL, NULL)
# f(unlist(pars))
f <- create_lik(parTab)
f(parTab$values)
devtools::install_github("jameshay218/lazymcmc")
library(lazymcmc)
before = Sys.time()
mcmcPars <- c("iterations"=25000,popt=0.44,opt_freq=1000,thin=10,burnin=0,adaptive_period=4000,save_block=100)
Rprof(tmp <- tempfile())
res <- run_MCMC(parTab,data,mcmcPars,"test",create_lik,NULL,NULL,0.1)
Rprof()
summaryRprof(tmp)
after = Sys.time()
after - before
chain <- read.csv(res$file)
plot(coda::as.mcmc(chain))
# ____ _____ ____ _ _ _ _____ ____ ___ _____ __ __ ____ __ __ ____
# | _ \| ____/ ___|| | | | | |_ _/ ___| / _ \| ___| | \/ |/ ___| \/ |/ ___|
# | |_) | _| \___ \| | | | | | | \___ \ | | | | |_ | |\/| | | | |\/| | |
# | _ <| |___ ___) | |_| | |___| | ___) | | |_| | _| | | | | |___| | | | |___
# |_| \_\_____|____/ \___/|_____|_| |____/ \___/|_| |_| |_|\____|_| |_|\____|
batch_folder = "0806"
mcmc_results = read.csv("0106_many_chains_actual_MCMC_top_fit_5_repetition_1_univariate_chain.csv")
plot(mcmc_results[,"lnlike"])
plot(as.mcmc(mcmc_results))
#from 100 onwards?
# mcmc_results_without_burnin = mcmc_results[c((length(mcmc_results[,1])-1999):length(mcmc_results[,1])),]
mcmc_results_without_burnin = mcmc_results[c(99:length(mcmc_results[,1])),]
mcmc_results_without_burnin = mcmc_results_without_burnin[seq(1, length(mcmc_results_without_burnin[,1]), 20),]
# test_set = mcmc_results_without_burnin[45,]
# names(test_set)
epi_start = 1986
# epi_end = 2030
epi_end = 2035
# setup -------------------------------------------------------------------
par_seq = c("c_comm", "c_noncomm")
condom_seq = c("fc_y_comm", "fc_y_noncomm", "n_y_comm", "n_y_noncomm")
groups_seq = c("ProFSW", "LowFSW", "GPF", "FormerFSW", "Client", "GPM", "VirginF", "VirginM", "FormerFSWoutside")
years_seq = seq(1985, 2016)
time <- seq(epi_start, epi_end, length.out = epi_end - epi_start + 1)
# mcmc_results_without_burnin = mcmc_results_without_burnin[c(1:3),]
parameters = list()
for(j in 1:length(mcmc_results_without_burnin[,1]))
{
ranges_post = ranges
for(i in 1:length(ranges[,1])) {
ranges_post[i,1] = as.numeric(mcmc_results_without_burnin[j,][rownames(ranges)[i]])
ranges_post[i,2] = ranges_post[i,1]
}
parameters[[j]] = cotonou::lhs_parameters(number_simulations, set_pars = best_set, Ncat = 9, time = time,
ranges = ranges_post, par_seq = par_seq, condom_seq = condom_seq, groups_seq = groups_seq, years_seq = years_seq)
}
parameters = lapply(parameters, function(x) x[[1]])
# ################ CEA ----------------------------------------------------
CEA_outputs = unique(c(
"prev_non_ben_fsw_1993",
"prev_non_ben_fsw_2015",
"testpar","pfFSW", "prop_FSW_I0_1", "prop_FSW_I0_2", "prop_FSW_I0_3", "prop_FSW_I0_4", "prop_FSW_I0_5",
"gamma32_without_supp",
"gamma33_without_supp",
"gamma34_without_supp",
"alpha33_without_supp",
"alpha34_without_supp",
"alpha35_without_supp",
"gamma32",
"gamma33",
"gamma34",
"alpha33",
"alpha34",
"alpha35",
"cost_Initiation_of_ART_study_FSW",
"cost_Initiation_of_ART_government_FSW",
"cost_1_year_of_ART_study_FSW",
"cost_1_year_of_ART_government_FSW",
"cost_Initiation_ART_rest_of_population",
"cost_1_year_of_ART_rest_of_population",
"cost_FSW_initiation_ART_Patient_costs",
"cost_FSW_1_year_ART_Patient_costs",
"cost_Initiation_of_PrEP_study",
"cost_1_year_PrEP_perfect_adherence_study",
"cost_1_year_PrEP_intermediate_adherence_study",
"cost_1_year_PrEP_non_adherence_study",
"cost_Initiation_of_PrEP_government",
"cost_1_year_PrEP_perfect_adherence_government",
"cost_1_year_PrEP_intermediate_adherence_government",
"cost_1_year_PrEP_non_adherence_government",
"cost_PREP_initiation_Patient_costs",
"cost_PREP_1_year_ART_Patient_costs",
"W0", "W1", "W2", "W3", "Number_DALY_W1","Number_DALY_W2", "Number_DALY_W3", "FSW_On_PrEP_all_cats", "PrEPinitiations", "PrEPinitiations1a",
"PrEPinitiations1b", "PrEPinitiations1c", "pc_susceptible_FSW_On_PrEP", "pc_all_FSW_On_PrEP", "Number_Susceptibles",
"HIV_positive_On_ART", "HIV_positive_Diagnosed_Off_ART", "Primary_Off_ART",
"CD4_above_500_Off_ART", "CD4_350_500_Off_ART", "CD4_200_350_Off_ART", "CD4_below_200_Off_ART", "cumuDeaths_On_ART", "HIV_positive", "ec", "cumuARTinitiations","cumuARTREinitiations", "rate_leave_pro_FSW","tau_intervention",
"testing_prob", "tau", "N", "S0", "S1a", "S1b", "S1c", "S1d", "I01", "I11", "I02", "I03", "I04",
"I05", "I22", "I23", "I24", "I25", "I32", "I33", "I34", "I35", "I42", "I43", "I44", "I45", "prev",
"frac_N", "Ntot", "epsilon", "rate_leave_client", "alphaItot", "prev_FSW", "prev_LowFSW", "prev_client",
"prev_men", "prev_women", "c_comm_balanced", "c_noncomm_balanced", "who_believe_comm", "ART_coverage_FSW",
"ART_coverage_men", "ART_coverage_women", "ART_coverage_all", "rho", "n_comm", "n_noncomm", "fc_comm",
"fc_noncomm", "N", "cumuHIVDeaths", "lambda_sum_0", "lambda_sum_1a", "lambda_sum_1b", "lambda_sum_1c",
"lambda_sum_1d", "S0", "S1a", "S1b", "S1c", "S1d", "OnPrEP1a", "OnPrEP1b",
"OnPrEP1c", "ART_eligible_CD4_above_500", "ART_eligible_CD4_350_500","ART_eligible_CD4_200_349","ART_eligible_CD4_below_200",
"cumuAllDeaths", "cumuHIVDeaths", "cumuARTinitiations", "cumuARTREinitiations",
"OnPrEP", "ART_sex_ratio", "pc_S1b", "pc_S1a", "pc_S1c", "cumuInf",
"intervention_ART_increase", "testing_prob", "rho_intervention",
"ART_eligible_CD4_above_500", "ART_eligible_CD4_350_500", "ART_eligible_CD4_200_349",
"ART_eligible_CD4_below_200", "new_people_in_group_FSW_only", "rate_move_out", "rate_move_in",
"FSW_out", "FSW_in", "zeta", "tau", "prep_offering_rate", "intervention_testing_increase", "sigma",
"PrEPOnOff", "prev", "frac_N", "Ntot", "epsilon", "rate_leave_client", "alphaItot", "prev_FSW",
"prev_LowFSW", "prev_client", "prev_men", "prev_women", "c_comm_balanced", "c_noncomm_balanced",
"who_believe_comm", "ART_coverage_FSW", "ART_coverage_men", "ART_coverage_women", "ART_coverage_all",
"rho", "n_comm", "n_noncomm", "fc_comm", "fc_noncomm", "N", "cumuHIVDeaths", "lambda_0", "lambda_1a",
"lambda_1b", "lambda_1c", "lambda_1d"))
# ################ CEA ----------------------------------------------------
########################################################################################## CEA
res_best_runs = lapply(parameters, cotonou::return_outputs, cotonou::main_model, time = time, outputs = CEA_outputs)
# res_best_runs[[1]]$lambda_sum_0[which(time == 2013)]
lapply(res_best_runs, function(x) x$lambda_sum_0[which(time == 2013)])
# res_best_runs[[1]] = res
# ignore these ######################################
# frac_ProFSW_F = data.frame(time, t(apply(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$frac_N), function(x) (x[,1]/(x[,1] + x[,2] + x[,3] + x[,4] + x[,7])))), 2, cotonou::quantile_95)))
# frac_ProFSW_F = data.frame(time, t(apply(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$frac_N/x$frac_F), function(x) x[,1])), 2, cotonou::quantile_95)))
annual_client_volume_pro_FSW = data.frame(time, t(do.call(rbind, lapply(res_best_runs, function(x) {x$c_comm[,1]}))))
colnames(annual_client_volume_pro_FSW) = c("time", as.character(seq(1, (length(annual_client_volume_pro_FSW[1,])-1))))
annual_client_volume_pro_FSW_melted = reshape2::melt(annual_client_volume_pro_FSW, id.vars = c("time"))
N_Pro_FSW = data.frame(time, t(do.call(rbind, lapply(lapply(res_best_runs, function(x) {x$N}), function(x) {return(x[,c(1)])}))))
colnames(N_Pro_FSW) = c("time", as.character(seq(1, (length(N_Pro_FSW[1,])-1))))
N_Pro_FSW_melted = reshape2::melt(N_Pro_FSW, id.vars = c("time"))
colnames(N_Pro_FSW_melted) = c("time", "run", "value")
N_Low_FSW = data.frame(time, t(do.call(rbind, lapply(lapply(res_best_runs, function(x) {x$N}), function(x) {return(x[,c(2)])}))))
colnames(N_Low_FSW) = c("time", as.character(seq(1, (length(N_Low_FSW[1,])-1))))
N_Low_FSW_melted = reshape2::melt(N_Low_FSW, id.vars = c("time"))
colnames(N_Low_FSW_melted) = c("time", "run", "value")
N_Client = data.frame(time, t(do.call(rbind, lapply(lapply(res_best_runs, function(x) {x$N}), function(x) {return(x[,c(5)])}))))
colnames(N_Client) = c("time", as.character(seq(1, (length(N_Client[1,])-1))))
N_Client_melted = reshape2::melt(N_Client, id.vars = c("time"))
colnames(N_Client_melted) = c("time", "run", "value")
Fraction_F = do.call(rbind, lapply(lapply(lapply(res_best_runs, function(x) {x$N}), function(x) {return(x[,c(1, 2, 3, 4, 7)])}), rowSums)) /
(do.call(rbind, lapply(lapply(lapply(res_best_runs, function(x) {x$N}), function(x) {return(x[,c(5, 6)])}), rowSums)) + do.call(rbind, lapply(lapply(lapply(res_best_runs, function(x) {x$N}), function(x) {return(x[,c(1, 2, 3, 4, 7)])}), rowSums)))
lambda_sum_0_ProFSW = data.frame(t(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$lambda_sum_0), function(x) x[,1]))))
lambda_sum_0_LowFSW = data.frame(t(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$lambda_sum_0), function(x) x[,2]))))
lambda_sum_0_GPF = data.frame(t(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$lambda_sum_0), function(x) x[,3]))))
lambda_sum_0_FormerFSW = data.frame(t(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$lambda_sum_0), function(x) x[,4]))))
lambda_sum_0_Client = data.frame(t(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$lambda_sum_0), function(x) x[,5]))))
lambda_sum_0_GPM = data.frame(t(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$lambda_sum_0), function(x) x[,6]))))
lambda_sum_0_Virgin_Female = data.frame(t(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$lambda_sum_0), function(x) x[,7]))))
lambda_sum_0_Virgin_Male = data.frame(t(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$lambda_sum_0), function(x) x[,8]))))
lambda_sum_0_Former_FSW_Outside = data.frame(t(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$lambda_sum_0), function(x) x[,9]))))
lambda_sum_0_indiv = rbind(lambda_sum_0_ProFSW, lambda_sum_0_LowFSW, lambda_sum_0_GPF,
lambda_sum_0_FormerFSW, lambda_sum_0_Client, lambda_sum_0_GPM)
lambda_sum_0_indiv = data.frame(time, rep(c("Pro FSW", "Low-level FSW", "GPF", "Former FSW in Cotonou", "Clients", "GPM"), each = length(time)), lambda_sum_0_indiv)
colnames(lambda_sum_0_indiv) = c("time", "variable", as.character(seq(1, length(lambda_sum_0_ProFSW[1,]))))
lambda_sum_0_indiv_melted = reshape2::melt(lambda_sum_0_indiv, id.vars = c("time", "variable"))
colnames(lambda_sum_0_indiv_melted) = c("time", "variable", "run", "value")
lambda_sum_0_indiv_melted$variable = factor(lambda_sum_0_indiv_melted$variable, levels = c("Pro FSW", "Low-level FSW", "GPF", "Former FSW in Cotonou", "Clients", "GPM"))
frac_ProFSW = data.frame(time, t(apply(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$frac_N*100), function(x) x[,1])), 2, cotonou::quantile_95)))
frac_LowFSW = data.frame(time, t(apply(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$frac_N*100), function(x) x[,2])), 2, cotonou::quantile_95)))
frac_GPF = data.frame(time, t(apply(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$frac_N*100), function(x) x[,3])), 2, cotonou::quantile_95)))
frac_FormerFSW = data.frame(time, t(apply(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$frac_N*100), function(x) x[,4])), 2, cotonou::quantile_95)))
frac_Client = data.frame(time, t(apply(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$frac_N*100), function(x) x[,5])), 2, cotonou::quantile_95)))
frac_GPM = data.frame(time, t(apply(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$frac_N*100), function(x) x[,6])), 2, cotonou::quantile_95)))
frac_Virgin_Female = data.frame(time, t(apply(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$frac_N*100), function(x) x[,7])), 2, cotonou::quantile_95)))
frac_Virgin_Male = data.frame(time, t(apply(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$frac_N*100), function(x) x[,8])), 2, cotonou::quantile_95)))
frac_Former_FSW_Outside = data.frame(time, t(apply(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$frac_N*100), function(x) x[,9])), 2, cotonou::quantile_95)))
frac_Active_FSW = data.frame(time, t(apply(do.call(rbind, lapply(res_best_runs, function(x) {100*(x$frac_N[,1] + x$frac_N[,2])})), 2, cotonou::quantile_95)))
Ratio_Low_Pro = data.frame(time, t(apply(do.call(rbind, lapply(res_best_runs, function(x) {x$frac_N[,2]/ x$frac_N[,1]})), 2, cotonou::quantile_95)))
frac = rbind(frac_ProFSW, frac_LowFSW, frac_GPF, frac_FormerFSW, frac_Client, frac_GPM, frac_Virgin_Female, frac_Virgin_Male, frac_Former_FSW_Outside, frac_Active_FSW, Ratio_Low_Pro)
frac = data.frame(frac, group = rep(c("Pro FSW", "Low-level FSW", "GPF", "Former FSW in Cotonou", "Clients", "GPM", "Virgin female", "Virgin male", "Former FSW outside Cotonou", "Active FSW", "Low Pro Ratio"), each = length(time)))
colnames(frac) = c("time", "Lower", "Median", "Upper", "variable")
frac$variable = factor(frac$variable, levels = c("Pro FSW", "Low-level FSW", "GPF", "Former FSW in Cotonou", "Clients", "GPM", "Virgin female", "Virgin male", "Former FSW outside Cotonou", "Active FSW", "Low Pro Ratio"))
prev_FSW = t(apply(do.call(rbind, lapply(res_best_runs, function(x) x$prev_FSW)), 2, cotonou::quantile_95))
prev_LowFSW = t(apply(do.call(rbind, lapply(res_best_runs, function(x) x$prev_LowFSW)), 2, cotonou::quantile_95))
prev_client = t(apply(do.call(rbind, lapply(res_best_runs, function(x) x$prev_client)), 2, cotonou::quantile_95))
prev_women = t(apply(do.call(rbind, lapply(res_best_runs, function(x) x$prev_women)), 2, cotonou::quantile_95))
prev_men = t(apply(do.call(rbind, lapply(res_best_runs, function(x) x$prev_men)), 2, cotonou::quantile_95))
prev = rbind(prev_FSW, prev_LowFSW, prev_client, prev_women, prev_men)
prev = data.frame(time, prev, rep(c("Pro FSW", "Low-level FSW", "Clients", "Women", "Men"), each = length(time)))
colnames(prev) = c("time", "Lower", "Median", "Upper", "variable")
prev$variable = factor(prev$variable, levels = c("Pro FSW", "Low-level FSW", "Clients", "Women", "Men"))
prev_FSW_indiv = t(do.call(rbind, lapply(res_best_runs, function(x) x$prev_FSW)))
prev_LowFSW_indiv = t(do.call(rbind, lapply(res_best_runs, function(x) x$prev_LowFSW)))
prev_client_indiv = t(do.call(rbind, lapply(res_best_runs, function(x) x$prev_client)))
prev_women_indiv = t(do.call(rbind, lapply(res_best_runs, function(x) x$prev_women)))
prev_men_indiv = t(do.call(rbind, lapply(res_best_runs, function(x) x$prev_men)))
prev_indiv = rbind(prev_FSW_indiv, prev_LowFSW_indiv, prev_client_indiv, prev_women_indiv, prev_men_indiv)
prev_indiv = data.frame(time, rep(c("Pro FSW", "Low-level FSW", "Clients", "Women", "Men"), each = length(time)), prev_indiv)
colnames(prev_indiv) = c("time", "variable", as.character(seq(1, length(prev_FSW_indiv[1,]))))
prev_indiv_melted = reshape2::melt(prev_indiv, id.vars = c("time", "variable"))
colnames(prev_indiv_melted) = c("time", "variable", "run", "value")
prev_indiv_melted$variable = factor(prev_indiv_melted$variable, levels = c("Pro FSW", "Low-level FSW", "Clients", "Women", "Men"))
Ntot = data.frame(time, t(apply(do.call(rbind, lapply(res_best_runs, function(x) x$Ntot)), 2, cotonou::quantile_95)))
colnames(Ntot) = c("time", "Lower", "Median", "Upper")
ART_coverage_women = t(apply(do.call(rbind, lapply(res_best_runs, function(x) x$ART_coverage_women)), 2, cotonou::quantile_95))
ART_coverage_men = t(apply(do.call(rbind, lapply(res_best_runs, function(x) x$ART_coverage_men)), 2, cotonou::quantile_95))
ART_coverage_FSW = t(apply(do.call(rbind, lapply(res_best_runs, function(x) x$ART_coverage_FSW)), 2, cotonou::quantile_95))
ART_coverage_all = t(apply(do.call(rbind, lapply(res_best_runs, function(x) x$ART_coverage_all)), 2, cotonou::quantile_95))
ART_coverage = rbind(ART_coverage_women, ART_coverage_men, ART_coverage_FSW, ART_coverage_all)
ART_coverage = data.frame(time, ART_coverage, rep(c("Women", "Men", "Pro FSW", "All"), each = length(time)))
colnames(ART_coverage) = c("time", "Lower", "Median", "Upper", "variable")
ART_coverage$variable = factor(ART_coverage$variable, levels = c("Pro FSW", "Women", "Men", "All"))
ART_coverage = ART_coverage[ART_coverage$variable == "All" | ART_coverage$variable == "Pro FSW",]
ART_FSW_indiv = t(do.call(rbind, lapply(res_best_runs, function(x) x$ART_coverage_FSW)))
ART_women_indiv = t(do.call(rbind, lapply(res_best_runs, function(x) x$ART_coverage_women)))
ART_men_indiv = t(do.call(rbind, lapply(res_best_runs, function(x) x$ART_coverage_men)))
ART_all_indiv = t(do.call(rbind, lapply(res_best_runs, function(x) x$ART_coverage_all)))
ART_indiv = rbind(ART_FSW_indiv, ART_all_indiv)
ART_indiv = data.frame(time, rep(c("Pro FSW", "All"), each = length(time)), ART_indiv)
colnames(ART_indiv) = c("time", "variable", as.character(seq(1, length(ART_FSW_indiv[1,]))))
ART_indiv_melted = reshape2::melt(ART_indiv, id.vars = c("time", "variable"))
colnames(ART_indiv_melted) = c("time", "variable", "run", "value")
ART_indiv_melted$variable = factor(ART_indiv_melted$variable, levels = c("Pro FSW", "All"))
# N of FSW on ART
N_ART_FSW_indiv = t(do.call(rbind, lapply(res_best_runs, function(x) x$HIV_positive_On_ART[,1])))
N_ART_FSW_indiv = data.frame(time, rep(c("Pro FSW"), each = length(time)), N_ART_FSW_indiv)
colnames(N_ART_FSW_indiv) = c("time", "variable", as.character(seq(1, length(N_ART_FSW_indiv[1,])-2)))
N_ART_FSW_indiv_melted = reshape2::melt(N_ART_FSW_indiv, id.vars = c("time", "variable"))
colnames(N_ART_FSW_indiv_melted) = c("time", "variable", "run", "value")
# N on ART and N off ART diagnosed
Diagnosed_Off_ART_FSW = data.frame(time, t(do.call(rbind, lapply(lapply(res_best_runs, function(x) {x$I22 + x$I23 + x$I24 + x$I25}), function(x) x[,1]))))
Diagnosed_On_ART_FSW = data.frame(time, t(do.call(rbind, lapply(lapply(res_best_runs, function(x) {x$I32 + x$I33 + x$I34 + x$I35}), function(x) x[,1]))))
Diagnosed_Dropout_ART_FSW = data.frame(time, t(do.call(rbind, lapply(lapply(res_best_runs, function(x) {x$I42 + x$I43 + x$I44 + x$I45}), function(x) x[,1]))))
Diagnosed_ALL_FSW = data.frame(time, t(do.call(rbind, lapply(lapply(res_best_runs, function(x) {x$I22 + x$I23 + x$I24 + x$I25 + x$I32 + x$I33 + x$I34 + x$I35 + x$I42 + x$I43 + x$I44 + x$I45}), function(x) x[,1]))))
Diagnosed_FSW = data.frame(rbind(Diagnosed_Off_ART_FSW, Diagnosed_On_ART_FSW, Diagnosed_Dropout_ART_FSW, Diagnosed_ALL_FSW),
rep(c("Diagnosed Off ART", "Diagnosed On ART", "Dropout", "All diagnosed"), each = length(time)))
colnames(Diagnosed_FSW) = c("time", as.character(seq(1, (length(Diagnosed_Off_ART_FSW[1,])-1))), "variable")
Diagnosed_FSW_melted = reshape2::melt(Diagnosed_FSW, id.vars = c("time", "variable"))
Diagnosed_FSW_melted$group = "FSW"
colnames(Diagnosed_FSW_melted) = c("time", "variable", "run", "value", "group")
Diagnosed_Off_ART_All = data.frame(time, t(do.call(rbind, lapply(lapply(lapply(res_best_runs, function(x) {x$I22 + x$I23 + x$I24 + x$I25}), function(x) {return(x[,c(1, 2, 3, 4, 5, 6, 7, 8)])}), rowSums))))
Diagnosed_On_ART_All = data.frame(time, t(do.call(rbind, lapply(lapply(lapply(res_best_runs, function(x) {x$I32 + x$I33 + x$I34 + x$I35}), function(x) {return(x[,c(1, 2, 3, 4, 5, 6, 7, 8)])}), rowSums))))
Diagnosed_Dropout_ART_All = data.frame(time, t(do.call(rbind, lapply(lapply(lapply(res_best_runs, function(x) {x$I42 + x$I43 + x$I44 + x$I45}), function(x) {return(x[,c(1, 2, 3, 4, 5, 6, 7, 8)])}), rowSums))))
Diagnosed_All = data.frame(rbind(Diagnosed_Off_ART_All, Diagnosed_On_ART_All, Diagnosed_Dropout_ART_All),
rep(c("Diagnosed Off ART", "Diagnosed On ART", "Dropout"), each = length(time)))
colnames(Diagnosed_All) = c("time", as.character(seq(1, (length(Diagnosed_Off_ART_All[1,])-1))), "variable")
Diagnosed_All_melted = reshape2::melt(Diagnosed_All, id.vars = c("time", "variable"))
Diagnosed_All_melted$group = "All"
colnames(Diagnosed_All_melted) = c("time", "variable", "run", "value", "group")
Diagnosed = rbind(Diagnosed_FSW_melted, Diagnosed_All_melted)
Diagnosed_Off_ART_Men = data.frame(time, t(do.call(rbind, lapply(lapply(lapply(res_best_runs, function(x) {x$I22 + x$I23 + x$I24 + x$I25}), function(x) {return(x[,c(5, 6)])}), rowSums))))
Diagnosed_On_ART_Men = data.frame(time, t(do.call(rbind, lapply(lapply(lapply(res_best_runs, function(x) {x$I32 + x$I33 + x$I34 + x$I35}), function(x) {return(x[,c(5, 6)])}), rowSums))))
Diagnosed_Dropout_ART_Men = data.frame(time, t(do.call(rbind, lapply(lapply(lapply(res_best_runs, function(x) {x$I42 + x$I43 + x$I44 + x$I45}), function(x) {return(x[,c(5, 6)])}), rowSums))))
Diagnosed_Men = data.frame(rbind(Diagnosed_Off_ART_Men, Diagnosed_On_ART_Men, Diagnosed_Dropout_ART_Men),
rep(c("Diagnosed Off ART", "Diagnosed On ART", "Dropout"), each = length(time)))
colnames(Diagnosed_Men) = c("time", as.character(seq(1, (length(Diagnosed_Off_ART_Men[1,])-1))), "variable")
Diagnosed_Men_melted = reshape2::melt(Diagnosed_Men, id.vars = c("time", "variable"))
Diagnosed_Men_melted$group = "Men"
colnames(Diagnosed_Men_melted) = c("time", "variable", "run", "value", "group")
Diagnosed_Off_ART_Women = data.frame(time, t(do.call(rbind, lapply(lapply(lapply(res_best_runs, function(x) {x$I22 + x$I23 + x$I24 + x$I25}), function(x) {return(x[,c(1, 2, 3, 4)])}), rowSums))))
Diagnosed_On_ART_Women = data.frame(time, t(do.call(rbind, lapply(lapply(lapply(res_best_runs, function(x) {x$I32 + x$I33 + x$I34 + x$I35}), function(x) {return(x[,c(1, 2, 3, 4)])}), rowSums))))
Diagnosed_Dropout_ART_Women = data.frame(time, t(do.call(rbind, lapply(lapply(lapply(res_best_runs, function(x) {x$I42 + x$I43 + x$I44 + x$I45}), function(x) {return(x[,c(1, 2, 3, 4)])}), rowSums))))
Diagnosed_Women = data.frame(rbind(Diagnosed_Off_ART_Women, Diagnosed_On_ART_Women, Diagnosed_Dropout_ART_Women),
rep(c("Diagnosed Off ART", "Diagnosed On ART", "Dropout"), each = length(time)))
colnames(Diagnosed_Women) = c("time", as.character(seq(1, (length(Diagnosed_Off_ART_Women[1,])-1))), "variable")
Diagnosed_Women_melted = reshape2::melt(Diagnosed_Women, id.vars = c("time", "variable"))
Diagnosed_Women_melted$group = "Women"
colnames(Diagnosed_Women_melted) = c("time", "variable", "run", "value", "group")
Diagnosed_Women_Men = rbind(Diagnosed_Women_melted, Diagnosed_Men_melted)
Diagnosed_Women_Men_ratio = data.frame(Diagnosed_Women_melted[,c("time", "variable", "run")],
value = Diagnosed_Women_melted$value/Diagnosed_Men_melted$value)
HIV_deaths = data.frame(time[-1], t(do.call(rbind, lapply(lapply(lapply(res_best_runs, function(x) x$cumuHIVDeaths), rowSums), diff) )))
colnames(HIV_deaths) = c("time", as.character(seq(1, (length(HIV_deaths[1,])-1))))
HIV_deaths_melted = reshape2::melt(HIV_deaths, id.vars = c("time"))
pc_OfWomen_ProFSW = data.frame(time, t(100*do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,1]))/
(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,1])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,2])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,3])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,4]))+
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,7]))
)))
pc_OfWomen_LowFSW = data.frame(time, t(100*do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,2]))/
(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,1])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,2])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,3])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,4]))+
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,7])))))
pc_OfWomen_Active_FSW = data.frame(time, t(100*(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,1])) + do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,2])))/
(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,1])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,2])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,3])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,4]))+
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,7])))))
pc_OfWomen_GPF = data.frame(time, t(100*do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,3]))/
(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,1])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,2])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,3])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,4]))+
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,7])))))
pc_OfWomen_FormerFSW = data.frame(time, t(100*do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,4]))/
(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,1])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,2])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,3])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,4]))+
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,7])))))
pc_OfMen_Client = data.frame(time, t(100*do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,5]))/
(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,5])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,6]))+
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,8])))))
pc_OfMen_GPM = data.frame(time, t(100*do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,6]))/
(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,5])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,6]))+
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,8])))))
pc_OfWomen_VF = data.frame(time, t(100*do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,7]))/
(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,1])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,2])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,3])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,4]))+
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,7])))))
pc_OfMen_VM = data.frame(time, t(100*do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,8]))/
(do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,5])) +
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,6]))+
do.call(rbind, lapply(lapply(res_best_runs, function(x) x$N), function(x) x[,8])))))
Ratio_Low_Pro = data.frame(time, t(do.call(rbind, lapply(res_best_runs, function(x) {x$frac_N[,2]/ x$frac_N[,1]}))))
frac_by_gender = rbind(pc_OfWomen_ProFSW,
pc_OfWomen_LowFSW,
pc_OfWomen_GPF,
pc_OfWomen_FormerFSW,
pc_OfMen_Client,
pc_OfMen_GPM,
pc_OfWomen_VF,
pc_OfMen_VM,
pc_OfWomen_Active_FSW,
Ratio_Low_Pro
)
frac_by_gender = data.frame(frac_by_gender, group = rep(c("Pro FSW", "Low-level FSW", "GPF", "Former FSW in Cotonou", "Clients", "GPM", "Virgin female", "Virgin male", "Active FSW", "Low Pro Ratio"), each = length(time)))
colnames(frac_by_gender) = c("time",as.character(seq(1, (length(pc_OfWomen_ProFSW[1,])-1))), "variable")
frac_by_gender$variable = factor(frac_by_gender$variable, levels = c("Pro FSW", "Low-level FSW", "GPF", "Former FSW in Cotonou", "Clients", "GPM", "Virgin female", "Virgin male", "Active FSW", "Low Pro Ratio"))
frac_by_gender_melted = reshape2::melt(frac_by_gender, id.vars = c("time", "variable"))
colnames(frac_by_gender_melted) = c("time", "variable", "run", "value")
colnames(pc_OfWomen_ProFSW) = c("time", as.character(seq(1, (length(pc_OfWomen_ProFSW[1,])-1))))
colnames(pc_OfWomen_LowFSW) = c("time", as.character(seq(1, (length(pc_OfWomen_LowFSW[1,])-1))))
colnames(pc_OfWomen_Active_FSW) = c("time", as.character(seq(1, (length(pc_OfWomen_Active_FSW[1,])-1))))
colnames(pc_OfMen_Client) = c("time", as.character(seq(1, (length(pc_OfMen_Client[1,])-1))))
colnames(pc_OfWomen_VF) = c("time", as.character(seq(1, (length(pc_OfWomen_VF[1,])-1))))
colnames(pc_OfMen_VM) = c("time", as.character(seq(1, (length(pc_OfMen_VM[1,])-1))))
pc_OfWomen_ProFSW_melted = reshape2::melt(pc_OfWomen_ProFSW, id.vars = c("time"))
pc_OfWomen_LowFSW_melted = reshape2::melt(pc_OfWomen_LowFSW, id.vars = c("time"))
pc_OfWomen_Active_FSW_melted = reshape2::melt(pc_OfWomen_Active_FSW, id.vars = c("time"))
pc_OfMen_Client_melted = reshape2::melt(pc_OfMen_Client, id.vars = c("time"))
pc_OfWomen_VF_melted = reshape2::melt(pc_OfWomen_VF, id.vars = c("time"))
pc_OfMen_VM_melted = reshape2::melt(pc_OfMen_VM, id.vars = c("time"))
condom_Pro_FSW_comm = t((do.call(rbind, lapply(res_best_runs, function(x) {x$fc_comm[,1,][,5]}))))
condom_Pro_FSW_noncomm = t((do.call(rbind, lapply(res_best_runs, function(x) {x$fc_noncomm[,1,][,5]}))))
condom_Pro_FSW = data.frame(time, rbind(condom_Pro_FSW_comm, condom_Pro_FSW_noncomm), rep(c("Commercial", "Non commercial"), each = length(time)))
colnames(condom_Pro_FSW) = c("time", as.character(seq(1, (length(condom_Pro_FSW[1,])-2))), "variable")
condom_Pro_FSW_melted = reshape2::melt(condom_Pro_FSW, id.vars = c("time", "variable"))
colnames(condom_Pro_FSW_melted) = c("time", "variable", "run", "value")
condom_Low_FSW_comm = t((do.call(rbind, lapply(res_best_runs, function(x) {x$fc_comm[,2,][,5]}))))
condom_Low_FSW_noncomm = t((do.call(rbind, lapply(res_best_runs, function(x) {x$fc_noncomm[,2,][,5]}))))
condom_Low_FSW = data.frame(time, rbind(condom_Low_FSW_comm, condom_Low_FSW_noncomm), rep(c("Commercial", "Non commercial"), each = length(time)))
colnames(condom_Low_FSW) = c("time", as.character(seq(1, (length(condom_Low_FSW[1,])-2))), "variable")
condom_Low_FSW_melted = reshape2::melt(condom_Low_FSW, id.vars = c("time", "variable"))
colnames(condom_Low_FSW_melted) = c("time", "variable", "run", "value")
condom_GPF_noncomm = data.frame(time, t((do.call(rbind, lapply(res_best_runs, function(x) {x$fc_noncomm[,3,][,6]})))))
colnames(condom_GPF_noncomm) = c("time", as.character(seq(1, (length(condom_GPF_noncomm[1,])-1))))
condom_GPF_noncomm_melted = reshape2::melt(condom_GPF_noncomm, id.vars = c("time"))
condom_GPM_noncomm = data.frame(time, t((do.call(rbind, lapply(res_best_runs, function(x) {x$fc_noncomm[,6,][,3]})))))
colnames(condom_GPM_noncomm) = c("time", as.character(seq(1, (length(condom_GPM_noncomm[1,])-1))))
condom_GPM_noncomm_melted = reshape2::melt(condom_GPM_noncomm, id.vars = c("time"))
testing_rate_ratio_F_M = data.frame(time, t((do.call(rbind, lapply(lapply(res_best_runs, function(x) {x$testing_prob}), function(x) {return(x[,c(3)])}))/do.call(rbind, lapply(lapply(res_best_runs, function(x) {x$testing_prob}), function(x) {return(x[,c(6)])})))))
colnames(testing_rate_ratio_F_M) = c("time", as.character(seq(1, (length(testing_rate_ratio_F_M[1,])-1))))
testing_rate_ratio_F_M_melted = reshape2::melt(testing_rate_ratio_F_M, id.vars = c("time"))
colnames(testing_rate_ratio_F_M_melted) = c("time", "run", "value")
Diagnosed_Women_Men_ratio_On_ART = Diagnosed_Women_Men_ratio[Diagnosed_Women_Men_ratio$variable == "Diagnosed On ART",]
testing_rate_ratio_F_M_Women_Men_ratio_On_ART = data.frame(x=testing_rate_ratio_F_M_melted$value, y=Diagnosed_Women_Men_ratio_On_ART$value)
pfFSW = data.frame(time, t(do.call(rbind, lapply(res_best_runs, function(x) {x$pfFSW[,1]}))))
colnames(pfFSW) = c("time", as.character(seq(1, (length(pfFSW[1,])-1))))
pfFSW_melted = reshape2::melt(pfFSW, id.vars = c("time"))
colnames(pfFSW_melted) = c("time", "run", "value")
ART_inits = data.frame(time = time[-length(time)], t(do.call(rbind,lapply(res_best_runs, function(x) {diff(rowSums(x$cumuARTinitiations))}))))
colnames(ART_inits) = c("time", as.character(seq(1, (length(ART_inits[1,])-1))))
ART_inits_melted = reshape2::melt(ART_inits, id.vars = c("time"))
colnames(ART_inits_melted) = c("time", "run", "value")
# Ratio female:male of rate of initiation per infected = 8.0/4.68 = 1.7
# to compare to the data, I am taking the number of initiations that years as a fraction of the total nubmer of PLHIV that year
# DENOMINTOR = ALL POSTIIVES
female_ART_init_rate_from_all_HIV_pos = data.frame(time = time[-length(time)], t(do.call(rbind,lapply(res_best_runs, function(x) {
a = diff(rowSums(x$cumuARTinitiations[,c(1,2,3,4)])) # number of female ART initiations per year
b = rowSums(x$HIV_positive[,c(1,2,3,4)])[-length(time)] # number of HIV positive females at the beginning of the year
return(a/b)
}))))
male_ART_init_rate_from_all_HIV_pos = data.frame(time = time[-length(time)], t(do.call(rbind,lapply(res_best_runs, function(x) {
c = diff(rowSums(x$cumuARTinitiations[,c(5,6)])) # number of male ART initiations per year
d = rowSums(x$HIV_positive[,c(5,6)])[-length(time)] # number of HIV positive males at the beginning of the year
return(c/d)
}))))
# DENOMINTOR = ALL POSTIIVES NOT ON ART
female_ART_init_rate_from_all_HIV_pos_NOT_ON_ART = data.frame(time = time[-length(time)], t(do.call(rbind,lapply(res_best_runs, function(x) {
a = diff(rowSums(x$cumuARTinitiations[,c(1,2,3,4)])) # number of female ART initiations per year
b = rowSums(x$HIV_positive[,c(1,2,3,4)] - x$HIV_positive_On_ART[,c(1,2,3,4)])[-length(time)] # number of HIV positive females not on ART at the beginning of the year
return(a/b)
}))))
# female_ART_init_rate_from_all_HIV_pos_NOT_ON_ART
male_ART_init_rate_from_all_HIV_pos_NOT_ON_ART = data.frame(time = time[-length(time)], t(do.call(rbind,lapply(res_best_runs, function(x) {
a = diff(rowSums(x$cumuARTinitiations[,c(5,6)])) # number of female ART initiations per year
b = rowSums(x$HIV_positive[,c(5,6)] - x$HIV_positive_On_ART[,c(5,6)])[-length(time)] # number of HIV positive females not on ART at the beginning of the year
return(a/b)
}))))
# male_ART_init_rate_from_all_HIV_pos_NOT_ON_ART
# FEMALE TO MALE RATIO OF INITIATION RATES WHERE DENOMINATOR IS ART POSITIVE NOT ON ART
FtM_ratio_ART_initiation_rates_from_HIVpos_not_on_ART = data.frame(time = time[-length(time)], t(do.call(rbind, lapply(res_best_runs, function(x) {
a = diff(rowSums(x$cumuARTinitiations[,c(1,2,3,4)])) # number of female ART initiations per year
b = rowSums(x$HIV_positive[,c(1,2,3,4)] - x$HIV_positive_On_ART[,c(1,2,3,4)])[-length(time)] # number of HIV positive females not on ART at the beginning of the year
c = diff(rowSums(x$cumuARTinitiations[,c(5,6)])) # number of female ART initiations per year
d = rowSums(x$HIV_positive[,c(5,6)] - x$HIV_positive_On_ART[,c(5,6)])[-length(time)] # number of HIV positive females not on ART at the beginning of the year
return((a/b)/(c/d))
}))))
colnames(FtM_ratio_ART_initiation_rates_from_HIVpos_not_on_ART) = c("time", as.character(seq(1, (length(FtM_ratio_ART_initiation_rates_from_HIVpos_not_on_ART[1,])-1))))
FtM_ratio_ART_initiation_rates_from_HIVpos_not_on_ART_melted = reshape2::melt(FtM_ratio_ART_initiation_rates_from_HIVpos_not_on_ART, id.vars = c("time"))
ART_inits_ratio_F_over_M = data.frame(time = time[-length(time)], t(do.call(rbind,lapply(res_best_runs, function(x) {diff(rowSums(x$cumuARTinitiations[,c(1,2,3,4)])) /
diff(rowSums(x$cumuARTinitiations[,c(5,6)]))}))))
ART_inits_ratio_F_over_M_melted = reshape2::melt(ART_inits_ratio_F_over_M, id.vars = c("time"))
ART_inits_cumu = data.frame(time = time, t(do.call(rbind,lapply(res_best_runs, function(x) {rowSums(x$cumuARTinitiations)}))))
colnames(ART_inits_cumu) = c("time", as.character(seq(1, (length(ART_inits_cumu[1,])-1))))
ART_inits_cumu_melted = reshape2::melt(ART_inits_cumu, id.vars = c("time"))
colnames(ART_inits_cumu_melted) = c("time", "run", "value")
ART_REinits_cumu = data.frame(time = time, t(do.call(rbind,lapply(res_best_runs, function(x) {rowSums(x$cumuARTREinitiations)}))))
colnames(ART_REinits_cumu) = c("time", as.character(seq(1, (length(ART_REinits_cumu[1,])-1))))
ART_REinits_cumu_melted = reshape2::melt(ART_REinits_cumu, id.vars = c("time"))
colnames(ART_REinits_cumu_melted) = c("time", "run", "value")
ART_REinits = data.frame(time = time[-length(time)], t(do.call(rbind,lapply(res_best_runs, function(x) {diff(rowSums(x$cumuARTREinitiations))}))))
colnames(ART_REinits) = c("time", as.character(seq(1, (length(ART_REinits[1,])-1))))
ART_REinits_melted = reshape2::melt(ART_REinits, id.vars = c("time"))
colnames(ART_REinits_melted) = c("time", "run", "value")
frac_ART_inits_REinits = data.frame(time = time[-length(time)])
for(i in 2:(length(res_best_runs)+1))
{
frac_ART_inits_REinits[,i] = ART_REinits[,i]/(ART_inits[,i] + ART_REinits[,i]) * 100
}
colnames(frac_ART_inits_REinits) = c("time", as.character(seq(1, (length(frac_ART_inits_REinits[1,])-1))))
frac_ART_inits_REinits_melted = reshape2::melt(frac_ART_inits_REinits, id.vars = c("time"))
colnames(frac_ART_inits_REinits_melted) = c("time", "run", "value")
ART_init_rate_from_all_HIV_pos_NOT_ON_ART_by_sex = data.frame(
Women_2014 = unlist(female_ART_init_rate_from_all_HIV_pos_NOT_ON_ART[which(time == 2014),c(2:(length(res_best_runs)+1))]),
Women_2015 = unlist(female_ART_init_rate_from_all_HIV_pos_NOT_ON_ART[which(time == 2015),c(2:(length(res_best_runs)+1))]),
Men_2014 = unlist(male_ART_init_rate_from_all_HIV_pos_NOT_ON_ART[which(time == 2014),c(2:(length(res_best_runs)+1))]),
Men_2015 = unlist(male_ART_init_rate_from_all_HIV_pos_NOT_ON_ART[which(time == 2015),c(2:(length(res_best_runs)+1))])
)
ART_init_rate_from_all_HIV_pos_NOT_ON_ART_by_sex_melted = reshape2::melt(ART_init_rate_from_all_HIV_pos_NOT_ON_ART_by_sex)
levels(ART_init_rate_from_all_HIV_pos_NOT_ON_ART_by_sex_melted) = c("Women_2014","Women_2015", "Men_2014","Men_2015")
###################
prev_indiv_ribbon = data.frame(time, rep(c("Pro FSW", "Low-level FSW", "Clients", "Women", "Men"), each = length(time)), rbind(t(apply(prev_FSW_indiv, 1, function(x) quantile(x, c(0.05, 0.95)))),
t(apply(prev_LowFSW_indiv, 1, function(x) quantile(x, c(0.05, 0.95)))),
t(apply(prev_client_indiv, 1, function(x) quantile(x, c(0.05, 0.95)))),
t(apply(prev_women_indiv, 1, function(x) quantile(x, c(0.05, 0.95)))),
t(apply(prev_men_indiv, 1, function(x) quantile(x, c(0.05, 0.95))))))
colnames(prev_indiv_ribbon) = c("time", "variable", "Lower", "Upper")
# prev_indiv_ribbon_melted = reshape2::melt(prev_indiv_ribbon,id.vars = c("time", "variable"))
N_ART_FSW_indiv_ribbon = data.frame(time, "Pro FSW", t(apply(N_ART_FSW_indiv[, -c(1, 2)], 1, function(x) quantile(x, c(0.05, 0.95)))))
colnames(N_ART_FSW_indiv_ribbon) = c("time", "variable", "Lower", "Upper")
ART_indiv_ribbon = data.frame(time,
rep(c("Pro FSW", "All"), each = length(time)),
rbind(t(apply(ART_FSW_indiv, 1, function(x) quantile(x, c(0.05, 0.95)))),
t(apply(ART_all_indiv, 1, function(x) quantile(x, c(0.05, 0.95))))
))
colnames(ART_indiv_ribbon) = c("time", "variable", "Lower", "Upper")
require(ggplot2)
# want plots on number of initiations of ART, dropout, reinitiations! compare to data points from the reports!
# need to do difference between years for cumulative art inits etc
prev_points = data.frame(time = c(1986, 1987, 1988, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015,
1998, 2002, 2005, 2008, 2012, 2015,
1998, 2008, 2011,
1998, 2008, 2011,
2012, 2015),
variable = c(rep("Pro FSW", 11),
rep("Clients", 6),
rep("Women", 3),
rep("Men", 3),
rep("Low-level FSW", 2)),
value = c(3.3, 8.2, 19.2, 53.3, 48.7, 40.6, 38.9, 34.8, 29.3, 27.4, 18.7,
100*0.084, 9, 6.9, 5.8, 100*0.028, 100*0.016,
100*0.035, 100*0.04, 2.2,
100*0.033, 100*0.02, 1.6,
100*0.084, 100*0.043),
lower = c(3.3, 8.2, 19.2, 48.02, 43.02, 36.58, 31.97, 30.42, 24.93, 23.01, 15.71,
100*0.05898524, 100*0.068218538, 100*0.04293149, 100*0.034772131, 100*0.012660836, 100*0.006039259,
100*0.024181624, 100*0.030073668, 100*0.012980254,
100*0.022857312, 100*0.012427931, 100*0.007517563,
# 100*0.091838441, 100*0.026704897),
100*0.055700329, 100*0.024043597),
upper = c(3.3, 8.2, 19.2, 58.48, 54.42, 44.67, 46.27, 39.38, 33.88, 32.23, 22.01,
100*0.11561791, 100*0.115608811, 100*0.105215792, 100*0.090216628, 100*0.051602442, 100*0.035338436,
100*0.047726245, 100*0.052817187, 100*0.035296286,
100*0.047183668, 100*0.029774338, 100*0.028546718,
100*0.120857355, 100*0.069311506))
prev_points_all = prev_points
prev_points = prev_points[-c(1,2,3),]
prev_axes = data.frame(variable = c(rep("Pro FSW", 2),
rep("Clients", 2),
rep("Women", 2),
rep("Men", 2),
rep("Low-level FSW", 2)),
time = c(rep(c(1986, 2025), 5)),
value = c(0, 70, 0, 70, 0, 15, 0, 15, 0, 70)
)
prev_points_80s = prev_points_all[c(1,2,3),]
prev_points[prev_points$time == "2015", "lower"][1] = 13.79
# plot prevalence in each group indiv runs
# g1=ggplot() + geom_line(data = prev_indiv_melted, aes(x = time, y = value, factor = variable, factor = run), alpha = 0.3) + theme_bw() + facet_wrap(~variable, scales = "free") + labs(y = "prevalence (%)") +
# geom_point(data = prev_points, aes(x = time, y = value))+ geom_errorbar(data = prev_points, aes(x = time, ymin = lower, ymax = upper))+
# geom_point(data = prev_points_80s, aes(x = time, y = value), colour = "red")+
# geom_blank(data = prev_axes, aes(x = time, y = value))+
# theme(text = element_text(size=20))
#
# g1
# levels(prev_indiv_ribbon$variable) <- factor(prev_indiv_ribbon$variable, levels(prev_indiv_ribbon$variable)[c(4,2,1,5,3)])
# c("Pro FSW", "Low-level FSW", "Clients", "Men", "Women")
g1=ggplot() +
# geom_line(data = prev_indiv_melted, aes(x = time, y = value, factor = variable, factor = run), alpha = 0.3) +
theme_bw() + facet_wrap(~variable, scales = "free") + labs(y = "prevalence (%)") +
geom_point(data = prev_points, aes(x = time, y = value))+
geom_errorbar(data = prev_points, aes(x = time, ymin = lower, ymax = upper))+
geom_point(data = prev_points_80s, aes(x = time, y = value), colour = "red")+
geom_blank(data = prev_axes, aes(x = time, y = value))+
theme(text = element_text(size=20)) +
geom_ribbon(data = prev_indiv_ribbon, aes(x = time, ymin = Lower, ymax = Upper), alpha = 0.3)
# geom_polygon(data = prev_indiv_melted, aes(x = time, y = value, factor = variable, factor = run))
g1
ART_data_points_with_numbers_FSW = data.frame(time = c(2012, 2014, 2015, 2016, 2017),
Lower = c(27, 34, 42, 53, 73),
Upper = c(37, 46, 56, 83, 107),
variable = rep("Pro FSW", 5),
datatype = c("data", "data", "data", "counterfactual", "counterfactual"))
ART_data_points_with_numbers_FSW_points = data.frame(time = c(2012, 2014, 2015, 2016, 2017),
value = c(32, 40, 49, 53, 73),
variable = rep("Pro FSW", 5),
datatype = c("data", "data", "data", "counterfactual", "counterfactual"))
g7=ggplot() +
# geom_line(data = N_ART_FSW_indiv_melted, aes(x = time, y = value, factor = variable, factor = run), alpha = 1) +
ggtitle("Number of professional FSW on ART") +theme_bw() +
geom_errorbar(data = ART_data_points_with_numbers_FSW, aes(x = time, ymin = Lower, ymax = Upper, col = datatype)) + labs(x = "Year") +
scale_x_continuous(breaks = seq(2000, 2018, 1), limits = c(2000,2018))+
labs(y = "")+
geom_point(data = ART_data_points_with_numbers_FSW_points, aes(x = time, y = value, col = datatype))+
theme(text = element_text(size=20),plot.title = element_text(hjust = 0.5))+
theme(text = element_text(size=24),
legend.text=element_text(size=18),
legend.key.size = unit(1.3, 'lines'))+ annotate("text", x = 2008.8, y = 70, size = 7,col = "black", label = "Data on numbers of FSW file active \n from Luc +- 15% \n Fitting to 2015 estimate only")+
geom_ribbon(data = N_ART_FSW_indiv_ribbon, aes(x = time, ymin = Lower, ymax = Upper), alpha = 0.3)
g7
frac_by_gender_discard_points_no_FSW_LB = data.frame(variable = c("Pro FSW", "Clients", "Virgin female", "Virgin male", "Active FSW", "Low Pro Ratio"),
min = c(0, 0.074, 0.07896475, 0.07039551, 0.0048, .01),
max = c(0.0143/2, 0.3 , 0.2, 0.17, 0.0143, .05))
# plot fraction in each group
g2=ggplot(frac_by_gender_melted) + geom_line(aes(x = time, y = value, factor = run)) +
theme_bw() + labs(y = "Percent in each group in their respective gender (%)") + facet_wrap(~variable, scales = "free") +
# geom_point(data = frac_N_data_points, aes(x = time, y = point), size = I(2), color = "red", shape = 15) +
geom_hline(data = frac_by_gender_discard_points_no_FSW_LB, aes(yintercept = 100*min), size = I(0.5), color = "red", linetype = 1, alpha = 0.7) +
geom_hline(data = frac_by_gender_discard_points_no_FSW_LB, aes(yintercept = 100*max), size = I(0.5), color = "red", linetype = 1, alpha = 0.7)+
theme(text = element_text(size=20))
g2
g3=ggplot() + geom_line(data = pc_OfWomen_ProFSW_melted, aes(x = time, y = value, factor = variable), alpha = 1) + theme_bw() + labs(y = "% of all women 15-59 that are pro FSW")+
theme(text = element_text(size=20)) + geom_hline(yintercept = 0.24,col = "orange", linetype="dashed") + geom_hline(yintercept = 0.72, linetype="dashed",col = "red")+
theme(text = element_text(size=24),
legend.text=element_text(size=18),
legend.key.size = unit(1.3, 'lines'))+ annotate("text", x = 2002.8, y = 0.7, size = 7,col = "red", label = "Fitted to below this line")+
theme(text = element_text(size=24),
legend.text=element_text(size=18),
legend.key.size = unit(1.3, 'lines'))+ annotate("text", x = 2003.8, y = 0.3, size = 7,col = "orange", label = "Original estimate lower bound ")
g3
# g3b=ggplot() + geom_line(data = pc_OfWomen_ProFSW_melted, aes(x = time, y = value, factor = variable), alpha = 1) + theme_bw() + labs(y = "% of all women 15-59 that are pro FSW")+
# theme(text = element_text(size=20)) + geom_hline(yintercept = 0.24,col = "orange", linetype="dashed") + geom_hline(yintercept = 0.72, linetype="dashed",col = "red")+
# theme(text = element_text(size=24),
# legend.text=element_text(size=18),
# legend.key.size = unit(1.3, 'lines'))+ annotate("text", x = 2002.8, y = 0.7, size = 7,col = "red", label = "Fitted to below this line")+
# theme(text = element_text(size=24),
# legend.text=element_text(size=18),
# legend.key.size = unit(1.3, 'lines'))+ annotate("text", x = 2003.8, y = 0.3, size = 7,col = "orange", label = "Original estimate lower bound ")
# g3b
g4=ggplot() + geom_line(data = pc_OfWomen_LowFSW_melted, aes(x = time, y = value, factor = variable), alpha = 0.3) + theme_bw() + labs(y = "% of all women 15-59 that are low level FSW")+
theme(text = element_text(size=20))
g4b = ggplot() + geom_line(data = pc_OfWomen_VF_melted, aes(x = time, y = value, factor = variable), alpha = 0.3) + theme_bw() + labs(y = "% of all women 15-59 that are virgins")+
theme(text = element_text(size=20))+ geom_hline(yintercept = c(7.9, 20), col = "red") +
theme(text = element_text(size=24), legend.text=element_text(size=18),
legend.key.size = unit(1.3, 'lines'))+ annotate("text", x = 2010, y = 19, size = 7,col = "red", label = "Fitting")
g4b
g5=ggplot() + geom_line(data = pc_OfMen_Client_melted, aes(x = time, y = value, factor = variable), alpha = 1) + theme_bw() + labs(y = "% of all men 15-59 that are clients")+
theme(text = element_text(size=20)) + geom_hline(yintercept = c(7.4, 30), col = "red") +
theme(text = element_text(size=24), legend.text=element_text(size=18),
legend.key.size = unit(1.3, 'lines'))+ annotate("text", x = 2003.8, y = 28, size = 7,col = "red", label = "Fitting")
g5
g5b = ggplot() + geom_line(data = pc_OfMen_VM_melted, aes(x = time, y = value, factor = variable), alpha = 0.3) + theme_bw() + labs(y = "% of all men 15-59 that are virgins")+
theme(text = element_text(size=20))+ geom_hline(yintercept = c(7, 17), col = "red")+
theme(text = element_text(size=24), legend.text=element_text(size=18),
legend.key.size = unit(1.3, 'lines'))+ annotate("text", x = 2010, y = 15, size = 7,col = "red", label = "Fitting")
g5b
# # plot prevalence in each group
# ggplot() + geom_line(data = prev, aes(x = time, y = Median))+ geom_ribbon(data = prev, aes(x = time, ymin = Lower, ymax = Upper), alpha = 0.5) + theme_bw() + facet_wrap(~variable, scales = "free") + labs(y = "prevalence (%)") +
# geom_point(data = prev_points, aes(x = time, y = value))+ geom_errorbar(data = prev_points, aes(x = time, ymin = lower, ymax = upper)) +
# geom_point(data = prev_points_80s, aes(x = time, y = value), colour = "red")+
# geom_blank(data = prev_axes, aes(x = time, y = value))
#
# plot total population size
g6=ggplot(Ntot) + geom_line(aes(x = time, y = Median)) + geom_ribbon(aes(x = time, ymin = Lower, ymax = Upper), alpha = 0.5) +
theme_bw() + labs(y = "Total population size of Grand Cotonou") +
geom_point(data = Ntot_data_points, aes(x = time, y = point, color = colour), size = I(2), shape = 15) + geom_errorbar(data = Ntot_data_points, aes(x = time, ymax = upper, ymin = lower, color = colour), width = 2)+
theme(text = element_text(size=20))
# ggplot(ART_coverage) +
# geom_line(aes(x = time, y = Median))+ geom_ribbon(aes(x = time, ymin = Lower, ymax = Upper), alpha = 0.5) + theme_bw() +
# facet_wrap(~variable) + labs(y = "ART coverage ") +
# geom_errorbar(data = ART_data_points, aes(x = time, ymin = Lower, ymax = Upper), colour = "darkred")
ART_data_points = data.frame(time = c(2011, 2012, 2013, 2014, 2015, 2016, 2017,
2012, 2014, 2015, 2016, 2017
),
Lower = c(0.33, 0.42, 0.44, 0.39, 0.44, 0.51, 0.57,
0.09, 0.14, 0.2, 0.22, 0.3
),
Upper = c(0.52, 0.66, 0.69, 0.61, 0.69, 0.8, 0.8,
0.13, 0.2, 0.28, 0.3, 0.42
),
variable = c("All", "All", "All", "All", "All", "All", "All",
"Pro FSW", "Pro FSW", "Pro FSW", "Pro FSW", "Pro FSW"))
ART_data_points$Data = c(rep("UNAIDS (fitted 2011 & 2017 only)", 7), rep("File active FSW (X-validation)", 5))
art_cov = data.frame(time = 2017, variable = "All", Upper = 0.91, Lower = 0.6, Data = "Audit report 2017 (x-validation)")
ART_data_points = rbind(ART_data_points, art_cov)
ART_text = data.frame(time = 2010, variable = "Pro FSW", value = 50)
g8=ggplot() +
# geom_line(data = ART_indiv_melted, aes(x = time, y = value*100, factor = variable, factor = run), alpha = 0.3) +
theme_bw() + facet_wrap(~variable, scales = "free") + labs(y = "ART coverage (%)") +
geom_errorbar(data = ART_data_points, aes(x = time, ymin = Lower*100, ymax = Upper*100, col = Data), size = I(1))+
theme(text = element_text(size=20)) + theme(legend.position = "top") + geom_text(data = ART_text, aes(x = time, y = value, factor = variable),
label = "2016 & 2017 are counterfactuals \n NOT including FSW on TasP/PrEP", size = 5)+
geom_ribbon(data = ART_indiv_ribbon, aes(x = time, ymin = Lower*100, ymax = Upper*100), alpha = 0.3)
g8
g9=ggplot() + geom_line(data = lambda_sum_0_indiv_melted[lambda_sum_0_indiv_melted$time > 2009,], aes(x = time, y = value*100, factor = variable, factor = run), alpha = 0.3) + theme_bw() + facet_wrap(~variable, scales = "free") + labs(y = "Force of infection on susceptibles off PrEP (%)")+
theme(text = element_text(size=20))
Diagnosed_FSW_melted$col = Diagnosed_FSW_melted$variable
Diagnosed_FSW_melted$col = ifelse(Diagnosed_FSW_melted$col == "Diagnosed Off ART", "#7fc97f",
ifelse(Diagnosed_FSW_melted$col == "Diagnosed On ART", "#beaed4",
ifelse(Diagnosed_FSW_melted$col == "Dropout", "#fdc086", 0)))
# 2012 FSW IBBA ever tested 0.75 - 0.83
# EQS mapping 2012 N = 889 - 1391
# 2012 prevalence = 27.4% [0.23 - 0.322]
# lower N diagnosed = 0.75 * 889 * 0.23 =
# upper N diagnosed = 0.83 * 1391 * 0.322 =
# bounds of infected ever tested = 304 * 0.75, 304 * 0.83 = 228 - 252
pc_diagnosed = data.frame(time = 2012, Lower = 153, Upper = 372, variable = "Pro FSW")
g10=ggplot() + geom_line(data = Diagnosed_FSW_melted, aes(x = time, y = value, factor = variable, factor = run, colour = variable), alpha = 1) + theme_bw() + labs(y = "Numbers of pro FSW in Grand Cotonou")+
scale_x_continuous(breaks = seq(min(time), max(time), 4) ) +
scale_y_continuous(breaks = seq(min(Diagnosed$value), #max(Diagnosed_FSW_melted$value), 20) )+
380, 40) )+
theme(text = element_text(size=20), legend.position = "top")+
geom_errorbar(data = ART_data_points_with_numbers_FSW, size=1,aes(x = time, ymin = Lower, ymax = Upper), col ="#4daf4a") +
labs(x = "Year")+
geom_errorbar(data = pc_diagnosed, size=1, aes(x = time, ymin = Lower, ymax = Upper), col ="#e41a1c") +
labs(x = "Year")+ scale_colour_brewer(type = "qual", palette = 6) +
theme(text = element_text(size=24),
legend.text=element_text(size=18))+ annotate("text", x = 2003.8, y = 222, size = 7,col = "#e41a1c", label = "Number of FSW * \n prevalence * \n % ever tested \n (cross-validation)")
g10
# #
# ggplot() + geom_line(data = Diagnosed_All_melted, aes(x = time, y = value, factor = variable, factor = run, colour = variable), alpha = 1) + theme_bw() + labs(y = "Numbers of all people in Grand Cotonou")+
# scale_x_continuous(breaks = seq(min(time), max(time)+1, 4) ) +
# scale_y_continuous(breaks = seq(min(Diagnosed$value), max(Diagnosed$value), 1000) ) +
# geom_vline(xintercept=2005, col = "black", alpha = 0.5, linetype = "dashed") +
# geom_text(aes(x=2005, label="CD4 < 200", y=-400))+
# geom_vline(xintercept=2012, col = "black", alpha = 0.5, linetype = "dashed") +
# geom_text(aes(x=2012, label="CD4 < 350", y=-800))+
# geom_vline(xintercept=2015, col = "black", alpha = 0.5, linetype = "dashed") +
# geom_text(aes(x=2015, label="CD4 < 500", y=-400))+
# geom_vline(xintercept=2016, col = "black", alpha = 0.5, linetype = "dashed") +
# geom_text(aes(x=2016, label="Everyone", y=-800))
# #
# I2x= Diagnosed_All_melted[Diagnosed_All_melted$variable == "Diagnosed Off ART",]
# I3x= Diagnosed_All_melted[Diagnosed_All_melted$variable == "Diagnosed On ART",]
# I4x= Diagnosed_All_melted[Diagnosed_All_melted$variable == "Dropout",]
I2xF= Diagnosed_Women_melted[Diagnosed_Women_melted$variable == "Diagnosed Off ART",]
I3xF= Diagnosed_Women_melted[Diagnosed_Women_melted$variable == "Diagnosed On ART",]
I4xF= Diagnosed_Women_melted[Diagnosed_Women_melted$variable == "Dropout",]
I2xM= Diagnosed_Men_melted[Diagnosed_Men_melted$variable == "Diagnosed Off ART",]
I3xM= Diagnosed_Men_melted[Diagnosed_Men_melted$variable == "Diagnosed On ART",]
I4xM= Diagnosed_Men_melted[Diagnosed_Men_melted$variable == "Dropout",]
diag_off_art_auditF = data.frame(time = 2017, Lower = 195, Upper = 209, variable = "Women")
diag_off_art_auditM = data.frame(time = 2017, Lower = 100, Upper = 109, variable = "Men")
g10a = ggplot() + geom_line(data = I2xF, aes(x = time, y = value, factor = variable, factor = run), colour = "black", alpha = 1) + theme_bw() +
labs(y = "")+
scale_x_continuous(breaks = seq(min(time), max(time)+1, 4) ) +
scale_y_continuous(breaks = seq(min(Diagnosed$value), max(Diagnosed$value), 1000) ) +
geom_vline(xintercept=2005, col = "black", alpha = 0.5, linetype = "dashed") +
geom_text(aes(x=2005, label="CD4 < 200", y=-400))+
geom_vline(xintercept=2012, col = "black", alpha = 0.5, linetype = "dashed") +
geom_text(aes(x=2012, label="CD4 < 350", y=-800))+
geom_vline(xintercept=2015, col = "black", alpha = 0.5, linetype = "dashed") +
geom_text(aes(x=2015, label="CD4 < 500", y=-400))+
geom_vline(xintercept=2016, col = "black", alpha = 0.5, linetype = "dashed") +
theme(text = element_text(size=20), legend.position = "none",plot.title = element_text(hjust = 0.5))+
geom_errorbar(data = diag_off_art_auditF, size=1, aes(x = time, ymin = Lower, ymax = Upper), col ="#e41a1c") +
ggtitle("Numbers of all women in Grand Cotonou diagnosed NOT on ART compared to data in red (Pre-ART audit report 2017)")
g10a
g10b = ggplot() + geom_line(data = I2xM, aes(x = time, y = value, factor = variable, factor = run), colour = "black", alpha = 1) + theme_bw() +
labs(y = "")+
scale_x_continuous(breaks = seq(min(time), max(time)+1, 4) ) +
scale_y_continuous(breaks = seq(min(Diagnosed$value), max(Diagnosed$value), 1000) ) +
geom_vline(xintercept=2005, col = "black", alpha = 0.5, linetype = "dashed") +
geom_text(aes(x=2005, label="CD4 < 200", y=-400))+
geom_vline(xintercept=2012, col = "black", alpha = 0.5, linetype = "dashed") +
geom_text(aes(x=2012, label="CD4 < 350", y=-800))+
geom_vline(xintercept=2015, col = "black", alpha = 0.5, linetype = "dashed") +
geom_text(aes(x=2015, label="CD4 < 500", y=-400))+
geom_vline(xintercept=2016, col = "black", alpha = 0.5, linetype = "dashed") +
theme(text = element_text(size=20), legend.position = "none",plot.title = element_text(hjust = 0.5))+
geom_errorbar(data = diag_off_art_auditM, size=1, aes(x = time, ymin = Lower, ymax = Upper), col ="#e41a1c") +
ggtitle("Numbers of all men in Grand Cotonou diagnosed NOT on ART compared to data in red (Pre-ART audit report 2017)")
g10b
I2xF_2017 = data.frame(
# Women_2016=I2xF[I2xF$time == 2016, "value"],
Women_2017=I2xF[I2xF$time == 2017, "value"],
# Men_2016=I2xM[I2xM$time == 2016, "value"],
Men_2017=I2xM[I2xM$time == 2017, "value"]
)
I2xF_2017_melted = reshape2::melt(I2xF_2017)
g10c = ggplot(I2xF_2017_melted, aes(x = variable, y=value)) + geom_boxplot() + theme_bw() +
geom_errorbar(aes(x = "Women_2017", ymax = 201, ymin = 215), col = "red", width = I(0.2))+
geom_errorbar(aes(x = "Men_2017", ymax = 108, ymin = 117), col = "red", width = I(0.2))+
theme(text = element_text(size=20), legend.position = "none",plot.title = element_text(hjust = 0.5)) +
labs(x = "Sex", y = "")+ ggtitle("Across the difference runs, the number of women and men in Grand Cotonou \n diagnosed but NOT on ART in 2017 vs the data in red for 2017")
g10c
I3xF_2017 = data.frame(Women=I3xF[I3xF$time == 2017, "value"], Men=I3xM[I3xM$time == 2017, "value"])
I3xF_2017_melted = reshape2::melt(I3xF_2017)
g10d = ggplot(I3xF_2017_melted, aes(x = variable, y=value)) + geom_boxplot() + theme_bw() +
geom_errorbar(aes(x = "Women", ymax = 6155, ymin = 13050), col = "red", width = I(0.2))+
geom_errorbar(aes(x = "Men", ymax = 2369, ymin = 5223), col = "red", width = I(0.2))+
theme(text = element_text(size=20), legend.position = "none",plot.title = element_text(hjust = 0.5)) +
labs(x = "Sex", y = "")+
ggtitle("Number of women and men in Grand Cotonou diagnosed \n ON ART in 2017 vs the data in red")
g10d
I4xF_2015 = data.frame(Women=I4xF[I4xF$time == 2015, "value"], Men=I4xM[I4xM$time == 2015, "value"])
I4xF_2017_melted = data.frame(variable = "All", value = rowSums(I4xF_2015))
g12=ggplot() + geom_line(data = Diagnosed_Women_Men_ratio[Diagnosed_Women_Men_ratio$variable == "Diagnosed On ART",], aes(x = time, y = value, factor = variable, factor = run), alpha = 1) +
theme_bw() + labs(y = "")+
ggtitle("Ratio of Women:Men on ART compared to data (2017) in red from audit 2017")+
scale_x_continuous(breaks = seq(min(time), (max(time)+1), 4) ) +
scale_y_continuous(breaks = seq(1, 5, 0.2) )+
theme(text = element_text(size=20), legend.position = "none",plot.title = element_text(hjust = 0.5)) +
geom_point(aes(x = 2017, y = 2.5), col = "red", size = I(03))+
# geom_hline(yintercept = c(1,2), col = "orange") +
theme(text = element_text(size=24),
legend.text=element_text(size=18))#+ annotate("text", x = 1994, y = 2.5, size = 7,col = "orange")
g12
g10e = ggplot(I4xF_2017_melted, aes(x = variable, y=value)) + geom_boxplot() + theme_bw() +
# geom_errorbar(aes(x = "Women", ymax = 6155, ymin = 13050), col = "red", width = I(0.2))+
geom_errorbar(aes(x = "All", ymax = 56, ymin = 56), col = "red", width = I(0.2))+
theme(text = element_text(size=20), legend.position = "none",plot.title = element_text(hjust = 0.5)) +
labs(x = "Sex", y = "")+
ggtitle("Number of women and men in Grand Cotonou in dropout \n from ART in 2017 vs the data in red (Lower bound)")
g10e
g10f = ggplot() + geom_line(data = ART_inits_cumu_melted, aes(x = time, y = value, factor = run), alpha = 1) + theme_bw() +
labs(y = "")+
theme(text = element_text(size=20), legend.position = "none",plot.title = element_text(hjust = 0.5))+
ggtitle("Cumulative ART initiations of all groups combined")
g10f
g10g = ggplot() + geom_line(data = ART_inits_melted, aes(x = time, y = value, factor = run), alpha = 1) + theme_bw() +
labs(y = "")+
scale_x_continuous(breaks = seq(min(time), max(time), 2) ) +
theme(text = element_text(size=20), legend.position = "none",plot.title = element_text(hjust = 0.5)) +
geom_errorbar(aes(x = 2015, ymax = 2051, ymin = 1318), col = "red", width = I(0.5), size = I(02))+
ggtitle("Numbers of ART initiations each year of all groups combined vs. the data in red")
g10g
g10h=ggplot() + geom_line(data = ART_inits_ratio_F_over_M_melted, aes(x = time, y = value, factor = variable), alpha = 1) +
theme_bw() + labs(y = "")+
scale_x_continuous(breaks = seq(min(time), (max(time)+3), 2) ) +
scale_y_continuous(breaks = seq(1, 5, 0.2) )+
theme(text = element_text(size=20),plot.title = element_text(hjust = 0.5)) +
geom_point(aes(x = 2015, y = 2.4), col = "red", size = I(03))+
ggtitle("Ratio of ART initiations women:Men, \n compared with data from audit 2017 in red")
g10h
g10i=ggplot() + geom_line(data = FtM_ratio_ART_initiation_rates_from_HIVpos_not_on_ART_melted, aes(x = time, y = value, factor = variable), alpha = 1) +
theme_bw() + labs(y = "")+
scale_x_continuous(breaks = seq(min(time), (max(time)+3), 2) ) +
scale_y_continuous(breaks = seq(1, 5, 0.2) )+
theme(text = element_text(size=20),plot.title = element_text(hjust = 0.5)) +
geom_point(aes(x = 2015, y = 1.7), col = "red", size = I(03))+
ggtitle("Ratio of ART initiation rate from HIV+ not on ART Women:Men, \n compared with data from audit 2017 in red")
g10i
g10j = ggplot(ART_init_rate_from_all_HIV_pos_NOT_ON_ART_by_sex_melted, aes(x = variable, y = value)) + geom_boxplot() + theme_bw() +
theme(text = element_text(size=20),plot.title = element_text(hjust = 0.5)) +
labs(y = "")+
geom_point(aes(x = "Women_2015", y = 0.169), col = "red", size = 3)+
geom_point(aes(x = "Men_2015", y = 0.1), col = "red", size = 3)+
ggtitle("ART initiation rate from HIV+ not on ART Women & Men in 2014 and 2015, \n compared with data for 2015 in red")
g10j
# Ratio female:male of rate of initiation per infected = 8.0/4.68 = 1.7
g10k = ggplot() + geom_line(data = ART_REinits_cumu_melted, aes(x = time, y = value, factor = run), alpha = 1) + theme_bw() +
labs(y = "")+
scale_x_continuous(breaks = seq(min(time), max(time), 2) ) +
theme(text = element_text(size=20), legend.position = "none",plot.title = element_text(hjust = 0.5)) +
# geom_errorbar(aes(x = 2015, ymax = 2051, ymin = 1318), col = "red", width = I(0.5), size = I(02))+
ggtitle("Cumulative numbers of ART RE initiations each year of all groups combined")
g10k
g10l = ggplot() + geom_line(data = ART_REinits_melted, aes(x = time, y = value, factor = run), alpha = 1) + theme_bw() +
labs(y = "")+
scale_x_continuous(breaks = seq(min(time), max(time), 2) ) +
theme(text = element_text(size=20), legend.position = "none",plot.title = element_text(hjust = 0.5)) +
geom_point(aes(x = 2015, y = 213), size = 5, col = "red") +
ggtitle("Numbers of ART RE initiations each year of all groups combined vs. the data in red") +
geom_text(col = "red",aes(x=2002, label="Data from monitoring reports for 2015 \n for Littoral + Atlantique + Oueme", y=300), size = 8)
g10l
g10m = ggplot() + geom_line(data = frac_ART_inits_REinits_melted[frac_ART_inits_REinits_melted$time > 1999,], aes(x = time, y = value, factor = run), alpha = 1) + theme_bw() +
labs(y = "")+
scale_x_continuous(breaks = seq(min(time), max(time), 2) ) +
theme(text = element_text(size=20), legend.position = "none",plot.title = element_text(hjust = 0.5)) +
geom_point(aes(x = 2013.5, y = 8), size = 5, col = "red") +
geom_point(aes(x = 2014.5, y = 7), size = 5, col = "red") +
geom_point(aes(x = 2015, y = 7), size = 5, col = "red") +
geom_point(aes(x = 2015.5, y = 9), size = 5, col = "red") +
geom_point(aes(x = 2016, y = 12), size = 5, col = "red") +
geom_point(aes(x = 2017, y = 11), size = 5, col = "red") +
ggtitle("Percent of ART initiations which are RE initiations (%) vs. the data in red") +
geom_text(col = "red",aes(x=2005, label="Data from monitoring reports by semester for BENIN", y=30), size = 8)
g10m
# ggplot() + geom_line(data = I2x, aes(x = time, y = value, factor = variable, factor = run, colour = variable), alpha = 1) + theme_bw() + labs(y = "Numbers of all people in Grand Cotonou")+
# scale_x_continuous(breaks = seq(min(time), max(time)+1, 4) ) +
# scale_y_continuous(breaks = seq(min(Diagnosed$value), max(Diagnosed$value), 1000) ) +
# geom_vline(xintercept=2005, col = "black", alpha = 0.5, linetype = "dashed") +
# geom_text(aes(x=2005, label="CD4 < 200", y=-400))+
# geom_vline(xintercept=2012, col = "black", alpha = 0.5, linetype = "dashed") +
# geom_text(aes(x=2012, label="CD4 < 350", y=-800))+
# geom_vline(xintercept=2015, col = "black", alpha = 0.5, linetype = "dashed") +
# geom_text(aes(x=2015, label="CD4 < 500", y=-400))+
# geom_vline(xintercept=2016, col = "black", alpha = 0.5, linetype = "dashed")
#
#
# ggplot() + geom_line(data = Diagnosed_Women_Men, aes(x = time, y = value, factor = variable, factor = run, colour = variable, linetype = group), alpha = 1) + theme_bw() + labs(y = "Numbers of all people in Grand Cotonou")+
# scale_x_continuous(breaks = seq(min(time), max(time), 4) ) +
# scale_y_continuous(breaks = seq(min(Diagnosed$value), max(Diagnosed$value), 1000) )
#
# ggplot() + geom_line(data = Diagnosed_Women_Men_ratio, aes(x = time, y = value, factor = variable, factor = run, colour = variable), alpha = 1) +
# theme_bw() + labs(y = "Ratio of Women:Men")+
# scale_x_continuous(breaks = seq(min(time), max(time), 4) ) +
# scale_y_continuous(breaks = seq(1, 5, 0.2) )
# ggplot() + geom_line(data = Diagnosed, aes(x = time, y = value, factor = variable, factor = run, colour = group), alpha = 1) + theme_bw() + facet_wrap(~variable, scales = "fixed") + labs(y = "Numbers of all people in Grand Cotonou") +
# scale_x_continuous(breaks = seq(min(time), max(time), 4) ) +
# scale_y_continuous(breaks = seq(min(Diagnosed$value), max(Diagnosed$value), 1000) )
#
# #
g11=ggplot() + geom_line(data = Diagnosed_Women_Men_ratio[Diagnosed_Women_Men_ratio$variable == "Diagnosed Off ART",], aes(x = time, y = value, factor = variable, factor = run, colour = variable), alpha = 1) +
theme_bw() + labs(y = "Ratio of Women:Men")+
scale_x_continuous(breaks = seq(min(time), max(time), 4) ) +
scale_y_continuous(breaks = seq(1, 5, 0.2) )+
theme(text = element_text(size=20))
# cross validation to ratio female:male on ART from audit report 2017
HIV_deaths_aidsinfo = data.frame(time = seq(1990, 2016),
lower = c(100,100,100,100,100,100,100,100,200,200,200,500,500,500,500,500,500,500,500,200,200,200,200,200,200,200,100),
upper = c(300,300,300,300,300,400,400,500,600,900,1200,1500,1500,1500,1500,2000,2000,2000,2000,900,900,800,1300,1300,1300,1300,700))
g13=ggplot() + geom_line(data = HIV_deaths_melted, aes(x = time, y = value, factor = variable), alpha = 1) + theme_bw() +
labs(y = "HIV deaths") +
geom_errorbar(data = HIV_deaths_aidsinfo, aes(x = time, ymin = lower, ymax = upper), col = "orange")+
theme(text = element_text(size=20)) +
theme(text = element_text(size=24),
legend.text=element_text(size=18))+ annotate("text", x = 1992, y = 1700, size = 7,col = "orange",
label = "Cross validation")
g13
# EQS 2012 mapping 889 – 1391
EQS_mapping = data.frame(time = 2012, lower = 889, upper = 1391)
g14=ggplot() + geom_line(data = N_Pro_FSW_melted, aes(x = time, y = value, factor = run), alpha = 1) + theme_bw() +
labs(y = "Number of Pro FSW") +
theme(text = element_text(size=20)) +
geom_errorbar(data = EQS_mapping, aes(x = time, ymin = lower, ymax = upper), col = "red")+
theme(text = element_text(size=24),
legend.text=element_text(size=18))+ annotate("text", x = 1992, y = 1200, size = 7,col = "red",
label = "Fitting to EQS mapping,\n CI given by Michel")
g14
N_low_level = data.frame(time = 2012, lower = 2000, upper = 3000)
g15=ggplot() + geom_line(data = N_Low_FSW_melted, aes(x = time, y = value, factor = run), alpha = 1) + theme_bw() +
labs(y = "Number of Low level FSW") +
theme(text = element_text(size=20)) +
geom_errorbar(data = N_low_level, aes(x = time, ymin = lower, ymax = upper), col = "orange")+
theme(text = element_text(size=24),
legend.text=element_text(size=18))+ annotate("text", x = 1992, y = 2750, size = 7,col = "orange",
label = "Cross validating")
low_pro_ratio = N_Low_FSW_melted
low_pro_ratio$value = N_Low_FSW_melted$value / N_Pro_FSW_melted$value
g22 = ggplot() + geom_line(data = low_pro_ratio, aes(x = time, y = value, factor = run), alpha = 1) + theme_bw() +
labs(y = "Ratio of low level FSW to pro FSW") +
theme(text = element_text(size=20)) +
theme(text = element_text(size=24),
legend.text=element_text(size=18))+ annotate("text", x = 1992, y = 4.5, size = 7,col = "red",
label = "Fitting") +
geom_hline(yintercept = c(1, 5), col = "red")
# lowndes 2002 19970
g16=ggplot() + geom_line(data = N_Client_melted, aes(x = time, y = value, factor = run), alpha = 1) + theme_bw() +
labs(y = "Number of Clients")
g17=ggplot() + geom_line(data = condom_Pro_FSW_melted, aes(x = time, y = value, colour = variable, factor = run), alpha = 1) +
theme_bw() + labs(y = "Condom use of Pro FSW")+
theme(text = element_text(size=20))
g18=ggplot() + geom_line(data = condom_Low_FSW_melted, aes(x = time, y = value, colour = variable, factor = run), alpha = 1) +
theme_bw() + labs(y = "Condom use of Low level FSW")+
theme(text = element_text(size=20))
g19=ggplot() + geom_line(data = condom_GPF_noncomm_melted, aes(x = time, y = value, factor = variable), alpha = 1) +
theme_bw() + labs(y = "Condom use of GPF")+
theme(text = element_text(size=20))
# ggplot() + geom_line(data = condom_GPM_noncomm_melted, aes(x = time, y = value, factor = variable), alpha = 1) +
# theme_bw() + labs(y = "Condom use of GPM")
#
g20=ggplot() + geom_line(data = annual_client_volume_pro_FSW_melted, aes(x = time, y = value, factor = variable), alpha = 1) +
theme_bw() + labs(y = "Annual client volume for Pro FSW")+
theme(text = element_text(size=20))
g21=ggplot() + geom_line(data = pfFSW_melted, aes(x = time, y = value*100, factor = run), alpha = 1) +
theme_bw() + labs(y = "HIV prevalence of incoming FSWs from neighbouring countries(%)")+
theme(text = element_text(size=20))
plotlist = list(g1, g3, g4,g4b, g5,g5b, g6, g7, g8, g9, g10,
g10a, g10b, g10c, g10d,g12, g10e, g10f, g10g, g10h, g10i, g10j,
g13, g14, g15, g16, g17, g18, g19, g20, g21, g22)
# outputfolder = "0103"
for(i in 1:length(plotlist))
{
jpeg(filename=paste0("C:\\Users\\eg1012\\Google Drive\\Imperial\\Cotonou\\Models\\Transmission Model\\Outputs/",
batch_folder, "/","all", i,".jpg"), width = 1000, height = 600)
print(plotlist[i][[1]])
dev.off()
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.