# 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)
# 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_ext_ent_full, inc_transm_ext_20to30, inc_transm_ext_ent_ge30,
newinf_transm_full, newinf_transm_20to30, newinf_transm_ge30
)
## # A tibble: 42 x 9
## prevalence prev_20to30 prev_ge30 inc_transm_ext_ent… inc_transm_ext_20…
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 31.3 27.9 56.1 5.14 5.42
## 2 31.1 27.9 55.4 4.71 4.82
## 3 30.0 26.6 54.0 4.10 4.12
## 4 29.3 25.9 53.2 3.66 3.75
## 5 28.3 24.6 52.4 3.29 3.26
## 6 28.0 24.3 51.4 2.87 2.87
## 7 31.7 28.5 55.9 4.89 5.22
## 8 31.0 27.7 55.0 4.94 5.30
## 9 29.9 26.4 54.1 4.04 4.20
## 10 29.4 26.0 53.1 3.66 3.80
## # ... with 32 more rows, and 4 more variables:
## # inc_transm_ext_ent_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
# set input
ranges = c(20, 29, 30, 34)
# call function
dt <- summarize_directory(loc, ranges=ranges)
# obtain tenth year prevalence
prev.dt <- dt[["prevs"]]
prev.dt[, param_set := rep(1:42, each=100)]
prev.dt[, seed := rep(rep(1:10, each=10), 42)]
prev.dt.10 <- prev.dt[ycat == 10]
avg.prev.10yr <- prev.dt.10[, .(avg = mean(mean),
avg_1=mean(mean_1),
avg_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 <- dt[["incs"]]
inc.dt[, param_set := rep(1:42, each=100)]
inc.dt[, seed := rep(rep(1:10, each=10), 42)]
inc.dt.10 <- inc.dt[ycat == 10]
avg.inc.10yr <- inc.dt.10[, .(avg = mean(full_inc_mean),
avg_1=mean(full_inc_mean_1),
avg_2=mean(full_inc_mean_2)),
by=param_set]
write.csv(avg.inc.10yr, file="out/finer-ret-scale-simple-fix-tenyearinc.csv")
# 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)
avg.inc.10yr[(instances %/% 10)+1,] #outcomes
## param_set avg avg_1 avg_2
## 1: 1 2.928766 3.170921 3.517576
## 2: 2 2.219319 2.198376 3.340121
## 3: 3 3.362988 3.532704 4.678183
## 4: 4 3.888863 4.162984 5.151426
## 5: 5 3.940419 4.134094 5.574615
## 6: 6 3.268182 3.393609 4.566469
avg.prev.10yr[(instances %/% 10)+1,]#outcomes
## param_set avg avg_1 avg_2
## 1: 1 24.87569 21.49115 46.81848
## 2: 2 23.17132 19.60936 44.36740
## 3: 3 26.24690 22.65603 49.21755
## 4: 4 28.78796 25.25648 52.66046
## 5: 5 28.31240 24.99845 51.52118
## 6: 6 26.73589 23.21845 49.94654
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)
avg.inc.10yr[(instances %/% 10),] #outcomes
## param_set avg avg_1 avg_2
## 1: 6 3.268182 3.393609 4.566469
## 2: 7 2.902677 2.981638 4.290509
## 3: 8 2.577205 2.741957 3.561334
## 4: 9 2.263576 2.312849 2.997994
## 5: 10 4.337166 4.678085 5.742401
## 6: 11 3.783215 4.112595 4.836742
avg.prev.10yr[(instances %/% 10),]
## param_set avg avg_1 avg_2
## 1: 6 26.73589 23.21845 49.94654
## 2: 7 25.63116 21.70287 49.25254
## 3: 8 24.61189 21.11406 46.65316
## 4: 9 22.64407 18.92724 44.16184
## 5: 10 29.51878 26.41256 53.00337
## 6: 11 28.52639 25.15322 52.34238
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)
avg.inc.10yr[(instances %/% 10),] #outcomes
## param_set avg avg_1 avg_2
## 1: 36 2.231206 2.220044 3.483239
## 2: 37 2.678126 2.839922 3.555495
## 3: 38 2.682054 2.859131 3.450525
## 4: 39 4.181784 4.551016 5.182648
## 5: 40 4.072229 4.442913 5.053291
## 6: 41 3.367206 3.573004 4.329220
avg.prev.10yr[(instances %/% 10),]
## param_set avg avg_1 avg_2
## 1: 36 22.96361 19.28122 44.77001
## 2: 37 23.35776 19.60337 45.42309
## 3: 38 24.31546 20.84455 45.96893
## 4: 39 29.73143 26.36205 53.79677
## 5: 40 28.40354 25.26861 51.52285
## 6: 41 26.72329 23.03008 50.46884
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))
avg.inc.10yr[(instances %/% 10)+1,] #outcomes
## param_set avg avg_1 avg_2
## 1: 1 2.928766 3.170921 3.517576
## 2: 7 2.902677 2.981638 4.290509
## 3: 13 3.268258 3.404504 4.543762
## 4: 19 3.491043 3.756467 4.414642
## 5: 25 3.674073 3.828298 5.273357
## 6: 31 3.763120 3.976919 5.089353
## 7: 37 2.678126 2.839922 3.555495
avg.prev.10yr[(instances %/% 10)+1,]
## param_set avg avg_1 avg_2
## 1: 1 24.87569 21.49115 46.81848
## 2: 7 25.63116 21.70287 49.25254
## 3: 13 26.30037 22.63343 49.31533
## 4: 19 27.20402 23.63585 50.58056
## 5: 25 28.06193 24.49705 51.97699
## 6: 31 28.44314 24.95468 52.34155
## 7: 37 23.35776 19.60337 45.42309
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 avg_1 avg_2
## 1: 2 2.219319 2.198376 3.340121
## 2: 8 2.577205 2.741957 3.561334
## 3: 14 2.766974 2.949722 3.655929
## 4: 20 3.181791 3.376253 4.350116
## 5: 26 3.262191 3.384837 4.824910
## 6: 32 3.519112 3.617153 4.878617
## 7: 38 2.682054 2.859131 3.450525
avg.prev.10yr[(instances %/% 10)+1,]
## param_set avg avg_1 avg_2
## 1: 2 23.17132 19.60936 44.36740
## 2: 8 24.61189 21.11406 46.65316
## 3: 14 25.02445 21.09733 48.56232
## 4: 20 25.94901 22.25740 49.24862
## 5: 26 26.75122 23.05620 50.42007
## 6: 32 27.68454 24.14497 51.11163
## 7: 38 24.31546 20.84455 45.96893
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)
avg.inc.10yr[(instances %/% 10)+1,] #outcomes
## param_set avg avg_1 avg_2
## 1: 4 3.888863 4.162984 5.151426
## 2: 10 4.337166 4.678085 5.742401
## 3: 16 2.623057 2.884943 3.339083
## 4: 22 2.239443 2.328124 3.053183
## 5: 28 2.663192 2.772312 3.814459
## 6: 34 2.879779 3.059467 3.667090
## 7: 40 4.072229 4.442913 5.053291
avg.prev.10yr[(instances %/% 10)+1,]
## param_set avg avg_1 avg_2
## 1: 4 28.78796 25.25648 52.66046
## 2: 10 29.51878 26.41256 53.00337
## 3: 16 24.12519 20.68156 45.90032
## 4: 22 23.17840 19.24572 45.66207
## 5: 28 24.63218 20.92702 47.27975
## 6: 34 25.23030 21.62914 47.63809
## 7: 40 28.40354 25.26861 51.52285
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)
avg.inc.10yr[(instances %/% 10)+1,] #outcomes
## param_set avg avg_1 avg_2
## 1: 6 3.268182 3.393609 4.566469
## 2: 12 3.590351 3.834485 5.090767
## 3: 18 3.811055 4.052969 4.986550
## 4: 24 4.023896 4.373673 5.410518
## 5: 30 4.347269 4.546415 6.075718
## 6: 36 2.231206 2.220044 3.483239
## 7: 42 3.040515 3.273993 3.923342
avg.prev.10yr[(instances %/% 10)+1,]
## param_set avg avg_1 avg_2
## 1: 6 26.73589 23.21845 49.94654
## 2: 12 27.34727 23.80040 50.85419
## 3: 18 28.23326 24.71487 51.87023
## 4: 24 28.61777 25.14014 52.90467
## 5: 30 29.25896 25.99466 52.79065
## 6: 36 22.96361 19.28122 44.77001
## 7: 42 25.43980 21.83194 47.76844
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%")
# Population size summary ----
responses_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
responses_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="age-20-30-34-w-newinf.RData")
#