create_patients() now samples patient age and gender in addition to m…

…utation status

DESCRIPTION

@@ -23,7 +23,8 @@ Imports:
     MASS,
     R6,
     stats,
-    survival
+    survival,
+    truncnorm
 Remotes: InnovationValueInitiative/hesim    
 Suggests:
     covr,

R/data.R

@@ -15,6 +15,31 @@
 #' print(treatments)
 "treatments"
 
+#' Patient age
+#'
+#' Distribution of patient ages.
+#'
+#' @format A \code{data.table} with columns:
+#' \describe{
+#'   \item{age_cat}{Age category.}
+#'   \item{age_lower}{Lower bound of age category.}
+#'   \item{age_upper}{Upper bound of age category.}
+#'   \item{age_mid}{Midpoint of age category.}
+#'   \item{prop}{Proportion of patients in age category.}
+#' }
+#' Also contains the attributes:
+#' \describe{
+#' \item{mean}{Mean age.}
+#' \item{sd}{Standard deviation of age.}
+#' } 
+#' @source {\url{https://seer.cancer.gov/archive/csr/1975_2010/results_merged/topic_med_age.pdf}}
+#'
+#' @examples 
+#' print(age_dist)
+#' attr(age_dist, "mean")
+#' attr(age_dist, "sd")
+"age_dist"
+
 #' Multi-state NMA parameters
 #'
 #' Posterior distributions of the regression coefficients of the continuous 

R/model-setup.R

@@ -157,6 +157,9 @@ create_trans_mat <- function(object){
 #' 
 #' Create a data table of patients to model.
 #' @param n Number of patients to model.
+#' @param female_prop The proportion of patients that are female.
+#' @param age_mean Mean age. Based on sources cited in \code{\link{age_dist}}.
+#' @param age_sd Standard deviation of age. Based on sources cited in \code{\link{age_dist}}.
 #' @param mutation_prob The probability of a T790M mutation. The default value
 #' is based on Table 3 from the article by Ma et al. cited below.
 #' @examples
@@ -169,16 +172,29 @@ create_trans_mat <- function(object){
 #' \describe{
 #' \item{patient_id}{An integer from 1 to \code{n} denoting a unique patient.}
 #' \item{mutation}{1 if a patient has a T790M mutation and 0 otherwise.}
+#' \item{female}{1 if a patient is female and 0 otherwise.}
 #' }
 #' 
 #' @export
-create_patients <- function(n, mutation_prob = .52){
+create_patients <- function(n, female_prop = .45, 
+                            age_mean = 70.39, age_sd = 11.68,
+                            mutation_prob = .52){
   patient_id <- 1:n
-  n_mutations <- round(mutation_prob * n)
-  mutation <- c(rep(0, n - n_mutations),
-                rep(1, n_mutations))
+
+  # Age
+  age <- truncnorm::rtruncnorm(n, a = 0, b = 100, mean = age_mean, sd = age_sd)
+
+  # Gender
+  female <- stats::rbinom(n, 1, female_prop)
+
+  # Mutations
+  mutation <- stats::rbinom(n, 1, mutation_prob)
+
+  # Create dataset
   object <- data.table(patient_id = patient_id,
-                    mutation = mutation)
+                       female = female,
+                       age = age,
+                       mutation = mutation)
   setattr(object, "class", c("patients", "data.table", "data.frame"))
   return(object)
 }

data-raw/Make-data.R

@@ -5,6 +5,7 @@
 rm(list = ls())
 unlink("../data/*") # deletes all files in data directory
 source("treatments.R")
+source("age.R")
 source("params_mstate_nma.R")
 source("adverse_events.R")
 source("params_ae_nma.R")

data-raw/age-distribution.csv

@@ -0,0 +1 @@
+age_bot,age_top,percentage0,19,020,34,0.0235,44,145,54,7.455,64,21.765,74,33.475,84,26.885,Inf,9.4

data-raw/age.R

@@ -0,0 +1,51 @@
+rm(list = ls())
+library("data.table")
+library("ggplot2")
+theme_set(theme_bw())
+
+# Age distribution table
+# Source: https://seer.cancer.gov/archive/csr/1975_2010/results_merged/topic_med_age.pdf
+age_dist <- fread("age-distribution.csv")
+age_dist[, age_top := ifelse(is.infinite(age_top), 100, age_top)]
+age_dist[, age_cat := paste0(age_bot, " - ", age_top)]
+age_dist[, age_cat := ifelse(age_bot == 85, "85+", age_cat)]
+age_dist[, age_cat := factor(age_cat)]
+age_dist[, prop := percentage/100]
+age_dist[, age_mid := (age_top - age_bot)/2 + age_bot]
+age_mean <- weighted.mean(x = age_dist$age_mid, w = age_dist$prop)
+print(paste0("Mean age: ", age_mean))
+age_var <-  sum(age_dist$prop * (age_dist$age_mid - age_mean)^2)
+age_sd <- sqrt(age_var)
+attr(age_dist, "mean") <- age_mean
+attr(age_dist, "sd") <- age_sd
+
+# A barplot to examine the distribution
+p <- ggplot(age_dist, aes(x = age_cat, y = prop)) +
+  geom_bar(stat = "identity", fill = "lightblue", color = "blue") + xlab("Age category") +
+  ylab("Proportion")
+ggsave("figs/age-barplot.pdf", p, width = 5, height = 7)
+
+# Can we replicate the distribution using a normal distribution?
+sim_data <- vector(mode = "list", length = nrow(age_dist))
+for (i in 1:length(sim_data)){
+  sim_data[[i]] <- runif(1000 * age_dist[i]$prop,
+                        age_dist[i]$age_bot,
+                        age_dist[i]$age_top)
+}
+sim_data <- unlist(sim_data)
+sim_data <- data.table(age = c(sim_data))
+
+p <- ggplot(sim_data, aes(x = age)) +
+  geom_histogram(aes(y = ..density..), fill = "lightblue", color = "blue") +
+  stat_function(data = age_dist,
+                mapping = aes(x = age_mid),
+                fun = dnorm,
+                args = list(mean = age_mean, sd = age_sd),
+                col = "red")
+ggsave("figs/age-density.pdf", p, width = 5, height = 7)
+
+
+# Save 
+setcolorder(age_dist, c("age_cat", "age_bot", "age_top", "age_mid", "percentage", "prop"))
+age_dist[, percentage := NULL]
+save(age_dist, file = "../data/age_dist.rda", compress = "bzip2")