diff --git a/results/accra/six_by_five_scenarios_16.Rds b/results/accra/six_by_five_scenarios_16.Rds deleted file mode 100644 index 81e1670a..00000000 Binary files a/results/accra/six_by_five_scenarios_16.Rds and /dev/null differ diff --git a/results/accra/six_by_one_scenarios_1024.Rds b/results/accra/six_by_one_scenarios_1024.Rds deleted file mode 100644 index a8c36d02..00000000 Binary files a/results/accra/six_by_one_scenarios_1024.Rds and /dev/null differ diff --git a/results/accra/six_by_one_scenarios_2048.Rds b/results/accra/six_by_one_scenarios_2048.Rds deleted file mode 100644 index 8359ea29..00000000 Binary files a/results/accra/six_by_one_scenarios_2048.Rds and /dev/null differ diff --git a/results/accra/six_by_one_scenarios_4096.Rds b/results/accra/six_by_one_scenarios_4096.Rds deleted file mode 100644 index 4374c351..00000000 Binary files a/results/accra/six_by_one_scenarios_4096.Rds and /dev/null differ diff --git a/results/accra/six_by_one_scenarios_8192.Rds b/results/accra/six_by_one_scenarios_8192.Rds deleted file mode 100644 index 73f9afe2..00000000 Binary files a/results/accra/six_by_one_scenarios_8192.Rds and /dev/null differ diff --git a/results/dose_response/AP/cvd_ihd25.pdf b/results/dose_response/AP/cvd_ihd25.pdf deleted file mode 100644 index 0808a2c8..00000000 Binary files a/results/dose_response/AP/cvd_ihd25.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_ihd30.pdf b/results/dose_response/AP/cvd_ihd30.pdf deleted file mode 100644 index f54f2081..00000000 Binary files a/results/dose_response/AP/cvd_ihd30.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_ihd35.pdf b/results/dose_response/AP/cvd_ihd35.pdf deleted file mode 100644 index 45ad1792..00000000 Binary files a/results/dose_response/AP/cvd_ihd35.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_ihd40.pdf b/results/dose_response/AP/cvd_ihd40.pdf deleted file mode 100644 index 12552c5d..00000000 Binary files a/results/dose_response/AP/cvd_ihd40.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_ihd45.pdf b/results/dose_response/AP/cvd_ihd45.pdf deleted file mode 100644 index c4191d8c..00000000 Binary files a/results/dose_response/AP/cvd_ihd45.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_ihd50.pdf b/results/dose_response/AP/cvd_ihd50.pdf deleted file mode 100644 index a2a89b28..00000000 Binary files a/results/dose_response/AP/cvd_ihd50.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_ihd55.pdf b/results/dose_response/AP/cvd_ihd55.pdf deleted file mode 100644 index cb0d23ae..00000000 Binary files a/results/dose_response/AP/cvd_ihd55.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_ihd60.pdf b/results/dose_response/AP/cvd_ihd60.pdf deleted file mode 100644 index e13b9611..00000000 Binary files a/results/dose_response/AP/cvd_ihd60.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_ihd65.pdf b/results/dose_response/AP/cvd_ihd65.pdf deleted file mode 100644 index a44a9ef0..00000000 Binary files a/results/dose_response/AP/cvd_ihd65.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_ihd70.pdf b/results/dose_response/AP/cvd_ihd70.pdf deleted file mode 100644 index 15f8d297..00000000 Binary files a/results/dose_response/AP/cvd_ihd70.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_ihd75.pdf b/results/dose_response/AP/cvd_ihd75.pdf deleted file mode 100644 index e3682508..00000000 Binary files a/results/dose_response/AP/cvd_ihd75.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_ihd80.pdf b/results/dose_response/AP/cvd_ihd80.pdf deleted file mode 100644 index 81cf6a6a..00000000 Binary files a/results/dose_response/AP/cvd_ihd80.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_ihd85.pdf b/results/dose_response/AP/cvd_ihd85.pdf deleted file mode 100644 index 2237d90b..00000000 Binary files a/results/dose_response/AP/cvd_ihd85.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_ihd90.pdf b/results/dose_response/AP/cvd_ihd90.pdf deleted file mode 100644 index 9fbe5615..00000000 Binary files a/results/dose_response/AP/cvd_ihd90.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_ihd95.pdf b/results/dose_response/AP/cvd_ihd95.pdf deleted file mode 100644 index 77546860..00000000 Binary files a/results/dose_response/AP/cvd_ihd95.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_ihdDR-sim.pdf b/results/dose_response/AP/cvd_ihdDR-sim.pdf deleted file mode 100644 index 059e3c0e..00000000 Binary files a/results/dose_response/AP/cvd_ihdDR-sim.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_ihdDR.pdf b/results/dose_response/AP/cvd_ihdDR.pdf deleted file mode 100644 index 7b866080..00000000 Binary files a/results/dose_response/AP/cvd_ihdDR.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_stroke25.pdf b/results/dose_response/AP/cvd_stroke25.pdf deleted file mode 100644 index 08e0fc8a..00000000 Binary files a/results/dose_response/AP/cvd_stroke25.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_stroke30.pdf b/results/dose_response/AP/cvd_stroke30.pdf deleted file mode 100644 index 6d850c04..00000000 Binary files a/results/dose_response/AP/cvd_stroke30.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_stroke35.pdf b/results/dose_response/AP/cvd_stroke35.pdf deleted file mode 100644 index 4750d864..00000000 Binary files a/results/dose_response/AP/cvd_stroke35.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_stroke40.pdf b/results/dose_response/AP/cvd_stroke40.pdf deleted file mode 100644 index 2952476f..00000000 Binary files a/results/dose_response/AP/cvd_stroke40.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_stroke45.pdf b/results/dose_response/AP/cvd_stroke45.pdf deleted file mode 100644 index bce7f175..00000000 Binary files a/results/dose_response/AP/cvd_stroke45.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_stroke50.pdf b/results/dose_response/AP/cvd_stroke50.pdf deleted file mode 100644 index 422aebd5..00000000 Binary files a/results/dose_response/AP/cvd_stroke50.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_stroke55.pdf b/results/dose_response/AP/cvd_stroke55.pdf deleted file mode 100644 index 64aa74a8..00000000 Binary files a/results/dose_response/AP/cvd_stroke55.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_stroke60.pdf b/results/dose_response/AP/cvd_stroke60.pdf deleted file mode 100644 index 23939e26..00000000 Binary files a/results/dose_response/AP/cvd_stroke60.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_stroke65.pdf b/results/dose_response/AP/cvd_stroke65.pdf deleted file mode 100644 index 04f905f8..00000000 Binary files a/results/dose_response/AP/cvd_stroke65.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_stroke70.pdf b/results/dose_response/AP/cvd_stroke70.pdf deleted file mode 100644 index 7623f807..00000000 Binary files a/results/dose_response/AP/cvd_stroke70.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_stroke75.pdf b/results/dose_response/AP/cvd_stroke75.pdf deleted file mode 100644 index cc55af9f..00000000 Binary files a/results/dose_response/AP/cvd_stroke75.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_stroke80.pdf b/results/dose_response/AP/cvd_stroke80.pdf deleted file mode 100644 index 9ce165a1..00000000 Binary files a/results/dose_response/AP/cvd_stroke80.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_stroke85.pdf b/results/dose_response/AP/cvd_stroke85.pdf deleted file mode 100644 index 891ad53e..00000000 Binary files a/results/dose_response/AP/cvd_stroke85.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_stroke90.pdf b/results/dose_response/AP/cvd_stroke90.pdf deleted file mode 100644 index 739bc577..00000000 Binary files a/results/dose_response/AP/cvd_stroke90.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_stroke95.pdf b/results/dose_response/AP/cvd_stroke95.pdf deleted file mode 100644 index e639b233..00000000 Binary files a/results/dose_response/AP/cvd_stroke95.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_strokeDR-sim.pdf b/results/dose_response/AP/cvd_strokeDR-sim.pdf deleted file mode 100644 index b9cd15d3..00000000 Binary files a/results/dose_response/AP/cvd_strokeDR-sim.pdf and /dev/null differ diff --git a/results/dose_response/AP/cvd_strokeDR.pdf b/results/dose_response/AP/cvd_strokeDR.pdf deleted file mode 100644 index cbd7c8be..00000000 Binary files a/results/dose_response/AP/cvd_strokeDR.pdf and /dev/null differ diff --git a/results/dose_response/AP/lri99.pdf b/results/dose_response/AP/lri99.pdf deleted file mode 100644 index 9a96f41f..00000000 Binary files a/results/dose_response/AP/lri99.pdf and /dev/null differ diff --git a/results/dose_response/AP/lriDR-sim.pdf b/results/dose_response/AP/lriDR-sim.pdf deleted file mode 100644 index a0f96612..00000000 Binary files a/results/dose_response/AP/lriDR-sim.pdf and /dev/null differ diff --git a/results/dose_response/AP/lriDR.pdf b/results/dose_response/AP/lriDR.pdf deleted file mode 100644 index 9ce0b69c..00000000 Binary files a/results/dose_response/AP/lriDR.pdf and /dev/null differ diff --git a/results/dose_response/AP/neo_lung99.pdf b/results/dose_response/AP/neo_lung99.pdf deleted file mode 100644 index 67cca8d3..00000000 Binary files a/results/dose_response/AP/neo_lung99.pdf and /dev/null differ diff --git a/results/dose_response/AP/neo_lungDR-sim.pdf b/results/dose_response/AP/neo_lungDR-sim.pdf deleted file mode 100644 index 54ae95d0..00000000 Binary files a/results/dose_response/AP/neo_lungDR-sim.pdf and /dev/null differ diff --git a/results/dose_response/AP/neo_lungDR.pdf b/results/dose_response/AP/neo_lungDR.pdf deleted file mode 100644 index abd40a84..00000000 Binary files a/results/dose_response/AP/neo_lungDR.pdf and /dev/null differ diff --git a/results/dose_response/AP/resp_copd99.pdf b/results/dose_response/AP/resp_copd99.pdf deleted file mode 100644 index 2299032b..00000000 Binary files a/results/dose_response/AP/resp_copd99.pdf and /dev/null differ diff --git a/results/dose_response/AP/resp_copdDR-sim.pdf b/results/dose_response/AP/resp_copdDR-sim.pdf deleted file mode 100644 index 9dca736f..00000000 Binary files a/results/dose_response/AP/resp_copdDR-sim.pdf and /dev/null differ diff --git a/results/dose_response/AP/resp_copdDR.pdf b/results/dose_response/AP/resp_copdDR.pdf deleted file mode 100644 index aed0ca6e..00000000 Binary files a/results/dose_response/AP/resp_copdDR.pdf and /dev/null differ diff --git a/results/dose_response/PA/PA_dose_response.pdf b/results/dose_response/PA/PA_dose_response.pdf deleted file mode 100644 index 6e3a5180..00000000 Binary files a/results/dose_response/PA/PA_dose_response.pdf and /dev/null differ diff --git a/results/dose_response/PA/PA_dose_response_sample.pdf b/results/dose_response/PA/PA_dose_response_sample.pdf deleted file mode 100644 index af51daa5..00000000 Binary files a/results/dose_response/PA/PA_dose_response_sample.pdf and /dev/null differ diff --git a/results/injury_matrices/bangalore_whw.csv b/results/injury_matrices/bangalore_whw.csv deleted file mode 100644 index 519f5351..00000000 --- a/results/injury_matrices/bangalore_whw.csv +++ /dev/null @@ -1,12 +0,0 @@ -"","X","No.other.fixed.or.stationary.object","Pedestrian","Pedal.cycle","X2.3.Wheeled","Car.pick.up.van","Bus","Railway.train.railway.vehicle","non.motor.vehicle","unspecified","Non.collision","Unknown","Trucks" -"1","Pedestrian",0,0,0,119,87,39,0,0,7,0,59,62 -"2","Pedal Cycle",0,0,0,3,7,1,0,0,3,0,3,4 -"3","Motorcycle",29,3,1,43,47,59,0,0,10,0,8,102 -"4","3Wheeled",3,0,0,0,4,3,0,0,0,0,0,8 -"5","Car",9,0,0,3,1,2,0,0,0,0,2,7 -"6","Pick-up truck/van",0,0,0,0,0,2,0,0,0,0,0,0 -"7","Heavy transport",3,0,0,0,0,0,0,0,0,0,0,8 -"8","Bus",6,0,0,0,0,8,0,0,0,0,0,4 -"9","Other",0,0,0,0,0,0,0,0,9,0,3,0 -"10","Unknown",0,0,0,0,0,0,0,0,0,0,0,0 -"11","Railway",0,0,0,0,0,0,0,0,0,0,0,0 diff --git a/results/injury_matrices/delhi_bangalore_saopaulo_accra_whw_injury.xlsx b/results/injury_matrices/delhi_bangalore_saopaulo_accra_whw_injury.xlsx deleted file mode 100644 index fc356be0..00000000 Binary files a/results/injury_matrices/delhi_bangalore_saopaulo_accra_whw_injury.xlsx and /dev/null differ diff --git a/results/injury_matrices/delhi_whw.csv b/results/injury_matrices/delhi_whw.csv deleted file mode 100644 index 2460105d..00000000 --- a/results/injury_matrices/delhi_whw.csv +++ /dev/null @@ -1,12 +0,0 @@ -"","X.1","X","No.other.fixed.or.stationary.object","Pedestrian","Pedal.cycle","X2.3.Wheeled","Car.pick.up.van","Bus","Railway.train.railway.vehicle","non.motor.vehicle","unspecified","Non.collision","Unknown","Trucks" -"1",1,"Pedestrian",0,0,0,82,200,61,0,0,0,0,427,95 -"2",2,"Pedal Cycle",0,0,0,8,35,15,0,0,0,0,30,25 -"3",3,"Motorcycle",57,0,0,25,141,60,0,0,0,0,260,139 -"4",4,"3Wheeled",6,0,0,2,13,3,0,0,0,0,10,8 -"5",5,"Car",32,0,0,2,21,4,0,0,0,0,17,21 -"6",6,"Pick-up truck/van",9,0,0,0,6,2,0,0,0,0,4,4 -"7",7,"Heavy transport",15,0,0,1,11,1,0,0,0,0,5,14 -"8",8,"Bus",3,0,0,1,2,2,0,0,0,0,0,1 -"9",9,"Other",0,0,0,0,0,0,0,0,0,0,41,0 -"10",10,"Unknown",0,0,0,0,0,0,0,0,0,0,21,0 -"11",11,"Railway",0,0,0,0,0,0,0,0,0,0,0,0 diff --git a/results/injury_matrices/who_hit_who_accra.csv b/results/injury_matrices/who_hit_who_accra.csv deleted file mode 100644 index 5f8b0e0c..00000000 --- a/results/injury_matrices/who_hit_who_accra.csv +++ /dev/null @@ -1,12 +0,0 @@ -"","No other/fixed or stationary object","Pedestrian","Pedal cycle","2/3 Wheeled","Car/pick-up/van","Bus","Railway train/railway vehicle","non-motor vehicle","unspecified","Non-collision","Unknown","Trucks" -"Pedestrian",1,0,2,63,434,218,0,0,52,0,35,76 -"Pedal Cycle",0,0,0,2,24,11,0,0,1,0,2,10 -"Motorcycle",20,6,0,9,56,29,0,0,4,0,1,18 -"3Wheeled",0,0,0,0,0,0,0,0,0,0,0,0 -"Car",47,3,0,4,29,25,0,0,5,0,0,22 -"Pick-up truck/van",8,0,0,0,3,4,0,0,0,0,0,5 -"Heavy transport",22,0,0,0,3,4,0,0,0,0,0,4 -"Bus",33,1,1,1,8,14,0,0,3,0,0,13 -"Other",0,0,0,0,0,0,0,0,7,0,0,0 -"Unknown",0,0,0,0,0,0,0,0,0,0,0,0 -"Railway",0,0,0,0,0,0,0,0,0,0,0,0 diff --git a/results/injury_matrices/who_hit_who_sao_paulo.csv b/results/injury_matrices/who_hit_who_sao_paulo.csv deleted file mode 100644 index d3967027..00000000 --- a/results/injury_matrices/who_hit_who_sao_paulo.csv +++ /dev/null @@ -1,12 +0,0 @@ -"","No other/fixed or stationary object","Pedestrian","Pedal cycle","2/3 Wheeled","Car/pick-up/van","Bus","Railway train/railway vehicle","non-motor vehicle","unspecified","Non-collision","Unknown","Trucks" -"Pedestrian",14,0,4,459,1036,494,0,0,3,0,333,269 -"Pedal Cycle",11,0,48,12,60,57,0,0,0,0,0,36 -"Motorcycle",374,0,1,493,541,358,0,0,3,0,0,316 -"3Wheeled",0,0,0,0,0,0,0,0,0,0,0,0 -"Car",397,0,0,16,259,104,0,0,0,0,0,72 -"Pick-up truck/van",9,0,0,0,7,3,0,0,0,0,0,18 -"Heavy transport",21,0,0,16,19,28,0,0,0,0,0,40 -"Bus",3,0,0,1,1,3,0,0,0,0,0,0 -"Other",0,0,0,0,0,0,0,0,9,0,0,0 -"Unknown",0,0,0,0,0,0,0,0,0,0,0,0 -"Railway",0,0,0,0,0,0,0,0,0,0,0,0 diff --git a/tests/debugging_deletethisafteruse.R b/tests/debugging_deletethisafteruse.R deleted file mode 100644 index 9ada0d16..00000000 --- a/tests/debugging_deletethisafteruse.R +++ /dev/null @@ -1,58 +0,0 @@ -# predict(reg_model[[type]], -# newdata = remove_missing_levels(reg_model[[type]],injuries_list[[scen]][[type]]),type='response') -# Accra -# 1343 rows -# 1113 rows between 15 y 69 -# 1076 rows remove truck and other -# 928 whw -# 148 nov -# 966 rows after add_distance_column function -# 835 whw -# 131 nov -sum(injuries_for_model[[1]][[type]]$count) -summary(injuries_for_model[[1]][[type]]$count) -hist(injuries_for_model[[1]][[type]]$count) -table(injuries_for_model[[1]][[type]]$count) - -aux = test -aux_predict = predict(aux, newdata = injuries_for_model[[1]][[type]], - type = "response") - -sum(aux_predict) # Same number of trips. 835 in total -summary(aux_predict) -hist(aux_predict) -table(aux_predict) - -aux_predict_missinglabels = predict(aux, - newdata = remove_missing_levels(aux, -injuries_for_model[[1]][[type]]), - type = "response") - -sum(aux_predict_missinglabels) # Same number of trips. 835 in total -summary(aux_predict_missinglabels) -hist(aux_predict_missinglabels) -table(aux_predict_missinglabels) - -aux2 = trim_glm_object(test) -aux2_predict = predict(aux2, newdata = injuries_for_model[[1]][[type]], - type = "response") -sum(aux2_predict) # Same number of trips. 835 in total -summary(aux2_predict) -hist(aux2_predict) -table(aux2_predict) - -# Predict from injuries_list instead of injuries_for model -new_data = injuries_list$Baseline$whw %>% - mutate(injury_reporting_rate = 1, - weight = 10) # set weight to 10 years -aux_pred_list <- predict(aux, newdata = new_data, - type = "response") -sum(aux_pred_list) # Same number of trips. 835 in total - - -results <- data.frame(true_values = injuries_for_model[[1]][[type]]$count, - predicted = round(aux_predict,2) , - predicted_missing = round(aux_predict_missinglabels,2), - predicted2 = round(aux2_predict,2), - predicted_weigth = round(aux_pred_list, 2)) -View(results) diff --git a/tests/plot_city_distributions.R b/tests/plot_city_distributions.R deleted file mode 100644 index 70793867..00000000 --- a/tests/plot_city_distributions.R +++ /dev/null @@ -1,142 +0,0 @@ -lnorm_parameters <- list( - # lnorm parameters for CHRONIC_DISEASE_SCALAR - chronic_disease_scalar = list(accra=c(0,log(1.2)), - sao_paulo=c(0,log(1.2)), - delhi=c(0,log(1.2)), - bangalore=c(0,log(1.2))), - # lnorm parameters for PM_CONC_BASE - pm_concentration = list(accra=c(log(50),log(1.3)), - sao_paulo=c(log(20),log(1.3)), - delhi=c(log(122),log(1.3)), - bangalore=c(log(63),log(1.3))), - # lnorm parameters for BACKGROUND_PA_SCALAR - background_pa_scalar = list(accra=c(0,log(1.2)), - sao_paulo=c(0,log(1.2)), - delhi=c(0,log(1.2)), - bangalore=c(0,log(1.2))), - # lnorm parameters for BUS_WALK_TIME - bus_walk_time = list(accra=c(log(5),log(1.2)), - sao_paulo=c(log(5),log(1.2)), - delhi=c(log(5),log(1.2)), - bangalore=c(log(5),log(1.2))), - # lnorm parameters for MOTORCYCLE_TO_CAR_RATIO - mc_car_ratio = list(accra=c(-1.4,0.4), - sao_paulo=c(-1.4,0.4), - delhi=c(-1.4,0.4), - bangalore=c(-1.4,0.4)), - # lnorm parameters for TRUCK_TO_CAR_RATIO - truck_car_ratio = list(accra=c(-1.4,0.4), - sao_paulo=c(-1.4,0.4), - delhi=c(-1.4,0.4), - bangalore=c(-1.4,0.4)) -) - - -beta_parameters <- list( - # beta parameters for INJURY_REPORTING_RATE - injury_report_rate = list(accra=c(8,3), - sao_paulo=c(50,3), - delhi=c(50,3), - bangalore=c(50,3)), - # beta parameters for PM_TRANS_SHARE - pm_trans_share = list(accra=c(5,20), - sao_paulo=c(8,8), - delhi=c(8,8), - bangalore=c(8,8)), - # bus occupancy beta distribution - bus_to_passenger_ratio = list(accra=c(20,600), - sao_paulo=c(20,600), - delhi=c(20,600), - bangalore=c(20,600)), - # truck beta distribution - truck_to_car_ratio <- list(accra=c(3,10), - sao_paulo=c(3,10), - delhi=c(3,10), - bangalore=c(3,10)) -) - - -cols <- c('navyblue','hotpink','grey','darkorange') -{x11(height=6,width=6); par(mfrow=c(2,2),mar=c(5,5,2,2)) - x <- seq(0,1,length=100) - for(i in 1:length(beta_parameters)) { - toplot <- list() - for(j in 1:length(beta_parameters[[1]])) - toplot[[j]] <- dbeta(x,beta_parameters[[i]][[j]][1],beta_parameters[[i]][[j]][2]) - ymax <- max(unlist(toplot)) - for(j in 1:length(beta_parameters[[1]])) - if(j==1){ - plot(x,toplot[[j]],typ='l',lwd=2,col=cols[j],frame=F,xlab=names(beta_parameters)[i],ylab='Density',ylim=c(0,ymax)) - }else{ - lines(x,toplot[[j]],beta_parameters[[i]][[j]][2],typ='l',lwd=2,col=cols[j]) - } - legend(legend=names(beta_parameters[[1]]),x=0,y=ymax,col=cols,lwd=2,bty='n') - } -} -xmaxes <- c(3,200,3,10,1,1) -{x11(height=6,width=9); par(mfrow=c(2,3),mar=c(5,5,2,2)) - for(i in 1:length(lnorm_parameters)) { - x <- seq(0,xmaxes[i],length=100) - toplot <- list() - for(j in 1:length(lnorm_parameters[[1]])) - toplot[[j]] <- dlnorm(x,lnorm_parameters[[i]][[j]][1],lnorm_parameters[[i]][[j]][2]) - ymax <- max(unlist(toplot)) - for(j in 1:length(lnorm_parameters[[1]])) - if(j==1){ - plot(x,toplot[[j]],typ='l',lwd=2,col=cols[j],frame=F,xlab=names(lnorm_parameters)[i],ylab='Density',ylim=c(0,ymax),xlim=c(0,xmaxes[i])) - }else{ - lines(x,toplot[[j]],lnorm_parameters[[i]][[j]][2],typ='l',lwd=2,col=cols[j]) - } - legend(legend=names(lnorm_parameters[[1]]),x=0.5*xmaxes[i],y=ymax,col=cols,lwd=2,bty='n') - } -} - - -emission_parameters <- list(accra=list( - bus_driver=0.82, - car=0.228, - taxi=0.011, - motorcycle=0.011, - truck=0.859, - big_truck=0.711, - other=0.082 -), -sao_paulo=list(motorcycle=4, - car=4, - bus_driver=32, - big_truck=56, - truck=4), -delhi=list(motorcycle=1409, - auto_rickshaw=133, - car=2214, - bus_driver=644, - big_truck=4624, - truck=3337), -bangalore=list(motorcycle=1757, - auto_rickshaw=220, - car=4173, - bus_driver=1255, - big_truck=4455, - truck=703)) - -library(distr) -confidences <- c(0.9,0.5,0.2) -confidences <- (confidences - 0.6)*5 + confidences -n <- 1000 -for(i in 1:length(emission_parameters)){ - total <- sum(unlist(emission_parameters[[i]])) - d1 <- list() - for(k in 1:length(confidences)){ - confidence <- confidences[k] - distributions <- lapply(emission_parameters[[i]],function(x) Gammad(shape=x/total*100*(3^confidence),scale=1)) - samples <- sapply(distributions,function(x) r(x)(n)) - d1[[k]] <- apply(samples,2,function(x)density(x/rowSums(samples),from=0,to=1)) - } - x11(height=4,width=8); par(mfrow=c(2,4),mar=c(5,5,2,2)) - for(j in 1:length(d1[[1]])) - for(k in 1:length(d1)) - if(k==1) - plot(d1[[k]][[j]],frame=F,xlab=names(emission_parameters[[i]])[j],main='',ylab='Density',lwd=2,col=cols[k],ylim=c(0,max(d1[[k]][[j]]$y)/2)) - else lines(d1[[k]][[j]],col=cols[k],lwd=2) -} - diff --git a/tests/test_ap_dose_response.R b/tests/test_ap_dose_response.R deleted file mode 100644 index d3130cc2..00000000 --- a/tests/test_ap_dose_response.R +++ /dev/null @@ -1,122 +0,0 @@ -rm(list=ls()) -library(ithimr) - -global_path <- file.path(find.package('ithimr',lib.loc=.libPaths()), 'extdata/global/') -global_path <- paste0(global_path, "/") - -##reading GBD 2017 IER functions that were provided by Rick Burnett: this include Diabetes in addition to previous five disease end-points -DR_AP <- read.csv(paste0(global_path,"dose_response/drap/dose_response_2017.csv")) - -DISEASE_INVENTORY <- read.csv(paste0(global_path,"dose_response/disease_outcomes_lookup.csv")) -nsamples <- 100 -parameters <- ithim_setup_parameters(NSAMPLES=nsamples, - AP_DOSE_RESPONSE_QUANTILE=T) -for(i in 1:length(parameters)) - assign(names(parameters)[i],parameters[[i]],pos=1) -causes <- unique(DR_AP$cause_code) - -############## plot parameter samples -dr_ap_list <- list() -for ( j in c(1:nrow(DISEASE_INVENTORY))[DISEASE_INVENTORY$air_pollution == 1]){ - cause <- as.character(DISEASE_INVENTORY$ap_acronym[j]) - dr_ap <- subset(DR_AP,cause_code==cause) - ages <- unique(dr_ap$age_code) - dr_ap_list[[cause]] <- list() - print(cause) - print(ages) - for(age in ages){ - dr_ap_list[[cause]][[age]] <- list() - dr_ap_age <- subset(dr_ap,age_code==age) - alpha <- sapply(DR_AP_LIST,function(x)x[[cause]][[as.character(age)]]$alpha) - beta <- sapply(DR_AP_LIST,function(x)x[[cause]][[as.character(age)]]$beta) - gamma <- sapply(DR_AP_LIST,function(x)x[[cause]][[as.character(age)]]$gamma) - tmrel <- sapply(DR_AP_LIST,function(x)x[[cause]][[as.character(age)]]$tmrel) - dr_ap_list[[cause]][[age]]$alpha <- alpha - dr_ap_list[[cause]][[age]]$beta <- beta - dr_ap_list[[cause]][[age]]$gamma <- gamma - dr_ap_list[[cause]][[age]]$tmrel <- tmrel - #pdf(paste0('results/dose_response/AP/',cause,age,'.pdf')) - #x11(width=8); - par(mfrow=c(2,3),mar=c(4,4,1,1)) - plot(log(dr_ap_age$alpha),log(dr_ap_age$beta),main=paste0(cause,age),col='grey'); points(log(alpha),log(beta),col='hotpink',pch=18) - plot(log(dr_ap_age$alpha),(log(dr_ap_age$gamma)),main=paste0(cause,age),col='grey'); points(log(alpha),log(gamma),col='hotpink',pch=18) - plot(log(dr_ap_age$alpha),log(dr_ap_age$tmrel),main=paste0(cause,age),col='grey'); points(log(alpha),log(tmrel),col='hotpink',pch=18) - plot(log(dr_ap_age$beta),log(dr_ap_age$gamma),main=paste0(cause,age),col='grey'); points(log(beta),log(gamma),col='hotpink',pch=18) - plot(log(dr_ap_age$beta),log(dr_ap_age$tmrel),main=paste0(cause,age),col='grey'); points(log(beta),log(tmrel),col='hotpink',pch=18) - plot(log(dr_ap_age$gamma),log(dr_ap_age$tmrel),main=paste0(cause,age),col='grey'); points(log(gamma),log(tmrel),col='hotpink',pch=18) - #dev.off() - } -} - -# Save dr AP -write_rds(dr_ap_list, "inst/extdata/global/dose_response/drap/dr_ap_list.Rds") - - -######## plot curves -pm <- 1:300 -for ( j in 1:nrow(DISEASE_INVENTORY)) if (DISEASE_INVENTORY$air_pollution[j] == 1){ - cause <- as.character(DISEASE_INVENTORY$ap_acronym[j]) - dr_ap <- subset(DR_AP,cause_code==cause) - ages <- unique(dr_ap$age_code) - print(cause) - print(ages) - if(length(ages)>1){ - pdf(paste0('results/dose_response/AP/',cause,'DR-sim.pdf'),width=6,height=8); par(mfrow=c(1,1)) - layout.matrix <- matrix(1:length(ages),nrow=5,ncol=3,byrow=T) - graphics::layout(mat=layout.matrix,heights=c(4,3,3,3,4),widths=c(3.8,3,3)) - }else{ - pdf(paste0('results/dose_response/AP/',cause,'DR-sim.pdf')) - par(mfrow=c(1,1)) - } - for(age in ages){ - par(mar=c(ifelse(age>80,5,1),ifelse(age%in%c(99,25,40,55,70,85),5,1),ifelse(age%in%c(25:35,99),5,1),1)) - dr_ap_age <- subset(dr_ap,age_code==age) - alpha <- dr_ap_list[[cause]][[age]]$alpha - beta <- dr_ap_list[[cause]][[age]]$beta - gamma <- dr_ap_list[[cause]][[age]]$gamma - tmrel <- dr_ap_list[[cause]][[age]]$tmrel - for(i in 1:100){#length(alpha)){ - pmcurve <- ap_dose_response_curve(pm,alpha[i],beta[i],gamma[i],tmrel[i]) - pmcurve[is.na(pmcurve)] <- 1 - if(i==1) {plot(pm,pmcurve,typ='l',col=rgb((nsamples-i)/nsamples,0.2,i/nsamples,0.1),ylim=c(1,5),main=ifelse(age%in%c(30,99),as.character(DISEASE_INVENTORY$GBD_name)[j],''),cex.lab=1.5, - frame=F,xlab=ifelse(age>70,'Dose',''),ylab=ifelse(age%in%c(99,25,40,55,70,85),'Relative risk',''),xaxt='n',yaxt='n',cex.main=1.5) - axis(1,labels=ifelse(age>80,T,F),cex.axis=1.5) - axis(2,labels=ifelse(age%in%c(99,25,40,55,70,85),T,F),cex.axis=1.5) - if(length(ages)>1) mtext(side=3,text=paste0('Age ',age),line=0) - } - else lines(pm,pmcurve,col=rgb((nsamples-i)/nsamples,0.2,i/nsamples,0.1),cex.axis=1.5) - } - } - dev.off() - if(length(ages)>1){ - pdf(paste0('results/dose_response/AP/',cause,'DR.pdf'),width=6,height=8); par(mfrow=c(1,1)) - layout.matrix <- matrix(1:length(ages),nrow=5,ncol=3,byrow=T) - graphics::layout(mat=layout.matrix,heights=c(4,3,3,3,4),widths=c(3.8,3,3)) - }else{ - pdf(paste0('results/dose_response/AP/',cause,'DR.pdf')) - par(mfrow=c(1,1)) - } - for(age in ages){ - par(mar=c(ifelse(age>80,5,1),ifelse(age%in%c(99,25,40,55,70,85),5,1),ifelse(age%in%c(25:35,99),5,1),1)) - dr_ap_age <- subset(dr_ap,age_code==age) - alpha <- dr_ap_age$alpha - beta <- dr_ap_age$beta - gamma <- dr_ap_age$gamma - tmrel <- dr_ap_age$tmrel - for(i in 1:100){#length(alpha)){ - pmcurve <- ap_dose_response_curve(pm,alpha[i],beta[i],gamma[i],tmrel[i]) - pmcurve[is.na(pmcurve)] <- 1 - if(i==1) {plot(pm,pmcurve,typ='l',col=rgb((nsamples-i)/nsamples,0.2,i/nsamples,0.1),ylim=c(1,5),main=ifelse(age%in%c(30,99),as.character(DISEASE_INVENTORY$GBD_name)[j],''),cex.lab=1.5, - frame=F,xlab=ifelse(age>70,'Dose',''),ylab=ifelse(age%in%c(99,25,40,55,70,85),'Relative risk',''),xaxt='n',yaxt='n',cex.main=1.5) - axis(1,labels=ifelse(age>80,T,F),cex.axis=1.5) - axis(2,labels=ifelse(age%in%c(99,25,40,55,70,85),T,F),cex.axis=1.5) - if(length(ages)>1) mtext(side=3,text=paste0('Age ',age),line=0) - } - else lines(pm,pmcurve,col=rgb((nsamples-i)/nsamples,0.2,i/nsamples,0.1),cex.axis=1.5) - } - } - dev.off() -} - - - diff --git a/tests/test_pa_dose_response.R b/tests/test_pa_dose_response.R deleted file mode 100644 index c5fac876..00000000 --- a/tests/test_pa_dose_response.R +++ /dev/null @@ -1,142 +0,0 @@ -rm(list=ls()) -require(ithimr) - -global_path <- file.path(find.package('ithimr',lib.loc=.libPaths()), 'extdata/global/') -global_path <- paste0(global_path, "/") - -DISEASE_INVENTORY <- read.csv(paste0(global_path,"dose_response/disease_outcomes_lookup.csv")) - -list_of_files <- list.files(path = paste0(global_path,"dose_response/drpa/extdata/"), recursive = TRUE, pattern = "\\.csv$", full.names = TRUE) -for (i in 1:length(list_of_files)){ - assign(stringr::str_sub(basename(list_of_files[[i]]), end = -5), - readr::read_csv(list_of_files[[i]],col_types = cols()), - pos = 1) -} - -doses_vector <- as.numeric(0:50) -cols <- rainbow(sum(DISEASE_INVENTORY$physical_activity == 1)) -cause_indices <- which(DISEASE_INVENTORY$physical_activity == 1) - -############# median dose-response curves -PA_DOSE_RESPONSE_QUANTILE <<- 0.5 -#pdf('results/dose_response/PA_dose_response.pdf',width=11,height=11);par(mar=c(5,5,2,1)) -for (j in 1:length(cause_indices)){ - cause <- as.character(DISEASE_INVENTORY$pa_acronym[cause_indices[j]]) - outcome <- ifelse(DISEASE_INVENTORY$outcome[cause_indices[j]] == "deaths", "fatal", - "fatal-and-non-fatal") - # get name of disease - print(cause) - print(outcome) - return_vector <- drpa::dose_response(cause = cause, outcome_type = outcome, - dose = doses_vector) - - if(j==1){ - plot(doses_vector,return_vector$rr,type='l',ylim=0:1,main='',ylab='Relative risk',xlab='Dose', - frame=F,col=cols[j],xlim=range(doses_vector),cex.axis=1.5,cex.lab=1.5,cex.main=1.5,lwd=2) - }else{ - lines(doses_vector,return_vector$rr,col=cols[j],lwd=2) - } -} -legend(legend=DISEASE_INVENTORY$pa_acronym[DISEASE_INVENTORY$physical_activity==1],col=cols,bty='n',x=0,y=0.5,lwd=2) -#dev.off() - -############# sample dose-response curves -nsamples <- 10 -parameters <- ithim_setup_parameters(NSAMPLES=nsamples, - PA_DOSE_RESPONSE_QUANTILE=T) -#pdf('results/dose_response/PA_dose_response_sample.pdf',width=11,height=11); - -ls_dr <- list() -par(mar=c(5,5,2,1),mfrow=c(3,3)) -for ( j in 1:length(cause_indices)){ - cause <- as.character(DISEASE_INVENTORY$pa_acronym[cause_indices[j]]) - # get name of disease - print(cause) - for(seed in 1:nsamples){ - for(i in 1:length(parameters)) - assign(names(parameters)[i],parameters[[i]][[seed]],pos=1) - - outcome <- ifelse(DISEASE_INVENTORY$outcome[cause_indices[j]] == "deaths", "fatal", - "fatal-and-non-fatal") - return_vector <- drpa::dose_response(cause = cause, outcome_type = outcome, - dose = doses_vector, quantile = parameters[[i]][[seed]]) - - ls_dr[[cause]][[seed]] <- return_vector - - if(seed==1){ - plot(doses_vector,return_vector$rr,type='l',ylim=0:1,main=cause,ylab='Relative risk',xlab='Dose', - frame=F,col=cols[j],xlim=range(doses_vector),cex.axis=1.5,cex.lab=1.5,cex.main=1.5,lwd=1) - }else{ - lines(doses_vector,return_vector$rr,col=cols[j],lwd=1) - } - } -} -#dev.off() - -################ explanatory graph -library(diagram) -j <- 3 -ylow <- 0.6 -{x11(); par(mfrow=c(2,2),mar=c(2,5,1,1)) -pa_dn <- as.character(DISEASE_INVENTORY$pa_acronym[j]) -pa_n <- as.character(DISEASE_INVENTORY$acronym[j]) -outcome_type <- ifelse(pa_dn%in%c('lung_cancer','stroke'), 'incidence' , 'mortality') -dose <- doses_clean -if(pa_dn %in% c('total_cancer','coronary_heart_disease')) dose[dose>35] <- 35 -cause<-pa_dn -fname <- paste(cause, outcome_type, sep = "_") -lookup_table <- get(paste0(fname)) -lookup_df <- as.data.frame(lookup_table) -rr <- approx(x=lookup_df$dose,y=lookup_df$RR,xout=dose,yleft=1,yright=min(lookup_df$RR))$y -lb <- approx(x = lookup_df$dose,y = lookup_df$lb,xout = dose,yleft = 1,yright = min(lookup_df$lb))$y -ub <-approx(x = lookup_df$dose,y = lookup_df$ub,xout = dose,yleft = 1,yright = min(lookup_df$ub))$y -plot(doses_clean,rr,lwd=2,type='l',col=rgb(0.6,0.2,0.2), - xlim=range(doses_clean),ylim=c(ylow,1),ylab='Response',frame=F,cex.lab=1.5,cex.axis=1.5,xlab='') -lines(doses_clean,lb,lwd=2,type='l',lty=2,col=rgb(0.6,0.2,0.2,1)) -lines(doses_clean,ub,lwd=2,type='l',lty=2,col=rgb(0.6,0.2,0.2,1)) -text('Mean',x=40,y=0.8,col=rgb(0.6,0.2,0.2),cex=1.5) -text('Lower bound',x=35,y=0.66,col=rgb(0.6,0.2,0.2),cex=1.5) -text('Upper bound',x=35,y=0.93,col=rgb(0.6,0.2,0.2),cex=1.5) -par(mar=c(2,2,1,1)) -for(i in 1:500){ - x<-runif(1) - rr0 <- truncnorm::qtruncnorm(p=x, mean=rr, sd=(rr-lb)/1.96,a=0)#a=lb,b=ub)# - if(i==1){plot(doses_clean,rr0,type='l',col=rgb(0.6,0.2,0.2,0.1), - xlim=range(doses_clean),ylim=c(ylow,1),ylab='',frame=F,cex.lab=1.5,cex.axis=1.5,xlab='') - text('Samples',x=35,y=0.97,col=rgb(0.6,0.2,0.2),cex=1.5) - }else lines(doses_clean,rr0,col=rgb(0.6,0.2,0.2,0.1)) -} -quant <- 0.5 -par(mar=c(5,5,1,1)) -rr0 <- truncnorm::qtruncnorm(p=quant, mean=rr, sd=(rr-lb)/1.96,a=0) -plot(doses_clean,rr0,lwd=2,type='l',col='grey', - xlim=range(doses_clean),ylim=c(ylow,1),ylab='Response',frame=F,cex.lab=1.5,cex.axis=1.5,xlab='Dose') -doses <- c(5,25,45) -text('Median',x=10,y=0.97,col='grey',cex=1.5) -abline(h=ylow,lwd=2,lty=2,col='grey') -points(doses,rep(ylow,length(doses)),col=rgb(0.6,0.2,0.2),pch=18,cex=1.5) -points(doses,rr0[which(doses_clean%in%doses)],col=rgb(0.6,0.2,0.2),pch=15,cex=1.5) -for(dose in doses) arrows(x0=dose,y0=ylow,y1=rr0[which(doses_clean==dose)]-0.02,col=rgb(0.6,0.2,0.2),length=0.1,lty=1,lwd=2) -par(mar=c(5,2,1,1)) -cols <- c('hotpink','navyblue','darkorange2') -quants <- c(0.1,0.6,0.9) -for(i in 1:length(quants)){ - quant <- quants[i] - rr0 <- truncnorm::qtruncnorm(p=quant, mean=rr, sd=(rr-lb)/1.96,a=0) - if(i==1) {plot(doses_clean,rr0,lwd=2,type='l',col=cols[i], - xlim=range(doses_clean),ylim=c(ylow,1),ylab='',frame=F,cex.lab=1.5,cex.axis=1.5,xlab='Dose') - for(dose in doses) arrows(x0=dose,y0=ylow,y1=rr0[which(doses_clean==dose)]-0.02,col=cols[i],length=0.1,lty=1,lwd=2) - } else{ - lines(doses_clean,rr0,lwd=2,type='l',col=cols[i]) - for(dose in doses) - curvedarrow(from=c(dose,ylow),to=c(dose,rr0[which(doses_clean==dose)]-0.01),lty=1,lwd=2,lcol=cols[i], - arr.col=cols[i],arr.pos=0.65+0.05*i,curve=-70+30*i,arr.type='curved',arr.adj=0,endhead=T) - } - abline(h=ylow,lwd=2,lty=2,col='grey') - points(doses,rr0[which(doses_clean%in%doses)],col=cols[i],pch=15,cex=1.5) - - points(doses,rep(ylow,length(doses)),col=rgb(0.6,0.2,0.2),pch=18,cex=1.5) -} -} - - diff --git a/tests/testthat.R b/tests/testthat.R deleted file mode 100644 index 3b2e4e00..00000000 --- a/tests/testthat.R +++ /dev/null @@ -1,4 +0,0 @@ -library(testthat) -library(ithimr) - -test_check("ithimr") diff --git a/tests/testthat/accra_evppi_test.Rds b/tests/testthat/accra_evppi_test.Rds deleted file mode 100644 index 15f8960a..00000000 Binary files a/tests/testthat/accra_evppi_test.Rds and /dev/null differ diff --git a/tests/testthat/accra_results.Rds b/tests/testthat/accra_results.Rds deleted file mode 100644 index 5790cfa0..00000000 Binary files a/tests/testthat/accra_results.Rds and /dev/null differ diff --git a/tests/testthat/accra_walk_results.Rds b/tests/testthat/accra_walk_results.Rds deleted file mode 100644 index 44f39905..00000000 Binary files a/tests/testthat/accra_walk_results.Rds and /dev/null differ diff --git a/tests/testthat/test-ithimr.R b/tests/testthat/test-ithimr.R deleted file mode 100644 index 2a5eba54..00000000 --- a/tests/testthat/test-ithimr.R +++ /dev/null @@ -1,133 +0,0 @@ -context("test-ithimr") - -test_that("accra basic", { - # load saved result - accra_results <- readRDS('accra_results.Rds') - # generate new baseline accra results - ithim_object <- run_ithim_setup(REFERENCE_SCENARIO='Scenario 1') - ithim_object$outcomes <- run_ithim(ithim_object, seed = 1) - result_mat <- colSums(ithim_object$outcome$hb$ylls[,3:ncol(ithim_object$outcome$hb$ylls)]) - # test - expect_equal(result_mat,accra_results,tolerance=0.1) - -}) - -test_that("accra_test, walk scenario", { - # load saved result - accra_results <- readRDS('accra_walk_results.Rds') - # generate new baseline accra results - ithim_object <- run_ithim_setup( CITY='accra_test', - TEST_WALK_SCENARIO=T, - ADD_WALK_TO_BUS_TRIPS=F, - ADD_TRUCK_DRIVERS = F, - ADD_BUS_DRIVERS = F) - ithim_object$outcomes <- run_ithim(ithim_object, seed = 1) - result_mat <- colSums(ithim_object$outcome$hb$ylls[,3:ncol(ithim_object$outcome$hb$ylls)]) - # test - expect_equal(accra_results,result_mat,tolerance=0.1) - -}) - -test_that("accra uncertain parallel", { - # load saved result - #accra_evppi <- readRDS('accra_evppi_test.Rds') - # generate new accra evppi results - - ithim_object <- run_ithim_setup(REFERENCE_SCENARIO='Scenario 1', - NSAMPLES = 16, - BUS_WALK_TIME = c(log(5), log(1.2)), - MMET_CYCLING = c(log(5), log(1.2)), - MMET_WALKING = c(log(2.5), log(1.2)), - INJURY_REPORTING_RATE = c(8,3), - CHRONIC_DISEASE_SCALAR = c(0,log(1.2)), - PM_CONC_BASE = c(log(50),log(1.2)), - PM_TRANS_SHARE = c(5,20), - BACKGROUND_PA_SCALAR = c(0,log(1.2)), - MOTORCYCLE_TO_CAR_RATIO = c(-1.4,0.4), - DAY_TO_WEEK_TRAVEL_SCALAR = c(20,3), - INJURY_LINEARITY= c(log(1),log(1.2)), - CASUALTY_EXPONENT_FRACTION = c(8,8), - PA_DOSE_RESPONSE_QUANTILE = F, - AP_DOSE_RESPONSE_QUANTILE = F) - numcores <- detectCores() - ithim_object$outcomes <- mclapply(1:NSAMPLES, FUN = ithim_uncertainty, ithim_object = ithim_object,mc.cores = ifelse(Sys.info()[['sysname']] == "Windows", 1, numcores)) - ## calculate EVPPI - parameter_names <- names(ithim_object$parameters)[names(ithim_object$parameters)!="DR_AP_LIST"] - parameter_samples <- sapply(parameter_names,function(x)ithim_object$parameters[[x]]) - ## omit all-cause mortality - outcome <- t(sapply(ithim_object$outcomes, function(x) colSums(x$hb$deaths[,3:ncol(x$hb$deaths)]))) - NSCEN <- nrow(ithim_objects$scen_prop) - evppi <- matrix(0, ncol = NSCEN, nrow = ncol(parameter_samples)) - for(j in 1:(NSCEN)){ - y <- rowSums(outcome[,seq(NSCEN+j,ncol(outcome),by=NSCEN)]) - vary <- var(y) - for(i in 1:ncol(parameter_samples)){ - x <- parameter_samples[, i]; - model <- gam(y ~ s(x)) - evppi[i, j] <- (vary - mean((y - model$fitted) ^ 2)) / vary * 100 - } - } - colnames(evppi) <- c("base",paste0("scen", 1:NSCEN) )[c(1,3:6)] - rownames(evppi) <- colnames(parameter_samples) - # test - expect_equal(sum(evppi>0),60) -}) - -test_that("accra evppi", { - ## this takes a long time: only test if you really want to - test_evppi <- F - if(test_evppi){ - # load saved result - #accra_evppi <- readRDS('accra_evppi_test.Rds') - - # generate new accra evppi results - ithim_object <- run_ithim_setup(REFERENCE_SCENARIO = 'Scenario 1', - NSAMPLES = 1024, - BUS_WALK_TIME = c(log(5), log(1.2)), - MMET_CYCLING = c(log(5), log(1.2)), - MMET_WALKING = c(log(2.5), log(1.2)), - INJURY_REPORTING_RATE = c(8,3), - CHRONIC_DISEASE_SCALAR = c(0,log(1.2)), - PM_CONC_BASE = c(log(50),log(1.2)), - PM_TRANS_SHARE = c(5,20), - BACKGROUND_PA_SCALAR = c(0,log(1.2)), - MOTORCYCLE_TO_CAR_RATIO = c(-1.4,0.4), - PA_DOSE_RESPONSE_QUANTILE = T, - AP_DOSE_RESPONSE_QUANTILE = T) - numcores <- detectCores() - ithim_object$outcomes <- mclapply(1:NSAMPLES, FUN = ithim_uncertainty, ithim_object = ithim_object,mc.cores = ifelse(Sys.info()[['sysname']] == "Windows", 1, numcores)) - ## calculate EVPPI - parameter_names <- names(ithim_object$parameters)[names(ithim_object$parameters)!="DR_AP_LIST"] - parameter_samples <- sapply(parameter_names,function(x)ithim_object$parameters[[x]]) - ## omit all-cause mortality - outcome <- t(sapply(ithim_object$outcomes, function(x) colSums(x$hb$deaths[,3:ncol(x$hb$deaths)]))) - NSCEN <- nrow(ithim_objects$scen_prop) - evppi <- matrix(0, ncol = NSCEN, nrow = ncol(parameter_samples)) - for(j in 1:(NSCEN)){ - y <- rowSums(outcome[,seq(NSCEN+j,ncol(outcome),by=NSCEN)]) - vary <- var(y) - for(i in 1:ncol(parameter_samples)){ - x <- parameter_samples[, i]; - model <- gam(y ~ s(x)) - evppi[i, j] <- (vary - mean((y - model$fitted) ^ 2)) / vary * 100 - } - } - colnames(evppi) <- c("base",paste0("scen", 1:NSCEN) )[c(1,3:6)] - rownames(evppi) <- colnames(parameter_samples) - ## add four-dimensional EVPPI if AP_DOSE_RESPONSE is uncertain. - if("DR_AP_LIST"%in%names(ithim_object$parameters)&&NSAMPLES>=1024){ - AP_names <- sapply(names(ithim_object$parameters),function(x)length(strsplit(x,'AP_DOSE_RESPONSE_QUANTILE_ALPHA')[[1]])>1) - diseases <- sapply(names(ithim_object$parameters)[AP_names],function(x)strsplit(x,'AP_DOSE_RESPONSE_QUANTILE_ALPHA_')[[1]][2]) - evppi_for_AP <- mclapply(diseases, FUN = parallel_evppi_for_AP,parameter_samples,outcome,NSCEN, mc.cores = ifelse(Sys.info()[['sysname']] == "Windows", 1, numcores)) - names(evppi_for_AP) <- paste0('AP_DOSE_RESPONSE_QUANTILE_',diseases) - evppi <- rbind(evppi,do.call(rbind,evppi_for_AP)) - ## get rows to remove - keep_names <- sapply(rownames(evppi),function(x)!any(c('ALPHA','BETA','GAMMA','TMREL')%in%strsplit(x,'_')[[1]])) - evppi <- evppi[keep_names,] - } - # test - expect_equal(sum(evppi>0),95) - }else{ - expect_equal(0,0) - } -})