Ithim setup

functions/error_handling.R

@@ -0,0 +1,25 @@
+# Programmed by Javad Rahimipour Anaraki on 09/11/18
+# Ph.D. Candidate
+# Department of Computer Science
+# Memorial University of Newfoundland
+# jra066 [AT] mun [DOT] ca | www.cs.mun.ca/~jra066
+
+#   input: Error/Warning code, function name, and valiable name
+#  output: Error/Warning message
+
+error_handling <- function(code, function_name, variable_name) {
+  output <-
+    switch(
+      code,
+      paste0(
+        "In ",
+        function_name,
+        "(",
+        variable_name,
+        "): Default values should be changed!"
+      )
+    )
+
+  print(output)
+
+}

functions/ithim_setup_parameters.R

@@ -1,108 +1,229 @@
+# Modified by Ali Abbas
+# Research Associate
+# UKCRC Centre for Diet and Activity Research (CEDAR)
+# MRC Epidemiology Unit
+# University of Cambridge School of Clinical Medicine
+# aa797 [AT] medschl [DOT] cam [DOT] ac [DOT] uk | http://www.cedar.iph.cam.ac.uk/people/cdfs/ali-abbas/
+
+# Modified by Daniel Fuller
+# Canada Research Chair in Population Physical Activity
+# School of Human Kinetics and Recreation
+# Memorial University of Newfoundland
+# dfuller [AT] mun [DOT] ca | www.walkabilly.ca/home/
+
+# Modified by Javad Rahimipour Anaraki on 08/11/18
+# Ph.D. Candidate
+# Department of Computer Science
+# Memorial University of Newfoundland
+# jra066 [AT] mun [DOT] ca | www.cs.mun.ca/~jra066
+
+#   input:
+#  output:
+
 ithim_setup_parameters <- function(NSAMPLES = 1,
-                                   MEAN_BUS_WALK_TIME= 5,
+                                   MEAN_BUS_WALK_TIME = 5,
                                    MMET_CYCLING = 4.63,
                                    MMET_WALKING = 2.53,
-                                   PM_CONC_BASE = 50,  
+                                   PM_CONC_BASE = 50,
                                    PM_TRANS_SHARE = 0.225,
                                    PA_DOSE_RESPONSE_QUANTILE = F,
                                    AP_DOSE_RESPONSE_QUANTILE = F,
                                    BACKGROUND_PA_SCALAR = 1,
                                    SAFETY_SCALAR = 1,
-                                   CHRONIC_DISEASE_SCALAR = 1 ){
+                                   CHRONIC_DISEASE_SCALAR = 1) {
+  #=======================Function=========================
+  source("error_handling.R")
+
+
+  #=======================Warnings=========================
+  if (PM_CONC_BASE == 50 |
+      PM_TRANS_SHARE == 0.225)
+    error_handling(1, "ithim_setup_parameters", "PM_CONC_BASE, PM_TRANS_SHARE")
+
+
+  #=======================Variables========================
   ## PARAMETERS
   ##RJ parameters are assigned to the environment and so are set for every function. They are over-written when sample_parameters is called.
   parameters <- list()
-  if(length(MEAN_BUS_WALK_TIME) == 1 ) {
-    MEAN_BUS_WALK_TIME <<- MEAN_BUS_WALK_TIME
-  }else{
-    parameters$MEAN_BUS_WALK_TIME <- rlnorm(NSAMPLES,MEAN_BUS_WALK_TIME[1], MEAN_BUS_WALK_TIME[2])
-  }
-  if(length(MMET_CYCLING) == 1 ) {
-    MMET_CYCLING <<- MMET_CYCLING
-  }else{
-    parameters$MMET_CYCLING <- rlnorm(NSAMPLES,MMET_CYCLING[1], MMET_CYCLING[2])
-  }
-  if(length(MMET_WALKING) == 1 ) {
-    MMET_WALKING <<- MMET_WALKING
-  }else{
-    parameters$MMET_WALKING <- rlnorm(NSAMPLES,MMET_WALKING[1], MMET_WALKING[2])
-  }
-  if(length(PM_CONC_BASE) == 1 ) {
-    PM_CONC_BASE <<- PM_CONC_BASE
-  }else{
-    parameters$PM_CONC_BASE <- rlnorm(NSAMPLES,PM_CONC_BASE[1],PM_CONC_BASE[2])
-  }
-  if(length(PM_TRANS_SHARE) == 1 ) {
-    PM_TRANS_SHARE <<- PM_TRANS_SHARE
-  }else{
-    parameters$PM_TRANS_SHARE <- rbeta(NSAMPLES,PM_TRANS_SHARE[1],PM_TRANS_SHARE[2])
-  }
-  if(length(BACKGROUND_PA_SCALAR) == 1 ) {
-    BACKGROUND_PA_SCALAR <<- BACKGROUND_PA_SCALAR
-  }else{
-    parameters$BACKGROUND_PA_SCALAR <- rlnorm(NSAMPLES,BACKGROUND_PA_SCALAR[1],BACKGROUND_PA_SCALAR[2])
-  }
-  if(length(SAFETY_SCALAR) == 1 ) {
-    SAFETY_SCALAR <<- SAFETY_SCALAR
-  }else{
-    parameters$SAFETY_SCALAR <- rlnorm(NSAMPLES,SAFETY_SCALAR[1],SAFETY_SCALAR[2])
-  }
-  if(length(CHRONIC_DISEASE_SCALAR) == 1 ) {
-    CHRONIC_DISEASE_SCALAR <<- CHRONIC_DISEASE_SCALAR
-  }else{
-    parameters$CHRONIC_DISEASE_SCALAR <- rlnorm(NSAMPLES,CHRONIC_DISEASE_SCALAR[1],CHRONIC_DISEASE_SCALAR[2])
+  variables <- c("MEAN_BUS_WALK_TIME",
+                 "MMET_CYCLING",
+                 "MMET_WALKING",
+                 "PM_CONC_BASE",
+                 "PM_TRANS_SHARE",
+                 "BACKGROUND_PA_SCALAR",
+                 "SAFETY_SCALAR",
+                 "CHRONIC_DISEASE_SCALAR")
+
+  #==========================Main==========================
+  for (i in 1:length(variables)) {
+    name <- variables[i]
+    val <- get(variables[i])
+    if (length(val) == 1) {
+      assign(name, val, envir = .GlobalEnv)
+    } else {
+      parameters[[name]] <-
+        rlnorm(NSAMPLES, val[1], val[2])
+    }
   }
-  if(PA_DOSE_RESPONSE_QUANTILE == F ) {
+
+  # if (length(MEAN_BUS_WALK_TIME) == 1) {
+  #   MEAN_BUS_WALK_TIME <<- MEAN_BUS_WALK_TIME
+  # }
+  # else{
+  #   parameters$MEAN_BUS_WALK_TIME <-
+  #     rlnorm(NSAMPLES, MEAN_BUS_WALK_TIME[1], MEAN_BUS_WALK_TIME[2])
+  # }
+  # 
+  # 
+  # if (length(MMET_CYCLING) == 1) {
+  #   MMET_CYCLING <<- MMET_CYCLING
+  # } else{
+  #   parameters$MMET_CYCLING <-
+  #     rlnorm(NSAMPLES, MMET_CYCLING[1], MMET_CYCLING[2])
+  # }
+  # 
+  # if (length(MMET_WALKING) == 1) {
+  #   MMET_WALKING <<- MMET_WALKING
+  # } else{
+  #   parameters$MMET_WALKING <-
+  #     rlnorm(NSAMPLES, MMET_WALKING[1], MMET_WALKING[2])
+  # }
+  # 
+  # if (length(PM_CONC_BASE) == 1) {
+  #   PM_CONC_BASE <<- PM_CONC_BASE
+  # } else{
+  #   parameters$PM_CONC_BASE <-
+  #     rlnorm(NSAMPLES, PM_CONC_BASE[1], PM_CONC_BASE[2])
+  # }
+  # 
+  # if (length(PM_TRANS_SHARE) == 1) {
+  #   PM_TRANS_SHARE <<- PM_TRANS_SHARE
+  # } else{
+  #   parameters$PM_TRANS_SHARE <-
+  #     rbeta(NSAMPLES, PM_TRANS_SHARE[1], PM_TRANS_SHARE[2])
+  # }
+  # 
+  # if (length(BACKGROUND_PA_SCALAR) == 1) {
+  #   BACKGROUND_PA_SCALAR <<- BACKGROUND_PA_SCALAR
+  # } else{
+  #   parameters$BACKGROUND_PA_SCALAR <-
+  #     rlnorm(NSAMPLES, BACKGROUND_PA_SCALAR[1], BACKGROUND_PA_SCALAR[2])
+  # }
+  # 
+  # if (length(SAFETY_SCALAR) == 1) {
+  #   SAFETY_SCALAR <<- SAFETY_SCALAR
+  # } else{
+  #   parameters$SAFETY_SCALAR <-
+  #     rlnorm(NSAMPLES, SAFETY_SCALAR[1], SAFETY_SCALAR[2])
+  # }
+  # 
+  # if (length(CHRONIC_DISEASE_SCALAR) == 1) {
+  #   CHRONIC_DISEASE_SCALAR <<- CHRONIC_DISEASE_SCALAR
+  # } else{
+  #   parameters$CHRONIC_DISEASE_SCALAR <-
+  #     rlnorm(NSAMPLES,
+  #            CHRONIC_DISEASE_SCALAR[1],
+  #            CHRONIC_DISEASE_SCALAR[2])
+  # }
+
+  if (PA_DOSE_RESPONSE_QUANTILE == F) {
     PA_DOSE_RESPONSE_QUANTILE <<- PA_DOSE_RESPONSE_QUANTILE
-  }else{
-    pa_diseases <- subset(DISEASE_OUTCOMES,physical_activity==1)
+  } else{
+    pa_diseases <- subset(DISEASE_OUTCOMES, physical_activity == 1)
     dr_pa_list <- list()
-    for(disease in pa_diseases$pa_acronym)
-      parameters[[paste0('PA_DOSE_RESPONSE_QUANTILE_',disease)]] <- runif(NSAMPLES,0,1)
+    for (disease in pa_diseases$pa_acronym)
+      parameters[[paste0('PA_DOSE_RESPONSE_QUANTILE_', disease)]] <-
+      runif(NSAMPLES, 0, 1)
   }
-  if(AP_DOSE_RESPONSE_QUANTILE == F ) {
+
+  if (AP_DOSE_RESPONSE_QUANTILE == F) {
     AP_DOSE_RESPONSE_QUANTILE <<- AP_DOSE_RESPONSE_QUANTILE
     dr_ap_list <- list()
-    for ( j in 1:nrow(DISEASE_OUTCOMES)) if (DISEASE_OUTCOMES$air_pollution[j] == 1){ 
-      cause <- as.character(DISEASE_OUTCOMES$ap_acronym[j])
-      dr_ap <- subset(DR_AP,cause_code==cause)
-      dr_ap_list[[cause]] <- list()
-      for(age in unique(dr_ap$age_code)){
-        dr_ap_age <- subset(dr_ap,age_code==age)
-        dr_ap_list[[cause]][[as.character(age)]] <- data.frame(alpha=mean(dr_ap_age$alpha),beta=mean(dr_ap_age$beta),gamma=mean(dr_ap_age$gamma),tmrel=mean(dr_ap_age$tmrel))
+    for (j in 1:nrow(DISEASE_OUTCOMES))
+      if (DISEASE_OUTCOMES$air_pollution[j] == 1) {
+        cause <- as.character(DISEASE_OUTCOMES$ap_acronym[j])
+        dr_ap <- subset(DR_AP, cause_code == cause)
+        dr_ap_list[[cause]] <- list()
+        for (age in unique(dr_ap$age_code)) {
+          dr_ap_age <- subset(dr_ap, age_code == age)
+          dr_ap_list[[cause]][[as.character(age)]] <-
+            data.frame(
+              alpha = mean(dr_ap_age$alpha),
+              beta = mean(dr_ap_age$beta),
+              gamma = mean(dr_ap_age$gamma),
+              tmrel = mean(dr_ap_age$tmrel)
+            )
+        }
       }
-    }
     DR_AP_LIST <<- dr_ap_list
-  }else{
-    ap_diseases <- subset(DISEASE_OUTCOMES,air_pollution==1)
-    for(disease in ap_diseases$ap_acronym)
-      for(letter in c('ALPHA_','BETA_','GAMMA_','TMREL_'))
-        parameters[[paste0('AP_DOSE_RESPONSE_QUANTILE_',letter,disease)]] <- runif(NSAMPLES,0,1)
+  } else{
+    ap_diseases <- subset(DISEASE_OUTCOMES, air_pollution == 1)
+    for (disease in ap_diseases$ap_acronym)
+      for (letter in c('ALPHA_', 'BETA_', 'GAMMA_', 'TMREL_'))
+        parameters[[paste0('AP_DOSE_RESPONSE_QUANTILE_', letter, disease)]] <-
+          runif(NSAMPLES, 0, 1)
       dr_ap_list <- list()
-      for(disease in ap_diseases$ap_acronym){ 
-        dr_ap <- subset(DR_AP,cause_code==disease)
+      for (disease in ap_diseases$ap_acronym) {
+        dr_ap <- subset(DR_AP, cause_code == disease)
         dr_ap_list[[disease]] <- list()
-        for(age in unique(dr_ap$age_code)){
-          dr_ap_age <- subset(dr_ap,age_code==age)
+        for (age in unique(dr_ap$age_code)) {
+          dr_ap_age <- subset(dr_ap, age_code == age)
           # generate a value for alpha
-          alpha_val <- quantile(log(dr_ap_age$alpha),parameters[[paste0('AP_DOSE_RESPONSE_QUANTILE_ALPHA_',disease)]])
+          alpha_val <-
+            quantile(log(dr_ap_age$alpha), parameters[[paste0('AP_DOSE_RESPONSE_QUANTILE_ALPHA_', disease)]])
           # generate a value for beta given alpha
-          mod <- gam(log(beta)~ns(log(alpha),df=3),data=dr_ap_age)
-          pred_val <- predict(mod, newdata=data.frame(alpha=exp(alpha_val)),se.fit=T)
-          beta_val <- qnorm(parameters[[paste0('AP_DOSE_RESPONSE_QUANTILE_BETA_',disease)]],pred_val$fit,sqrt(mod$sig2))
+          mod <- gam(log(beta) ~ ns(log(alpha), df = 3), data = dr_ap_age)
+          pred_val <-
+            predict(mod,
+                    newdata = data.frame(alpha = exp(alpha_val)),
+                    se.fit = T)
+          beta_val <-
+            qnorm(parameters[[paste0('AP_DOSE_RESPONSE_QUANTILE_BETA_', disease)]], pred_val$fit, sqrt(mod$sig2))
           # generate a value for gamma given beta and alpha
-          mod <- gam(log(gamma)~ns(log(beta),df=3)+ns(log(alpha),df=3),data=dr_ap_age)
-          pred_val <- predict(mod, newdata=data.frame(alpha=exp(alpha_val),beta=exp(beta_val)),se.fit=T)
-          gamma_val <- qnorm(parameters[[paste0('AP_DOSE_RESPONSE_QUANTILE_GAMMA_',disease)]],pred_val$fit,sqrt(mod$sig2))
+          mod <-
+            gam(log(gamma) ~ ns(log(beta), df = 3) + ns(log(alpha), df = 3), data =
+                  dr_ap_age)
+          pred_val <-
+            predict(mod,
+                    newdata = data.frame(
+                      alpha = exp(alpha_val),
+                      beta = exp(beta_val)
+                    ),
+                    se.fit = T)
+          gamma_val <-
+            qnorm(parameters[[paste0('AP_DOSE_RESPONSE_QUANTILE_GAMMA_', disease)]], pred_val$fit, sqrt(mod$sig2))
           # generate a value for tmrel given alpha, beta and gamma
-          mod <- gam(log(tmrel)~ns(log(gamma),df=3)+ns(log(beta),df=3)+ns(log(alpha),df=3),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(parameters[[paste0('AP_DOSE_RESPONSE_QUANTILE_TMREL_',disease)]],pred_val$fit,sqrt(mod$sig2))
-          dr_ap_list[[disease]][[as.character(age)]] <- data.frame(alpha=exp(alpha_val),beta=exp(beta_val),gamma=exp(gamma_val),tmrel=exp(tmrel_val))
+          mod <-
+            gam(log(tmrel) ~ ns(log(gamma), df = 3) + ns(log(beta), df = 3) + ns(log(alpha), df =
+                                                                                   3),
+                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(parameters[[paste0('AP_DOSE_RESPONSE_QUANTILE_TMREL_', disease)]], pred_val$fit, sqrt(mod$sig2))
+          dr_ap_list[[disease]][[as.character(age)]] <-
+            data.frame(
+              alpha = exp(alpha_val),
+              beta = exp(beta_val),
+              gamma = exp(gamma_val),
+              tmrel = exp(tmrel_val)
+            )
         }
       }
       # 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,])))
+      parameters$DR_AP_LIST <-
+        lapply(1:NSAMPLES, function(x)
+          lapply(dr_ap_list, function(y)
+            lapply(y, function(z)
+              z[x, ])))
   }
   parameters
 }