scripts/lazymcmc_test_cluster.R
In Geidelberg/cotonou:
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 = NULL
#####################################################
ART_data_points_lazymcmc = data.frame(time = c(2017, 2017, 2015),
variable = c("Men", "Women","Pro FSW"),
x = c(3908, 9784, 49))
good_previous_fit = get(load("best_pars_combined_0904_top.Rdata"))
good_previous_fit_variedpars = good_previous_fit[[1]][rownames(ranges)]
good_previous_fit = get(load("best_pars_combined_0904_top.Rdata"))
good_previous_fit_variedpars = good_previous_fit[[1]][rownames(ranges)]
lazymcmc_for_cluster(good_previous_fit_variedpars = good_previous_fit_variedpars,
ranges = ranges,
best_set = best_set, time = time,
par_seq = par_seq, condom_seq = condom_seq,
groups_seq = groups_seq, years_seq = years_seq, outputs = outputs,
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,
iterations = 200,
popt = 0.44,
opt_freq = 100,
thin = 1,
burnin = 0,
adaptive_period = 0,
save_block = 100,
filename = "lol")
lazymcmc_for_cluster <- function(good_previous_fit_variedpars, ranges, best_set, time,
par_seq, condom_seq, groups_seq, years_seq, outputs,
prev_points,
frac_N_discard_points,
Ntot_data_points,
ART_data_points,
PrEP_fitting,
iterations,
popt,
opt_freq,
thin,
burnin,
adaptive_period,
save_block, filename) {
# CREATING PARTAB WHICH IS THE INPUT DATAFRAME
parTab_create = data.frame(
names = rownames(ranges),
values = ranges[,1],
fixed = 0,
steps = 0.1,
lower_bound = ranges[,1],
upper_bound = ranges[,2]
)
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
}
}
parTab_create[parTab_create$lower_bound ==parTab_create$upper_bound, "fixed"] = 1
rownames(parTab_create) = NULL
###
# create_lik <- function(parTab, ranges, best_set, time,
# par_seq, condom_seq, groups_seq, years_seq, outputs,
# prev_points,
# frac_N_discard_points,
# Ntot_data_points,
# ART_data_points,
# PrEP_fitting,
# PRIOR_FUNC,...)
create_lik <- function(parTab_create, data, PRIOR_FUNC,...)
{
par_names <- rownames(ranges)
# print(ranges)
## using the `...` bit.
likelihood_func <- function(pars){
# print(ranges)
ranges_new = cbind(pars, pars)
rownames(ranges_new) = par_names
parameters = cotonou::lhs_parameters(1, 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)
res = cotonou::return_outputs(parameters[[1]], cotonou::main_model, time = time, outputs = outputs)
lik = cotonou::likelihood_lazymcmc(res, time = time,
prev_points = prev_points,
frac_N_discard_points = frac_N_discard_points,
Ntot_data_points = Ntot_data_points,
ART_data_points = ART_data_points,
PrEP_fitting = PrEP_fitting)
return(lik)
}
return(likelihood_func)
}
f <- create_lik(parTab_create = parTab_create, data = NULL, PRIOR_FUNC = NULL, ranges, best_set, time,
par_seq, condom_seq, groups_seq, years_seq, outputs,
prev_points,
frac_N_discard_points,
Ntot_data_points,
ART_data_points,
PrEP_fitting)
if(is.numeric(f(parTab_create$values))) {
mcmcPars <- c("iterations"=iterations,popt=popt,opt_freq=opt_freq,thin=thin,burnin=burnin,adaptive_period=adaptive_period,save_block=save_block)
res <- lazymcmc::run_MCMC(parTab = parTab_create,
mcmcPars = mcmcPars,
filename = filename,
CREATE_POSTERIOR_FUNC = create_lik,
mvrPars = NULL, PRIOR_FUNC = NULL, OPT_TUNING = 0.1,
ranges = ranges,
best_set = best_set, time = time,
par_seq = par_seq, condom_seq = condom_seq,
groups_seq = groups_seq, years_seq = years_seq, outputs = outputs,
prev_points = prev_points,
frac_N_discard_points = frac_N_discard_points,
Ntot_data_points = Ntot_data_points,
ART_data_points = ART_data_points,
PrEP_fitting = PrEP_fitting)
} else {
res <- "failed"
}
return(res)
}
mcmcPars <- c("iterations"=25000,popt=0.44,opt_freq=1000,thin=10,burnin=0,adaptive_period=4000,save_block=100)
mcmcPars <- c("iterations"=25000,
popt=0.44,
opt_freq=1000,
thin=10,
burnin=0,
adaptive_period=4000,
save_block=100)
Geidelberg/cotonou documentation built on Aug. 29, 2020, 4:22 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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 = NULL
#####################################################
ART_data_points_lazymcmc = data.frame(time = c(2017, 2017, 2015),
variable = c("Men", "Women","Pro FSW"),
x = c(3908, 9784, 49))
good_previous_fit = get(load("best_pars_combined_0904_top.Rdata"))
good_previous_fit_variedpars = good_previous_fit[[1]][rownames(ranges)]
good_previous_fit = get(load("best_pars_combined_0904_top.Rdata"))
good_previous_fit_variedpars = good_previous_fit[[1]][rownames(ranges)]
lazymcmc_for_cluster(good_previous_fit_variedpars = good_previous_fit_variedpars,
ranges = ranges,
best_set = best_set, time = time,
par_seq = par_seq, condom_seq = condom_seq,
groups_seq = groups_seq, years_seq = years_seq, outputs = outputs,
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,
iterations = 200,
popt = 0.44,
opt_freq = 100,
thin = 1,
burnin = 0,
adaptive_period = 0,
save_block = 100,
filename = "lol")
lazymcmc_for_cluster <- function(good_previous_fit_variedpars, ranges, best_set, time,
par_seq, condom_seq, groups_seq, years_seq, outputs,
prev_points,
frac_N_discard_points,
Ntot_data_points,
ART_data_points,
PrEP_fitting,
iterations,
popt,
opt_freq,
thin,
burnin,
adaptive_period,
save_block, filename) {
# CREATING PARTAB WHICH IS THE INPUT DATAFRAME
parTab_create = data.frame(
names = rownames(ranges),
values = ranges[,1],
fixed = 0,
steps = 0.1,
lower_bound = ranges[,1],
upper_bound = ranges[,2]
)
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
}
}
parTab_create[parTab_create$lower_bound ==parTab_create$upper_bound, "fixed"] = 1
rownames(parTab_create) = NULL
###
# create_lik <- function(parTab, ranges, best_set, time,
# par_seq, condom_seq, groups_seq, years_seq, outputs,
# prev_points,
# frac_N_discard_points,
# Ntot_data_points,
# ART_data_points,
# PrEP_fitting,
# PRIOR_FUNC,...)
create_lik <- function(parTab_create, data, PRIOR_FUNC,...)
{
par_names <- rownames(ranges)
# print(ranges)
## using the `...` bit.
likelihood_func <- function(pars){
# print(ranges)
ranges_new = cbind(pars, pars)
rownames(ranges_new) = par_names
parameters = cotonou::lhs_parameters(1, 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)
res = cotonou::return_outputs(parameters[[1]], cotonou::main_model, time = time, outputs = outputs)
lik = cotonou::likelihood_lazymcmc(res, time = time,
prev_points = prev_points,
frac_N_discard_points = frac_N_discard_points,
Ntot_data_points = Ntot_data_points,
ART_data_points = ART_data_points,
PrEP_fitting = PrEP_fitting)
return(lik)
}
return(likelihood_func)
}
f <- create_lik(parTab_create = parTab_create, data = NULL, PRIOR_FUNC = NULL, ranges, best_set, time,
par_seq, condom_seq, groups_seq, years_seq, outputs,
prev_points,
frac_N_discard_points,
Ntot_data_points,
ART_data_points,
PrEP_fitting)
if(is.numeric(f(parTab_create$values))) {
mcmcPars <- c("iterations"=iterations,popt=popt,opt_freq=opt_freq,thin=thin,burnin=burnin,adaptive_period=adaptive_period,save_block=save_block)
res <- lazymcmc::run_MCMC(parTab = parTab_create,
mcmcPars = mcmcPars,
filename = filename,
CREATE_POSTERIOR_FUNC = create_lik,
mvrPars = NULL, PRIOR_FUNC = NULL, OPT_TUNING = 0.1,
ranges = ranges,
best_set = best_set, time = time,
par_seq = par_seq, condom_seq = condom_seq,
groups_seq = groups_seq, years_seq = years_seq, outputs = outputs,
prev_points = prev_points,
frac_N_discard_points = frac_N_discard_points,
Ntot_data_points = Ntot_data_points,
ART_data_points = ART_data_points,
PrEP_fitting = PrEP_fitting)
} else {
res <- "failed"
}
return(res)
}
mcmcPars <- c("iterations"=25000,popt=0.44,opt_freq=1000,thin=10,burnin=0,adaptive_period=4000,save_block=100)
mcmcPars <- c("iterations"=25000,
popt=0.44,
opt_freq=1000,
thin=10,
burnin=0,
adaptive_period=4000,
save_block=100)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.