-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathest_R_hosp.R~
153 lines (120 loc) · 5.32 KB
/
est_R_hosp.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
142
143
144
145
146
147
148
149
150
151
152
library(tidyverse)
library(readxl)
library(haven)
library(EpiEstim) # Install from github: mrc-ide/EpiEstim
library(bbmle)
library(cowplot)
library(ggthemes)
library(doParallel)
library(zoo)
source('cleaningFunctions.R')
source('fittingFunctions.R')
source('rtFunctions.R')
source('plottingFunctions.R')
source('imp_functions.R')
.args <- if(interactive){
c(7, #days
FALSE,
TRUE
}else{
commandArgs(trailingOnly =TRUE)
}
R_est_window <- as.numeric(.args[[1]])
analyse_only_deaths <- .args[[2]]
pubpriv_sector <- "ALL"
use_uncertain_SI <- .args[[3]
#Public or Private or ALL
max_onset_date <- as.Date('2022-04-30') # Update with new data
min_onset_date <- as.Date('2020-03-05') # Keep
trim_days_admissions <- 7
trim_days_deaths <- 7
# set max common (between-province) onset dates per sector, for constructing national time series
mcos <- c(as.Date('2022-04-28'), # public
as.Date('2022-04-28'), #private
as.Date('2022-04-28')) - #all
rep(trim_days_admissions, 3)
mcos_d <- c(as.Date('2022-04-20'), #public
as.Date('2022-04-20'), #private
as.Date('2022-04-20')) - #all
rep(trim_days_deaths, 3)
names(mcos) <- c("Public", "Private", "ALL")
names(mcos_d) <- c("Public", "Private", "ALL")
max_combo_onset_date_admit <- mcos[pubpriv_sector]# as.Date('2022-02-18') - trim_days_admissions
max_combo_onset_date_deaths <- mcos_d[pubpriv_sector] #as.Date('2022-02-09') - trim_days_deaths
n1_config = 25
n2_config = 25
use_cores <- 5
save_output <- TRUE
if(analyse_only_deaths){
plot_ylab <- "Lab-confirmed COVID-19 deaths"
file_suffix = sprintf("%s_deaths", str_to_title(pubpriv_sector))
}else{
plot_ylab <- "Lab-confirmed COVID-19 hospitalisations"
file_suffix = sprintf("%s_admit", str_to_title(pubpriv_sector))
}
max_onset2admit <- 21
max_onset2death <- 40
si_mean_gamma <- 6.63
si_sd_gamma <- 3.28
si_mean_std_gamma <- 0.51
si_std_std_gamma <- 0.27
set.seed(20200727) # Note that if running in parallel doesn't gaurantee reproducibility
registerDoParallel(use_cores) # Register cores
provs = c('za', 'gt', 'wc', 'ec', 'kzn', 'nw', 'mp', 'fs','lm', 'nc')
prov_names = c('ZA', 'Gauteng', 'Western Cape', 'Eastern Cape',
'KwaZulu-Natal', 'North West',
'Mpumalanga', 'Free State', 'Limpopo' , 'Northern Cape')
names(provs) = prov_names
for(province_name in prov_names){
prov_code <- provs[province_name]
o2r_fit_tmp <- readRDS(sprintf('./data/o2r_fit_%s_%s.RDS',)
ts_onset_tmp <- adjust_rcens_imputed_ts(
# combine_imp_get_ts_admit_par(n_imputations = n_imputations,
# n_splits = n_split,
# folder = imp_folder,
# filename_base = "admitted_cases_nb",
# prov = province_name,
# sect = pubpriv_sector,
# only_deaths = analyse_only_deaths,
# renew_calculation = recalc_ts)
) %>% filter(onset_date <= max(onset_date) - ifelse(analyse_only_deaths, trim_days_deaths, trim_days_admissions))
# print(max(ts_onset_tmp$onset_date))
# if(prov_code %in% c('gt','kzn')){
# ts_onset_tmp <- ts_onset_tmp %>% filter(onset_date <= as.Date('2021-07-10'))
# }
# ts_plot_tmp <- plot_ts_onset(ts_onset_tmp, qq = use_quant,
# y_lab = plot_ylab)
# rt_est(ts_onset_tmp, n1 = n1_config, n2 = n2_config, window_size = 5)
rt_est_tmp <- wrap_rt_est(ts_onset_tmp, n1 = n1_config, n2 = n2_config, window_size = R_est_window)
# rt_plot_tmp <- plot_rt(rt_est_tmp)
#rt_specific_date(rt_est_tmp, date = as.Date('2021-01-18'))
saveRDS(ts_onset_tmp, file = sprintf("./results/%s_ts_%s.RDS",file_suffix, prov_code))
saveRDS(rt_est_tmp, file = sprintf("./results/%s_rt_%s.RDS", file_suffix, prov_code))
readRDS(file = sprintf("./results/%s_ts_%s.RDS", file_suffix, prov_code))
}
alt_ts_za <- data.frame(onset_date = as.Date(character()), rep = numeric(), adj_cnt = numeric(), cnt = numeric(), prov = character())
#names(alt_ts_za) <- c("onset_date", "rep", "adj_cnt", "prov")
count_combo_provs = 0
for (province_code in provs){
if(province_code == "za") {next}
else{
tmp <-readRDS(sprintf("./results/%s_ts_%s.RDS", file_suffix, province_code)) %>% mutate(prov = province_code)
alt_ts_za <- rbind(alt_ts_za, tmp)
count_combo_provs = count_combo_provs + 1
}
}
if(count_combo_provs != 9){print("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! combo ts incomplete")}
alt_ts_za <- alt_ts_za %>%
group_by(onset_date, rep) %>%
summarise(adj_cnt = sum(adj_cnt), cnt = sum(cnt)) %>%
filter(if(analyse_only_deaths) onset_date <= max_combo_onset_date_deaths else
onset_date <= max_combo_onset_date_admit)
za_rt_alt <- wrap_rt_est(alt_ts_za, n1 = n1_config, n2 = n2_config, window_size = R_est_window)
saveRDS(alt_ts_za, file = sprintf('./results/%s_ts_za_alt.RDS', file_suffix))
saveRDS(za_rt_alt, file = sprintf('./results/%s_rt_za_alt.RDS', file_suffix))
# PERIOD ANALYSIS
for(province_code in provs){
print(province_code)
rt_per_tmp <- wrap_rt_est(readRDS(sprintf("./results/%s_ts_%s.RDS", file_suffix, ifelse(province_code == "za", "za_alt", province_code))), use_periods = T, browse = FALSE, n1 = n1_config, n2 = n2_config, window_size = 7)
saveRDS(rt_per_tmp, file = sprintf("./results/%s/Period_analyses/%s_per_%s.RDS", file_suffix, province_code))
}