Malignant Lung Cancer Prediction

1 Overview

1.1 Objectives

This is a part of multiphase research project to understand non-traditional environmental contributors to lung cancer and how modifications to existing screening policies could potentially change outcomes. Phase 1 of this research project looked at at more typical contributors to lung cancer, like the number of packs of cigarettes smocked in a year, and explored other potential factors and moderators, like socioeconomic status. Phase 2 is where the crux of this study begins. In this document we look at a subset of some characteristics of the sample to create features for our modeling in which we attempt to understand the most important contributors to confirmed cases of malignant lung cancer and try to predict what happens if we modify screening policies in an effort to decrease the rate of undetected malignancy.

1.2 Data Source and Study Population

The data used for analysis and modeling in this study were obtained from electronic medical records (EMR) from a large academic medical center in Chicago using a standardized form collecting data on demographics, smoking history, lung cancer screening eligibility,

1.3 Methods

Because of what is at stake with this project, especially in terms of equitable community health, we will find and use algorithms that allow for the greatest transparency while sacrificing as little prediction power as possible. Typically this would be tree based models. After exploring the characteristics of our sample establish a baseline of performance with existing lung cancer screening policies (which currently at 20 packs of cigarettes a year, otherwise known as a packyear in our data), and then vary the binning of packyear to maximize the detection of malignant lung cancer. As a sanity check, the subgroup statistics in the beginning can be used to compare against the results of the models.

2 Patient Sample Characteristics

2.1 Percentage of cases in the total sample based on smoking status

lung_emr_raw %>%
  group_by(smokingstatus) %>%
  summarise(n = n()) %>%
  mutate(Percentage = scales::percent(n/sum(n), accuracy = 0.01)) %>%
  kable(col.names = c("Smoking Status",
                      "Count",
                      "Percentage"),align = "lllll", 
        caption = "Smoking Status of Total Sample", 
        digits = 2, 
        format.args = list(big.mark = ",")) %>%
  kable_classic(font_size = 15,
                full_width = F,
                html_font = "Cambria") %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")

Smoking Status of Total Sample

Smoking Status	Count	Percentage
current	2,252	31.29%
former	1,886	26.20%
never	2,217	30.80%
NA	843	11.71%

2.2 Percentage of cases that were current, former, or never smokers by race

lung_emr_raw %>%
  group_by(raceethnic, smokingstatus) %>%
  summarise(n = n(), .groups = "drop_last") %>%
  mutate(Percentage = scales::percent(n/sum(n), accuracy = 0.01)) %>%
  kable(col.names = c("Race/Ethnicity",
                      "Smoking Status",
                      "Count",
                      "Percentage within Race"),
        align = "lllll", 
        caption = "Smoking Status of Total Sample by Race", 
        digits = 2, 
        format.args = list(big.mark = ",")) %>%
  kable_classic(font_size = 15,
                full_width = F,
                html_font = "Cambria") %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")

Smoking Status of Total Sample by Race

Race/Ethnicity	Smoking Status	Count	Percentage within Race
Black	current	1,668	36.09%
Black	former	1,182	25.57%
Black	never	1,268	27.43%
Black	NA	504	10.90%
Latinx	current	271	16.74%
Latinx	former	446	27.55%
Latinx	never	692	42.74%
Latinx	NA	210	12.97%
White	current	313	32.71%
White	former	258	26.96%
White	never	257	26.85%
White	NA	129	13.48%

lung_emr_raw %>%
  group_by(raceethnic, smokingstatus) %>%
  summarise(n = n(), .groups = "drop_last") %>%
  mutate(Percentage = scales::percent(n/sum(n))) %>%
  gt() %>%
  gt_plt_bar_pct(column = Percentage, scaled = TRUE, fill = "azure2") %>%
  cols_align("center", contains("Percentage")) %>%
  cols_width(3 ~ px(125))

2.3 Percentage of cases that were current, former, or never smokers by gender

lung_emr_raw %>%
  group_by(gender, smokingstatus) %>%
  summarise(n = n(), .groups = 'drop_last') %>%
  mutate(Percentage = scales::percent(n/sum(n), accuracy = 0.01)) %>%
  kable(col.names = c("Gender",
                      "Smoking Status",
                      "Count",
                      "Percentage within Gender"),
        align = "lllll", 
        caption = "Smoking Status of Total Sample by Gender", 
        digits = 2, 
        format.args = list(big.mark = ",")) %>%
  kable_classic(font_size = 15,
                full_width = F,
                html_font = "Cambria") %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")

Smoking Status of Total Sample by Gender

Gender	Smoking Status	Count	Percentage within Gender
FEMALE	current	1,135	28.07%
FEMALE	former	941	23.27%
FEMALE	never	1,549	38.30%
FEMALE	NA	419	10.36%
MALE	current	1,117	35.45%
MALE	former	944	29.96%
MALE	never	667	21.17%
MALE	NA	423	13.42%
UNKNOWN	former	1	33.33%
UNKNOWN	never	1	33.33%
UNKNOWN	NA	1	33.33%

2.4 Percentage of cases that were current, former, or never smokers by homicide exposure

lung_emr_raw %>%
  group_by(smokingstatus, homicidegtmean2) %>%
  summarise(n = n(), .groups = 'drop_last') %>%
  mutate(Percentage = scales::percent(n/sum(n), accuracy = 0.01)) %>%
  kable(col.names = c("Smoking Status",
                      "Homicide Rate Exposure",
                      "Count",
                      "Percentage within Smoking Group"),
        align = "lllll", 
        caption = "Smoking Status of Total Sample by Homicide Rate Exposure", 
        digits = 2, 
        format.args = list(big.mark = ",")) %>%
  kable_classic(font_size = 15,
                full_width = F,
                html_font = "Cambria") %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")

Smoking Status of Total Sample by Homicide Rate Exposure

Smoking Status	Homicide Rate Exposure	Count	Percentage within Smoking Group
current	<mean	1,308	58.08%
current	>=mean	938	41.65%
current	NA	6	0.27%
former	<mean	1,169	61.98%
former	>=mean	717	38.02%
never	<mean	1,473	66.44%
never	>=mean	741	33.42%
never	NA	3	0.14%
NA	<mean	554	65.72%
NA	>=mean	287	34.05%
NA	NA	2	0.24%

3 Descriptive Statistics on Patients Who Smoked

3.1 Mean pack years smoked for current and former smokers

lung_emr_raw %>%
  filter(smokingstatus != "never") %>%
  group_by(smokingstatus) %>%
  summarise(mean_packyear = mean(packyear, na.rm = T), .groups = 'drop') %>%
  kable(col.names = c("Smoking Status",
                      "Mean"),
        align = "lllll", 
        caption = "Mean Packyear of Smokers", 
        digits = 2, 
        format.args = list(big.mark = ",")) %>%
  kable_classic(font_size = 15,
                full_width = F,
                html_font = "Cambria") %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")

Mean Packyear of Smokers

Smoking Status	Mean
current	16.36
former	22.27

3.2 Mean pack years smoked for current and former smokers by Race

lung_emr_raw %>%
  filter(smokingstatus != "never") %>%
  group_by(smokingstatus, raceethnic) %>%
  summarise(mean_packyear = mean(packyear, na.rm = T), .groups = 'drop') %>%
  kable(col.names = c("Smoking Status",
                      "Race/Ethnicity",
                      "Mean"),
        align = "lllll", 
        caption = "Mean Packyear of Smokers by Race", 
        digits = 2, 
        format.args = list(big.mark = ",")) %>%
  kable_classic(font_size = 15,
                full_width = F,
                html_font = "Cambria") %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")

Mean Packyear of Smokers by Race

Smoking Status	Race/Ethnicity	Mean
current	Black	15.24
current	Latinx	17.30
current	White	21.48
former	Black	22.62
former	Latinx	16.95
former	White	29.07

3.3 Mean pack years smoked for current and former smokers by Gender

lung_emr_raw %>%
  filter(smokingstatus != "never") %>%
  group_by(smokingstatus, gender) %>%
  summarise(mean_packyear = mean(packyear, na.rm = T), .groups = 'drop') %>%
  kable(col.names = c("Smoking Status",
                      "Gender",
                      "Mean Packyears"),
        align = "lllll", 
        caption = "Mean Packyear of Smokers by Gender", 
        digits = 2, 
        format.args = list(big.mark = ",")) %>%
  kable_classic(font_size = 15,
                full_width = F,
                html_font = "Cambria") %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")

Mean Packyear of Smokers by Gender

Smoking Status	Gender	Mean Packyears
current	FEMALE	14.94
current	MALE	17.75
former	FEMALE	22.85
former	MALE	21.62
former	UNKNOWN	NaN

3.4 Mean pack years smoked for current and former smokers by homicide exposure

lung_emr_raw %>%
  filter(smokingstatus != "never") %>%
  group_by(smokingstatus, homicidegtmean2) %>%
  summarise(mean_packyear = mean(packyear, na.rm = T), .groups = 'drop') %>%
  kable(col.names = c("Smoking Status",
                      "Homicide Rate Exposure",
                      "Mean Packyears"),
        align = "lllll", 
        caption = "Mean Packyear of Smokers by Race", 
        digits = 2, 
        format.args = list(big.mark = ",")) %>%
  kable_classic(font_size = 15,
                full_width = F,
                html_font = "Cambria") %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")

Mean Packyear of Smokers by Race

Smoking Status	Homicide Rate Exposure	Mean Packyears
current	<mean	17.83
current	>=mean	14.11
current	NA	52.00
former	<mean	21.61
former	>=mean	23.28

4 LDCT Eligibility Breakdowns

Eligibility criteria: Aged 50 to 80, 20 pack year smoking history, and no prior history of lung cancer

lung_emr_ldct =
  lung_emr_raw %>%
  mutate(
    age_elig = case_when(age >= 50 & age <= 80 ~ "Yes",
                         .default = "No"),
    packyear_elig = case_when(packyear20 == 1 ~ "Yes",
                              .default = "No"),
    prior_elig = case_when(historyca == 0 | historyca == "#NULL!" ~ "Yes",
                           .default = "No")
  )

4.1 Eligibility by Race

race_age_elig = 
  lung_emr_ldct %>%
  group_by(raceethnic, age_elig) %>%
  summarise(n = n(), .groups = "drop_last") %>%
  mutate(Percentage = scales::percent(n/sum(n), accuracy = 0.01)) 

race_packyear_elig = 
  lung_emr_ldct %>%
  group_by(raceethnic, packyear_elig) %>%
  summarise(n = n(), .groups = 'drop_last') %>%
  mutate(Percentage = scales::percent(n/sum(n), accuracy = 0.01))

race_prior_elig = 
  lung_emr_ldct %>%
  group_by(raceethnic, prior_elig) %>%
  summarise(n = n(), .groups = 'drop_last') %>%
  mutate(Percentage = scales::percent(n/sum(n), accuracy = 0.01)) 

race_elig = 
  bind_cols(c(race_age_elig, race_packyear_elig, race_prior_elig)) %>%
  select(!matches("race.*..[2-9]")) %>%
  rename_with(~ gsub("\\.\\.\\.", " ", .x)) 

## New names:
## • `raceethnic` -> `raceethnic...1`
## • `n` -> `n...3`
## • `Percentage` -> `Percentage...4`
## • `raceethnic` -> `raceethnic...5`
## • `n` -> `n...7`
## • `Percentage` -> `Percentage...8`
## • `raceethnic` -> `raceethnic...9`
## • `n` -> `n...11`
## • `Percentage` -> `Percentage...12`

race_elig %>%
  kable(col.names = c("Race/Ethnicity",
                      "Age Eligibility",
                      "Count",
                      "Percent",
                      "Packyear Eligibility",
                      "Count",
                      "Percent",
                      "No Prior Diagnosis Eligibility",
                      "Count",
                      "Percent"),
        align = "lllll", 
        caption = "Mean Packyear of Smokers by Race", 
        digits = 2, 
        format.args = list(big.mark = ",")) %>%
  kable_classic(font_size = 15,
                full_width = F,
                html_font = "Cambria") %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")

Mean Packyear of Smokers by Race

Race/Ethnicity	Age Eligibility	Count	Percent	Packyear Eligibility	Count	Percent	No Prior Diagnosis Eligibility	Count	Percent
Black	No	1,046	22.63%	No	4,384	94.85%	No	4,075	88.17%
Black	Yes	3,576	77.37%	Yes	238	5.15%	Yes	547	11.83%
Latinx	No	446	27.55%	No	1,575	97.28%	No	1,469	90.74%
Latinx	Yes	1,173	72.45%	Yes	44	2.72%	Yes	150	9.26%
White	No	204	21.32%	No	891	93.10%	No	836	87.36%
White	Yes	753	78.68%	Yes	66	6.90%	Yes	121	12.64%

4.2 Eligibility by Gender

gender_age_elig = 
  lung_emr_ldct %>%
  group_by(gender, age_elig) %>%
  summarise(n = n(), .groups = "drop_last") %>%
  mutate(Percentage = scales::percent(n/sum(n), accuracy = 0.01)) 

gender_packyear_elig = 
  lung_emr_ldct %>%
  group_by(gender, packyear_elig) %>%
  summarise(n = n(), .groups = 'drop_last') %>%
  bind_rows(tibble(gender = c("UNKNOWN"), 
                   packyear_elig = c("Yes"), 
                   n = c(0))) %>%
  mutate(Percentage = scales::percent(n/sum(n), accuracy = 0.01))

gender_prior_elig = 
  lung_emr_ldct %>%
  group_by(gender, prior_elig) %>%
  summarise(n = n(), .groups = 'drop_last') %>%
  bind_rows(tibble(gender = c("UNKNOWN"), 
                   prior_elig = c("Yes"), 
                   n = c(0))) %>%
  mutate(Percentage = scales::percent(n/sum(n), accuracy = 0.01)) 

gender_elig = 
  bind_cols(c(gender_age_elig, gender_packyear_elig, gender_prior_elig)) %>%
  select(!matches("gender.*..[2-9]"))

## New names:
## • `gender` -> `gender...1`
## • `n` -> `n...3`
## • `Percentage` -> `Percentage...4`
## • `gender` -> `gender...5`
## • `n` -> `n...7`
## • `Percentage` -> `Percentage...8`
## • `gender` -> `gender...9`
## • `n` -> `n...11`
## • `Percentage` -> `Percentage...12`

gender_elig %>%
  kable(col.names = c("Gender",
                      "Age Eligibility",
                      "Count",
                      "Percent",
                      "Packyear Eligibility",
                      "Count",
                      "Percent",
                      "No Prior Diagnosis Eligibility",
                      "Count",
                      "Percent"),
        align = "lllll", 
        caption = "Mean Packyear of Smokers by Race", 
        digits = 2, 
        format.args = list(big.mark = ",")) %>%
  kable_classic(font_size = 15,
                full_width = F,
                html_font = "Cambria") %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")

Mean Packyear of Smokers by Race

Gender	Age Eligibility	Count	Percent	Packyear Eligibility	Count	Percent	No Prior Diagnosis Eligibility	Count	Percent
FEMALE	No	1,049	25.94%	No	3,877	95.87%	No	3,617	89.44%
FEMALE	Yes	2,995	74.06%	Yes	167	4.13%	Yes	427	10.56%
MALE	No	645	20.47%	No	2,970	94.26%	No	2,760	87.59%
MALE	Yes	2,506	79.53%	Yes	181	5.74%	Yes	391	12.41%
UNKNOWN	No	2	66.67%	No	3	100.00%	No	3	100.00%
UNKNOWN	Yes	1	33.33%	Yes	0	0.00%	Yes	0	0.00%

4.3 Eligibility by Race/Gender

racegender_age_elig = 
  lung_emr_ldct %>%
  group_by(gender, raceethnic, age_elig) %>%
  summarise(n = n(), .groups = "drop_last") %>%
  mutate(Percentage = scales::percent(n/sum(n), accuracy = 0.01)) %>%
  relocate(gender, .after = raceethnic) %>%
  arrange(raceethnic)

racegender_packyear_elig = 
  lung_emr_ldct %>%
  group_by(gender, raceethnic, packyear_elig) %>%
  summarise(n = n(), .groups = 'drop_last') %>%
  bind_rows(tibble(gender = c("UNKNOWN"), 
                   raceethnic = c("White"),
                   packyear_elig = c("No"), 
                   n = c(0))) %>%
  mutate(Percentage = scales::percent(n/sum(n), accuracy = 0.01)) %>%
  relocate(gender, .after = raceethnic) %>%
  arrange(raceethnic)

racegender_prior_elig = 
  lung_emr_ldct %>%
  group_by(gender, raceethnic, prior_elig) %>%
  summarise(n = n(), .groups = 'drop_last') %>%
  bind_rows(tibble(gender = c("UNKNOWN"), 
                   raceethnic = c("White"),
                   prior_elig = c("No"), 
                   n = c(0))) %>%
  mutate(Percentage = scales::percent(n/sum(n), accuracy = 0.01)) %>%
  relocate(gender, .after = raceethnic) %>%
  arrange(raceethnic)

racegender_elig = 
  bind_cols(c(racegender_age_elig, racegender_packyear_elig, racegender_prior_elig)) %>%
  select(!matches("(raceethnic|gender)...([3-9]|[0-9]{2})"))

## New names:
## • `raceethnic` -> `raceethnic...1`
## • `gender` -> `gender...2`
## • `n` -> `n...4`
## • `Percentage` -> `Percentage...5`
## • `raceethnic` -> `raceethnic...6`
## • `gender` -> `gender...7`
## • `n` -> `n...9`
## • `Percentage` -> `Percentage...10`
## • `raceethnic` -> `raceethnic...11`
## • `gender` -> `gender...12`
## • `n` -> `n...14`
## • `Percentage` -> `Percentage...15`

racegender_elig %>%
  kable(col.names = c("Race/Ethnicity",
                      "Gender",
                      "Age Eligibility",
                      "Count",
                      "Percent within Race/Gender",
                      "Packyear Eligibility",
                      "Count",
                      "Percent within Race/Gender",
                      "No Prior Diagnosis Eligibility",
                      "Count",
                      "Percent within Race/Gender"),
        align = "lllll", 
        caption = "Mean Packyear of Smokers by Race", 
        digits = 2, 
        format.args = list(big.mark = ",")) %>%
  kable_classic(font_size = 15,
                full_width = F,
                html_font = "Cambria") %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "left")

Mean Packyear of Smokers by Race

Race/Ethnicity	Gender	Age Eligibility	Count	Percent within Race/Gender	Packyear Eligibility	Count	Percent within Race/Gender	No Prior Diagnosis Eligibility	Count	Percent within Race/Gender
Black	FEMALE	No	711	25.32%	No	2,681	95.48%	No	2,500	89.03%
Black	FEMALE	Yes	2,097	74.68%	Yes	127	4.52%	Yes	308	10.97%
Black	MALE	No	334	18.43%	No	1,701	93.87%	No	1,573	86.81%
Black	MALE	Yes	1,478	81.57%	Yes	111	6.13%	Yes	239	13.19%
Black	UNKNOWN	No	1	50.00%	No	2	100.00%	No	2	100.00%
Black	UNKNOWN	Yes	1	50.00%	No	793	98.51%	No	737	91.55%
Latinx	FEMALE	No	241	29.94%	Yes	12	1.49%	Yes	68	8.45%
Latinx	FEMALE	Yes	564	70.06%	No	781	96.06%	No	731	89.91%
Latinx	MALE	No	204	25.09%	Yes	32	3.94%	Yes	82	10.09%
Latinx	MALE	Yes	609	74.91%	No	1	100.00%	No	1	100.00%
Latinx	UNKNOWN	No	1	100.00%	No	403	93.50%	No	380	88.17%
White	FEMALE	No	97	22.51%	Yes	28	6.50%	Yes	51	11.83%
White	FEMALE	Yes	334	77.49%	No	488	92.78%	No	456	86.69%
White	MALE	No	107	20.34%	Yes	38	7.22%	Yes	70	13.31%
White	MALE	Yes	419	79.66%	No	0	NA	No	0	NA

5 Modeling and Predicting Malignanc Lung Cancer Prevalence

In a predictive model including all patients who do not meet eligibility for LDCT, what is the overall predictive ability of following combined data to predict lung cancer diagnosis:

1. Age
2. Race/ethnicity
3. Gender
4. Smoking status (never vs former/current)
5. Packyear binned by 5 up to current packyear eligibility (see below)

levels(lung_emr$smokingstatus_packyear)

## [1] "[  0,  5) packs/year" "[  5, 10) packs/year" "[ 10, 15) packs/year"
## [4] "[ 15, 20) packs/year" "[ 20,250] packs/year" "never smoked"

5.1 Initial CHAID Model Predicting Malignant Lung Cancer

The CHAID model below is the result of optimizing over some of the following hyperparameters using 70% of the sample selected at random:

• alpha2: Level of significance used for merging of predictor categories (step 2).

• alpha3: If set to a positive value \(< 1\), level of significance used for the the splitting of former merged categories of the predictor (step 3). Otherwise, step 3 is omitted (the default).

• alpha4: Level of significance used for splitting of a node in the most significant predictor (step 5).

• minsplit: Number of observations in splitted response at which no further split is desired.

• minbucket: Minimum number of observations in terminal nodes.

• minprob: Mininimum frequency of observations in terminal nodes.

# custom function necessary for usage in caret train function
cgpCHAID <- list(label = "CGP CHAID",
                 library = "CHAID",
                 loop = NULL,
                 type = c("Classification"),
                 parameters = data.frame(parameter = c('alpha2',
                                                       'alpha3',
                                                       'alpha4',
                                                       'minsplit', 
                                                       'minbucket', 
                                                       'minprob', 
                                                       'maxheight'),
                                         class = rep('numeric', 7),
                                         label = c('alpha2',
                                                   'alpha3',
                                                   'alpha4', 
                                                   'Numb obs in response where no further split', 
                                                   "Minimum numb obs in terminal nodes", 
                                                   "Minimum freq of obs in terminal nodes.",
                                                   "Maximum height for the tree")
                 ),
                 grid = function(x, y, len = NULL, search = "grid") {
                   if(search == "grid") {
                     out <- data.frame(minsplit = c(20,30),
                                       minbucket = 7,
                                       minprob = c(0.05,0.01),
                                       maxheight = -1)
                   } else {
                     out <- data.frame(minsplit = c(20,30),
                                       minbucket = 7,
                                       minprob = c(0.05,0.01),
                                       maxheight = -1)
                   }
                   out
                 },
                 fit = function(x, y, wts, param, lev, last, classProbs, ...) {
                   dat <- if(is.data.frame(x)) x else as.data.frame(x)
                   dat$.outcome <- y
                   theDots <- list(...)
                   if(any(names(theDots) == "control")) {
                     theDots$control$alpha2 <- param$alpha2
                     theDots$control$alpha3 <- param$alpha3
                     theDots$control$alpha4 <- param$alpha4
                     theDots$control$minsplit <- param$minsplit
                     theDots$control$minbucket <- param$minbucket
                     theDots$control$minprob <- param$minprob
                     theDots$control$maxheight <- param$maxheight
                     ctl <- theDots$control
                     theDots$control <- NULL
                   } else ctl <- chaid_control(alpha2 <- param$alpha2,
                                               alpha3 <- param$alpha3,
                                               alpha4 <- param$alpha4,
                                               minsplit = param$minsplit,
                                               minbucket = param$minbucket,
                                               minprob = param$minprob,
                                               maxheight = param$maxheight)
                   ## pass in any model weights
                   if(!is.null(wts)) theDots$weights <- wts
                   modelArgs <- c(
                     list(
                       formula = as.formula(".outcome ~ ."),
                       data = dat,
                       control = ctl),
                     theDots)
                   out <- do.call(CHAID::chaid, modelArgs)
                   out
                 },
                 predict = function(modelFit, newdata, submodels = NULL) {
                   if(!is.data.frame(newdata)) newdata <- as.data.frame(newdata)
                   predict(modelFit, newdata)
                 },
                 prob = function(modelFit, newdata, submodels = NULL) {
                   if(!is.data.frame(newdata)) newdata <- as.data.frame(newdata)
                   predict(modelFit, newdata, type = "prob")
                 },
                 levels = function(x) x$obsLevels,
                 predictors = function(x, surrogate = TRUE, ...) {
                   predictors(terms(x))
                 },
                 tags = c('Tree-Based Model', "Implicit Feature Selection", "Two Class Only", "Accepts Case Weights"),
                 sort = function(x) x[order(-x$alpha2, -x$alpha3, -x$alpha4, 
                                            -x$minsplit, -x$minbucket, 
                                            -x$minprob, -x$maxheight),])

A visualization of the resulting optimized model is shown below:

lung_mod_fit_0 = readRDS("lung_model_0.rds")
lung_mod_fit_1 = readRDS("lung_model_1.rds")
lung_mod_fit_2 = readRDS("lung_model_2.rds")
lung_mod_fit_3 = readRDS("lung_model_3.rds")
lung_mod_fit_3_alt = readRDS("lung_model_3_alt.rds")

# plot(lung_mod_fit_0$finalModel, ep_args = list(justmin = 20), type='simple')
# plot(lung_mod_fit_0$finalModel, ep_args = list(justmin = 20))
# plot(lung_mod_fit_1$finalModel)
# plot(lung_mod_fit_2$finalModel)
plot(lung_mod_fit_3$finalModel, ep_args = list(justmin = 20), type='simple')

# plot(lung_mod_fit_3_alt$finalModel, ep_args = list(justmin = 20), type='simple')
# plot(lung_mod_fit_3_alt$finalModel, ep_args = list(justmin = 20))

lung_mod_var_imp = 
  # must specify namespace otherwise it will error
  bind_rows(partykit::varimp(lung_mod_fit_3_alt$finalModel)) %>%
  pivot_longer(cols = everything(),
               names_to = "variable",
               values_to = "importance") %>%
  mutate(importance_perc = round(importance/sum(importance), digits=3),
         variable = fct_reorder(variable, importance_perc, .desc = F))

lung_mod_var_imp_fig =
  plot_ly(
    lung_mod_var_imp,
    y = ~variable,
    x = ~importance_perc*100,
    text = ~scales::percent(importance_perc, accuracy = 0.01),
    type = 'bar',
    orientation = 'h'
  ) %>%
  layout(title = "Model Variable Importance",
         xaxis = list(title = "Importance Relative Percentage (%)"))

lung_mod_var_imp_fig

The predictions on the remaining \(30\%\) of the set, held out for testing, is shown below as a confusion matrix where the upper right represents False Negatives and lower left represents False Positives.

get_chaid_stats = function(mod_fit, 
                           file_name, 
                           predict_type='prob', 
                           test_set, 
                           caret=caret, 
                           confusionMatrix=confusionMatrix) {
  
  mod_prob = predict(mod_fit, test_set, type=predict_type)
  mod_pred = prediction(mod_prob[,2], test_set$malignanto)
  
  # find optimal cutoff probability to use to maximize accuracy
  mod_cutoff = optimalCutoff(test_set$malignanto, mod_prob[,1])[1]
  
  # column 1 of the model probabilities represent predictions of no malignant cancer
  mod_results = ifelse(mod_prob[,1] >= mod_cutoff, 0, 1)
  
  # create confusion matrix
  results = caret::confusionMatrix(as.factor(mod_results),
                                   as.factor(test_set$malignantoactuals))
  
  saveRDS(object = results, file = paste(toString(file_name), "results.rds", sep = "_"))
  
  show(results)
  
  auc_perf_1 = performance(mod_pred, "tpr", "fpr")
  plot(auc_perf_1, col = "blue", main = "ROC Curve")
  
  results
}

get_model_table =
  function(dataset, title, mod_results) {
    df1 =
      dataset %>%
      bind_cols(prediction = mod_results) %>%
      group_by(malignanto, ldct_elig) %>%
      summarise(ldct_n = n(), .groups = 'drop_last') %>%
      mutate(ldct_percentage = scales::percent(ldct_n / sum(ldct_n), accuracy = 0.01)) %>%
      # count() %>%
      arrange(desc(malignanto), ldct_elig)
    
    df2 =
      dataset %>%
      bind_cols(prediction = mod_results) %>%
      group_by(malignanto, prediction) %>%
      summarise(prediction_n = n(), .groups = 'drop_last') %>%
      mutate(prediction_percentage = scales::percent(prediction_n / sum(prediction_n), accuracy = 0.01)) %>%
      # count() %>%
      arrange(desc(malignanto), desc(prediction))
    
    combined = bind_cols(c(df1, df2), .name_repair = "unique") %>%
      select(!matches("[a-z]+...([3-9]|[0-9]{2})"))
    # show(combined)
    
    show(
      combined %>%
        kable(
          align = "lllll",
          caption = paste("Eligibility Versus Prediction for", title),
          digits = 2,
          format.args = list(big.mark = ",")
        ) %>%
        kable_classic(
          font_size = 15,
          full_width = F,
          html_font = "Cambria"
        ) %>%
        kable_styling(
          bootstrap_options = "striped",
          full_width = F,
          position = "left"
        )
    )
  }

plot_var_imp = function(model_fit, plot_title) {
  mod_var_imp =
    bind_rows(partykit::varimp(model_fit$finalModel)) %>%
    pivot_longer(cols = everything(),
                 names_to = "variable",
                 values_to = "importance") %>%
    mutate(
      importance_perc = importance / sum(importance),
      variable = fct_reorder(variable, importance_perc, .desc = F)
    )
  
  mod_var_imp_fig =
    plot_ly(
      mod_var_imp,
      y = ~ variable,
      x = ~ importance_perc * 100,
      text = ~ scales::percent(importance_perc, accuracy = 0.01),
      type = 'bar',
      orientation = 'h'
    ) %>%
    layout(title = plot_title,
           xaxis = list(title = "Importance Relative Percentage (%)"))
  
  show(mod_var_imp_fig)
  
  mod_var_imp_fig
}

lung_mod_fit_3_alt_results = get_chaid_stats(lung_mod_fit_3_alt$finalModel, "chaid_base_model_results", 'prob', lung_emr_test)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 1404  106
##          1  539  111
##                                           
##                Accuracy : 0.7014          
##                  95% CI : (0.6816, 0.7206)
##     No Information Rate : 0.8995          
##     P-Value [Acc > NIR] : 1               
##                                           
##                   Kappa : 0.1241          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.7226          
##             Specificity : 0.5115          
##          Pos Pred Value : 0.9298          
##          Neg Pred Value : 0.1708          
##              Prevalence : 0.8995          
##          Detection Rate : 0.6500          
##    Detection Prevalence : 0.6991          
##       Balanced Accuracy : 0.6171          
##                                           
##        'Positive' Class : 0               
## 

lung_mod_fit_3_alt_acc = scales::percent(lung_mod_fit_3_alt_results$overall[["Accuracy"]], accuracy = 0.01)

The resulting model has a lung_mod_fit_3_alt_acc accuracy with higher rate of false negatives than false positives.

5.2 CHAID Model with Homicide Rate Predicting Malignant Lung Cancer

Does adding exposure to neighborhood violence (Homicide rate > = Mean vs. < Mean) increase the predictive ability of the model?

a.  Expectation = Yes

Remodeling with the homicide rate per 100k variable, we get the following model:

lung_mod_hom_fit = readRDS("lung_model_hom.rds")
plot(lung_mod_hom_fit$finalModel, ep_args = list(justmin = 15), type='simple')

plot(lung_mod_hom_fit$finalModel, ep_args = list(justmin = 15))

Applying the model to the \(30\%\) held out test set results in the following:

lung_mod_hom_fit_results = get_chaid_stats(lung_mod_hom_fit$finalModel, 'chaid_hom_model_results', 'prob', lung_emr_test)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 1709  157
##          1  234   60
##                                          
##                Accuracy : 0.819          
##                  95% CI : (0.8021, 0.835)
##     No Information Rate : 0.8995         
##     P-Value [Acc > NIR] : 1.0000000      
##                                          
##                   Kappa : 0.1348         
##                                          
##  Mcnemar's Test P-Value : 0.0001213      
##                                          
##             Sensitivity : 0.8796         
##             Specificity : 0.2765         
##          Pos Pred Value : 0.9159         
##          Neg Pred Value : 0.2041         
##              Prevalence : 0.8995         
##          Detection Rate : 0.7912         
##    Detection Prevalence : 0.8639         
##       Balanced Accuracy : 0.5780         
##                                          
##        'Positive' Class : 0              
## 

lung_mod_hom_fit_acc = scales::percent(lung_mod_hom_fit_results$overall[["Accuracy"]], accuracy = 0.01)

This model has a lung_mod_hom_fit_acc accuracy with, again, a higher rate of false negatives than false positives. However, this model produces false negatives at lower rate than the previous model without the homicide rate included. Additionally there is an approximately lung_mod_hom_fit_acc - lung_mod_fit_3_alt_acc gain in accuracy.

lung_mod_hom_var_imp = 
  bind_rows(partykit::varimp(lung_mod_hom_fit$finalModel)) %>%
  pivot_longer(cols = everything(),
               names_to = "variable",
               values_to = "importance") %>%
  mutate(importance_perc = importance/sum(importance),
         variable = fct_reorder(variable, importance_perc, .desc = F))

lung_mod_hom_var_imp_fig =
  plot_ly(
    lung_mod_hom_var_imp,
    y = ~variable,
    x = ~importance_perc*100,
    text = ~scales::percent(importance_perc, accuracy = 0.01),
    type = 'bar',
    orientation = 'h'
  ) %>%
  layout(title = "Model with Homicide Rate Variable Importance",
         xaxis = list(title = "Importance Relative Percentage (%)"))

lung_mod_hom_var_imp_fig

5.3 CHAID Model with Narrower Bins

Using intervals of size 1 between 5 and 15, with ordinal categories as below:

levels(lung_emr$smokingstatus_packyear_1)

##  [1] "[  0,  5)"    "[  5,  6)"    "[  6,  7)"    "[  7,  8)"    "[  8,  9)"   
##  [6] "[  9, 10)"    "[ 10, 11)"    "[ 11, 12)"    "[ 12, 13)"    "[ 13, 14)"   
## [11] "[ 14, 15)"    "[ 15, 20)"    "[ 20,250]"    "never smoked"

tg_hom = 
  expand.grid(alpha2 = c(0.05),
              # alpha3 = c(0.05, 0.10, 0.20),
              alpha3 = -1,
              alpha4 = c(0.05, 0.10, 0.20),
              minsplit = c(10, 30),
              minbucket = c(10, 30),
              minprob = c(0.01, 0.1),
              maxheight = c(2, 3, 4)) %>% 
  filter(alpha2 > alpha3)

set.seed(825)

lung_mod_narrow_fit1 = train(
  x = lung_emr_train %>% select(all_of(x_cols_special_bin1)),
  y = lung_emr_train$malignanto,
  method = cgpCHAID,
  metric = 'ROC',
  tuneGrid = tg_hom,
  trControl = trainControl(classProbs = T,
                           summaryFunction = twoClassSummary)
)

saveRDS(lung_mod_narrow_fit1, "lung_model_chaid_narrow_bins1.rds")

lung_mod_narrow_fit1 = readRDS("lung_model_chaid_narrow_bins1.rds")
plot(lung_mod_narrow_fit1$finalModel, 
     ep_args = list(justmin = 20), 
     type='simple',
     main = "Lung Cancer Model With Narrower Bins on Total Sample")

print("Lung Cancer Model With Narrower Bins on Total Sample")

## [1] "Lung Cancer Model With Narrower Bins on Total Sample"

get_chaid_stats(
  lung_mod_narrow_fit1$finalModel,
  test_set = lung_emr_test,
  file_name = "chaid_narrow_bin1_results",
  predict_type = 'prob'
)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 1696  154
##          1  247   63
##                                           
##                Accuracy : 0.8144          
##                  95% CI : (0.7973, 0.8305)
##     No Information Rate : 0.8995          
##     P-Value [Acc > NIR] : 1               
##                                           
##                   Kappa : 0.1371          
##                                           
##  Mcnemar's Test P-Value : 4.343e-06       
##                                           
##             Sensitivity : 0.8729          
##             Specificity : 0.2903          
##          Pos Pred Value : 0.9168          
##          Neg Pred Value : 0.2032          
##              Prevalence : 0.8995          
##          Detection Rate : 0.7852          
##    Detection Prevalence : 0.8565          
##       Balanced Accuracy : 0.5816          
##                                           
##        'Positive' Class : 0               
## 

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 1696  154
##          1  247   63
##                                           
##                Accuracy : 0.8144          
##                  95% CI : (0.7973, 0.8305)
##     No Information Rate : 0.8995          
##     P-Value [Acc > NIR] : 1               
##                                           
##                   Kappa : 0.1371          
##                                           
##  Mcnemar's Test P-Value : 4.343e-06       
##                                           
##             Sensitivity : 0.8729          
##             Specificity : 0.2903          
##          Pos Pred Value : 0.9168          
##          Neg Pred Value : 0.2032          
##              Prevalence : 0.8995          
##          Detection Rate : 0.7852          
##    Detection Prevalence : 0.8565          
##       Balanced Accuracy : 0.5816          
##                                           
##        'Positive' Class : 0               
## 

Using intervals of size 5 between 0 and 10 and then size 1 between 10 and 15, with ordinal categories as below:

levels(lung_emr$smokingstatus_packyear_2)

##  [1] "[  0,  5)"    "[  5, 10)"    "[ 10, 11)"    "[ 11, 12)"    "[ 12, 13)"   
##  [6] "[ 13, 14)"    "[ 14, 15)"    "[ 15, 20)"    "[ 20,250]"    "never smoked"

tg_hom = 
  expand.grid(alpha2 = c(0.05),
              # alpha3 = c(0.05, 0.10, 0.20),
              alpha3 = -1,
              alpha4 = c(0.05, 0.10, 0.20),
              minsplit = c(10, 30),
              minbucket = c(10, 30),
              minprob = c(0.01, 0.1),
              maxheight = c(2, 3, 4)) %>% 
  filter(alpha2 > alpha3)

set.seed(825)

lung_mod_narrow_fit2 = train(
  x = lung_emr_train %>% select(all_of(x_cols_special_bin2)),
  y = lung_emr_train$malignanto,
  method = cgpCHAID,
  metric = 'ROC',
  tuneGrid = tg_hom,
  trControl = trainControl(classProbs = T,
                           summaryFunction = twoClassSummary)
)

saveRDS(lung_mod_narrow_fit2, "lung_model_chaid_narrow_bins2.rds")

lung_mod_narrow_fit2 = readRDS("lung_model_chaid_narrow_bins2.rds")
plot(lung_mod_narrow_fit2$finalModel, ep_args = list(justmin = 20), type='simple')

get_chaid_stats(
  lung_mod_narrow_fit2$finalModel,
  test_set = lung_emr_test,
  file_name = "chaid_narrow_bin2_results",
  predict_type = 'prob'
)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 1696  154
##          1  247   63
##                                           
##                Accuracy : 0.8144          
##                  95% CI : (0.7973, 0.8305)
##     No Information Rate : 0.8995          
##     P-Value [Acc > NIR] : 1               
##                                           
##                   Kappa : 0.1371          
##                                           
##  Mcnemar's Test P-Value : 4.343e-06       
##                                           
##             Sensitivity : 0.8729          
##             Specificity : 0.2903          
##          Pos Pred Value : 0.9168          
##          Neg Pred Value : 0.2032          
##              Prevalence : 0.8995          
##          Detection Rate : 0.7852          
##    Detection Prevalence : 0.8565          
##       Balanced Accuracy : 0.5816          
##                                           
##        'Positive' Class : 0               
## 

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 1696  154
##          1  247   63
##                                           
##                Accuracy : 0.8144          
##                  95% CI : (0.7973, 0.8305)
##     No Information Rate : 0.8995          
##     P-Value [Acc > NIR] : 1               
##                                           
##                   Kappa : 0.1371          
##                                           
##  Mcnemar's Test P-Value : 4.343e-06       
##                                           
##             Sensitivity : 0.8729          
##             Specificity : 0.2903          
##          Pos Pred Value : 0.9168          
##          Neg Pred Value : 0.2032          
##              Prevalence : 0.8995          
##          Detection Rate : 0.7852          
##    Detection Prevalence : 0.8565          
##       Balanced Accuracy : 0.5816          
##                                           
##        'Positive' Class : 0               
## 

5.4 CHAID Model on Subsets

chaid_subset_train = function(subset, name, cols = c(x_cols_combined, "homicidegtmean2")) {
  set.seed(88)
  lung_emr_sub_split = initial_split(subset, prop = 0.7)
  lung_emr_sub_train = training(lung_emr_sub_split)
  lung_emr_sub_test = testing(lung_emr_sub_split)
  
  tg_sub_hom = 
    expand.grid(alpha2 = c(0.05, 0.10),
                # alpha3 = c(0.05, 0.10, 0.20),
                alpha3 = -1,
                alpha4 = c(0.05, 0.10, 0.20),
                minsplit = c(10, 30),
                minbucket = c(10, 30),
                minprob = c(0.01, 0.1),
                maxheight = c(2, 3, 4)) %>% 
    filter(alpha2 > alpha3)
  
  lung_mod_hom_sub_fit = train(
    x = lung_emr_sub_train %>% select(cols),
    y = lung_emr_sub_train$malignanto,
    method = cgpCHAID,
    metric = 'ROC',
    tuneGrid = tg_sub_hom,
    trControl = trainControl(classProbs = T,
                             summaryFunction = twoClassSummary)
  )
  
  file_name = paste(name, ".rds", sep = "")
  
  show(paste("Saving model to ", file_name))
  
  saveRDS(object = lung_mod_hom_sub_fit, file = file_name)
  
  show(lung_mod_hom_sub_fit$bestTune)
  
  get_chaid_stats(lung_mod_hom_sub_fit$finalModel, 
                  test_set = lung_emr_sub_test, 
                  file_name = toString(name),
                  predict_type = 'prob')
}

5.4.1 Model on Race Subsets Using Size 5 Bins

chaid_by_race_model_fits =
  lapply(
    X = chaid_by_race_model_rds,
    FUN = function(file) {
      rds = readRDS(file)
      title = paste("Race/Ethnicity Subset with Bin Size 5:",
                    toTitleCase(str_extract(file, "_([a-z]+)(_results)?.rds", group = 1)))
      if (str_detect(file, "narrow", negate = T)) {
        if (str_detect(file, "results")) {
          print(title)
          show(rds)
        } else {
          plot(
            rds$finalModel,
            ep_args = list(justmin = 20),
            type = 'simple',
            main = title
          )
          plot_var_imp(
            model_fit = rds,
            plot_title = paste(title, "Model Variable Importance", sep = " ")
          )
        }
      }
    }
  )

## [1] "Race/Ethnicity Subset with Bin Size 5: Black"
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 974  86
##          1 263  64
##                                          
##                Accuracy : 0.7484         
##                  95% CI : (0.7247, 0.771)
##     No Information Rate : 0.8919         
##     P-Value [Acc > NIR] : 1              
##                                          
##                   Kappa : 0.141          
##                                          
##  Mcnemar's Test P-Value : <2e-16         
##                                          
##             Sensitivity : 0.7874         
##             Specificity : 0.4267         
##          Pos Pred Value : 0.9189         
##          Neg Pred Value : 0.1957         
##              Prevalence : 0.8919         
##          Detection Rate : 0.7022         
##    Detection Prevalence : 0.7642         
##       Balanced Accuracy : 0.6070         
##                                          
##        'Positive' Class : 0              
## 

## [1] "Race/Ethnicity Subset with Bin Size 5: Latinx"
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 448  30
##          1   8   0
##                                           
##                Accuracy : 0.9218          
##                  95% CI : (0.8943, 0.9441)
##     No Information Rate : 0.9383          
##     P-Value [Acc > NIR] : 0.9412190       
##                                           
##                   Kappa : -0.0267         
##                                           
##  Mcnemar's Test P-Value : 0.0006577       
##                                           
##             Sensitivity : 0.9825          
##             Specificity : 0.0000          
##          Pos Pred Value : 0.9372          
##          Neg Pred Value : 0.0000          
##              Prevalence : 0.9383          
##          Detection Rate : 0.9218          
##    Detection Prevalence : 0.9835          
##       Balanced Accuracy : 0.4912          
##                                           
##        'Positive' Class : 0               
## 

## [1] "Race/Ethnicity Subset with Bin Size 5: White"
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0  77   2
##          1 181  28
##                                           
##                Accuracy : 0.3646          
##                  95% CI : (0.3089, 0.4231)
##     No Information Rate : 0.8958          
##     P-Value [Acc > NIR] : 1               
##                                           
##                   Kappa : 0.0637          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.2984          
##             Specificity : 0.9333          
##          Pos Pred Value : 0.9747          
##          Neg Pred Value : 0.1340          
##              Prevalence : 0.8958          
##          Detection Rate : 0.2674          
##    Detection Prevalence : 0.2743          
##       Balanced Accuracy : 0.6159          
##                                           
##        'Positive' Class : 0               
## 

5.4.2 Model on Race Subsets Using Custom Sized Bins

apply(
  lung_emr %>% distinct(raceethnic),
  MARGIN = 1,
  FUN = function(race) {
    chaid_subset_train(
      lung_emr %>% filter(raceethnic == race),
      paste("lung_model_chaid",
            "narrow_bins",
            "race_subset", 
            toString(tolower(race)), 
            sep = "_"),
      c(x_cols_special_bin2, "homicidegtmean2")
    )
  }
)

chaid_by_race_narrow_bins_model_fits =
  lapply(
    X = chaid_by_race_narrow_bins_model_rds,
    FUN = function(file) {
      rds = readRDS(file)
      race = toTitleCase(str_extract(file, "_([a-z]+)(_results)?.rds", group = 1))
      title = paste("Race/Ethnicity Subset with Narrower Bins:",
                    toTitleCase(race))
      if (str_detect(file, "results")){
        print(title)
        show(rds)
      } else {
        mod_fit = rds$finalModel
        test_set = lung_emr %>% filter(raceethnic == race)
        
        mod_prob = predict(mod_fit, test_set, type="prob")
        mod_pred = prediction(mod_prob[,2], test_set$malignanto)
        
        # find optimal cutoff probability to use to maximize accuracy
        mod_cutoff = optimalCutoff(test_set$malignanto, mod_prob[,1])[1]
        
        # column 1 of the model probabilities represent predictions of no malignant cancer
        mod_results = ifelse(mod_prob[,1] >= mod_cutoff, 0, 1)
        
        get_model_table(
          lung_emr %>%
            filter(raceethnic == race) %>%
            select(
              c(
                x_cols_special_bin2,
                "homicidegtmean2",
                contains("ldct"),
                malignanto
              )
            ),
          paste(race, "Subset"),
          mod_results
        ) 
        
        plot(
          rds$finalModel,
          ep_args = list(justmin = 20),
          type = 'simple',
          main = title
        )
        plot_var_imp(
          model_fit = rds,
          plot_title = paste(title, "Model Variable Importance", sep = " ")
        )
      }
    }
  )

## New names:
## • `malignanto` -> `malignanto...1`
## • `malignanto` -> `malignanto...5`

## <table class=" lightable-classic table table-striped" style="font-size: 15px; font-family: Cambria; width: auto !important; margin-left: auto; margin-right: auto; width: auto !important; ">
## <caption style="font-size: initial !important;">Eligibility Versus Prediction for Black Subset</caption>
##  <thead>
##   <tr>
##    <th style="text-align:left;"> malignanto...1 </th>
##    <th style="text-align:left;"> ldct_elig </th>
##    <th style="text-align:left;"> ldct_n </th>
##    <th style="text-align:left;"> ldct_percentage </th>
##    <th style="text-align:left;"> prediction </th>
##    <th style="text-align:left;"> prediction_n </th>
##    <th style="text-align:left;"> prediction_percentage </th>
##   </tr>
##  </thead>
## <tbody>
##   <tr>
##    <td style="text-align:left;"> yes </td>
##    <td style="text-align:left;"> No </td>
##    <td style="text-align:left;"> 484 </td>
##    <td style="text-align:left;"> 95.09% </td>
##    <td style="text-align:left;"> 1 </td>
##    <td style="text-align:left;"> 202 </td>
##    <td style="text-align:left;"> 39.69% </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> yes </td>
##    <td style="text-align:left;"> Yes </td>
##    <td style="text-align:left;"> 25 </td>
##    <td style="text-align:left;"> 4.91% </td>
##    <td style="text-align:left;"> 0 </td>
##    <td style="text-align:left;"> 307 </td>
##    <td style="text-align:left;"> 60.31% </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> no </td>
##    <td style="text-align:left;"> No </td>
##    <td style="text-align:left;"> 3,938 </td>
##    <td style="text-align:left;"> 95.75% </td>
##    <td style="text-align:left;"> 1 </td>
##    <td style="text-align:left;"> 811 </td>
##    <td style="text-align:left;"> 19.72% </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> no </td>
##    <td style="text-align:left;"> Yes </td>
##    <td style="text-align:left;"> 175 </td>
##    <td style="text-align:left;"> 4.25% </td>
##    <td style="text-align:left;"> 0 </td>
##    <td style="text-align:left;"> 3,302 </td>
##    <td style="text-align:left;"> 80.28% </td>
##   </tr>
## </tbody>
## </table>

## [1] "Race/Ethnicity Subset with Narrower Bins: Black"
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 974  86
##          1 263  64
##                                          
##                Accuracy : 0.7484         
##                  95% CI : (0.7247, 0.771)
##     No Information Rate : 0.8919         
##     P-Value [Acc > NIR] : 1              
##                                          
##                   Kappa : 0.141          
##                                          
##  Mcnemar's Test P-Value : <2e-16         
##                                          
##             Sensitivity : 0.7874         
##             Specificity : 0.4267         
##          Pos Pred Value : 0.9189         
##          Neg Pred Value : 0.1957         
##              Prevalence : 0.8919         
##          Detection Rate : 0.7022         
##    Detection Prevalence : 0.7642         
##       Balanced Accuracy : 0.6070         
##                                          
##        'Positive' Class : 0              
## 

## New names:
## • `malignanto` -> `malignanto...1`
## • `malignanto` -> `malignanto...5`

## <table class=" lightable-classic table table-striped" style="font-size: 15px; font-family: Cambria; width: auto !important; margin-left: auto; margin-right: auto; width: auto !important; ">
## <caption style="font-size: initial !important;">Eligibility Versus Prediction for Latinx Subset</caption>
##  <thead>
##   <tr>
##    <th style="text-align:left;"> malignanto...1 </th>
##    <th style="text-align:left;"> ldct_elig </th>
##    <th style="text-align:left;"> ldct_n </th>
##    <th style="text-align:left;"> ldct_percentage </th>
##    <th style="text-align:left;"> prediction </th>
##    <th style="text-align:left;"> prediction_n </th>
##    <th style="text-align:left;"> prediction_percentage </th>
##   </tr>
##  </thead>
## <tbody>
##   <tr>
##    <td style="text-align:left;"> yes </td>
##    <td style="text-align:left;"> No </td>
##    <td style="text-align:left;"> 96 </td>
##    <td style="text-align:left;"> 96.97% </td>
##    <td style="text-align:left;"> 1 </td>
##    <td style="text-align:left;"> 8 </td>
##    <td style="text-align:left;"> 8.08% </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> yes </td>
##    <td style="text-align:left;"> Yes </td>
##    <td style="text-align:left;"> 3 </td>
##    <td style="text-align:left;"> 3.03% </td>
##    <td style="text-align:left;"> 0 </td>
##    <td style="text-align:left;"> 91 </td>
##    <td style="text-align:left;"> 91.92% </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> no </td>
##    <td style="text-align:left;"> No </td>
##    <td style="text-align:left;"> 1,490 </td>
##    <td style="text-align:left;"> 98.03% </td>
##    <td style="text-align:left;"> 1 </td>
##    <td style="text-align:left;"> 34 </td>
##    <td style="text-align:left;"> 2.24% </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> no </td>
##    <td style="text-align:left;"> Yes </td>
##    <td style="text-align:left;"> 30 </td>
##    <td style="text-align:left;"> 1.97% </td>
##    <td style="text-align:left;"> 0 </td>
##    <td style="text-align:left;"> 1,486 </td>
##    <td style="text-align:left;"> 97.76% </td>
##   </tr>
## </tbody>
## </table>

## [1] "Race/Ethnicity Subset with Narrower Bins: Latinx"
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 446  30
##          1  10   0
##                                           
##                Accuracy : 0.9177          
##                  95% CI : (0.8896, 0.9406)
##     No Information Rate : 0.9383          
##     P-Value [Acc > NIR] : 0.971979        
##                                           
##                   Kappa : -0.0318         
##                                           
##  Mcnemar's Test P-Value : 0.002663        
##                                           
##             Sensitivity : 0.9781          
##             Specificity : 0.0000          
##          Pos Pred Value : 0.9370          
##          Neg Pred Value : 0.0000          
##              Prevalence : 0.9383          
##          Detection Rate : 0.9177          
##    Detection Prevalence : 0.9794          
##       Balanced Accuracy : 0.4890          
##                                           
##        'Positive' Class : 0               
## 

## New names:
## • `malignanto` -> `malignanto...1`
## • `malignanto` -> `malignanto...5`

## <table class=" lightable-classic table table-striped" style="font-size: 15px; font-family: Cambria; width: auto !important; margin-left: auto; margin-right: auto; width: auto !important; ">
## <caption style="font-size: initial !important;">Eligibility Versus Prediction for White Subset</caption>
##  <thead>
##   <tr>
##    <th style="text-align:left;"> malignanto...1 </th>
##    <th style="text-align:left;"> ldct_elig </th>
##    <th style="text-align:left;"> ldct_n </th>
##    <th style="text-align:left;"> ldct_percentage </th>
##    <th style="text-align:left;"> prediction </th>
##    <th style="text-align:left;"> prediction_n </th>
##    <th style="text-align:left;"> prediction_percentage </th>
##   </tr>
##  </thead>
## <tbody>
##   <tr>
##    <td style="text-align:left;"> yes </td>
##    <td style="text-align:left;"> No </td>
##    <td style="text-align:left;"> 97 </td>
##    <td style="text-align:left;"> 97.98% </td>
##    <td style="text-align:left;"> 1 </td>
##    <td style="text-align:left;"> 64 </td>
##    <td style="text-align:left;"> 64.65% </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> yes </td>
##    <td style="text-align:left;"> Yes </td>
##    <td style="text-align:left;"> 2 </td>
##    <td style="text-align:left;"> 2.02% </td>
##    <td style="text-align:left;"> 0 </td>
##    <td style="text-align:left;"> 35 </td>
##    <td style="text-align:left;"> 35.35% </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> no </td>
##    <td style="text-align:left;"> No </td>
##    <td style="text-align:left;"> 812 </td>
##    <td style="text-align:left;"> 94.64% </td>
##    <td style="text-align:left;"> 1 </td>
##    <td style="text-align:left;"> 493 </td>
##    <td style="text-align:left;"> 57.46% </td>
##   </tr>
##   <tr>
##    <td style="text-align:left;"> no </td>
##    <td style="text-align:left;"> Yes </td>
##    <td style="text-align:left;"> 46 </td>
##    <td style="text-align:left;"> 5.36% </td>
##    <td style="text-align:left;"> 0 </td>
##    <td style="text-align:left;"> 365 </td>
##    <td style="text-align:left;"> 42.54% </td>
##   </tr>
## </tbody>
## </table>

## [1] "Race/Ethnicity Subset with Narrower Bins: White"
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 114  14
##          1 144  16
##                                           
##                Accuracy : 0.4514          
##                  95% CI : (0.3929, 0.5108)
##     No Information Rate : 0.8958          
##     P-Value [Acc > NIR] : 1               
##                                           
##                   Kappa : -0.0085         
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.4419          
##             Specificity : 0.5333          
##          Pos Pred Value : 0.8906          
##          Neg Pred Value : 0.1000          
##              Prevalence : 0.8958          
##          Detection Rate : 0.3958          
##    Detection Prevalence : 0.4444          
##       Balanced Accuracy : 0.4876          
##                                           
##        'Positive' Class : 0               
## 

race_mod_fits = c()
for (race in lung_emr %>% distinct(raceethnic)) {
  set.seed(8)
  lung_mod_hom_sub_fit =
    chaid_subset_train(
      lung_emr %>% filter(raceethnic == race),
      toString(tolower(race)),
      c(x_cols_combined, "homicidegtmean2")
    )
  append(race_mod_fits, lung_mod_hom_sub_fit)
}
length(race_mod_fits)

set.seed(825)

lung_emr_black_split = initial_split(
  lung_emr %>% filter(raceethnic == "Black"), 
  prop=0.7) 
lung_emr_black_train = training(lung_emr_black_split)
lung_emr_black_test = testing(lung_emr_black_split)

tg_hom = 
  expand.grid(alpha2 = c(0.01, 0.05, 0.10),
              # alpha3 = c(0.05, 0.10, 0.20),
              alpha3 = -1,
              alpha4 = c(0.05, 0.10, 0.20),
              minsplit = c(10, 30, 100),
              minbucket = c(10, 30, 100),
              minprob = c(0.01, 0.1),
              maxheight = c(2, 3, 4)) %>% 
  filter(alpha2 > alpha3)

x_cols_hom = c(x_cols_combined, "homicidegtmean2")

set.seed(825)

lung_mod_hom_black_fit = train(
  x = lung_emr_black_train %>% select(all_of(x_cols_hom)),
  y = lung_emr_black_train$malignanto,
  method = cgpCHAID,
  metric = 'ROC',
  tuneGrid = tg_hom,
  trControl = trainControl(classProbs = T,
                           summaryFunction = twoClassSummary)
)

saveRDS(lung_mod_hom_black_fit, "lung_model_hom_black.rds")

lung_mod_hom_black_fit$bestTune

5.5 Generalized Linear Model, Logistic Regression

Let’s compare the CHAID trees to a logistic regression model. Here we have an interaction term with smokingstatus and packyear since packyear would not be applicable without a “not never” status.

lung_emr_lr = 
  lung_emr %>%
  mutate(smokingstatus = as.factor(case_when(smokingstatus != "never" ~ 1,
                                   .default = 0)),
         packyear = case_when(is.na(packyear) & smokingstatus == 0 ~ 0,
                              .default = packyear))

lung_emr_lr_split = initial_split(lung_emr_lr, 
                               prop=0.7) 
lung_emr_lr_train = training(lung_emr_lr_split)
lung_emr_lr_test = testing(lung_emr_lr_split)

# show(summary(lung_emr_lr$smokingstatus))

lung_mod_glm = glm(
  malignanto ~ agecat + gender + raceethnic + smokingstatus * packyear + homicidegtmean2,
  data = lung_emr_lr_train,
  family = "binomial"
)

summary(lung_mod_glm)

## 
## Call:
## glm(formula = malignanto ~ agecat + gender + raceethnic + smokingstatus * 
##     packyear + homicidegtmean2, family = "binomial", data = lung_emr_lr_train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3859  -0.5568  -0.4072  -0.2726   3.1039  
## 
## Coefficients: (1 not defined because of singularities)
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              -3.388645   0.154065 -21.995  < 2e-16 ***
## agecat.L                  1.109462   0.162683   6.820 9.12e-12 ***
## agecat.Q                 -0.247364   0.121573  -2.035  0.04188 *  
## genderMALE               -0.002696   0.097786  -0.028  0.97801    
## genderUNKNOWN            -9.745831 221.541247  -0.044  0.96491    
## raceethnicLatinx         -0.461841   0.152911  -3.020  0.00253 ** 
## raceethnicWhite          -0.016659   0.156549  -0.106  0.91525    
## smokingstatus1            1.009121   0.140372   7.189 6.53e-13 ***
## packyear                  0.008041   0.003893   2.065  0.03888 *  
## homicidegtmean2.L         0.231492   0.076816   3.014  0.00258 ** 
## smokingstatus1:packyear         NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3276.4  on 5037  degrees of freedom
## Residual deviance: 3066.6  on 5028  degrees of freedom
## AIC: 3086.6
## 
## Number of Fisher Scoring iterations: 11

lung_mod_glm_prob = predict(lung_mod_glm, lung_emr_lr_test, type = 'link')
lung_mod_glm_pred = prediction(lung_mod_glm_prob, lung_emr_test$malignanto)

# find optimal cutoff probability to use to maximize accuracy
lung_mod_glm_cutoff = optimalCutoff(lung_emr_test$malignanto, lung_mod_glm_prob)[1]

# column 1 of the model probabilities represent predictions of no malignant cancer
lung_mod_glm_results = ifelse(lung_mod_glm_prob >= lung_mod_glm_cutoff, 0, 1)

# create confusion matrix
lung_mod_glm_conf = caret::confusionMatrix(as.factor(lung_mod_glm_results),
                                 as.factor(lung_emr_test$malignantoactuals))
show(lung_mod_glm_conf)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 1942  217
##          1    1    0
##                                           
##                Accuracy : 0.8991          
##                  95% CI : (0.8856, 0.9115)
##     No Information Rate : 0.8995          
##     P-Value [Acc > NIR] : 0.5465          
##                                           
##                   Kappa : -9e-04          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.9995          
##             Specificity : 0.0000          
##          Pos Pred Value : 0.8995          
##          Neg Pred Value : 0.0000          
##              Prevalence : 0.8995          
##          Detection Rate : 0.8991          
##    Detection Prevalence : 0.9995          
##       Balanced Accuracy : 0.4997          
##                                           
##        'Positive' Class : 0               
## 

lung_mod_glm_auc_perf = performance(lung_mod_glm_pred, "tpr", "fpr")
plot(lung_mod_glm_auc_perf, col = "blue", main = "ROC Curve")

Modeling only the part of the sample that smoked at any time in their lives produces the following model:

lung_emr_smokers_split = initial_split(
  lung_emr %>%  filter(smokingstatus != "never"), 
  prop=0.7
) 
lung_emr_smokers_train = training(lung_emr_smokers_split)
lung_emr_smokers_test = testing(lung_emr_smokers_split)

lung_mod_glm_smokers = glm(
  malignanto ~ agecat + gender + raceethnic + packyear + homicidegtmean2,
  data = lung_emr_smokers_train,
  family = "binomial"
)

summary(lung_mod_glm_smokers)

## 
## Call:
## glm(formula = malignanto ~ agecat + gender + raceethnic + packyear + 
##     homicidegtmean2, family = "binomial", data = lung_emr_smokers_train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.2919  -0.5593  -0.4731  -0.3285   2.6975  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -2.53318    0.12487 -20.287  < 2e-16 ***
## agecat.L            1.13337    0.18594   6.095 1.09e-09 ***
## agecat.Q           -0.20435    0.13757  -1.485 0.137417    
## genderMALE         -0.02743    0.10604  -0.259 0.795884    
## genderUNKNOWN     -10.61529  365.20689  -0.029 0.976812    
## raceethnicLatinx   -0.43241    0.17321  -2.496 0.012543 *  
## raceethnicWhite     0.17289    0.16132   1.072 0.283839    
## packyear            0.01537    0.00455   3.378 0.000729 ***
## homicidegtmean2.L   0.28555    0.08469   3.372 0.000747 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2588.6  on 3485  degrees of freedom
## Residual deviance: 2480.0  on 3477  degrees of freedom
## AIC: 2498
## 
## Number of Fisher Scoring iterations: 12

lung_mod_glm_smokers_prob = predict(lung_mod_glm_smokers, lung_emr_smokers_test, type = 'link')
lung_mod_glm_smokers_pred = prediction(lung_mod_glm_smokers_prob, lung_emr_smokers_test$malignanto)

# find optimal cutoff probability to use to maximize accuracy
lung_mod_glm_smokers_cutoff = optimalCutoff(lung_emr_smokers_test$malignanto, lung_mod_glm_smokers_prob)[1]

# column 1 of the model probabilities represent predictions of no malignant cancer
lung_mod_glm_smokers_results = ifelse(lung_mod_glm_smokers_prob >= lung_mod_glm_smokers_cutoff, 0, 1)

# create confusion matrix
lung_mod_glm_smokers_results = caret::confusionMatrix(as.factor(lung_mod_glm_smokers_results),
                                 as.factor(lung_emr_smokers_test$malignantoactuals))
show(lung_mod_glm_smokers_results)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 1303  191
##          1    1    0
##                                           
##                Accuracy : 0.8716          
##                  95% CI : (0.8535, 0.8881)
##     No Information Rate : 0.8722          
##     P-Value [Acc > NIR] : 0.55            
##                                           
##                   Kappa : -0.0013         
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.9992          
##             Specificity : 0.0000          
##          Pos Pred Value : 0.8722          
##          Neg Pred Value : 0.0000          
##              Prevalence : 0.8722          
##          Detection Rate : 0.8716          
##    Detection Prevalence : 0.9993          
##       Balanced Accuracy : 0.4996          
##                                           
##        'Positive' Class : 0               
## 

lung_mod_glm_smokers_auc_perf = performance(lung_mod_glm_smokers_pred, "tpr", "fpr")
plot(lung_mod_glm_smokers_auc_perf, col = "blue", main = "ROC Curve")

5.6 Other Tree Based Models

Supposing that are binning of the packyear variable into 5-packyear segments might be imprecise, a tree model that can take in continuous variables might produce different results in terms of what packyear split might produce the predictions on malignant lung cancer.

5.6.1 Gradient Boosted Trees

lung_emr_factors =
  lung_emr %>%  
  filter(smokingstatus != "never") %>%
  select(all_of(x_cols_packyear_hom) & !smokingstatus, malignanto) %>%
  select(where(is.factor))

lung_emr_factors1 = lung_emr_factors %>% select(!c(agecat, homicidegtmean2, malignanto))

lung_emr_dummy1 = dummyVars(" ~ .", data = lung_emr_factors1)

lung_emr_dummy_cols1 = predict(lung_emr_dummy1, newdata = lung_emr_factors1)

lung_emr_dummy =
  bind_cols(
    lung_emr_dummy_cols1,
    class2ind(lung_emr_factors$agecat),
    class2ind(lung_emr_factors$homicidegtmean2),
    lung_emr %>%  
      filter(smokingstatus != "never") %>%
      select(packyear),
    lung_emr_factors %>% 
      select(malignanto) %>%
      mutate(malignanto = as.factor(as.numeric(malignanto) - 1))
  )

str(lung_emr_dummy)

## 'data.frame':    4981 obs. of  13 variables:
##  $ gender.FEMALE    : num  1 1 1 1 1 0 0 1 1 1 ...
##  $ gender.MALE      : num  0 0 0 0 0 1 1 0 0 0 ...
##  $ gender.UNKNOWN   : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ raceethnic.Black : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ raceethnic.Latinx: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ raceethnic.White : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ 40-50            : num  0 0 0 0 1 0 1 0 0 0 ...
##  $ 50-60            : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ 60+              : num  1 1 1 1 0 1 0 1 1 1 ...
##  $ <mean            : num  1 0 0 1 0 1 0 0 0 1 ...
##  $ >=mean           : num  0 1 1 0 1 0 1 1 1 0 ...
##  $ packyear         : num  9.19 21 3.25 13.32 11.6 ...
##  $ malignanto       : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...

lung_emr_smokers_dummy_split = initial_split(
  lung_emr_dummy %>% select(malignanto, !malignanto), 
  prop=0.7
) 
lung_emr_smokers_dummy_train = training(lung_emr_smokers_dummy_split)
lung_emr_smokers_dummy_test = testing(lung_emr_smokers_dummy_split)

tg_xgbt = expand.grid(nrounds = c(50, 75, 100, 150, 200, 250),
                      max_depth = c(2, 3, 4, 5),
                      eta = c(0.05, 0.1, 0.15),
                      gamma = c(0),
                      colsample_bytree = c(0.7),
                      min_child_weight = c(5),
                      subsample = c(0.6))

set.seed(88)

lung_mod_xgbt_fit = train(
  formula(lung_emr_smokers_dummy_train),
  data = lung_emr_smokers_dummy_train,
  method = 'xgbTree',
  tuneGrid = tg_xgbt,
  verbosity = 0
)

saveRDS(lung_mod_xgbt_fit, "lung_mod_xgbt_fit.rds")

lung_mod_xgbt_fit = readRDS("lung_mod_xgbt_fit.rds")
xgb.plot.tree(model = lung_mod_xgbt_fit$finalModel, trees = 8:10)

## [19:06:42] WARNING: src/learner.cc:553: 
##   If you are loading a serialized model (like pickle in Python, RDS in R) generated by
##   older XGBoost, please export the model by calling `Booster.save_model` from that version
##   first, then load it back in current version. See:
## 
##     https://xgboost.readthedocs.io/en/latest/tutorials/saving_model.html
## 
##   for more details about differences between saving model and serializing.

# using DiagrammeR to render

lung_mod_xgbt_fit = readRDS("lung_mod_xgbt_fit.rds")
xgb_graph = xgb.plot.tree(model = lung_mod_xgbt_fit$finalModel, render = F)

## [19:06:42] WARNING: src/learner.cc:553: 
##   If you are loading a serialized model (like pickle in Python, RDS in R) generated by
##   older XGBoost, please export the model by calling `Booster.save_model` from that version
##   first, then load it back in current version. See:
## 
##     https://xgboost.readthedocs.io/en/latest/tutorials/saving_model.html
## 
##   for more details about differences between saving model and serializing.

xgb_graph %>% 
  select_edges_by_node_id(xgb_graph %>%
  select_nodes(conditions = data == "packyear") %>%
  get_node_df_ws() %>%
  select(id)) %>%
  get_edges()

##   [1] "1->2"     "2->4"     "6->7"     "7->9"     "11->12"   "12->14"  
##   [7] "13->16"   "18->19"   "19->21"   "23->24"   "24->26"   "28->29"  
##  [13] "29->31"   "30->33"   "35->36"   "36->38"   "37->40"   "42->43"  
##  [19] "43->45"   "44->47"   "49->50"   "50->52"   "51->54"   "56->57"  
##  [25] "57->59"   "61->62"   "62->64"   "66->67"   "67->69"   "71->72"  
##  [31] "72->74"   "73->76"   "78->79"   "79->81"   "80->83"   "85->86"  
##  [37] "86->88"   "87->90"   "92->93"   "93->95"   "94->97"   "99->100" 
##  [43] "100->102" "101->104" "106->107" "107->109" "108->111" "113->114"
##  [49] "114->116" "118->119" "119->121" "120->123" "125->126" "126->128"
##  [55] "127->130" "132->133" "133->135" "134->137" "139->140" "140->142"
##  [61] "141->144" "146->147" "147->149" "151->152" "152->154" "153->156"
##  [67] "158->159" "159->161" "160->163" "165->166" "166->168" "167->170"
##  [73] "172->173" "173->175" "177->178" "178->180" "179->182" "184->185"
##  [79] "185->187" "186->189" "191->192" "192->194" "193->196" "198->199"
##  [85] "199->201" "203->204" "204->206" "208->209" "209->211" "213->214"
##  [91] "214->216" "215->218" "220->221" "221->223" "222->225" "227->228"
##  [97] "228->230" "232->233" "233->235" "234->237" "239->240" "240->242"
## [103] "244->245" "245->247" "246->249" "251->252" "252->254" "256->257"
## [109] "257->259" "261->262" "264->265" "265->267" "269->270" "270->272"
## [115] "274->275" "275->277" "279->280" "280->282" "281->284" "286->287"
## [121] "287->289" "291->292" "292->294" "296->297" "297->299" "1->3"    
## [127] "2->5"     "6->8"     "7->10"    "11->13"   "12->15"   "13->17"  
## [133] "18->20"   "19->22"   "23->25"   "24->27"   "28->30"   "29->32"  
## [139] "30->34"   "35->37"   "36->39"   "37->41"   "42->44"   "43->46"  
## [145] "44->48"   "49->51"   "50->53"   "51->55"   "56->58"   "57->60"  
## [151] "61->63"   "62->65"   "66->68"   "67->70"   "71->73"   "72->75"  
## [157] "73->77"   "78->80"   "79->82"   "80->84"   "85->87"   "86->89"  
## [163] "87->91"   "92->94"   "93->96"   "94->98"   "99->101"  "100->103"
## [169] "101->105" "106->108" "107->110" "108->112" "113->115" "114->117"
## [175] "118->120" "119->122" "120->124" "125->127" "126->129" "127->131"
## [181] "132->134" "133->136" "134->138" "139->141" "140->143" "141->145"
## [187] "146->148" "147->150" "151->153" "152->155" "153->157" "158->160"
## [193] "159->162" "160->164" "165->167" "166->169" "167->171" "172->174"
## [199] "173->176" "177->179" "178->181" "179->183" "184->186" "185->188"
## [205] "186->190" "191->193" "192->195" "193->197" "198->200" "199->202"
## [211] "203->205" "204->207" "208->210" "209->212" "213->215" "214->217"
## [217] "215->219" "220->222" "221->224" "222->226" "227->229" "228->231"
## [223] "232->234" "233->236" "234->238" "239->241" "240->243" "244->246"
## [229] "245->248" "246->250" "251->253" "252->255" "256->258" "257->260"
## [235] "261->263" "264->266" "265->268" "269->271" "270->273" "274->276"
## [241] "275->278" "279->281" "280->283" "281->285" "286->288" "287->290"
## [247] "291->293" "292->295" "296->298" "297->300"