# PrEP retention and scale intervention for uptake = 5K
rm(list=ls())
#options(tibble.print_max = 50, tibble.print_min = 50, dplyr.width=Inf)
# Libraries ----
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(stringr)
# Data ----
full_res <- read.csv("../data/finer-ret-scale-simple-fix-final_results.csv", header=FALSE) #final results
loc <- "../../finer-ret-scale-simple-fix"
in1 <- read.table(paste0(loc,"/instance_1/output/parameters.txt"))
counts_dt <- read.csv(paste0(loc, "/instance_51/output/counts.csv"))
param_list <- read.csv("../../../data/prep-finer-retention-gradual-scale-interv.csv")
param_list_txt <- read.table("../../../data/prep-finer-retention-gradual-scale-interv.txt")
# Generate Response Dataframe ----
all_responses <- full_res[,c(1, 2, 4, 6, seq(17, 97, by=10), 98:109)]
names(all_responses) <- c("instance",
"prevalence", "prev20to30", "prevge30",
"inc_transm_full", "inc_transm_20to30", "inc_transm_ge30",
"inc_transm_ext_full", "inc_transm_ext_20to30", "inc_transm_ext_ge30",
"inc_transm_ext_ent_full", "inc_transm_ext_ent_20to30", "inc_transm_ext_ent_ge30",
"newinf_transm_full", "newinf_transm_20to30", "newinf_transm_ge30",
"newinf_transm_ext_full", "newinf_transm_ext_20to30", "newinf_transm_ext_ge30",
"newinf_transm_ext_ent_full", "newinf_transm_ext_ent_20to30", "newinf_transm_ext_ent_ge30",
"popsize_full", "popsize_20to30", "popsize_ge30"
)
all_responses$param_set <- rep(1:42, each=10)
all_responses$seed <- rep(1:10, 42)
dim(all_responses)
## [1] 420 27
# Compute means by seeds ----
by_param_set <- group_by(all_responses, param_set)
responses <- summarise(by_param_set,
prevalence = mean(prevalence),
prev_20to30 = mean(prev20to30),
prev_ge30 = mean(prevge30),
inc_transm_full = mean(inc_transm_full), #transm. olny
inc_transm_20to30 = mean(inc_transm_20to30),
inc_transm_ge30 = mean(inc_transm_ge30),
inc_transm_ext_full = mean(inc_transm_ext_full), #transm+ext.
inc_transm_ext_20to30 = mean(inc_transm_ext_20to30),
inc_transm_ext_ge30 = mean(inc_transm_ext_ge30),
inc_transm_ext_ent_full = mean(inc_transm_ext_ent_full),#transm+ext+entry
inc_transm_ext_ent_20to30 = mean(inc_transm_ext_ent_20to30),
inc_transm_ext_ent_ge30 = mean(inc_transm_ext_ent_ge30),
newinf_transm_full = mean(newinf_transm_full), #transm. olny
newinf_transm_20to30 = mean(newinf_transm_20to30),
newinf_transm_ge30 = mean(newinf_transm_ge30),
newinf_transm_ext_full = mean(newinf_transm_ext_full), #transm+ext.
newinf_transm_ext_20to30 = mean(newinf_transm_ext_20to30),
newinf_transm_ext_ge30 = mean(newinf_transm_ext_ge30),
newinf_transm_ext_ent_full = mean(newinf_transm_ext_ent_full),#transm+ext+entry
newinf_transm_ext_ent_20to30 = mean(newinf_transm_ext_ent_20to30),
newinf_transm_ext_ent_ge30 = mean(newinf_transm_ext_ent_ge30),
popsize_full = mean(popsize_full),
popsize_20to30 = mean(popsize_20to30),
popsize_ge30 = mean(popsize_ge30)
)
dim(responses)
## [1] 42 25
# All Responses ----
responses %>%
select(prevalence, prev_20to30, prev_ge30,
inc_transm_full, inc_transm_20to30, inc_transm_ge30,
newinf_transm_full, newinf_transm_20to30, newinf_transm_ge30
)
## # A tibble: 42 x 9
## prevalence prev_20to30 prev_ge30 inc_transm_full inc_transm_20to30
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 31.3 27.9 56.1 4.40 4.78
## 2 31.1 27.9 55.4 3.91 4.14
## 3 30.0 26.6 54.0 3.38 3.49
## 4 29.3 25.9 53.2 2.94 3.10
## 5 28.3 24.6 52.4 2.55 2.63
## 6 28.0 24.3 51.4 2.14 2.25
## 7 31.7 28.5 55.9 4.16 4.51
## 8 31.0 27.7 55.0 4.06 4.45
## 9 29.9 26.4 54.1 3.26 3.47
## 10 29.4 26.0 53.1 2.96 3.17
## # ... with 32 more rows, and 4 more variables: inc_transm_ge30 <dbl>,
## # newinf_transm_full <dbl>, newinf_transm_20to30 <dbl>,
## # newinf_transm_ge30 <dbl>
write.csv(responses, file="out/prep-finer-ret-scale-simple-fix-responses.csv")
# Select Unique Parameters From List ----
seq_of_ones <- seq(1, nrow(param_list), by=10)
uniq_param_list <- param_list[seq_of_ones,] #unique input parameters (without seed repetitions)
View(uniq_param_list) #unique input parameters (without seed repetitions)
write.csv(uniq_param_list, file="out/finer-ret-scale-simple-fix-param_list.csv")
# Summary of Results Regions ----
apply(responses, 2, summary)
## param_set prevalence prev_20to30 prev_ge30 inc_transm_full
## Min. 1.00 27.95009 24.28926 51.39962 2.035100
## 1st Qu. 11.25 29.12760 25.68321 53.02701 2.566475
## Median 21.50 30.05337 26.67738 54.13998 3.235050
## Mean 21.50 30.06730 26.65539 54.15937 3.206521
## 3rd Qu. 31.75 31.21506 27.86285 55.55633 3.925075
## Max. 42.00 32.03375 28.82924 56.69630 4.510200
## inc_transm_20to30 inc_transm_ge30 inc_transm_ext_full
## Min. 2.066500 2.651000 2.625500
## 1st Qu. 2.648025 3.433775 3.206675
## Median 3.445650 4.348300 3.885050
## Mean 3.388971 4.336986 3.847636
## 3rd Qu. 4.191900 5.083100 4.557300
## Max. 4.905100 6.437700 5.193900
## inc_transm_ext_20to30 inc_transm_ext_ge30 inc_transm_ext_ent_full
## Min. 2.680700 3.249900 2.706900
## 1st Qu. 3.288225 4.117975 3.303300
## Median 4.022300 5.033150 3.978750
## Mean 4.031764 4.970486 3.933443
## 3rd Qu. 4.796200 5.750975 4.646925
## Max. 5.597100 7.020000 5.283500
## inc_transm_ext_ent_20to30 inc_transm_ext_ent_ge30
## Min. 2.680700 3.249900
## 1st Qu. 3.288225 4.117975
## Median 4.022300 5.033150
## Mean 4.031764 4.970486
## 3rd Qu. 4.796200 5.750975
## Max. 5.597100 7.020000
## newinf_transm_full newinf_transm_20to30 newinf_transm_ge30
## Min. 133.790 93.9100 29.65000
## 1st Qu. 149.825 105.7175 33.20000
## Median 167.960 118.9550 36.47500
## Mean 167.221 118.7255 36.07524
## 3rd Qu. 186.190 132.2825 39.71000
## Max. 199.300 143.0900 41.65000
## newinf_transm_ext_full newinf_transm_ext_20to30
## Min. 160.0500 111.1400
## 1st Qu. 175.8575 123.1725
## Median 193.0250 135.8650
## Mean 192.6633 135.6864
## 3rd Qu. 211.4050 149.5300
## Max. 223.9300 159.3800
## newinf_transm_ext_ge30 newinf_transm_ext_ent_full
## Min. 33.37000 163.5700
## 1st Qu. 37.15500 179.0125
## Median 40.35500 196.7250
## Mean 40.01071 196.2769
## 3rd Qu. 43.65250 215.2725
## Max. 45.68000 227.5300
## newinf_transm_ext_ent_20to30 newinf_transm_ext_ent_ge30
## Min. 111.1400 33.37000
## 1st Qu. 123.1725 37.15500
## Median 135.8650 40.35500
## Mean 135.6864 40.01071
## 3rd Qu. 149.5300 43.65250
## Max. 159.3800 45.68000
## popsize_full popsize_20to30 popsize_ge30
## Min. 5551.000 3537.300 1288.200
## 1st Qu. 5575.325 3554.600 1290.900
## Median 5589.800 3568.050 1292.150
## Mean 5587.700 3566.862 1292.438
## 3rd Qu. 5601.200 3575.875 1294.050
## Max. 5619.100 3597.600 1296.200
# Annual prevalence and incidence ----
# source function
source("../../../R/ph_summarize_functions.R")
## Loading required package: data.table
##
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
##
## between, first, last
## Loading required package: gtools
# set input
ranges = c(20, 29, 30, 34)
# call function
dt <- summarize_directory(loc, ranges=ranges)
# dt_w_instance <-
# lapply(dt, function (x)
# mutate(x, instance = as.numeric(str_extract(expid, "(\\d)+")))
# )
# prev.dt <- dt_w_instance[["prevs"]]
# inc.dt <- dt_w_instance[["incs"]]
prev.dt <- dt[["prevs"]] %>%
mutate(ycat = as.numeric(ycat))
inc.dt <- dt[["incs"]] %>%
mutate(ycat = as.numeric(ycat))
prev.dt <- as.data.table(prev.dt)
inc.dt <- as.data.table(inc.dt)
# obtain tenth year prevalence
#prev.dt <- as.data.table(arrange(prev.dt, instance))
prev.dt[, param_set := rep(1:42, each=100)]
prev.dt[, seed := rep(rep(1:10, each=10), 42)]
avg.prev.across.years <-
prev.dt[, .(avg_prev_across_years = mean(mean),
avg_prev_across_year_age1 = mean(mean_1),
avg_prev_across_year_age2 = mean(mean_2)),
by=param_set]
prev.dt.10 <- prev.dt[ycat == 10]
avg.prev.10yr <- prev.dt.10[, .(avg_prev = mean(mean),
avg_prev_1=mean(mean_1),
avg_prev_2=mean(mean_2)),
by=param_set]
write.csv(avg.prev.10yr, file="out/finer-ret-scale-simple-fix-tenyearprev.csv")
# obtain tenth year incidence
#inc.dt <- as.data.table(arrange(inc.dt, instance))
inc.dt[, param_set := rep(1:42, each=100)]
inc.dt[, seed := rep(rep(1:10, each=10), 42)]
avg.inc.across.years <-
inc.dt[, .(avg_inc_across_years = mean(full_inc_mean),
avg_inc_across_year_age1 = mean(full_inc_mean_1),
avg_inc_across_year_age2 = mean(full_inc_mean_2)),
by=param_set]
inc.dt.10 <- inc.dt[ycat == 10]
avg.inc.10yr <- inc.dt.10[, .(avg_inc = mean(full_inc_mean),
avg_inc_age1=mean(full_inc_mean_1),
avg_inc_age2=mean(full_inc_mean_2)),
by=param_set]
write.csv(avg.inc.10yr, file="out/finer-ret-scale-simple-fix-tenyearinc.csv")
tenthyr.outcomes <- round(left_join(avg.prev.10yr, avg.inc.10yr, by ="param_set"), 2)
avg.across.years <- round(left_join(avg.prev.across.years, avg.inc.across.years,
by ="param_set"), 2)
head(tenthyr.outcomes)
## param_set avg_prev avg_prev_1 avg_prev_2 avg_inc avg_inc_age1
## 1 1 29.20 25.91 53.20 4.40 4.78
## 2 2 28.73 25.52 52.11 3.91 4.14
## 3 3 26.73 23.17 49.84 3.38 3.49
## 4 4 25.15 21.58 47.45 2.94 3.10
## 5 5 23.01 19.40 44.49 2.55 2.63
## 6 6 21.98 18.34 42.67 2.14 2.25
## avg_inc_age2
## 1 5.38
## 2 5.61
## 3 5.03
## 4 3.86
## 5 3.59
## 6 2.70
head(avg.across.years)
## param_set avg_prev_across_years avg_prev_across_year_age1
## 1 1 31.30 27.92
## 2 2 31.14 27.92
## 3 3 30.03 26.65
## 4 4 29.31 25.87
## 5 5 28.33 24.65
## 6 6 27.95 24.34
## avg_prev_across_year_age2 avg_inc_across_years avg_inc_across_year_age1
## 1 56.07 5.02 5.35
## 2 55.39 4.84 5.15
## 3 53.96 4.36 4.63
## 4 53.20 4.02 4.31
## 5 52.40 3.64 3.81
## 6 51.44 3.38 3.58
## avg_inc_across_year_age2
## 1 7.22
## 2 6.83
## 3 6.14
## 4 5.61
## 5 5.41
## 6 4.77
head(responses)
## # A tibble: 6 x 25
## param_set prevalence prev_20to30 prev_ge30 inc_transm_full
## <int> <dbl> <dbl> <dbl> <dbl>
## 1 1 31.3 27.9 56.1 4.40
## 2 2 31.1 27.9 55.4 3.91
## 3 3 30.0 26.6 54.0 3.38
## 4 4 29.3 25.9 53.2 2.94
## 5 5 28.3 24.6 52.4 2.55
## 6 6 28.0 24.3 51.4 2.14
## # ... with 20 more variables: inc_transm_20to30 <dbl>,
## # inc_transm_ge30 <dbl>, inc_transm_ext_full <dbl>,
## # inc_transm_ext_20to30 <dbl>, inc_transm_ext_ge30 <dbl>,
## # inc_transm_ext_ent_full <dbl>, inc_transm_ext_ent_20to30 <dbl>,
## # inc_transm_ext_ent_ge30 <dbl>, newinf_transm_full <dbl>,
## # newinf_transm_20to30 <dbl>, newinf_transm_ge30 <dbl>,
## # newinf_transm_ext_full <dbl>, newinf_transm_ext_20to30 <dbl>,
## # newinf_transm_ext_ge30 <dbl>, newinf_transm_ext_ent_full <dbl>,
## # newinf_transm_ext_ent_20to30 <dbl>, newinf_transm_ext_ent_ge30 <dbl>,
## # popsize_full <dbl>, popsize_20to30 <dbl>, popsize_ge30 <dbl>
tenthyr.outcomes <- apply(tenthyr.outcomes, c(1,2), function (x) round(x, 2))
# PrEP summary ----
qplot(data=counts_dt[-1,], x=tick/365, y=on_prep/uninfected, geom="line")
# pull out data with constant retention, increasing uptake
## example: retention = 180 days
instances <- seq(1, 51, by=10)
tenthyr.outcomes[(instances %/% 10)+1,] #10th year outcomes
## param_set avg_prev avg_prev_1 avg_prev_2 avg_inc avg_inc_age1
## [1,] 1 29.20 25.91 53.20 4.40 4.78
## [2,] 2 28.73 25.52 52.11 3.91 4.14
## [3,] 3 26.73 23.17 49.84 3.38 3.49
## [4,] 4 25.15 21.58 47.45 2.94 3.10
## [5,] 5 23.01 19.40 44.49 2.55 2.63
## [6,] 6 21.98 18.34 42.67 2.14 2.25
## avg_inc_age2
## [1,] 5.38
## [2,] 5.61
## [3,] 5.03
## [4,] 3.86
## [5,] 3.59
## [6,] 2.70
avg.across.years[(instances %/% 10)+1,]
## param_set avg_prev_across_years avg_prev_across_year_age1
## 1 1 31.30 27.92
## 2 2 31.14 27.92
## 3 3 30.03 26.65
## 4 4 29.31 25.87
## 5 5 28.33 24.65
## 6 6 27.95 24.34
## avg_prev_across_year_age2 avg_inc_across_years avg_inc_across_year_age1
## 1 56.07 5.02 5.35
## 2 55.39 4.84 5.15
## 3 53.96 4.36 4.63
## 4 53.20 4.02 4.31
## 5 52.40 3.64 3.81
## 6 51.44 3.38 3.58
## avg_inc_across_year_age2
## 1 7.22
## 2 6.83
## 3 6.14
## 4 5.61
## 5 5.41
## 6 4.77
uniq_param_list[(instances %/% 10)+1,4:7]#inputs
## prep.yearly.increment.lt prep.yearly.increment.gte
## 1 0.0000 0.0000
## 11 0.0126 0.0126
## 21 0.0326 0.0326
## 31 0.0526 0.0526
## 41 0.0726 0.0726
## 51 0.0926 0.0926
## prep.mean.days.usage.lt prep.mean.days.usage.gte
## 1 180 180
## 11 180 180
## 21 180 180
## 31 180 180
## 41 180 180
## 51 180 180
ret180 <- lapply(instances, function(x)
read.csv(paste0(loc, "/instance_", x, "/output/counts.csv"))
)
ret180 <- lapply(1:length(ret180), function (x) ret180[[x]][-1,])
names(ret180) <- c("13.7%-control", "20%", "30%", "40%", "50%", "60%")
ggplot(bind_rows(ret180, .id="f.prep.uptake"),
aes(x=tick/365, y=on_prep/uninfected, color=f.prep.uptake))+
geom_line()+
ggtitle("Retention 180 days")
## example: retention = 270 days
instances <- seq(61, 111, by=10)
tenthyr.outcomes[(instances %/% 10)+1,] #tenth year outcomes
## param_set avg_prev avg_prev_1 avg_prev_2 avg_inc avg_inc_age1
## [1,] 7 29.84 26.50 54.02 4.16 4.51
## [2,] 8 28.51 25.31 51.70 4.06 4.45
## [3,] 9 26.32 22.60 49.80 3.26 3.47
## [4,] 10 25.30 21.67 47.68 2.96 3.17
## [5,] 11 23.22 19.87 44.33 2.36 2.47
## [6,] 12 22.02 18.27 43.03 2.08 2.15
## avg_inc_age2
## [1,] 5.28
## [2,] 5.08
## [3,] 4.15
## [4,] 3.86
## [5,] 3.04
## [6,] 2.91
avg.across.years[(instances %/% 10)+1,] #avg outcomes
## param_set avg_prev_across_years avg_prev_across_year_age1
## 7 7 31.68 28.50
## 8 8 30.96 27.69
## 9 9 29.90 26.44
## 10 10 29.40 26.01
## 11 11 28.52 25.03
## 12 12 27.96 24.39
## avg_prev_across_year_age2 avg_inc_across_years avg_inc_across_year_age1
## 7 55.95 5.07 5.48
## 8 55.04 4.77 5.08
## 9 54.13 4.33 4.56
## 10 53.14 4.05 4.31
## 11 52.26 3.67 3.89
## 12 51.40 3.37 3.54
## avg_inc_across_year_age2
## 7 6.86
## 8 6.72
## 9 6.28
## 10 5.65
## 11 5.18
## 12 4.85
uniq_param_list[(instances %/% 10)+1,4:7]#inputs
## prep.yearly.increment.lt prep.yearly.increment.gte
## 61 0.0000 0.0000
## 71 0.0126 0.0126
## 81 0.0326 0.0326
## 91 0.0526 0.0526
## 101 0.0726 0.0726
## 111 0.0926 0.0926
## prep.mean.days.usage.lt prep.mean.days.usage.gte
## 61 270 270
## 71 270 270
## 81 270 270
## 91 270 270
## 101 270 270
## 111 270 270
ret270 <- lapply(instances, function(x)
read.csv(paste0(loc, "/instance_", x, "/output/counts.csv"))
)
ret270 <- lapply(1:length(ret270), function (x) ret270[[x]][-1,])
names(ret270) <- c("13.7%-control", "20%", "30%", "40%", "50%", "60%")
ggplot(bind_rows(ret270, .id="f.prep.uptake"),
aes(x=tick/365, y=on_prep/uninfected, color=f.prep.uptake))+
geom_line()+
ggtitle("Retention 270 days")
## example: retention = 730 days
instances <- seq(361, 411, by=10)
tenthyr.outcomes[(instances %/% 10)+1,] #tenth year outcomes
## param_set avg_prev avg_prev_1 avg_prev_2 avg_inc avg_inc_age1
## [1,] 37 29.74 26.47 53.59 4.50 4.68
## [2,] 38 28.93 25.56 52.62 3.96 4.24
## [3,] 39 27.61 23.88 51.53 3.35 3.47
## [4,] 40 26.15 22.38 49.61 3.07 3.16
## [5,] 41 23.86 20.07 46.08 2.62 2.72
## [6,] 42 23.00 19.22 45.01 2.11 2.17
## avg_inc_age2
## [1,] 6.44
## [2,] 5.08
## [3,] 4.73
## [4,] 4.20
## [5,] 3.38
## [6,] 2.95
avg.across.years[(instances %/% 10)+1,] #avg outcomes
## param_set avg_prev_across_years avg_prev_across_year_age1
## 37 37 31.80 28.51
## 38 38 31.44 28.16
## 39 39 30.84 27.44
## 40 40 30.32 26.97
## 41 41 29.22 25.74
## 42 42 28.80 25.27
## avg_prev_across_year_age2 avg_inc_across_years avg_inc_across_year_age1
## 37 56.33 5.14 5.48
## 38 55.77 4.92 5.26
## 39 55.13 4.59 4.88
## 40 54.37 4.30 4.60
## 41 53.04 3.84 4.09
## 42 52.61 3.62 3.84
## avg_inc_across_year_age2
## 37 7.25
## 38 6.90
## 39 6.48
## 40 5.78
## 41 5.31
## 42 5.05
uniq_param_list[(instances %/% 10)+1,4:7]#inputs
## prep.yearly.increment.lt prep.yearly.increment.gte
## 361 0.0000 0.0000
## 371 0.0126 0.0126
## 381 0.0326 0.0326
## 391 0.0526 0.0526
## 401 0.0726 0.0726
## 411 0.0926 0.0926
## prep.mean.days.usage.lt prep.mean.days.usage.gte
## 361 730 730
## 371 730 730
## 381 730 730
## 391 730 730
## 401 730 730
## 411 730 730
ret730 <- lapply(instances, function(x)
read.csv(paste0(loc, "/instance_", x, "/output/counts.csv"))
)
ret730 <- lapply(1:length(ret730), function (x) ret730[[x]][-1,])
names(ret730) <- c("13.7%-control", "20%", "30%", "40%", "50%", "60%")
ggplot(bind_rows(ret730, .id="f.prep.uptake"),
aes(x=tick/365, y=on_prep/uninfected, color=f.prep.uptake))+
geom_line()+
ggtitle("Retention 730 days")
## uptake = control level, vary retention
instances <- c(1, seq(61, 361, by=60))
tenthyr.outcomes[(instances %/% 10)+1,] #tenth year outcomes
## param_set avg_prev avg_prev_1 avg_prev_2 avg_inc avg_inc_age1
## [1,] 1 29.20 25.91 53.20 4.40 4.78
## [2,] 7 29.84 26.50 54.02 4.16 4.51
## [3,] 13 29.59 26.21 53.64 4.27 4.53
## [4,] 19 30.17 27.06 54.09 4.51 4.91
## [5,] 25 30.10 27.06 54.05 4.26 4.57
## [6,] 31 29.63 26.19 53.78 4.20 4.51
## [7,] 37 29.74 26.47 53.59 4.50 4.68
## avg_inc_age2
## [1,] 5.38
## [2,] 5.28
## [3,] 5.89
## [4,] 5.85
## [5,] 6.01
## [6,] 5.35
## [7,] 6.44
avg.across.years[(instances %/% 10)+1,] #avg outcomes
## param_set avg_prev_across_years avg_prev_across_year_age1
## 1 1 31.30 27.92
## 7 7 31.68 28.50
## 13 13 31.66 28.41
## 19 19 31.86 28.67
## 25 25 32.03 28.83
## 31 31 31.65 28.39
## 37 37 31.80 28.51
## avg_prev_across_year_age2 avg_inc_across_years avg_inc_across_year_age1
## 1 56.07 5.02 5.35
## 7 55.95 5.07 5.48
## 13 56.07 5.05 5.44
## 19 56.23 5.18 5.56
## 25 56.70 5.22 5.63
## 31 56.05 5.05 5.42
## 37 56.33 5.14 5.48
## avg_inc_across_year_age2
## 1 7.22
## 7 6.86
## 13 6.96
## 19 7.19
## 25 7.20
## 31 7.02
## 37 7.25
uniq_param_list[(instances %/% 10)+1,4:7]#inputs
## prep.yearly.increment.lt prep.yearly.increment.gte
## 1 0 0
## 61 0 0
## 121 0 0
## 181 0 0
## 241 0 0
## 301 0 0
## 361 0 0
## prep.mean.days.usage.lt prep.mean.days.usage.gte
## 1 180 180
## 61 270 270
## 121 365 365
## 181 450 450
## 241 540 540
## 301 630 630
## 361 730 730
uptake_control <- lapply(instances, function(x)
read.csv(paste0(loc, "/instance_", x, "/output/counts.csv"))
)
uptake_control <- lapply(1:length(uptake_control), function (x) uptake_control[[x]][-1,])
names(uptake_control) <- c("180d-control", "270d", "365d", "450d", "540d", "630d", "730d")
ggplot(bind_rows(uptake_control, .id="retention"),
aes(x=tick/365, y=on_prep/uninfected, color=retention))+
geom_line()+
ggtitle("Uptake: Control=13.7%")
## uptake = final 20%, vary retention
instances <- seq(11, 420, by=60)
avg.inc.10yr[(instances %/% 10)+1,] #outcomes
## param_set avg_inc avg_inc_age1 avg_inc_age2
## 1: 2 3.905369 4.137651 5.611401
## 2: 8 4.063552 4.451783 5.083601
## 3: 14 3.876223 4.110918 5.265175
## 4: 20 3.758785 4.049055 4.921404
## 5: 26 3.938484 4.232780 5.038376
## 6: 32 3.931684 4.209932 5.348477
## 7: 38 3.961347 4.240930 5.081770
avg.prev.10yr[(instances %/% 10)+1,]
## param_set avg_prev avg_prev_1 avg_prev_2
## 1: 2 28.72903 25.52474 52.10605
## 2: 8 28.51085 25.30748 51.70097
## 3: 14 28.60673 25.17921 52.42787
## 4: 20 28.60078 25.06677 52.74354
## 5: 26 28.52012 25.08827 52.03066
## 6: 32 28.99623 25.60659 53.14828
## 7: 38 28.92503 25.55694 52.62377
uptake_f_20pct <- lapply(instances, function(x)
read.csv(paste0(loc, "/instance_", x, "/output/counts.csv"))
)
uptake_f_20pct <- lapply(1:length(uptake_f_20pct), function (x) uptake_f_20pct[[x]][-1,])
names(uptake_f_20pct) <- c("180d-control", "270d", "365d", "450d", "540d", "630d", "730d")
ggplot(bind_rows(uptake_f_20pct, .id="retention"),
aes(x=tick/365, y=on_prep/uninfected, color=retention))+
geom_line()+
ggtitle("Uptake: Final=20%")
## uptake = final 40%, vary retention
instances <- seq(31, 420, by=60)
tenthyr.outcomes[(instances %/% 10)+1,] #tenth year outcomes
## param_set avg_prev avg_prev_1 avg_prev_2 avg_inc avg_inc_age1
## [1,] 4 25.15 21.58 47.45 2.94 3.10
## [2,] 10 25.30 21.67 47.68 2.96 3.17
## [3,] 16 25.36 21.33 49.12 2.81 2.89
## [4,] 22 25.83 21.97 49.23 3.02 3.18
## [5,] 28 25.89 22.23 49.03 3.11 3.28
## [6,] 34 25.85 22.15 48.75 3.21 3.43
## [7,] 40 26.15 22.38 49.61 3.07 3.16
## avg_inc_age2
## [1,] 3.86
## [2,] 3.86
## [3,] 4.17
## [4,] 3.98
## [5,] 4.36
## [6,] 4.49
## [7,] 4.20
avg.across.years[(instances %/% 10)+1,] #avg outcomes
## param_set avg_prev_across_years avg_prev_across_year_age1
## 4 4 29.31 25.87
## 10 10 29.40 26.01
## 16 16 29.79 26.25
## 22 22 29.90 26.54
## 28 28 30.08 26.71
## 34 34 29.95 26.50
## 40 40 30.32 26.97
## avg_prev_across_year_age2 avg_inc_across_years avg_inc_across_year_age1
## 4 53.20 4.02 4.31
## 10 53.14 4.05 4.31
## 16 54.12 4.12 4.42
## 22 53.80 4.17 4.47
## 28 54.15 4.25 4.59
## 34 53.89 4.24 4.48
## 40 54.37 4.30 4.60
## avg_inc_across_year_age2
## 4 5.61
## 10 5.65
## 16 5.76
## 22 5.78
## 28 5.76
## 34 6.01
## 40 5.78
uniq_param_list[(instances %/% 10)+1,4:7]#inputs
## prep.yearly.increment.lt prep.yearly.increment.gte
## 31 0.0526 0.0526
## 91 0.0526 0.0526
## 151 0.0526 0.0526
## 211 0.0526 0.0526
## 271 0.0526 0.0526
## 331 0.0526 0.0526
## 391 0.0526 0.0526
## prep.mean.days.usage.lt prep.mean.days.usage.gte
## 31 180 180
## 91 270 270
## 151 365 365
## 211 450 450
## 271 540 540
## 331 630 630
## 391 730 730
uptake_f_40pct <- lapply(instances, function(x)
read.csv(paste0(loc, "/instance_", x, "/output/counts.csv"))
)
uptake_f_40pct <- lapply(1:length(uptake_f_40pct), function (x) uptake_f_40pct[[x]][-1,])
names(uptake_f_40pct) <- c("180d-control", "270d", "365d", "450d", "540d", "630d", "730d")
ggplot(bind_rows(uptake_f_40pct, .id="retention"),
aes(x=tick/365, y=on_prep/uninfected, color=retention))+
geom_line()+
ggtitle("Uptake: Final=40%")
## uptake = final 60%, vary retention
instances <- seq(51, 420, by=60)
tenthyr.outcomes[(instances %/% 10)+1,] #tenth year outcomes
## param_set avg_prev avg_prev_1 avg_prev_2 avg_inc avg_inc_age1
## [1,] 6 21.98 18.34 42.67 2.14 2.25
## [2,] 12 22.02 18.27 43.03 2.08 2.15
## [3,] 18 21.89 18.22 42.90 2.07 2.07
## [4,] 24 22.11 18.26 43.98 2.09 2.18
## [5,] 30 22.35 18.53 44.16 2.06 2.11
## [6,] 36 22.54 18.75 44.41 2.04 2.10
## [7,] 42 23.00 19.22 45.01 2.11 2.17
## avg_inc_age2
## [1,] 2.70
## [2,] 2.91
## [3,] 2.97
## [4,] 2.84
## [5,] 2.82
## [6,] 2.65
## [7,] 2.95
avg.across.years[(instances %/% 10)+1,] #avg outcomes
## param_set avg_prev_across_years avg_prev_across_year_age1
## 6 6 27.95 24.34
## 12 12 27.96 24.39
## 18 18 27.96 24.29
## 24 24 28.29 24.68
## 30 30 28.51 24.94
## 36 36 28.54 25.05
## 42 42 28.80 25.27
## avg_prev_across_year_age2 avg_inc_across_years avg_inc_across_year_age1
## 6 51.44 3.38 3.58
## 12 51.40 3.37 3.54
## 18 51.73 3.41 3.55
## 24 52.08 3.48 3.68
## 30 52.19 3.51 3.71
## 36 52.17 3.54 3.75
## 42 52.61 3.62 3.84
## avg_inc_across_year_age2
## 6 4.77
## 12 4.85
## 18 5.13
## 24 4.95
## 30 4.93
## 36 4.94
## 42 5.05
uniq_param_list[(instances %/% 10)+1,4:7]#inputs
## prep.yearly.increment.lt prep.yearly.increment.gte
## 51 0.0926 0.0926
## 111 0.0926 0.0926
## 171 0.0926 0.0926
## 231 0.0926 0.0926
## 291 0.0926 0.0926
## 351 0.0926 0.0926
## 411 0.0926 0.0926
## prep.mean.days.usage.lt prep.mean.days.usage.gte
## 51 180 180
## 111 270 270
## 171 365 365
## 231 450 450
## 291 540 540
## 351 630 630
## 411 730 730
uptake_f_60pct <- lapply(instances, function(x)
read.csv(paste0(loc, "/instance_", x, "/output/counts.csv"))
)
uptake_f_60pct <- lapply(1:length(uptake_f_60pct), function (x) uptake_f_60pct[[x]][-1,])
names(uptake_f_60pct) <- c("180d-control", "270d", "365d", "450d", "540d", "630d", "730d")
ggplot(bind_rows(uptake_f_60pct, .id="retention"),
aes(x=tick/365, y=on_prep/uninfected, color=retention))+
geom_line()+
ggtitle("Uptake: Final=60%")
# Plot Annualized Incidences ----
#avg.inc.each.year <- inc.dt
inc.dt[, c("avg_inc_each_year",
"avg_inc_each_age1",
"avg_inc_year_age2") := list(
mean(full_inc_mean),
mean(full_inc_mean_1),
mean(full_inc_mean_2)),
by=param_set]
final_uptake <- c("13.7%-control", "20%",
"30%", "40%", "50%", "60%")
retention <- c("180d-control", "270d", "365d",
"450d", "540d", "630d", "730d")
summ_prev_dt <- prev.dt %>%
group_by(param_set, ycat) %>%
summarize(avg_annual_prev = mean(mean),
sd_annual_prev = sd(mean))
summ_prev_dt$final_uptake = rep(final_uptake, each = 10)
summ_prev_dt$retention <- rep(retention, each=420/7)
summ_inc_dt <- inc.dt %>%
group_by(param_set, ycat) %>%
summarize(avg_annual_inc = mean(full_inc_mean),
sd_annual_inc = sd(full_inc_mean))
summ_inc_dt$final_uptake = rep(final_uptake, each = 10)
summ_inc_dt$retention <- rep(retention, each=420/7)
p <-
ggplot(filter(summ_prev_dt, factor(param_set) %in% as.factor(1:6)),
aes(x=ycat, y=avg_annual_prev, color=final_uptake))+
geom_line()+geom_point()+
ggtitle("Average PrEP Retention = 6 months")+
xlab("year of intervention")+
ylab("prevalence (%)")+
scale_y_continuous(limit=c(10,35))+
scale_x_continuous(breaks = 1:10)
p+
geom_errorbar(data=p[["data"]],
aes(x=ycat, ymin=(avg_annual_prev-sd_annual_prev),
ymax=avg_annual_prev+sd_annual_prev),
width = 0.15)
p1 <-
ggplot(filter(summ_inc_dt, factor(param_set) %in% as.factor(1:6)),
aes(x=ycat, y=avg_annual_inc, color=final_uptake))+
geom_line()+geom_point()+
ggtitle("Average PrEP Retention = 6 months")+
xlab("year of intervention")+
ylab("incidence (per 100 py)")+
scale_y_continuous(limit=c(0,8))+
scale_x_continuous(breaks = 1:10)
p1+
geom_errorbar(data=p1[["data"]],
aes(x=ycat, ymin=(avg_annual_inc-sd_annual_inc),
ymax=avg_annual_inc+sd_annual_inc),
width = 0.15)
control.uptake.instances <- seq(1, 42, 6)
control.uptake.inc.data <-
filter(summ_inc_dt,factor(param_set) %in% control.uptake.instances)
control.uptake.prev.data <-
filter(summ_prev_dt,factor(param_set) %in% control.uptake.instances)
p2 <-
ggplot(control.uptake.prev.data,
aes(x=ycat, y=avg_annual_prev, color=retention))+
geom_line()+geom_point()+
ggtitle("PrEP Uptake = 13.7% (Control)")+
xlab("year of intervention")+
ylab("prevalence (%)")+
scale_y_continuous(limit=c(10,35))+
scale_x_continuous(breaks = 1:10)
p2+
geom_errorbar(data=p2[["data"]],
aes(x=ycat, ymin=(avg_annual_prev-sd_annual_prev),
ymax=avg_annual_prev+sd_annual_prev),
width = 0.15)
p3 <-
ggplot(control.uptake.inc.data,
aes(x=ycat, y=avg_annual_inc, color=retention))+
geom_line()+geom_point()+
ggtitle("PrEP Uptake = 13.7% (Control)")+
xlab("year of intervention")+
ylab("incidence (per 100 py)")+
scale_y_continuous(limit=c(0,8))+
scale_x_continuous(breaks = 1:10)
p3+
geom_errorbar(data=p3[["data"]],
aes(x=ycat, ymin=(avg_annual_inc-sd_annual_inc),
ymax=avg_annual_inc+sd_annual_inc),
width = 0.15)
## Retention at 30% final uptake
f.30pct.uptake.instances <- seq(3, 42, 6)
f.30pct.uptake.inc.data <-
filter(summ_inc_dt,factor(param_set) %in% f.30pct.uptake.instances)
f.30pct.uptake.prev.data <-
filter(summ_prev_dt,factor(param_set) %in% f.30pct.uptake.instances)
p4 <-
ggplot(f.30pct.uptake.prev.data,
aes(x=ycat, y=avg_annual_prev, color=retention))+
geom_line()+geom_point()+
ggtitle("Final PrEP Uptake = 30%")+
xlab("year of intervention")+
ylab("prevalence (%)")+
scale_y_continuous(limit=c(10,35))+
scale_x_continuous(breaks = 1:10)
p4 +
geom_errorbar(data=p4[["data"]],
aes(x=ycat, ymin=(avg_annual_prev-sd_annual_prev),
ymax=avg_annual_prev+sd_annual_prev),
width = 0.15)
p5 <-
ggplot(f.30pct.uptake.inc.data,
aes(x=ycat, y=avg_annual_inc, color=retention))+
geom_line()+geom_point()+
ggtitle("Final PrEP Uptake = 30%")+
xlab("year of intervention")+
ylab("incidence (per 100 py)")+
scale_y_continuous(limit=c(0,8))+
scale_x_continuous(breaks = 1:10)
p5+
geom_errorbar(data=p5[["data"]],
aes(x=ycat, ymin=(avg_annual_inc-sd_annual_inc),
ymax=avg_annual_inc+sd_annual_inc),
width = 0.15)
## Viewing p3[["data"]] shows that in the tenth year, control incidence (retention = 6 m, row 10 of p3 data),
## is less than the incidence for the case where retention = 2 yrs (row 70 of p3 data)
# Population size summary ----
avg.across.years_popsize <-
responses %>%
select(popsize_full, popsize_20to30, popsize_ge30)
# Compute annual growth rate
starting.popsize <- counts_dt$vertex_count[2]
#starting.popsize.20to30 <- counts_dt$vertex_count_18[2]
#starting.popsize.ge30 <- counts_dt$vertex_count_B[2]
n.ts <- nrow(counts_dt) - 1
n.yrs <- n.ts/365
growth.rate <- ((5600/5000)^(1/n.yrs) - 1)*100 #example
avg.across.years_popsize %>%
mutate(growthrate_allages = ((popsize_full/starting.popsize)^(1/n.yrs) - 1)*100
#growthrate_20to30 = ((popsize_20to30/starting.popsize.20to30)^(1/n.yrs) - 1)*100,
#growthrate_ge30 = ((popsize_ge30/starting.popsize.ge30)^(1/n.yrs) - 1)*100
)
## # A tibble: 42 x 4
## popsize_full popsize_20to30 popsize_ge30 growthrate_allages
## <dbl> <dbl> <dbl> <dbl>
## 1 5619 3579 1294 -0.0884
## 2 5615 3596 1290 -0.0950
## 3 5594 3568 1296 -0.133
## 4 5609 3594 1290 -0.107
## 5 5608 3585 1290 -0.109
## 6 5569 3552 1295 -0.179
## 7 5589 3572 1290 -0.142
## 8 5595 3569 1296 -0.132
## 9 5604 3586 1291 -0.116
## 10 5589 3571 1293 -0.141
## # ... with 32 more rows
# # Network data ----
#
# main_net <- readRDS(paste0(loc, "main_network_50000.RDS"))
# cas_net <- readRDS(paste0(loc, "casual_network_50000.RDS"))
# main_net
# cas_net
#
# # Check simulated vs target statistics ----
#
# ## mean degree
# network.edgecount(main_net)*2/network.size(main_net);
# network.edgecount(cas_net)*2/network.size(cas_net); 687*2/3000
#
# ## degree distribution
# degreedist(main_net)/network.size(main_net); c(1698, 1236, 54)/3000
# degreedist(cas_net)/network.size(cas_net); c(1767, 966, 204)/3000
#
# ## partnership duration
# ## see `compute-partnership-durations.R` in
# ##`/project/khanna7/Projects/midway2/BARS/transmission_model/r/testing` for example
# # pt_data <- fread(paste0(loc, "partnership_events.csv"))
# # pt_data_main <- pt_data[network_type == 0] #main network
# # pt_data_casual <- pt_data[network_type == 1] #casual network
# #
# # pt_data_main[,.(counts = .N),by=.(p1,p2),][,unique(counts)] # 1 2 4 3
# # pt_data_main[,.(counts = .N),by=.(p1,p2),][,.N,by=counts]
# # pt_data_casual[,.(counts = .N),by=.(p1,p2),][,unique(counts)] # 1 2 3 4 6
# # pt_data_casual[,.(counts = .N),by=.(p1,p2),][,.N,by=counts]
# #
# # res_main_2 <- pt_data_main[,.(tick, counts = .N),by=.(p1,p2),][counts == 2,.(dur = abs(diff(
# # tick))), by=.(p1,p2)][,dur]
# # res_main_4 <- pt_data_main[,.(tick, counts = .N),by=.(p1,p2),][counts == 4,.(dur = abs(diff(
# # tick))), by=.(p1,p2)][,dur]
# # res_casual_2 <- pt_data_casual[,.(tick, counts = .N),by=.(p1,p2),][counts == 2,.(dur = abs(diff(
# # tick))), by=.(p1,p2)][,dur]
# # res_casual_4 <- pt_data_casual[,.(tick, counts = .N),by=.(p1,p2),][counts == 4,.(dur = abs(diff(
# # tick))), by=.(p1,p2)][,dur]
# # res_casual_6 <- pt_data_casual[,.(tick, counts = .N),by=.(p1,p2),][counts == 6,.(dur = abs(diff(
# # tick))), by=.(p1,p2)][,dur]
# #
# # summary(c(res_main_2, res_main_4))
# # summary(c(res_casual_2, res_casual_4, res_casual_6))
# #
# #
# # ## mean age diff
# # main_el <- as.matrix(as.edgelist(main_net))
# # cas_el <- as.matrix(as.edgelist(cas_net))
# # mean(abs((main_net%v%"age")[main_el[,1]] - (main_net%v%"age")[main_el[,2]])); 2.91
# # mean(abs((cas_net%v%"age")[cas_el[,1]] - (cas_net%v%"age")[cas_el[,2]])); 3.14
#
# ## population distribution
# length(which(main_net%v%"age"<26))/network.size(main_net)
# length(which(main_net%v%"age" >= 26))/network.size(main_net)
#
# ## for one instance, add other forms of incidence
# summarize_inc <- function(filename="counts.csv"){
#
# counts <- read.csv(filename)
# counts <- counts[-1,]
# counts$total_inc <- counts$infected_via_transmission+counts$infected_at_entry+counts$infected_externally
# counts_yr_chunks <- split(counts,
# ceiling(seq_along(counts$tick)/365))
#
# mean_inc_total <- lapply(counts_yr_chunks, function (x)
# mean(x[2:nrow(x), "total_inc"] / x[1:(nrow(x) - 1), "uninfected"])
# )
# mean_inc_total_ten_yrs <- tail(mean_inc_total, 10)
# mean_inc_total_ten_yrs <- unlist(mean_inc_total_ten_yrs)*365*100
#
# mean_inc <- lapply(counts_yr_chunks, function (x)
# mean(x[2:nrow(x), "infected_via_transmission"] / x[1:(nrow(x) - 1), "uninfected"])
# )
# mean_inc_ten_yrs <- tail(mean_inc, 10)
# mean_inc_ten_yrs <- unlist(mean_inc_ten_yrs)*365*100
# full_result <- c(round(mean_inc_ten_yrs[10], 3), round(mean_inc_total_ten_yrs[10], 3))
#
# gte26_result <- round(mean_inc_ten_yrs, 3)
#
# return(c(full_result))
# }
#
# summarize_inc(paste0(loc, "counts.csv"))
#
# # Save data ----
#
save.image(file="out/finer-ret-scale-simple-fix.RData")
#