-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathload_params.R
141 lines (136 loc) · 8.19 KB
/
load_params.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
# functions
# rm(list=ls());
# currentdir_path=dirname(rstudioapi::getSourceEditorContext()$path); setwd(currentdir_path)
x1 <- c("here","tidyverse","deSolve","gtools","lubridate","Rcpp",
"RcppRoll","pracma","lhs","ISOweek","tidybayes")
# "rstudioapi","devtools","wpp2019",
# ,"epiR"
x2 <- x1 %in% row.names(installed.packages())
if (any(x2 == FALSE)) { install.packages(x1[!x2],repos="http://cran.us.r-project.org") }
# Load all packages (unless already loaded) # as.Date <- zoo::as.Date
lapply(x1,library,character.only=TRUE)
select <- dplyr::select; # row_number <- dplyr::row_number; summarise <- dplyr::summarise
source("fcns/essential_fcns.R")
# retired functions in `source('fcns/RSV_model_functions.R')`
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###
# set plotting theme
standard_theme=theme(plot.title=element_text(hjust=0.5,size=16),
axis.text.x=element_text(size=9,angle=90,vjust=1/2),
axis.text.y=element_text(size=9),axis.title=element_text(size=14),
panel.grid.minor=element_blank())
# panel.grid=element_line(linetype="solid",colour="black",size=0.1),
# text=element_text(family="Calibri") # some R versions/OS will throw errors if font set to Calibri
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###
# SET PARAMETERS --------------------------------------------------------
### ### ### ### ### ### ### ###
# constant parameters
# selected country
country_sel="United Kingdom"
# time resolution (in days)
elem_time_step=1
# population data
standard_age_groups <- fun_cntr_agestr(country_sel,i_year="2020",seq(0,75,5),c(seq(4,74,5),99))
popul_struct=fcn_cntr_fullpop(n_year="2020",country_sel)
# RSV age groups (population data from wpp2019)
rsv_age_groups <- fun_rsv_agegroups(standard_age_groups,popul_struct,
rsv_age_groups_low=c(0,0.5,1,1.5, 2,3,4, 5,15, 45, 65),
rsv_age_group_sizes=c(rep(0.4,4),rep(0.9,3), 9, 29, 19, 34))
# rsv_age_groups$value=rsv_age_groups$value*67e6/sum(rsv_age_groups$value)
ons_2020_midyear_estimates_uk <- read_csv(here::here("repo_data/ons_2020_midyear_estimates_uk.csv")) %>%
mutate(age_num=as.numeric(gsub("\\+","",age)))
low_inds<-findInterval(rsv_age_groups$age_low,ons_2020_midyear_estimates_uk$age_num)
high_inds <- findInterval(rsv_age_groups$age_low+rsv_age_groups$duration-0.1,
ons_2020_midyear_estimates_uk$age_num)
rsv_age_groups$value <- unlist(lapply(1:length(low_inds),
function(x) sum(ons_2020_midyear_estimates_uk$value[low_inds[x]:high_inds[x]])*ifelse(
rsv_age_groups$duration[x]<1,rsv_age_groups$duration[x],1) ))
# DEATHS (2019: 530841 deaths [England and Wales!])
# "uk_death_rate_byage_rsv_agegroups.csv" is for 1000 population!
# read_csv("data/uk_death_rate_byage_rsv_agegroups.csv")
# I adjusted the age-specific deaths rates to get a stationary population close to the 2019 age struct
uk_death_rate=c(rep(1e-5,2)*3,rep(1e-6,5),rep(0,2),1e-6,1.79e-4)
# number of age groups, reinfections and variables (S,I,R)
n_age=nrow(rsv_age_groups); varname_list=c('S','I','R'); n_compartment=length(varname_list); n_inf=3
dim_sys=n_age*n_compartment*n_inf; n_days_year=365
# BIRTH RATE (births into S_1_1)
daily_births=2314; birth_rates=matrix(c(daily_births,rep(0,dim_sys-1)),dim_sys,1)
# we want population to be stationary (at 2019 or 2020 value), so deaths = births
if (!any(grepl("death",colnames(rsv_age_groups)))){
rsv_age_groups <- rsv_age_groups %>%
mutate(deaths_per_person_per_day=uk_death_rate,
stationary_popul=fcn_calc_stat_popul(rsv_age_groups,rsv_age_groups$duration,
daily_births,uk_death_rate,output_type=""))
}
# query variables: fun_sub2ind(1:3,11,"R",varname_list,n_age,n_inf)
# force of infection terms
# linear indices of the I & S variables
l_inf_susc=fun_inf_susc_index_lists(n_age,n_inf,varname_list); inf_vars_inds=l_inf_susc[[1]]
susc_vars_inds=l_inf_susc[[2]]
# CONTACT MATRIX
# contact matrix from covidm ("home","work","school","other")
# cm_path="~/Desktop/research/models/epid_models/covid_model/lmic_model/covidm/"
# if UK -> England's contact matrix # check: cm_parameters_SEI3R(cm_uk_locations("UK", 1))$pop[[1]]$matrices
# list_contmatrs=fun_covidm_contactmatrix(country_sel,currentdir_path,cm_path=cm_path)
# make matrix reciprocal
# C_m_polymod=Reduce('+',list_contmatrs) # fun_recipr_contmatr(Reduce('+',list_contmatrs),
# age_group_sizes=standard_age_groups$values)
C_m_polymod<-readRDS(here::here("repo_data/UK_contact_matrix_sum.RDS"))
# modify contact matrix to correspond to our age groups
C_m_merged_nonrecipr=fun_create_red_C_m(C_m_polymod,rsv_age_groups,
orig_age_groups_duration=standard_age_groups$duration,orig_age_groups_sizes=standard_age_groups$values)
# make it reciprocal for the larger group
C_m=fun_recipr_contmatr(C_m_merged_nonrecipr,age_group_sizes=rsv_age_groups$stationary_popul)
# bc of reinfections we need to input contact matrix repeatedly:
# the force of infection is the same for the 1st, 2nd, 3rd infections,
# but in the full equation they are multiplied by a different
# susceptibility parameter (delta1>delta2>delta3).
# Therefore we repeat the contact matrix's rows three times,
# corresponding to 1st, 2nd, 3rd infxs of each age group
# Also the columns of the matrix are normalised by the age group sizes,
# as they will multiply the infection terms by age groups
contmatr_rowvector = t(do.call(cbind,
lapply(1:nrow(C_m), function(x){diag(C_m[x,]) %*% matrix(1,n_age,n_inf)}))) /
matrix(rep(rsv_age_groups$stationary_popul,n_age*n_inf),nrow=n_age*n_inf,byrow=T)
# matrix(1,nrow=33,ncol = 11)/matrix(rep(1:11,33),nrow=33,byrow=T)
every_nth = function(n) { return(function(x) {x[c(TRUE, rep(FALSE, n-1))]}) }
# build kinetic matrix (=linear terms of recovery, immunity waning and aging)
# WANING (immunity) terms: R_i_j -> S_min(i+1,n_inf)_j
omega=1/350 # 1/runif(1,60,200)
# RECOVERY
rho=1/7 # 1/rho=rweibull(1, shape=4.1,scale=8.3)
# KINETIC MATRIX (aging terms need to be scaled by duration of age groups!)
K_m=fun_K_m_sirs_multiage(dim_sys,n_age,n_inf,n_compartment,rho,omega,varname_list,
agegroup_durations=rsv_age_groups$duration)
# agegroup indices for maternal immunity
mat_imm_flag <- TRUE; mat_imm_inds<-list(fun_sub2ind(i_inf=1,j_age=1,"R",c("S","I","R"),n_age,3),
fun_sub2ind(i_inf=c(1,2,3),j_age=9,"R",c("S","I","R"),n_age,3),
fun_sub2ind(i_inf=c(1,2,3),j_age=9,"S",c("S","I","R"),n_age,3))
# load estimated attack rates
message("load attack rate estimates")
attack_rate_tol_factor = 2
estim_attack_rates <- read_csv(here("repo_data/kenya_attack_rates.csv")) %>%
mutate(min_est=1e2*sympt_attack_rate/attack_rate_tol_factor,
max_est=1e2*attack_rate_tol_factor*sympt_attack_rate)
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###
# load hospitalisation data: `hosp_probabilities` contains
# hospitalisation per infection probabilities to convert cases to hospitalisations
# the source for these estimates is Hodgson 2021
message("load hospitalisation per infection estimates")
hosp_probabilities <- read_csv("repo_data/hosp_probabilities.csv")
# aggregate into our age groups
hosp_probabilities_broad_age =
hosp_probabilities %>% mutate(size=rsv_age_groups$value,
broad_age=findInterval(agegroup,c(2,4,7,10)+1)+1) %>%
group_by(broad_age) %>%
# this is not entirely correct, because it merges some age groups
summarise(prob_hosp_per_infection=sum(prob_hosp_per_infection*size/sum(size)),
agegroup=unique(broad_age) ) %>%
select(agegroup,prob_hosp_per_infection)
# for the manuscript we need larger fonts
manuscript_large_font_theme <- theme(
strip.text=element_text(size=18),legend.title=element_text(size=16),
legend.text=element_text(size=15),legend.position="top",axis.text.x=element_text(size=14),
axis.text.y=element_text(size=12),axis.title.x=element_text(size=16),
axis.title.y=element_text(size=18))
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###
message("end of load_params.R")