From 1d8d1281be715c268848a825bb7c84ba6b5dd393 Mon Sep 17 00:00:00 2001
From: Ali Abbas <m.aliabbas@gmail.com>
Date: Wed, 20 Oct 2021 16:32:02 +0100
Subject: [PATCH 1/2] Adjust code to read dr for AP for both constant and
 uncertain modes

---
 R/ithim_setup_parameters.R | 64 +++++++-------------------------------
 1 file changed, 12 insertions(+), 52 deletions(-)

diff --git a/R/ithim_setup_parameters.R b/R/ithim_setup_parameters.R
index 024243cb..54218d20 100644
--- a/R/ithim_setup_parameters.R
+++ b/R/ithim_setup_parameters.R
@@ -163,60 +163,20 @@ ithim_setup_parameters <- function(NSAMPLES = 1,
   }
   
   #### AP DOSE RESPONSE
-  AP_DOSE_RESPONSE_QUANTILE <<- AP_DOSE_RESPONSE_QUANTILE
-  ## shortcut: use saved median values
-  if(!AP_DOSE_RESPONSE_QUANTILE){
-    global_path <- file.path(find.package('ithimr',lib.loc=.libPaths()), 'extdata/global/')
-    global_path <- paste0(global_path, "/")
-    DR_AP_LIST <<- readRDS(paste0(global_path,"dose_response/drap/dr_ap_list.Rds"))
-  }else{
+  if(AP_DOSE_RESPONSE_QUANTILE == T ) {
+    ap_diseases <- subset(DISEASE_INVENTORY, air_pollution==1)
     dr_ap_list <- list()
-    ap_diseases <- subset(DISEASE_INVENTORY,air_pollution==1)
-    ap_parameters <- list()
-    for(disease in ap_diseases$ap_acronym){ 
-      for(letter in c('ALPHA_','BETA_','GAMMA_','TMREL_')){
-        if(AP_DOSE_RESPONSE_QUANTILE){
-          ap_parameters[[paste0('AP_DOSE_RESPONSE_QUANTILE_',letter,disease)]] <- runif(NSAMPLES,0,1)
-          parameters[[paste0('AP_DOSE_RESPONSE_QUANTILE_',letter,disease)]] <- ap_parameters[[paste0('AP_DOSE_RESPONSE_QUANTILE_',letter,disease)]]
-        } else {
-          ap_parameters[[paste0('AP_DOSE_RESPONSE_QUANTILE_',letter,disease)]] <- 0.5
-        }
-      }
-      dr_ap <- subset(DR_AP,cause_code==disease)
-      dr_ap_list[[disease]] <- list()
-      quant1 <- ap_parameters[[paste0('AP_DOSE_RESPONSE_QUANTILE_GAMMA_',disease)]]
-      quant2 <- ap_parameters[[paste0('AP_DOSE_RESPONSE_QUANTILE_BETA_',disease)]]
-      quant3 <- ap_parameters[[paste0('AP_DOSE_RESPONSE_QUANTILE_ALPHA_',disease)]]
-      quant4 <- ap_parameters[[paste0('AP_DOSE_RESPONSE_QUANTILE_TMREL_',disease)]]
-      for(age in unique(dr_ap$age_code)){
-        dr_ap_age <- subset(dr_ap,age_code==age)
-        #######################################
-        lbeta <- log(dr_ap_age$beta)
-        lgamma <- log(dr_ap_age$gamma)
-        gamma_val <- quantile(density(lgamma),quant1)
-        beta_val <- c()
-        for(i in 1:ifelse(AP_DOSE_RESPONSE_QUANTILE,NSAMPLES,1)){
-          den <- kde2d(lgamma,lbeta,n=c(1,100),h=0.2,lims=c(gamma_val[i],gamma_val[i],min(lbeta)-1,max(lbeta)+1))
-          beta_val[i] <- approx(x=cumsum(den$z)/sum(den$z),y=den$y,xout=quant2[i])$y
-        }
-        mod <- gam(log(alpha)~te(log(gamma),log(beta)),data=dr_ap_age)
-        pred_val <- predict(mod, newdata=data.frame(beta=exp(beta_val),gamma=exp(gamma_val)),se.fit=T)
-        alpha_val <- qnorm(quant3,pred_val$fit,sqrt(mod$sig2))
-        # generate a value for tmrel given alpha, beta and gamma
-        mod <- gam(log(tmrel)~ns(log(gamma),df=8)+ns(log(beta),df=8)+ns(log(alpha),df=8),data=dr_ap_age)
-        pred_val <- predict(mod, newdata=data.frame(alpha=exp(alpha_val),beta=exp(beta_val),gamma=exp(gamma_val)),se.fit=T)
-        tmrel_val <- qnorm(quant4,pred_val$fit,sqrt(mod$sig2))
-        dr_ap_list[[disease]][[age]] <- data.frame(alpha=exp(alpha_val),beta=exp(beta_val),gamma=exp(gamma_val),tmrel=exp(tmrel_val))
-      }
-      if(AP_DOSE_RESPONSE_QUANTILE){ 
-        # turn list inside out, so it's indexed first by sample
-        parameters$DR_AP_LIST <- lapply(1:NSAMPLES,function(x)lapply(dr_ap_list,function(y) lapply(y,function(z)z[x,])))
-      }else{
-        DR_AP_LIST <<- dr_ap_list
-      }
-    }
-    
+    for(disease in ap_diseases$ap_acronym)
+      parameters[[paste0('AP_DOSE_RESPONSE_QUANTILE_',disease)]] <- runif(NSAMPLES,0,1)
   }
+  
+  # Read dr_ap_list.Rds which contains parameter values by diseases and ages, for both constant and uncertain mode
+  global_path <- file.path(find.package('ithimr',lib.loc=.libPaths()), 'extdata/global/')
+  global_path <- paste0(global_path, "/")
+  DR_AP_LIST <<- readRDS(paste0(global_path,"dose_response/drap/dr_ap_list.Rds"))
+  
+  parameters$DR_AP_LIST <- DR_AP_LIST
+  
   parameters
 }
 

From 0d75a8260a93a1846eb134df6e34312adaaf8610 Mon Sep 17 00:00:00 2001
From: Ali Abbas <m.aliabbas@gmail.com>
Date: Wed, 20 Oct 2021 17:01:28 +0100
Subject: [PATCH 2/2] Initialize parameters only for uncertain run

---
 R/ithim_setup_parameters.R | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/R/ithim_setup_parameters.R b/R/ithim_setup_parameters.R
index 54218d20..87f823f0 100644
--- a/R/ithim_setup_parameters.R
+++ b/R/ithim_setup_parameters.R
@@ -175,7 +175,8 @@ ithim_setup_parameters <- function(NSAMPLES = 1,
   global_path <- paste0(global_path, "/")
   DR_AP_LIST <<- readRDS(paste0(global_path,"dose_response/drap/dr_ap_list.Rds"))
   
-  parameters$DR_AP_LIST <- DR_AP_LIST
+  if(AP_DOSE_RESPONSE_QUANTILE)
+    parameters$DR_AP_LIST <- DR_AP_LIST
   
   parameters
 }