Merge pull request #102 from KateJohnson/master

Pathway for overdiagnosis added

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: epicR
 Title: R Package for Evaluation Platform in COPD
-Version: 0.26.2
+Version: 0.26.3
 Authors@R: c(
     person("Mohsen", "Sadatsafavi", email = "[email protected]", role = c("aut", "cph")),
     person("Amin", "Adibi", email = "[email protected]", role = c("aut", "cre")))

NAMESPACE

@@ -46,6 +46,7 @@ export(validate_exacerbation)
 export(validate_gpvisits)
 export(validate_lung_function)
 export(validate_mortality)
+export(validate_overdiagnosis)
 export(validate_payoffs)
 export(validate_population)
 export(validate_smoking)

R/input.R

@@ -428,6 +428,17 @@ init_input <- function() {
   input_ref$diagnosis$logit_p_diagnosis_by_sex <- "Kate's regression on CanCOLD, provided on 2019-05-29"
   input$diagnosis$p_hosp_diagnosis <- 0.5
 
+  input_help$diagnosis$logit_p_overdiagnosis_by_sex <- "Probability of being overdiagnosed for non-COPD subjects"
+  input$diagnosis$logit_p_overdiagnosis_by_sex <- cbind(male=c(intercept=-5.2169, age=0.0025, smoking=0.6911, gpvisits=0.0075,
+                                                               cough=0.7264, phlegm=0.7956, wheeze=0.66, dyspnea=0.8798,
+                                                               case_detection=0),
+                                                    female=c(intercept=-5.2169+0.2597, age=0.0025, smoking=0.6911, gpvisits=0.0075,
+                                                             cough=0.7264, phlegm=0.7956, wheeze=0.66, dyspnea=0.8798,
+                                                             case_detection=0))
+  input_ref$diagnosis$logit_p_overdiagnosis_by_sex <- "Kate's regression on CanCOLD, provided on 2019-07-16"
+  input$diagnosis$p_correct_overdiagnosis <- 0.9
+
+
   # Case detection;
 
   input_help$diagnosis$p_case_detection <- "Probability of having case detection given an undiagnosed patient meets the selection criteria"

R/validation.R

@@ -1325,6 +1325,78 @@ validate_treatment<- function(n_sim = 1e+04) {
   cat(mean(exac.diff$Severe),"more severe exacerbations, and\n")
   cat(mean(exac.diff$VerySevere),"more very severe exacerbations per year.\n")
 
+  ###
+  cat("\n")
+  cat("Now, set all COPD patients to diagnosed, then undiagnosed, and compare the exacerbation rates.\n")
+  cat("\n")
+
+  init_session(settings = settings)
+
+  input_nd <- model_input$values
+
+  input_nd$diagnosis$logit_p_diagnosis_by_sex <- cbind(male=c(intercept=-100, age=-0.0324, smoking=0.3711, fev1=-0.8032,
+                                                              gpvisits=0.0087, cough=0.208, phlegm=0.4088, wheeze=0.0321, dyspnea=0.722,
+                                                              case_detection=0),
+                                                       female=c(intercept=-100, age=-0.0324, smoking=0.3711, fev1=-0.8032,
+                                                                gpvisits=0.0087, cough=0.208, phlegm=0.4088, wheeze=0.0321, dyspnea=0.722,
+                                                                case_detection=0))
+
+  res <- run(input = input_nd)
+  if (res < 0)
+    stop("Execution stopped.\n")
+
+  output_ex_nd <- Cget_output_ex()
+
+  exac_rate_nodiag <- rowSums(output_ex_nd$n_exac_by_ctime_severity)/rowSums(output_ex_nd$n_COPD_by_ctime_sex)
+
+  terminate_session()
+
+  ###
+
+  init_session(settings = settings)
+
+  input_d <- model_input$values
+
+  input_d$diagnosis$logit_p_diagnosis_by_sex <- cbind(male=c(intercept=100, age=-0.0324, smoking=0.3711, fev1=-0.8032,
+                                                              gpvisits=0.0087, cough=0.208, phlegm=0.4088, wheeze=0.0321, dyspnea=0.722,
+                                                              case_detection=0),
+                                                       female=c(intercept=100, age=-0.0324, smoking=0.3711, fev1=-0.8032,
+                                                                gpvisits=0.0087, cough=0.208, phlegm=0.4088, wheeze=0.0321, dyspnea=0.722,
+                                                                case_detection=0))
+
+  res <- run(input = input_d)
+  if (res < 0)
+    stop("Execution stopped.\n")
+
+  inputs_d <- Cget_inputs()
+  output_ex_d <- Cget_output_ex()
+
+  exac_rate_diag <- rowSums(output_ex_d$n_exac_by_ctime_severity)/rowSums(output_ex_d$n_COPD_by_ctime_sex)
+
+  ##
+  cat("Annual exacerbation rate:\n")
+  cat("\n")
+
+  trt_effect<- data.frame(Year=1:inputs_d$global_parameters$time_horizon,
+                          Diagnosed = exac_rate_diag,
+                          Undiagnosed = exac_rate_nodiag)
+
+  trt_effect$Delta <- trt_effect$Undiagnosed - trt_effect$Diagnosed
+
+  print(trt_effect)
+
+  cat("\n")
+  cat("Treatment reduces the rate of exacerbations by a mean of:", mean(trt_effect$Delta),"\n")
+
+  # plot
+  trt.plot <- tidyr::gather(data=trt_effect, key="Diagnosis", value="Rate", Diagnosed:Undiagnosed)
+
+  trt.plotted <- ggplot2::ggplot(trt.plot, aes(x=Year, y=Rate, col=Diagnosis)) +
+                            geom_line() + geom_point() + expand_limits(y = 0) +
+                            theme_bw() + ylab("Annual exacerbation rate") + xlab("Years")
+
+  plot(trt.plotted)
+
   terminate_session()
 }
 
@@ -1503,3 +1575,53 @@ test_case_detection <- function(n_sim = 1e+04, p_of_CD=0.1, min_age=40, min_pack
 }
 
 
+#' Returns results of validation tests for overdiagnosis
+#' @param n_sim number of agents
+#' @return validation test results
+#' @export
+validate_overdiagnosis <- function(n_sim = 1e+04) {
+  cat("Let's take a look at overdiagnosis\n")
+  petoc()
+
+  settings <- default_settings
+  settings$record_mode <- record_mode["record_mode_none"]
+  settings$agent_stack_size <- 0
+  settings$n_base_agents <- n_sim
+  settings$event_stack_size <- 0
+  init_session(settings = settings)
+
+  input <- model_input$values
+
+  res <- run(input = input)
+  if (res < 0)
+    stop("Execution stopped.\n")
+
+  inputs <- Cget_inputs()
+  output_ex <- Cget_output_ex()
+
+  cat("Here are the proportion of non-COPD subjects overdiagnosed over model time: \n")
+
+  overdiag <- data.frame(Year=1:inputs$global_parameters$time_horizon,
+                     NonCOPD=output_ex$n_COPD_by_ctime_severity[,1],
+                     Overdiagnosed=rowSums(output_ex$n_Overdiagnosed_by_ctime_sex))
+
+  overdiag$Proportion <- overdiag$Overdiagnosed/overdiag$NonCOPD
+
+  print(overdiag)
+
+  cat("The average proportion overdiagnosed from year", round(length(overdiag$Proportion)/2,0), "to", length(overdiag$Proportion), "is",
+      mean(overdiag$Proportion[(round(length(overdiag$Proportion)/2,0)):(length(overdiag$Proportion))]),"\n")
+
+  overdiag.plot <- tidyr::gather(data=overdiag, key="Variable", value="Number", c(NonCOPD, Overdiagnosed))
+
+  overdiag.plotted <- ggplot2::ggplot(overdiag.plot, aes(x=Year, y=Number, col=Variable)) +
+    geom_line() + geom_point() + expand_limits(y = 0) +
+    theme_bw() + ylab("Number of non-COPD subjects") + xlab("Years")
+
+  plot(overdiag.plotted)
+
+  cat("\n")
+
+  terminate_session()
+
+  }

src/model.WIP.cpp

@@ -574,6 +574,8 @@ struct input
   {
     double logit_p_diagnosis_by_sex[10][2];
     double p_hosp_diagnosis;
+    double logit_p_overdiagnosis_by_sex[9][2];
+    double p_correct_overdiagnosis;
     double p_case_detection;
     double min_cd_age;
     double min_cd_pack_years;
@@ -744,6 +746,8 @@ List Cget_inputs()
     Rcpp::Named("diagnosis")=Rcpp::List::create(
     Rcpp::Named("logit_p_diagnosis_by_sex")=AS_MATRIX_DOUBLE(input.diagnosis.logit_p_diagnosis_by_sex),
     Rcpp::Named("p_hosp_diagnosis")=input.diagnosis.p_hosp_diagnosis,
+    Rcpp::Named("logit_p_overdiagnosis_by_sex")=AS_MATRIX_DOUBLE(input.diagnosis.logit_p_overdiagnosis_by_sex),
+    Rcpp::Named("p_correct_overdiagnosis")=input.diagnosis.p_correct_overdiagnosis,
     Rcpp::Named("p_case_detection")=input.diagnosis.p_case_detection,
     Rcpp::Named("min_cd_age")=input.diagnosis.min_cd_age,
     Rcpp::Named("min_cd_pack_years")=input.diagnosis.min_cd_pack_years,
@@ -867,6 +871,8 @@ int Cset_input_var(std::string name, NumericVector value)
 
   if(name=="diagnosis$logit_p_diagnosis_by_sex") READ_R_MATRIX(value,input.diagnosis.logit_p_diagnosis_by_sex);
   if(name=="diagnosis$p_hosp_diagnosis") {input.diagnosis.p_hosp_diagnosis=value[0]; return(0);};
+  if(name=="diagnosis$logit_p_overdiagnosis_by_sex") READ_R_MATRIX(value,input.diagnosis.logit_p_overdiagnosis_by_sex);
+  if(name=="diagnosis$p_correct_overdiagnosis") {input.diagnosis.p_correct_overdiagnosis=value[0]; return(0);};
   if(name=="diagnosis$p_case_detection") {input.diagnosis.p_case_detection=value[0]; return(0);};
   if(name=="diagnosis$min_cd_age") {input.diagnosis.min_cd_age=value[0]; return(0);};
   if(name=="diagnosis$min_cd_pack_years") {input.diagnosis.min_cd_pack_years=value[0]; return(0);};
@@ -1009,11 +1015,12 @@ struct agent
   bool dyspnea;
   bool wheeze;
 
-  double gpvisits;
+  int gpvisits;
   double tmp_gpvisits_rate;
   int diagnosis;
   double p_hosp_diagnosis;
-  double case_detection;
+  double p_correct_overdiagnosis;
+  int case_detection;
   double p_case_detection;
   double min_cd_age;
   double min_cd_pack_years;
@@ -1411,6 +1418,7 @@ struct output_ex
   int n_COPD_by_ctime_severity[100][5]; //no COPD to GOLD 4;
   int n_COPD_by_age_sex[111][2];
   int n_Diagnosed_by_ctime_sex[1000][2];
+  int n_Overdiagnosed_by_ctime_sex[1000][2];
   int n_Diagnosed_by_ctime_severity[1000][5];
   int cumul_time_by_ctime_GOLD[100][5];
 #endif
@@ -1522,6 +1530,7 @@ List Cget_output_ex()
     out["n_COPD_by_ctime_severity"]=AS_MATRIX_INT_SIZE(output_ex.n_COPD_by_ctime_severity,input.global_parameters.time_horizon),
     out["n_COPD_by_age_sex"]=AS_MATRIX_INT(output_ex.n_COPD_by_age_sex),
     out["n_Diagnosed_by_ctime_sex"]=AS_MATRIX_INT_SIZE(output_ex.n_Diagnosed_by_ctime_sex,input.global_parameters.time_horizon),
+    out["n_Overdiagnosed_by_ctime_sex"]=AS_MATRIX_INT_SIZE(output_ex.n_Overdiagnosed_by_ctime_sex,input.global_parameters.time_horizon),
     out["n_Diagnosed_by_ctime_severity"]=AS_MATRIX_INT_SIZE(output_ex.n_Diagnosed_by_ctime_severity,input.global_parameters.time_horizon),
     out("cumul_time_by_ctime_GOLD")=AS_MATRIX_INT_SIZE(output_ex.cumul_time_by_ctime_GOLD,input.global_parameters.time_horizon),
 #endif
@@ -1630,7 +1639,8 @@ void update_output_ex(agent *ag)
       output_ex.n_COPD_by_ctime_age[time][age-1]+=((*ag).gold>0)*1;
       output_ex.n_COPD_by_ctime_severity[time][((*ag).gold)]+=1;
       output_ex.n_COPD_by_age_sex[age-1][(*ag).sex]+=1;
-      output_ex.n_Diagnosed_by_ctime_sex[time][(*ag).sex]+=((*ag).diagnosis>0)*1;
+      if((*ag).gold>0) output_ex.n_Diagnosed_by_ctime_sex[time][(*ag).sex]+=((*ag).diagnosis>0)*1;
+      if((*ag).gold==0) output_ex.n_Overdiagnosed_by_ctime_sex[time][(*ag).sex]+=((*ag).diagnosis>0)*1;
       if((*ag).gold>0) output_ex.n_Diagnosed_by_ctime_severity[time][(*ag).gold]+=((*ag).diagnosis>0)*1;
       if((*ag).local_time>0) output_ex.cumul_time_by_ctime_GOLD [time][((*ag).gold)]+=1;
 
@@ -1829,16 +1839,16 @@ double apply_case_detection(agent *ag)
  double update_diagnosis(agent *ag)
 {
 
-   if((*ag).diagnosis>0) return(0);
-
   double p_diagnosis = 0;
 
   if ((*ag).gpvisits!=0) {
 
   apply_case_detection(ag);
 
-  if((*ag).gold!=0)
-  {
+  if((*ag).gold!=0)  {
+
+    if((*ag).diagnosis>0) return(0);
+
     p_diagnosis = exp(input.diagnosis.logit_p_diagnosis_by_sex[0][(*ag).sex] +
       input.diagnosis.logit_p_diagnosis_by_sex[1][(*ag).sex]*((*ag).local_time+(*ag).age_at_creation) +
       input.diagnosis.logit_p_diagnosis_by_sex[2][(*ag).sex]*((*ag).smoking_status) +
@@ -1857,10 +1867,52 @@ double apply_case_detection(agent *ag)
             (*ag).medication_status=MED_CLASS_LAMA;
             medication_LPT(ag);
           }
-        }
+
+  } else {
+
+    double correct_overdiagnosis;
+    correct_overdiagnosis = input.diagnosis.p_correct_overdiagnosis;
+
+      if((*ag).diagnosis>0) {
+
+        if(rand_unif() < correct_overdiagnosis) {
+
+        (*ag).diagnosis = 0;
+
+         return(0);
       }
+
+    } else {
+
+      double p_overdiagnosis = 0;
+
+      p_overdiagnosis = exp(input.diagnosis.logit_p_overdiagnosis_by_sex[0][(*ag).sex] +
+        input.diagnosis.logit_p_overdiagnosis_by_sex[1][(*ag).sex]*((*ag).local_time+(*ag).age_at_creation) +
+        input.diagnosis.logit_p_overdiagnosis_by_sex[2][(*ag).sex]*((*ag).smoking_status) +
+        input.diagnosis.logit_p_overdiagnosis_by_sex[3][(*ag).sex]*((*ag).gpvisits) +
+        input.diagnosis.logit_p_overdiagnosis_by_sex[4][(*ag).sex]*((*ag).cough) +
+        input.diagnosis.logit_p_overdiagnosis_by_sex[5][(*ag).sex]*((*ag).phlegm) +
+        input.diagnosis.logit_p_overdiagnosis_by_sex[6][(*ag).sex]*((*ag).wheeze) +
+        input.diagnosis.logit_p_overdiagnosis_by_sex[7][(*ag).sex]*((*ag).dyspnea) +
+        input.diagnosis.logit_p_overdiagnosis_by_sex[8][(*ag).sex]*((*ag).case_detection));
+
+      p_overdiagnosis = p_overdiagnosis / (1 + p_overdiagnosis);
+
+      if (rand_unif() < p_overdiagnosis) {
+        (*ag).diagnosis = 1;
+
+      } else {
+
+        (*ag).diagnosis = 0;
+
+          }
+        }
+    }
+  }
   return(0);
- }
+}
+
+
 
 ////
 

tests/testthat/testDiagnosis.R

@@ -0,0 +1,36 @@
+library(testthat)
+library(epicR)
+library(dplyr)
+
+context("Diagnosis Tests")
+test_that("Proportion diagnosed in GOLD 1 and 2 is no greater than 40%, and Proporsed increased by Gold stage", {
+  init_session()
+  run()
+  inputs <- Cget_inputs()
+  output_ex <- Cget_output_ex()
+
+  prop <- data.frame(Year=1:inputs$global_parameters$time_horizon,
+                     output_ex$n_Diagnosed_by_ctime_severity/output_ex$n_COPD_by_ctime_severity)[,c(1,3,4,5,6)]
+  names(prop) <- c("Year","GOLD1","GOLD2","GOLD3","GOLD4")
+  prop <- prop[-1,]
+
+
+  #Total proportion of diagnosis of Gold 1 and Gold 2 
+  Gold_I_meanprop <- mean(prop$GOLD1[1:nrow(prop)])
+  Gold_II_meanprop <-mean(prop$GOLD2[1:nrow(prop)])
+
+
+  # Test total proportion on average years
+  expect_lt (Gold_I_meanprop,0.4)
+  expect_lt (Gold_II_meanprop,0.4)
+
+  # Test proportion year by year
+  Gold_I_prop <- prop %>% filter(GOLD1>0.4)
+  expect_equal (nrow(Gold_I_prop),0)
+
+  #Test diagnosis proportion increases by GOLD stages
+  Gold_meanprop= data.frame(colMeans(prop[,2:5]))
+  all(diff(Gold_meanprop) > 0)
+
+terminate_session()
+})