Skip to content

Commit fa7df73

Browse files
robj411robj411
committed
merging
1 parent a1b77c1 commit fa7df73

2 files changed

Lines changed: 24 additions & 25 deletions

File tree

‎code/scenarios/accra/master_script_rj.R‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -85,12 +85,14 @@ ithim_object <- run_ithim_setup(NSAMPLES = 64,
8585
MMET_CYCLING = c(log(5), log(1.2)),
8686
PM_CONC_BASE = c(log(50), log(1.2)),
8787
PM_TRANS_SHARE = c(5, 5),
88+
INJURY_REPORTING_RATE = c(5, 5),
8889
MC_TO_CAR_RATIO = c(-1.4,0.4),
89-
PA_DOSE_RESPONSE_QUANTILE = T,
90-
AP_DOSE_RESPONSE_QUANTILE = T)
90+
PA_DOSE_RESPONSE_QUANTILE = F,
91+
AP_DOSE_RESPONSE_QUANTILE = F)
9192

9293
numcores <- detectCores()
9394
ithim_object$outcomes <- mclapply(1:NSAMPLES, FUN = ithim_uncertainty, ithim_object = ithim_object, mc.cores = ifelse(Sys.info()[['sysname']] == "Windows", 1, numcores))
95+
for(i in 1:NSAMPLES) print(length(ithim_object$outcomes[[i]]))
9496

9597
plot(ithim_object$parameters$MC_TO_CAR_RATIO,sapply(ithim_object$outcomes,function(x)sum(x$hb$deaths[,40])))
9698

‎ithim_r_functions.R‎

Lines changed: 20 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -181,8 +181,6 @@ ithim_setup_parameters <- function(NSAMPLES = 1,
181181
for(age in unique(dr_ap$age_code)){
182182
dr_ap_age <- subset(dr_ap,age_code==age)
183183
#######################################
184-
##RJ I recommend the following as a better approximation to the distribution but it is currently v e r y slow
185-
## so I leave it here commented out until we want to develop it and/or use it
186184
lbeta <- log(dr_ap_age$beta)
187185
lgamma <- log(dr_ap_age$gamma)
188186
gamma_val <- quantile(density(lgamma),quant1)
@@ -393,7 +391,7 @@ edit_accra_trips <- function(raw_trip_set){
393391
total_mc_distance <- total_car_distance*VEHICLE_INVENTORY$distance_ratio_to_car[VEHICLE_INVENTORY$trip_mode==new_mode]
394392
mc_duration <- total_mc_distance/VEHICLE_INVENTORY$speed[VEHICLE_INVENTORY$trip_mode==new_mode]*60
395393
residual_mc_duration <- mc_duration - sum(subset(raw_trip_set,trip_mode==new_mode)$trip_duration)
396-
duration_range <- 15:100
394+
#duration_range <- 15:100
397395
nTrips <- 1
398396
nPeople <- 20#round(residual_mc_duration/nTrips/mean(duration_range))
399397
duration <- residual_mc_duration/nPeople
@@ -636,8 +634,7 @@ create_all_scenarios <- function(trip_set){
636634
target_new_trips <- round(0.1 * tt) - sum(rdr$trip_mode=='Motorcycle')
637635
mcycle_trips_sample <- create_scenario(rdr, scen_name = 'Scenario 3', source_modes = source_modes,
638636
combined_modes = T, target_modes = target_modes,
639-
source_distance_cats = DIST_CAT, source_trips = c(round(0.1 * tt) -
640-
nrow(filter(rdr, trip_mode == 'Motorcycle'))))
637+
source_distance_cats = DIST_CAT, source_trips = target_new_trips)
641638
# Update selected rows for mode and duration
642639
rdr$trip_mode[match(mcycle_trips_sample$row_id,rdr$row_id)] <- mcycle_trips_sample$trip_mode
643640
rdr$trip_duration[match(mcycle_trips_sample$row_id,rdr$row_id)] <- mcycle_trips_sample$trip_duration
@@ -904,8 +901,19 @@ distances_for_injury_function <- function(trip_scen_sets){
904901
}
905902
}
906903

907-
908-
list(relative_distances=relative_distances,scen_dist=scen_dist,true_distances=true_distances,injuries_list=injuries_list)
904+
##!! move to distance calculation. store predictions and scale by distance.
905+
reg_model <- list()
906+
##TODO write formulae without prior knowledge of column names
907+
##TODO use all ages. ns.
908+
##TODO different formulae for whw and noov
909+
for(type in c('whw','noov')){
910+
injuries_list[[1]][[type]]$injury_reporting_rate <- 1
911+
suppressWarnings(reg_model[[type]] <- glm(count~cas_mode*strike_mode+cas_age+cas_gender,data=injuries_list[[1]][[type]],family='poisson',
912+
offset=0.5*log(cas_distance)+0.5*log(strike_distance)-log(injury_reporting_rate)))
913+
for(scen in SCEN)
914+
injuries_list[[scen]][[type]] <- subset(injuries_list[[scen]][[type]],year==2016)
915+
}
916+
list(relative_distances=relative_distances,scen_dist=scen_dist,true_distances=true_distances,injuries_list=injuries_list,reg_model=reg_model)
909917
}
910918

911919
######################################################################
@@ -1128,11 +1136,8 @@ PA_dose_response <- function (cause, outcome_type, dose, confidence_intervals =
11281136
stop('Incidence does not exist for all_cause')
11291137
}
11301138
fname <- paste(cause, outcome_type, sep = "_")
1131-
lookup_table <- get(paste0(fname))
1139+
lookup_table <- get(fname)
11321140
lookup_df <- setDT(lookup_table)
1133-
#pert_75 <- stringr::str_sub(basename(list_of_files[[1]]), end = -5)
1134-
##RJ previously:
1135-
## cond <- ifelse(use_75_pert, abs(lookup_table$dose - dose), which.min(abs(lookup_table$dose - dose)))
11361141
rr <- approx(x=lookup_df$dose,y=lookup_df$RR,xout=dose,yleft=1,yright=min(lookup_df$RR))$y
11371142
if (confidence_intervals || PA_DOSE_RESPONSE_QUANTILE==T) {
11381143
lb <-
@@ -1154,7 +1159,7 @@ PA_dose_response <- function (cause, outcome_type, dose, confidence_intervals =
11541159
}
11551160
if (PA_DOSE_RESPONSE_QUANTILE==T){
11561161
#rr <- truncnorm::qtruncnorm(get(paste0('PA_DOSE_RESPONSE_QUANTILE_',cause)), rr, sd=rr-lb,a=0, b=1)
1157-
rr <- qnorm(get(paste0('PA_DOSE_RESPONSE_QUANTILE_',cause)), rr, sd=(ub-lb)/1.96)
1162+
rr <- qnorm(get(paste0('PA_DOSE_RESPONSE_QUANTILE_',cause)), mean=rr, sd=(ub-lb)/1.96)
11581163
rr[rr<0] <- 0
11591164
}
11601165
if (confidence_intervals) {
@@ -1186,22 +1191,14 @@ combined_rr_pa_pa <- function(ind_pa,ind_ap){
11861191
ind_ap_pa
11871192
}
11881193

1189-
injuries_function_2 <- function(true_distances,injuries_list){
1190-
##!! move to distance calculation. store predictions and scale by distance.
1191-
reg_model <- list()
1192-
##TODO write formulae without prior knowledge of column names
1193-
##TODO use all ages. ns.
1194-
##TODO different formulae for whw and noov
1195-
for(type in c('whw','noov'))
1196-
reg_model[[type]] <- glm(count~cas_mode+strike_mode+cas_age+cas_gender,data=injuries_list[[1]][[type]],family='poisson',
1197-
offset=0.5*log(cas_distance)+0.5*log(strike_distance))
1194+
injuries_function_2 <- function(true_distances,injuries_list,reg_model){
11981195
##
11991196
## For predictive uncertainty, we could sample a number from the predicted distribution
12001197
injuries <- true_distances
12011198
injuries$Bus_driver <- 0
12021199
for(scen in SCEN){
12031200
for(type in c('whw','noov')){
1204-
injuries_list[[scen]][[type]] <- subset(injuries_list[[scen]][[type]],year==2016)
1201+
injuries_list[[scen]][[type]]$injury_reporting_rate <- INJURY_REPORTING_RATE
12051202
injuries_list[[scen]][[type]]$pred <- predict(reg_model[[type]],newdata = injuries_list[[scen]][[type]],type='response')
12061203
}
12071204
for(injured_mode in unique(injuries_list[[1]]$whw$cas_mode))
@@ -1444,7 +1441,7 @@ run_ithim <- function(ithim_object,seed=1){
14441441
for(i in 1:length(inj_distances))
14451442
assign(names(inj_distances)[i],inj_distances[[i]])
14461443
(injuries <- injuries_function(relative_distances,scen_dist))
1447-
#system.time(injuries <- injuries_function_2(true_distances,injuries_list))
1444+
(injuries <- injuries_function_2(true_distances,injuries_list,reg_model))
14481445
(deaths_yll_injuries <- injury_death_to_yll(injuries))
14491446
ref_injuries <- deaths_yll_injuries$ref_injuries
14501447
############################

0 commit comments

Comments
 (0)