MDC_PredictiveModeling.Rmd

---
title: "Multidisciplinary breast program in New York City - Predictive Modeling" 
author: "Anni Liu"
date: "`r format(Sys.Date(), '%B %d, %Y')`"
knit: knitautomator::knit_filename 
output:
  word_document:
    fig_caption: no
    fig_height: 4
    fig_width: 4
    highlight: null
    toc: yes
    reference_docx: MDCAsian_manuscript_style_V0.docx
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  options(scipen = 999, digits = 4),
  cache = TRUE,
  error = FALSE,
  message = FALSE,
  warning = FALSE,
  tidy.opts = list(width.cutoff = 60),
  tidy = TRUE,
  fig.width = 12,
  fig.height = 8
)
```


```{r include=FALSE}
## Attach libraries and functions
easypackages::libraries("tidyverse", "readxl", "gtsummary", "flextable", "stringi", "fst", "broom", "MASS") |> suppressPackageStartupMessages()

"%_%" <- function(m, n) paste0(m, "_", n)
"%0%" <- function(m, n) paste0(m, n)
"%+%" <- function(m) paste0(m, collapse = "+")
alter.glm.rownames <- function(model) {
  var.names <- names(coef(model))
  for (v in 1:length(var.names)) {
    i <- var.names[v]
    if (i == "(Intercept)") {
      var.names[v] <- "Intercept"
    } else if (i == "Age.Dx") {
      var.names[v] <- "Age at Diagnosis"
    } else if (grepl("*.[a-z][A-Z]", i)) {
      if (grepl("^Race.A", i)) {
        var.names[v] <- substring(i, regexpr("*.[a-z][A-Z]", i) + 2) %0% " (vs Non-Asian)"
        if (grepl(":(.*)", i)) {
          var.names[v] <- gsub(pattern = "Insure.type|MMGD.2", replacement = "", x = substring(i, regexpr("*.[a-z][A-Z]", i) + 2))
        }
      } else if (grepl("^RaceE", i)) {
        var.names[v] <- substring(i, regexpr("*.[0-9][A-Z]", i) + 2) %0% " (vs Non-Hispanic White)"
        if (grepl(":(.*)", i)) { # *.:Insure.type|*.:MMGD.2
          var.names[v] <- gsub(pattern = "Insure.type|MMGD.2", replacement = "", x = substring(i, regexpr("*.[0-9][A-Z]", i) + 2))
        }
      } else if (grepl("^Ins", i)) {
        var.names[v] <- substring(i, regexpr("*.[a-z][A-Z]", i) + 2) %0% " (vs Private)"
      } else if (grepl("^Trip", i)) {
        var.names[v] <- "Triple Negative Yes (vs Triple Negative No)"
      } else if (grepl("^Sub", i)) {
        var.names[v] <- substring(i, regexpr("*.[a-z][A-Z]", i) + 2) %0% " (vs HR+, HER2-)"
      }
    } else if (i == "MMGSDNo") {
      var.names[v] <- "Mammo Screen-Detected No (vs Mammo Screen-Detected Yes)"
    } else if (i == "MMGOYes") {
      var.names[v] <- "Mammo-Occult Yes (vs Mammo-Occult No)"
    } else if (i == "MRISDYes") {
      var.names[v] <- "MRI Screen-Detected Yes (vs MRI Screen-Detected No)"
    } else if (grepl("*.[0-9][A-Z]", i)) {
        if (grepl("^Hos", i)) {
          var.names[v] <- "NYPQ (vs WCM)"
        } else if (grepl("^MMGD", i)) {
          var.names[v] <- substring(i, regexpr("*.[0-9][A-Z]", i) + 2) %0% " (vs Heterogenously dense)"
        } else {
          var.names[v] <- "HER2 Positive (vs HER2 Negative)" }
    } else if (i == "ERNegative") {
      var.names[v] <- "ER Negative (vs ER Positive)"
    } else if (i == "PRNegative") {
      var.names[v] <- "PR Negative (vs PR Positive)"
    }
  }
  return(var.names)
}
```


# Data preparation
Notice that the categorical levels `Self Pay` and `None` are removed, because each of two levels has only 1 patient.
```{r load_data, include=FALSE}
data.full <- read_fst(path = "../data/derived/2023May14_dat_MDCAsian.RData")

data.full$Early.stage <- factor(data.full$Early.stage, levels = c("Non-early", "Early"), labels = c("No", "Yes"))
data.full$TripleNegative <- with(data.full, ifelse(Subtype.2 == "Receptor positive", "No", "Yes")) |> factor(levels = c("No", "Yes"))
data.full$HER2 <- factor(data.full$HER2, levels = c("Negative", "Positive"))

# Combine Southeast Asian and Asian Other of few people showing early stage to one category
# Characteristics  Early      Non-early
# South Asian	     72 (2.8%)	46 (4.9%)
# Southeast Asian	 71 (2.7%)	24 (2.6%)
data.full$RaceEthnicity.4 <- with(data.full, case_when(RaceEthnicity.4 == "Southeast Asian" ~ "Asian Other", TRUE ~ RaceEthnicity.4)) |> factor(levels = c("Non-Hispanic White", "East Asian", "South Asian", "Asian Other", "Non-Hispanic Black", "Hispanic", "Other"))

# Drop patients with self pay or no insurance
data.full <- data.full[-which(data.full$Insure.type %in% c("Self Pay", "None")), ]

response <- "Early.stage"
predictors <- c("Age.Dx", "Race.Asian", "RaceEthnicity.4", "Hospital.2", "BMI", "Insure.type", "MMGSD", "MMGO", "MRISD", "MMGD.2", "ER", "PR", "HER2", "Subtype", "TripleNegative")
```


## Variable definitions
* `Age.Dx`: Age at Diagnosis (in years)
* `Race.Asian`: Non-Asian (reference group); East Asian; South Asian; Southeast Asian; Asian Other
* `RaceEthnicity.4`: Non-Hispanic White (reference group); East Asian; South Asian; Asian Other (combine Southeast Asian and Asian Other in `Race.Asian`) 
* `Hospital.2`: WCM (reference group); NYPQ
* `BMI`: Body Mass Index ($kg/m^2$)
* `Insure.type`: Private (reference group); Medicaid; Medicare
* `MMGSD`: Mammo Screen-Detected Status: Yes (reference group); No
* `MMGO`: Mammo-Occult Status: No (reference group); Yes
* `MRISD`: MRI Screen-Detected Status: No (reference group); Yes
* `MMGD.2`: Mammographic Density: Heterogenously dense (reference group); Almost entirely fatty; Scattered areas of fibrogladular density; Extremely dense 
* `ER`: Estrogen Receptor Status: ER Positive (reference group); ER Negative
* `PR`: Progesterone Receptor Status: PR Positive (reference group); PR Negative
* `HER2`: Human Epidermal Growth Factor Receptor 2 Status: HER2 Negative (reference group); HER2 Positive
* `TripleNegative`: Triple Negative Status: No (reference group); Yes
* `Subtype`: Tumor Subtype: HR+, HER2- (reference group); HER2+; Triple negative  


# Multivariable analysis on breast cancer stage at presentation

```{r include=FALSE}
covariates <- predictors[-which(predictors %in% c("RaceEthnicity.4", "BMI", "PR", "HER2", "Subtype", "TripleNegative"))] 
data.work <- na.omit(data.full[c("StudyID", response, covariates)])
size.1.1 <- nrow(data.work)
```
## Analysis set 1.1 - `r "%+%"(covariates)`
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, Race.Asian, Hospital.2, Insure.type, MMGSD, MMGO, MRISD, MMGD.2, ER
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


## Backward stepwise selection for 1.1 
```{r include=FALSE, warning=FALSE}
glm.fit.back <- glm.fit |> stepAIC(direction = "backward", trace = FALSE)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit.back)), 
                           lower.limit = exp(confint(glm.fit.back))[, 1], 
                           upper.limit = exp(confint(glm.fit.back))[, 2], 
                           p.value = summary(glm.fit.back)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column

## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit.back)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                            align = "lcccc")
```


```{r include=FALSE}
covariates <- predictors[-which(predictors %in% c("Race.Asian", "BMI", "PR", "HER2", "Subtype", "TripleNegative"))] 
data.work <- na.omit(data.full[c("StudyID", response, covariates)])
size.1.2 <- nrow(data.work)
```
## Analysis set 1.2 - `r "%+%"(covariates)`
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, RaceEthnicity.4, Hospital.2, Insure.type, MMGSD, MMGO, MRISD, MMGD.2, ER
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


## Backward stepwise selection for 1.2
```{r include=FALSE, warning=FALSE}
glm.fit.back <- glm.fit |> stepAIC(direction = "backward", trace = FALSE)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit.back)), 
                           lower.limit = exp(confint(glm.fit.back))[, 1], 
                           upper.limit = exp(confint(glm.fit.back))[, 2], 
                           p.value = summary(glm.fit.back)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL

## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit.back)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


```{r include=FALSE}
data.temp <- filter(data.full, RaceEthnicity.1 == "Chinese")
covariates0 <- predictors[-which(predictors %in% c("Race.Asian", "RaceEthnicity.4", "BMI", "PR", "HER2", "Subtype", "TripleNegative"))] 
data.work0 <- na.omit(data.temp[c("StudyID", response, covariates0)])

data.work0$MMGO |> table()
data.work0$MRISD |> table()
data.work0$MMGD.2 |> table()

covariates <- covariates0[-which(covariates0 %in% c("MMGO", "MRISD"))] 
data.work <- na.omit(data.temp[c("StudyID", response, covariates)]) |> filter(MMGD.2 != "Almost entirely fatty")
size.1.3 <- nrow(data.work)
```
## Analysis set 1.3 - `r "%+%"(covariates)` in Chinese patients
Notice that the binary variables `Mammo-Occult` and `MRI Screen-Detected`, and the categorical level `Almost entirely fatty` are removed, because only 1 Chinese patient has Mammo-Occcult Yes, 6 Chinese patients have MRI Screen-Detected Yes, and 3 Chinese patient have almost entirely fatty in Chinese patients with complete data in `r "%+%"(covariates0)`. These small samples generate unstable estimates.
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, Hospital.2, Insure.type, MMGSD, MRISD, MMGD.2, ER
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


```{r include=FALSE}
data.temp <- filter(data.full, Race.Asian == "South Asian")
covariates0 <- predictors[-which(predictors %in% c("Race.Asian", "RaceEthnicity.4", "BMI", "PR", "HER2", "Subtype", "TripleNegative"))]
data.work0 <- na.omit(data.temp[c("StudyID", response, covariates0)])

data.work0$MMGO |> table()
data.work0$MRISD |> table()
data.work0$MMGD.2 |> table()

covariates <- covariates0[-which(covariates0 %in% c("MMGO", "MRISD"))] 
data.work <- na.omit(data.temp[c("StudyID", response, covariates)])
size.1.4 <- nrow(data.work)
```
## Analysis set 1.4 - `r "%+%"(covariates)` in South Asian patients
Notice that the binary variables `Mammo-Occult` and `MRI Screen-Detected` are removed, because only 1 South Asian patient has Mammo-Occcult Yes and 1 South Asian patient has MRI Screen-Detected Yes in South Asian patients with complete data in `r "%+%"(covariates0)`. These small samples generate unstable estimates.
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, Hospital.2, Insure.type, MMGSD, MRISD, MMGD.2, ER
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


```{r include=FALSE}
data.temp <- filter(data.full, RaceEthnicity.4 == "Non-Hispanic Black")
covariates0 <- predictors[-which(predictors %in% c("Race.Asian", "RaceEthnicity.4", "BMI", "PR", "HER2", "Subtype", "TripleNegative"))]
data.work0 <- na.omit(data.temp[c("StudyID", response, covariates0)])

data.work0$MMGO |> table()
data.work0$MRISD |> table()
data.work0$MMGD.2 |> table()

covariates <- covariates0[-which(covariates0 %in% c("MMGO", "MRISD"))] 
data.work <- na.omit(data.temp[c("StudyID", response, covariates)])
size.1.5 <- nrow(data.work)
```
## Analysis set 1.5 - `r "%+%"(covariates)` in Non-Hispanic Black patients
Notice that the binary variables `Mammo-Occult` and `MRI Screen-Detected` are removed, because only 2 Non-Hispanic Black patient has Mammo-Occcult Yes and 3 Non-Hispanic Black patients have MRI Screen-Detected Yes, among Non-Hispanic Black patients with complete data in `r "%+%"(covariates0)`. These small samples generate unstable estimates.
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, Hospital.2, Insure.type, MMGSD, MRISD, MMGD.2, ER
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


```{r include=FALSE}
covariates <- predictors[-which(predictors %in% c("RaceEthnicity.4", "PR", "HER2", "Subtype", "TripleNegative"))] 
data.work <- na.omit(data.full[c("StudyID", response, covariates)])
size.2.1 <- nrow(data.work)
```
## Analysis set 2.1 - `r "%+%"(covariates)`
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, Race.Asian, Hospital.2, Insure.type, MMGSD, MMGO, MRISD, MMGD.2, ER, BMI
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


## Backward stepwise selection for 2.1
```{r include=FALSE, warning=FALSE}
glm.fit.back <- glm.fit |> stepAIC(direction = "backward", trace = FALSE)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit.back)), 
                           lower.limit = exp(confint(glm.fit.back))[, 1], 
                           upper.limit = exp(confint(glm.fit.back))[, 2], 
                           p.value = summary(glm.fit.back)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL

## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit.back)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


```{r include=FALSE}
covariates <- predictors[-which(predictors %in% c("RaceEthnicity.4", "PR", "HER2", "Subtype", "TripleNegative"))] 
data.work <- na.omit(data.full[c("StudyID", response, covariates)])

table(data.work$Race.Asian, data.work$Insure.type)

data.work <- filter(data.work, Race.Asian != "Asian Other")

size.2.2 <- nrow(data.work)
```
## Analysis set 2.2 - `r "%+%"(c(covariates, "Race.Asian*Insure.type"))`
Notice that the categorical level `Asian Other` is removed, because only 2 Asian Other patients have Medicaid, and 6 Asian Other patients have Medicare, in all patients with complete data in `r "%+%"(covariates)`. These small samples generate unstable estimates.
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, Race.Asian, Hospital.2, Insure.type, MMGSD, MMGO, MRISD, MMGD.2, ER, BMI, Race.Asian*Insure.type
glm.fit <- glm("%0%"(response%0%"~", "%+%"(c(covariates, "Race.Asian*Insure.type"))), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


## Backward stepwise selection for 2.2
```{r include=FALSE, warning=FALSE}
glm.fit.back <- glm.fit |> stepAIC(direction = "backward", trace = FALSE)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit.back)), 
                           lower.limit = exp(confint(glm.fit.back))[, 1], 
                           upper.limit = exp(confint(glm.fit.back))[, 2], 
                           p.value = summary(glm.fit.back)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL

## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit.back)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```



```{r include=FALSE}
covariates <- predictors[-which(predictors %in% c("Race.Asian", "PR", "HER2", "Subtype", "TripleNegative"))] 
data.work <- na.omit(data.full[c("StudyID", response, covariates)])
size.2.3 <- nrow(data.work)
```
## Analysis set 2.3 - `r "%+%"(covariates)`
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, RaceEthnicity.4, Hospital.2, Insure.type, MMGSD, MMGO, MRISD, MMGD.2, ER, BMI
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```

## Backward stepwise selection for 2.3
```{r include=FALSE, warning=FALSE}
glm.fit.back <- glm.fit |> stepAIC(direction = "backward", trace = FALSE)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit.back)), 
                           lower.limit = exp(confint(glm.fit.back))[, 1], 
                           upper.limit = exp(confint(glm.fit.back))[, 2], 
                           p.value = summary(glm.fit.back)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL

## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit.back)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


```{r include=FALSE}
covariates <- predictors[-which(predictors %in% c("Race.Asian", "PR", "HER2", "Subtype", "TripleNegative"))] 
data.work <- na.omit(data.full[c("StudyID", response, covariates)])

table(data.work$RaceEthnicity.4, data.work$Insure.type)

size.2.4 <- nrow(data.work)
```
## Analysis set 2.4 - `r "%+%"(c(covariates, "RaceEthnicity.4*Insure.type"))`
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, RaceEthnicity.4, Hospital.2, Insure.type, MMGSD, MMGO, MRISD, MMGD.2, ER, BMI, RaceEthnicity.4*Insure.type
glm.fit <- glm("%0%"(response%0%"~", "%+%"(c(covariates, "RaceEthnicity.4*Insure.type"))), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


## Backward stepwise selection for 2.4
```{r include=FALSE, warning=FALSE}
glm.fit.back <- glm.fit |> stepAIC(direction = "backward", trace = FALSE)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit.back)), 
                           lower.limit = exp(confint(glm.fit.back))[, 1], 
                           upper.limit = exp(confint(glm.fit.back))[, 2], 
                           p.value = summary(glm.fit.back)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL

## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit.back)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


```{r include=FALSE}
covariates <- predictors[-which(predictors %in% c("Race.Asian", "PR", "HER2", "Subtype", "TripleNegative"))] 
data.work <- na.omit(data.full[c("StudyID", response, covariates)])

table(data.work$RaceEthnicity.4, data.work$Insure.type)

table(data.work$RaceEthnicity.4, data.work$MMGD.2)

data.work <- filter(data.work, MMGD.2 != "Almost entirely fatty")

size.2.5 <- nrow(data.work)
```
## Analysis set 2.5 - `r "%+%"(c(covariates, "RaceEthnicity.4*Insure.type", "RaceEthnicity.4*MMGD.2"))`
Notice that the categorical level `Almost entirely fatty` is removed, because only 4 East Asian patients, 3 South Asian patients, 0 Asian Other patient, 4 Hispanic patients, and 2 Other patients have Almost entirely fatty, in all patients with complete data in `r "%+%"(covariates)`. These small samples generate unstable estimates.
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, RaceEthnicity.4, Hospital.2, Insure.type, MMGSD, MMGO, MRISD, MMGD.2, ER, BMI, RaceEthnicity.4*Insure.type, RaceEthnicity.4*MMGD.2
glm.fit <- glm("%0%"(response%0%"~", "%+%"(c(covariates, "RaceEthnicity.4*Insure.type", "RaceEthnicity.4*MMGD.2"))), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


## Backward stepwise selection for 2.5
```{r include=FALSE, warning=FALSE}
glm.fit.back <- glm.fit |> stepAIC(direction = "backward", trace = FALSE)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit.back)), 
                           lower.limit = exp(confint(glm.fit.back))[, 1], 
                           upper.limit = exp(confint(glm.fit.back))[, 2], 
                           p.value = summary(glm.fit.back)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL

## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit.back)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```



```{r include=FALSE}
data.temp <- filter(data.full, RaceEthnicity.1 == "Chinese")
covariates0 <- predictors[-which(predictors %in% c("Race.Asian", "RaceEthnicity.4", "PR", "HER2", "Subtype", "TripleNegative"))] 
data.work0 <- na.omit(data.temp[c("StudyID", response, covariates0)])

data.work0$MMGO |> table()
data.work0$MRISD |> table()
data.work0$MMGD.2 |> table()

covariates <- covariates0[-which(covariates0 %in% c("MMGO", "MRISD"))] 
data.work <- na.omit(data.temp[c("StudyID", response, covariates)]) |> filter(MMGD.2 != "Almost entirely fatty")
size.2.6 <- nrow(data.work)
```
## Analysis set 2.6 - `r "%+%"(covariates)` in Chinese patients
Notice that the binary variables `Mammo-Occult` and `MRI Screen-Detected`, and the categorical level `Almost entirely fatty` are removed, because only 1 Chinese patient has Mammo-Occcult Yes, 6 Chinese patients have MRI Screen-Detected Yes, and 3 Chinese patients have almost entirely fatty, in Chinese patients with complete data in `r "%+%"(covariates0)`. These small samples generate unstable estimates.
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, Hospital.2, Insure.type, MMGSD, MRISD, MMGD.2, ER
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


```{r include=FALSE}
data.temp <- filter(data.full, Race.Asian == "South Asian")
covariates0 <- predictors[-which(predictors %in% c("Race.Asian", "RaceEthnicity.4", "PR", "HER2", "Subtype", "TripleNegative"))] 
data.work0 <- na.omit(data.temp[c("StudyID", response, covariates0)])

data.work0$MMGO |> table()
data.work0$MRISD |> table()
data.work0$MMGD.2 |> table()

covariates <- covariates0[-which(covariates0 %in% c("MMGO", "MRISD"))] 
data.work <- na.omit(data.temp[c("StudyID", response, covariates)])
size.2.7 <- nrow(data.work)
```
## Analysis set 2.7 - `r "%+%"(covariates)` in South Asian patients
Notice that the binary variables `Mammo-Occult` and `MRI Screen-Detected` are removed, because only 1 South Asian patient has Mammo-Occcult Yes, and 1 South Asian patient has MRI Screen-Detected Yes, among South Asian patients with complete data in `r "%+%"(covariates0)`. These small samples generate unstable estimates.
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, Hospital.2, Insure.type, MMGSD, MRISD, MMGD.2, ER
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


```{r include=FALSE}
data.temp <- filter(data.full, RaceEthnicity.4 == "Non-Hispanic Black")
covariates0 <- predictors[-which(predictors %in% c("Race.Asian", "RaceEthnicity.4", "PR", "HER2", "Subtype", "TripleNegative"))] 
data.work0 <- na.omit(data.temp[c("StudyID", response, covariates0)])

data.work0$MMGO |> table()
data.work0$MRISD |> table()
data.work0$MMGD.2 |> table()

covariates <- covariates0[-which(covariates0 %in% c("MMGO", "MRISD"))] 
data.work <- na.omit(data.temp[c("StudyID", response, covariates)])
size.2.8 <- nrow(data.work)
```
## Analysis set 2.8 - `r "%+%"(covariates)` in Non-Hispanic Black patients
Notice that the binary variables `Mammo-Occult` and `MRI Screen-Detected` are removed, because only 1 Non-Hispanic Black patient has Mammo-Occcult Yes, and 2 Non-Hispanic Black patients have MRI Screen-Detected Yes, in Non-Hispanic Black patients with complete data in `r "%+%"(covariates0)`. These small samples generate unstable estimates.
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, Hospital.2, Insure.type, MMGSD, MRISD, MMGD.2, ER
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


```{r include=FALSE}
covariates <- predictors[-which(predictors %in% c("RaceEthnicity.4", "Subtype"))] 
data.work <- na.omit(data.full[c("StudyID", response, covariates)])
size.3.1 <- nrow(data.work)
```
## Analysis set 3.1 - `r "%+%"(covariates)`
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, Race.Asian, Hospital.2, Insure.type, MMGSD, MMGO, MRISD, MMGD.2, ER, BMI, PR, HER2, TripleNegative
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


## Backward stepwise selection for 3.1
```{r include=FALSE, warning=FALSE}
glm.fit.back <- glm.fit |> stepAIC(direction = "backward", trace = FALSE)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit.back)), 
                           lower.limit = exp(confint(glm.fit.back))[, 1], 
                           upper.limit = exp(confint(glm.fit.back))[, 2], 
                           p.value = summary(glm.fit.back)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL

## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit.back)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


```{r include=FALSE}
covariates <- predictors[-which(predictors %in% c("RaceEthnicity.4", "Subtype"))] 
data.work <- na.omit(data.full[c("StudyID", response, covariates)])

table(data.work$Race.Asian, data.work$Insure.type)

data.work <- filter(data.work, Race.Asian != "Asian Other")

size.3.2 <- nrow(data.work)
```
## Analysis set 3.2 - `r "%+%"(c(covariates, "Race.Asian*Insure.type"))`
Notice that the categorical level `Asian Other` is removed, because only 1 Asian Other patient has Medicaid, and 3 Asian Other patients have Medicare, in all patients with complete data in `r "%+%"(covariates)`. These small samples generate unstable estimates.
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, Race.Asian, Hospital.2, Insure.type, MMGSD, MMGO, MRISD, MMGD.2, ER, BMI, Race.Asian*Insure.type
glm.fit <- glm("%0%"(response%0%"~", "%+%"(c(covariates, "Race.Asian*Insure.type"))), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


## Backward stepwise selection for 3.2
```{r include=FALSE, warning=FALSE}
glm.fit.back <- glm.fit |> stepAIC(direction = "backward", trace = FALSE)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit.back)), 
                           lower.limit = exp(confint(glm.fit.back))[, 1], 
                           upper.limit = exp(confint(glm.fit.back))[, 2], 
                           p.value = summary(glm.fit.back)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL

## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit.back)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


```{r include=FALSE}
covariates <- predictors[-which(predictors %in% c("Race.Asian", "Subtype"))] 
data.work <- na.omit(data.full[c("StudyID", response, covariates)])
size.3.3 <- nrow(data.work)
```
## Analysis set 3.3 - `r "%+%"(covariates)`
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, RaceEthnicity.4, Hospital.2, Insure.type, MMGSD, MMGO, MRISD, MMGD.2, ER, BMI, PR, HER2, TripleNegative
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


## Backward stepwise selection for 3.3
```{r include=FALSE, warning=FALSE}
glm.fit.back <- glm.fit |> stepAIC(direction = "backward", trace = FALSE)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit.back)), 
                           lower.limit = exp(confint(glm.fit.back))[, 1], 
                           upper.limit = exp(confint(glm.fit.back))[, 2], 
                           p.value = summary(glm.fit.back)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL

## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit.back)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


```{r include=FALSE}
covariates <- predictors[-which(predictors %in% c("Race.Asian", "Subtype"))] 
data.work <- na.omit(data.full[c("StudyID", response, covariates)])

table(data.work$RaceEthnicity.4, data.work$Insure.type)

size.3.4 <- nrow(data.work)
```
## Analysis set 3.4 - `r "%+%"(c(covariates, "RaceEthnicity.4*Insure.type"))`
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, RaceEthnicity.4, Hospital.2, Insure.type, MMGSD, MMGO, MRISD, MMGD.2, ER, BMI, RaceEthnicity.4*Insure.type
glm.fit <- glm("%0%"(response%0%"~", "%+%"(c(covariates, "RaceEthnicity.4*Insure.type"))), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


## Backward stepwise selection for 3.4
```{r include=FALSE, warning=FALSE}
glm.fit.back <- glm.fit |> stepAIC(direction = "backward", trace = FALSE)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit.back)), 
                           lower.limit = exp(confint(glm.fit.back))[, 1], 
                           upper.limit = exp(confint(glm.fit.back))[, 2], 
                           p.value = summary(glm.fit.back)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL

## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit.back)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


```{r include=FALSE}
data.temp <- filter(data.full, RaceEthnicity.1 == "Chinese")
covariates0 <- predictors[-which(predictors %in% c("Race.Asian", "RaceEthnicity.4", "Subtype"))] 
data.work0 <- na.omit(data.temp[c("StudyID", response, covariates0)])

data.work0$MMGO |> table()
data.work0$MRISD |> table()
data.work0$MMGD.2 |> table()

# all.fac <- sapply(names(data.work0)[sapply(data.work0, is.factor)], function(x) with(data.work0, table(get(x), useNA = "ifany")))
# View(all.fac)

covariates <- covariates0[-which(covariates0 %in% c("MMGO", "MRISD"))] 
data.work <- na.omit(data.temp[c("StudyID", response, covariates)]) |> filter(MMGD.2 != "Almost entirely fatty")
size.3.5 <- nrow(data.work)
```
## Analysis set 3.5 - `r "%+%"(covariates)` in Chinese patient
Notice that the binary variables `Mammo-Occult` and `MRI Screen-Detected`, and the categorical level `Almost entirely fatty` are removed, because only 1 Chinese patient has Mammo-Occcult Yes, 2 Chinese patients have MRI Screen-Detected Yes, and 3 Chinese patients have almost entirely fatty in Chinese patients with complete data in `r "%+%"(covariates0)`. These small samples generate unstable estimates.
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, RaceEthnicity.4, Hospital.2, Insure.type, MMGSD, MMGO, MRISD, MMGD.2, ER, BMI, PR, HER2, TripleNegative
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


```{r include=FALSE}
data.temp <- filter(data.full, Race.Asian == "South Asian")
covariates0 <- predictors[-which(predictors %in% c("Race.Asian", "RaceEthnicity.4", "Subtype"))] 
data.work0 <- na.omit(data.temp[c("StudyID", response, covariates0)])

data.work0$MMGO |> table()
data.work0$MRISD |> table()
data.work0$MMGD.2 |> table()

covariates <- covariates0[-which(covariates0 %in% c("MMGO", "MRISD"))] 
data.work <- na.omit(data.temp[c("StudyID", response, covariates)])
size.3.6 <- nrow(data.work)
```
## Analysis set 3.6 - `r "%+%"(covariates)` in South Asian patients
Notice that the binary variables `Mammo-Occult` and `MRI Screen-Detected` are removed, because only 1 South Asian patient has Mammo-Occcult Yes, and 1 South Asian patient has MRI Screen-Detected Yes, in South Asian patients with complete data in `r "%+%"(covariates0)`. These small samples generate unstable estimates.
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, RaceEthnicity.4, Hospital.2, Insure.type, MMGSD, MMGO, MRISD, MMGD.2, ER, BMI, PR, HER2, TripleNegative
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


```{r include=FALSE}
data.temp <- filter(data.full, RaceEthnicity.4 == "Non-Hispanic Black")
covariates0 <- predictors[-which(predictors %in% c("Race.Asian", "RaceEthnicity.4", "Subtype"))] 
data.work0 <- na.omit(data.temp[c("StudyID", response, covariates0)])

data.work0$MMGO |> table()
data.work0$MRISD |> table()
data.work0$MMGD.2 |> table()

# all.fac <- sapply(names(data.work0)[sapply(data.work0, is.factor)], function(x) with(data.work0, table(get(x), useNA = "ifany")))
# View(all.fac)

covariates <- covariates0[-which(covariates0 %in% c("MMGO", "MRISD"))] 
data.work <- na.omit(data.temp[c("StudyID", response, covariates)])
size.3.7 <- nrow(data.work)
```
## Analysis set 3.7 - `r "%+%"(covariates)` in Non-Hispanic Black patients
Notice that the binary variables `Mammo-Occult` and `MRI Screen-Detected` are removed, because only 1 Non-Hispanic Black patient has Mammo-Occcult Yes, and 1 Non-Hispanic Black patient has MRI Screen-Detected Yes, among Non-Hispanic Black patients with complete data in `r "%+%"(covariates0)`. These small samples generate unstable estimates.
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, RaceEthnicity.4, Hospital.2, Insure.type, MMGSD, MMGO, MRISD, MMGD.2, ER, BMI, PR, HER2, TripleNegative
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


```{r include=FALSE}
covariates <- predictors[-which(predictors %in% c("RaceEthnicity.4", "PR", "HER2", "ER", "TripleNegative"))] 
data.work <- na.omit(data.full[c("StudyID", response, covariates)])
size.3.8 <- nrow(data.work)
```
## Analysis set 3.8 - `r "%+%"(covariates)`
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, Race.Asian, Hospital.2, Insure.type, MMGSD, MMGO, MRISD, MMGD.2, Subtype, BMI
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


## Backward stepwise selection for 3.8
```{r include=FALSE, warning=FALSE}
glm.fit.back <- glm.fit |> stepAIC(direction = "backward", trace = FALSE)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit.back)), 
                           lower.limit = exp(confint(glm.fit.back))[, 1], 
                           upper.limit = exp(confint(glm.fit.back))[, 2], 
                           p.value = summary(glm.fit.back)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL

## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit.back)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```



```{r include=FALSE}
covariates <- predictors[-which(predictors %in% c("RaceEthnicity.4", "PR", "HER2", "ER", "TripleNegative"))] 
data.work <- na.omit(data.full[c("StudyID", response, covariates)])

table(data.work$Race.Asian, data.work$Insure.type)

data.work <- filter(data.work, Race.Asian != "Asian Other")

size.3.9 <- nrow(data.work)
```
## Analysis set 3.9 - `r "%+%"(c(covariates, "Race.Asian*Insure.type"))`
Notice that the categorical level `Asian Other` is removed, because only 1 Asian Other patient has Medicaid, and 3 Asian Other patients have Medicare, in all patients with complete data in `r "%+%"(covariates)`. These small samples generate unstable estimates.
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, Race.Asian, Hospital.2, Insure.type, MMGSD, MMGO, MRISD, MMGD.2, Subtype, BMI, Race.Asian*Insure.type
glm.fit <- glm("%0%"(response%0%"~", "%+%"(c(covariates, "Race.Asian*Insure.type"))), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


## Backward stepwise selection for 3.9
```{r include=FALSE, warning=FALSE}
glm.fit.back <- glm.fit |> stepAIC(direction = "backward", trace = FALSE)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit.back)), 
                           lower.limit = exp(confint(glm.fit.back))[, 1], 
                           upper.limit = exp(confint(glm.fit.back))[, 2], 
                           p.value = summary(glm.fit.back)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL

## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit.back)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```



```{r include=FALSE}
covariates <- predictors[-which(predictors %in% c("Race.Asian", "PR", "HER2", "ER", "TripleNegative"))] 
data.work <- na.omit(data.full[c("StudyID", response, covariates)])
size.3.10 <- nrow(data.work)
```
## Analysis set 3.10 - `r "%+%"(covariates)`
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, RaceEthnicity.4, Hospital.2, Insure.type, MMGSD, MMGO, MRISD, MMGD.2, Subtype, BMI
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


## Backward stepwise selection for 3.10
```{r include=FALSE, warning=FALSE}
glm.fit.back <- glm.fit |> stepAIC(direction = "backward", trace = FALSE)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit.back)), 
                           lower.limit = exp(confint(glm.fit.back))[, 1], 
                           upper.limit = exp(confint(glm.fit.back))[, 2], 
                           p.value = summary(glm.fit.back)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL

## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit.back)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```



```{r include=FALSE}
covariates <- predictors[-which(predictors %in% c("Race.Asian", "PR", "HER2", "ER", "TripleNegative"))] 
data.work <- na.omit(data.full[c("StudyID", response, covariates)])

table(data.work$RaceEthnicity.4, data.work$Insure.type)

size.3.11 <- nrow(data.work)
```
## Analysis set 3.11 - `r "%+%"(c(covariates, "RaceEthnicity.4*Insure.type"))`
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, RaceEthnicity.4, Hospital.2, Insure.type, MMGSD, MMGO, MRISD, MMGD.2, ER, BMI, RaceEthnicity.4*Insure.type
glm.fit <- glm("%0%"(response%0%"~", "%+%"(c(covariates, "RaceEthnicity.4*Insure.type"))), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```



## Backward stepwise selection for 3.11
```{r include=FALSE, warning=FALSE}
glm.fit.back <- glm.fit |> stepAIC(direction = "backward", trace = FALSE)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit.back)), 
                           lower.limit = exp(confint(glm.fit.back))[, 1], 
                           upper.limit = exp(confint(glm.fit.back))[, 2], 
                           p.value = summary(glm.fit.back)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL

## Add a new column that formats the default variable name outputted by glm
result.multi$object <- alter.glm.rownames(glm.fit.back)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```



```{r include=FALSE}
data.temp <- filter(data.full, RaceEthnicity.1 == "Chinese")
covariates0 <- predictors[-which(predictors %in% c("Race.Asian", "RaceEthnicity.4", "PR", "HER2", "ER", "TripleNegative"))] 
data.work0 <- na.omit(data.temp[c("StudyID", response, covariates0)])

data.work0$MMGO |> table()
data.work0$MRISD |> table()
data.work0$MMGD.2 |> table()

covariates <- covariates0[-which(covariates0 %in% c("MMGO", "MRISD"))] 
data.work <- na.omit(data.temp[c("StudyID", response, covariates)]) |> filter(MMGD.2 != "Almost entirely fatty")
size.3.12 <- nrow(data.work)
```
## Analysis set 3.12 - `r "%+%"(covariates)` in Chinese patients
Notice that the binary variables `Mammo-Occult` and `MRI Screen-Detected`, and the categorical level `Almost entirely fatty` are removed, because only 1 Chinese patient has Mammo-Occcult Yes, 2 Chinese patients have MRI Screen-Detected Yes, and 3 Chinese patients have almost entirely fatty, in Chinese patients with complete data in `r "%+%"(covariates0)`. These small samples generate unstable estimates.
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, Hospital.2, Insure.type, MMGSD, MRISD, MMGD.2, BMI, Subtype
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


```{r include=FALSE}
data.temp <- filter(data.full, Race.Asian == "South Asian")
covariates0 <- predictors[-which(predictors %in% c("Race.Asian", "RaceEthnicity.4", "PR", "HER2", "ER", "TripleNegative"))] 
data.work0 <- na.omit(data.temp[c("StudyID", response, covariates0)])

data.work0$MMGO |> table()
data.work0$MRISD |> table()
data.work0$MMGD.2 |> table()

covariates <- covariates0[-which(covariates0 %in% c("MMGO", "MRISD"))] 
data.work <- na.omit(data.temp[c("StudyID", response, covariates)])
size.3.13 <- nrow(data.work)
```
## Analysis set 3.13 - `r "%+%"(covariates)` in South Asian patients
Notice that the binary variables `Mammo-Occult` and `MRI Screen-Detected` are removed, because only 1 South Asian patient has Mammo-Occcult Yes, and 1 South Asian patient has MRI Screen-Detected Yes, among South Asian patients with complete data in `r "%+%"(covariates0)`. These small samples generate unstable estimates.
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, Hospital.2, Insure.type, MMGSD, MRISD, MMGD.2, BMI, Subtype
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


```{r include=FALSE}
data.temp <- filter(data.full, RaceEthnicity.4 == "Non-Hispanic Black")
covariates0 <- predictors[-which(predictors %in% c("Race.Asian", "RaceEthnicity.4", "PR", "HER2", "ER", "TripleNegative"))] 
data.work0 <- na.omit(data.temp[c("StudyID", response, covariates0)])

data.work0$MMGO |> table()
data.work0$MRISD |> table()
data.work0$MMGD.2 |> table()

covariates <- covariates0[-which(covariates0 %in% c("MMGO", "MRISD"))] 
data.work <- na.omit(data.temp[c("StudyID", response, covariates)])
size.3.14 <- nrow(data.work)
```
## Analysis set 3.14 - `r "%+%"(covariates)` in Non-Hispanic Black patients
Notice that the binary variables `Mammo-Occult` and `MRI Screen-Detected` are removed, because only 1 Non-Hispanic Black patient has Mammo-Occcult Yes, and 1 Non-Hispanic Black patient has MRI Screen-Detected Yes, in Non-Hispanic Black patients with complete data in `r "%+%"(covariates0)`. These small samples generate unstable estimates.
```{r include=FALSE, warning=FALSE}
# Run multivariate logistic models using Early.stage, Age.Dx, Hospital.2, Insure.type, MMGSD, MRISD, MMGD.2, BMI, Subtype
glm.fit <- glm("%0%"(response%0%"~", "%+%"(covariates)), data = data.work, family = binomial)

# Output odds ratio and its 95% CI
result.multi <- data.frame(OR = exp(coef(glm.fit)), 
                           lower.limit = exp(confint(glm.fit))[, 1], 
                           upper.limit = exp(confint(glm.fit))[, 2], 
                           p.value = summary(glm.fit)$coef[, 4]) 

# Format the output
## Change the row names
rownames(result.multi) <- NULL
## Add a new column
result.multi$object <- alter.glm.rownames(glm.fit)

result.fit <- result.multi %>%
  dplyr::select(object, everything()) %>%
  mutate(across(.cols = c(OR, lower.limit, upper.limit), .fns = round, 2)) %>% 
  mutate(p.value = round(p.value, 3)) %>% 
  mutate(p.value = ifelse(p.value < 0.001, "< 0.001", p.value)) %>%
  filter(object != "Intercept")
```


```{r echo=FALSE, results='asis'}
result.fit |> knitr::kable(col.names = c("Covariate", "Adjusted Odds Ratio", "95% CI Lower Limit", "95% CI Upper Limit", "P-value"), # Add column names
                           align = "lcccc")
```


# Sample size for each analysis set
```{r echo=FALSE, results='asis'}
sample.size <- mget(c("size.1."%0%1:5, "size.2."%0%1:8, "size.3."%0%1:14)) |> unlist()
names(sample.size) <- NULL
data.frame(
  `Analysis set` = c("1."%0%1:5, "2."%0%1:8, "3."%0%1:14),
  `Sample Size` = sample.size
  ) |> knitr::kable()
```


# Session information
```{r echo=FALSE}
sessionInfo()
```