---
title: "Predictive Modeling of Blood Pressure Categories: Integrating Demographic and Dietary Factors for Personalized Management"
authors:
  - name: Li Yuan
    department: Department of Statistics
    affiliation: University of Michigan, Ann Arbor
    location: Ann Arbor, MI 48109
    email: leeyuan@umich.edu
abstract: |
  This study delves into predictive modeling of blood pressure levels, focusing on the United States, addressing the global health concern of hypertension. Mainly utilizing demographic and dietary data from the Centers for Disease Control and Prevention (CDC) National Health and Nutrition Examination Survey (NHANES) 2017-2018, aims to craft personalized management strategies. Drawing on research emphasizing the multifaceted determinants of hypertension, we leverage the multinomial regression model with lasso regularization as a baseline. Furthermore, the study advances to the eXtreme Gradient Boosting (XGBoost) algorithm, achieving a better performance than multinomial regression. Evaluation metrics include accuracy and Area Under the ROC Curve (ROC-AUC) in a 10-fold cross validation framework. The study provides possible personal blood pressure management solution.
keywords:
  - blood pressure
  - multinomial regression
  - lasso
  - tree
  - gradient boosting
  - xgboost
  - R
bibliography: references.bib
biblio-style: unsrt
output: rticles::arxiv_article
header-includes:
   - \usepackage{amsmath}
   - \usepackage{graphicx}
   - \usepackage{subcaption}
---

```{r setup, include=FALSE}
rm(list = ls())
knitr::opts_chunk$set(echo = TRUE)
set.seed(88)
library(caret)
library(dplyr)
library(GGally)
library(ggpubr)
library(glmnet)
library(grid)
library(gridExtra)
library(haven)
library(knitr)
library(patchwork)
library(pROC)
library(tidyr)
library(xgboost)
knitr::opts_chunk$set(fig.pos = "ht", out.extra = "")
```

# Introduction

Hypertension, a prevalent and frequently asymptomatic health condition, persists as a global health concern, significantly contributing to the burden of cardiovascular diseases [@forouzanfar_2016_global]. Against the backdrop of recent advancements in data science and machine learning, this paper initiates a meticulous investigation into the predictive modeling of blood pressure categories, specifically concentrating on the United States. The study underscores the importance of integrating demographic and dietary information to formulate personalized management strategies tailored to the unique health landscape of the United States.

Contemporary research on hypertension highlights the multifaceted nature of its determinants, necessitating a comprehensive approach to prediction and management. A seminal study by @iqbal_2021_demographic emphasizes the significance of demographic factors in predicting hypertension prevalence, underscoring the need for nuanced models that account for individual characteristics. Additionally, the work of @johnson_2009_dietary advocates for a personalized approach, accentuating the substantial influence of dietary habits on blood pressure regulation.

The dataset utilized in this study originates from the Centers for Disease Control and Prevention (CDC) National Health and Nutrition Examination Survey (NHANES), focusing specifically on the years 2017-2018. Our dataset encompasses a comprehensive array of variables, including blood pressure measurements, demographic profiles, nutrient intakes, diabetes indicators, and health insurance details. By integrating these diverse factors, the study aims to delve into the intricate relationships between various determinants and blood pressure outcomes, ultimately seeking to develop a predictive model tailored to the specific context of blood pressure management in the United States.

As a baseline, we employ multinomial regression with lasso regularization and 10-fold cross-validation. This choice is motivated by the desire for a robust baseline that reflects the complexities of medical datasets, particularly in disease prediction scenarios. In further pursuit, inspired by @islam_2023_predicting in their research on hypertension modeling in Ethiopia, we leverage the eXtreme Gradient Boosting (XGB) algorithm to surpass the multinomial regression with lasso regularization baseline. The adoption of these methodologies aligns with the evolving landscape of predictive modeling in hypertension research. With a focus on predicting blood pressure categories (normal, elevated, and high), we evaluate these models based on their test accuracy and test Area under the ROC Curve (ROC-AUC).

As we embark on this study, we draw from a rich tapestry of existing research to contribute novel insights into the predictive modeling of blood pressure categories. Our aim extends beyond advancing the technical aspects of machine learning applications; we seek to provide practical and personalized strategies for hypertension management, aligning with the evolving landscape of precision medicine in cardiovascular health.

# Data

According to the @centersfordiseasecontrolandprevention_2023_20192020, NHANES field operations were suspended in March 2020 because of COVID-19. Consequently, data collection for the NHANES 2019-2020 cycle was incomplete, rendering it non-nationally representative. In response to the disruption caused by the COVID-19 pandemic, we only use those data collected in the 2017-2018 cycle to ensure the study's relevance and generalizability to the U.S. civilian non-institutionalized population. 

NHANES employs a complex, multistage probability design for sampling the civilian, noninstitutionalized population in the U.S. In 2017-2018, 16,211 persons were selected from 30 survey locations, with 9,254 completing interviews and 8,704 undergoing examinations. Each participant has a unique identification number `SEQN`. To ensure representation, materials were translated into various languages, and cultural competency training was provided to staff [@cdc_2017_2018_nhanes].

## Overview

In the context of this study, the most important data we selected is the examination data of blood pressure (`BPX_J`), which "provides data for three consecutive blood pressure (BP) measurements and other methodological measurements to obtain an accurate BP. Heart rate or pulse, depending on age, are also reported" [@centersfordiseasecontrolandprevention_2020_20172018]. This data contains 4 readings of systolic blood pressure and 4 readings of diastolic blood pressure for each participant. In order to create a response variable about blood pressure level (`BPXLEVEL`), we first average the 4 readings of systolic blood pressure and diastolic blood pressure of each participant respectively. Then we follow the definition of normal, elevated, and hypertension provided by @whelton_2018_2017 to divide our average systolic blood pressure and diastolic blood pressure into three blood pressure levels shown in table \ref{tab:blood_pressure}.

```{r data-clean1, include=FALSE}
findNonNAColumns = function(data) {
  return(colnames(data)[apply(data, 2, function(col) all(!is.na(col)))])
}
# Import 2017 - 2018 blood pressure data
# Doc: https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/BPX_J.htm
BPX_J = read_xpt("https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/BPX_J.XPT")
BPX_J = BPX_J %>%
  mutate(BPXSYAVG = rowMeans(select(., starts_with("BPXSY")), na.rm = TRUE),
         BPXDIAVG = rowMeans(select(., starts_with("BPXDI")), na.rm = TRUE)) %>%
  filter(complete.cases(BPXSYAVG, BPXDIAVG)) %>%
  mutate(BPXLEVEL = case_when(
    # Normal blood pressure
    BPXSYAVG < 120 & BPXDIAVG < 80 ~ 0,
    # Elevated blood pressure
    (BPXSYAVG >= 120 & BPXSYAVG < 130) & (BPXDIAVG < 80) ~ 1,
    # High blood pressure
    (BPXSYAVG >= 130) | (BPXDIAVG >= 80) ~ 2
  ))
columns_no_na = findNonNAColumns(BPX_J)
BPX_J = BPX_J[c("BPXLEVEL", "SEQN", "BPACSZ", "BPXPLS",
                "BPXPTY", "BPXSYAVG", "BPXDIAVG")]
BPX_J = BPX_J %>%
  mutate(
    BPXLEVEL = as.factor(BPXLEVEL),
    BPACSZ = as.factor(BPACSZ),
    BPXPTY = as.factor(BPXPTY)
  )
BPXAVG = BPX_J[,c(6, 7)]
BPX_J = BPX_J[,-c(6, 7)]
```

\begin{table}[ht]
  \caption{Blood Pressure Levels Divided by Systolic and Diastolic Blood Pressure}
  \centering
  \begin{tabular}{llll}
    \toprule
    Blood Pressure Levels & Systolic Blood Pressure &  & Diastolic Blood Pressure \\
    \midrule
    Normal (\texttt{BPXLEVEL} = 0)              & $<$ 120 mmHg            & and & $<$ 80 mmHg          \\
    Elevated (\texttt{BPXLEVEL} = 1)            & 120-129 mmHg            & and & $<$ 80 mmHg          \\
    Hypertension (\texttt{BPXLEVEL} = 2)        & $\geq$ 130 mmHg         & or  & $\geq$ 80 mmHg        \\
    \bottomrule
  \end{tabular}
  \label{tab:blood_pressure}
\end{table}

After obtaining blood pressure levels (`BPXLEVEL`), we integrated additional data from NHANES, specifically Demographic Variables and Sample Weights (`DEMO_J`), which offers individual, family, and household-level information [@centersfordiseasecontrolandprevention_2020_20172018DEMO_J]. Additionally, we incorporated Dietary Interview - Total Nutrient Intakes, First Day (`DR1TOT_J`), containing detailed dietary intake information from NHANES participants. This data is used to estimate the types and amounts of foods and beverages consumed during the 24-hour period before the interview, along with estimating intakes of energy, nutrients, and other food components [@centersfordiseasecontrolandprevention_2020_20172018DR1TOT_J].

Furthermore, we included Diabetes (`DIQ_J`), providing personal interview data on diabetes, prediabetes, use of insulin or oral hypoglycemic medications, and diabetic retinopathy [@centersfordiseasecontrolandprevention_2020_20172018DIQ_J]. Lastly, Health Insurance (`HIQ_J`) data was merged, offering respondent-level interview information on insurance coverage, type of coverage, coverage of prescription drugs, and uninsured status during the past 12 months [@centersfordiseasecontrolandprevention_2020_20172018HIQ_J]. 

By merging these data based on participants' unique identification number `SEQN`, selecting relevant predictors, and removing some missing values, we got a curated data frame with 6,125 observations and 80 variables. Description of these 80 variables are adapted from @cdc_search_variable and shown in table \ref{tab:variable_list1}.

```{r data-clean2, include=FALSE}
# Import 2017 - 2018 demographic data
# Doc: https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/DEMO_J.htm
DEMO_J = read_xpt("https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/DEMO_J.XPT")
currentVar = colnames(BPX_J)
FULLDATA = BPX_J %>% left_join(DEMO_J, by = "SEQN")
columns_no_na1 = setdiff(findNonNAColumns(FULLDATA), currentVar)
FULLDATA = FULLDATA[c(colnames(BPX_J), "RIAGENDR", "RIDAGEYR", "RIDRETH3", 
                      "DMDHHSIZ", "DMDHHSZA", "DMDHHSZB", "DMDHHSZE")]
FULLDATA = FULLDATA %>%
  mutate(
    RIAGENDR = as.factor(RIAGENDR),
    RIDRETH3 = as.factor(RIDRETH3)
  )

# Import 2017 - 2018 Total Nutrient Intakes, First Day
# Doc: https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/DR1TOT_J.htm
DR1TOT_J = read_xpt("https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/DR1TOT_J.XPT")
FULLDATA = FULLDATA %>% left_join(DR1TOT_J, by = "SEQN")
FULLDATA = FULLDATA[, -c(13:42)]
FULLDATA = FULLDATA[, -c(79:149)]
FULLDATA = na.omit(FULLDATA)

# Import 2017 - 2018 Diabetes
# Doc: https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/DIQ_J.htm
currentVar = colnames(FULLDATA)
DIQ_J = read_xpt("https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/DIQ_J.XPT")
FULLDATA = FULLDATA %>% left_join(DIQ_J, by = "SEQN")
FULLDATA = FULLDATA[c(currentVar, "DIQ010", "DIQ050")]
FULLDATA = FULLDATA %>%
  mutate(
    DIQ010 = as.factor(DIQ010),
    DIQ050 = as.factor(DIQ050)
  )

# Import 2017 - 2018 Health Insurance
# Doc: https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/HIQ_J.htm
currentVar = colnames(FULLDATA)
HIQ_J = read_xpt("https://wwwn.cdc.gov/Nchs/Nhanes/2017-2018/HIQ_J.XPT")
FULLDATA = FULLDATA %>% left_join(HIQ_J, by = "SEQN")
FULLDATA = FULLDATA[c(currentVar, "HIQ011")]
FULLDATA = FULLDATA %>%
  mutate(
    HIQ011 = as.factor(HIQ011)
  )
# Prepare data for training and testing
FULLDATA = FULLDATA[,-2]
n.obs = nrow(FULLDATA)
index.train = sample(seq(n.obs), floor(n.obs * 0.8), replace = FALSE)
# Data frame train and test
train = FULLDATA[index.train, ]
test = FULLDATA[-index.train, ]
train.X = train[, -1]
train.Y = train$BPXLEVEL
test.X = test[, -1]
test.Y = test$BPXLEVEL
# Model matrix train and test
train.X.mm = model.matrix(~ . - 1, train.X)
test.X.mm = model.matrix(~ . - 1, test.X)
# Data frame train and test with dummy
train.X.dummy = cbind(train.Y, as.data.frame(train.X.mm))
colnames(train.X.dummy)[1] = "BPXLEVEL"
test.X.dummy = cbind(test.Y, as.data.frame(test.X.mm))
colnames(test.X.dummy)[1] = "BPXLEVEL"
```

\begin{table}[ht]
  \caption{Variable Names and Labels in the Curated Data Frame}
  \centering
  \begin{tabular}{p{1.45cm}p{4cm}p{1.55cm}|p{1.45cm}p{4cm}p{1.55cm}}
    \toprule
    Name & Label & Source & Name & Label & Source \\
    \midrule
    \texttt{BPXLEVEL} & Blood pressure levels & Derived & \texttt{BPACSZ} & Coded cuff size & \texttt{BPX\textunderscore J} \\
    \texttt{BPXPLS} & 60 sec. pulse & \texttt{BPX\textunderscore J} & \texttt{BPXPTY} & Pulse type & \texttt{BPX\textunderscore J} \\
    \texttt{RIAGENDR} & Gender & \texttt{DEMO\textunderscore J} & \texttt{RIDAGEYR} & Age in years at screening & \texttt{DEMO\textunderscore J} \\
    \texttt{RIDRETH3} & Race / Hispanic origin w/ NH Asian & \texttt{DEMO\textunderscore J} & \texttt{DMDHHSIZ} & Total number of people in the Household & \texttt{DEMO\textunderscore J} \\
    \texttt{DMDHHSZA} & Number of children 5 years or younger in HH & \texttt{DEMO\textunderscore J} & \texttt{DMDHHSZB} & Number of children 6-17 years old in HH & \texttt{DEMO\textunderscore J} \\
    \texttt{DMDHHSZE} & Number of adults 60 years or older in HH & \texttt{DEMO\textunderscore J} & \texttt{DR1TNUMF} & Number of foods / beverages reported & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TKCAL} & Energy (kcal) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TPROT} & Protein (gm) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TCARB} & Carbohydrate (gm) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TSUGR} & Total sugars (gm) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TFIBE} & Dietary fiber (gm) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TTFAT} & Total fat (gm) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TSFAT} & Total saturated fatty acids (gm) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TMFAT} & Total monounsaturated fatty acids (gm) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TPFAT} & Total polyunsaturated fatty acids (gm) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TCHOL} & Cholesterol (mg) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TATOC} & Vitamin E as alpha-tocopherol (mg) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TATOA} & Added alpha-tocopherol (Vitamin E) (mg) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TRET} & Retinol (mcg) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TVARA} & Vitamin A, RAE (mcg) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TACAR} & Alpha-carotene (mcg) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TBCAR} & Beta-carotene (mcg) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TCRYP} & Beta-cryptoxanthin (mcg) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TLYCO} & Lycopene (mcg) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TLZ} & Lutein + zeaxanthin (mcg) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TVB1} & Thiamin (Vitamin B1) (mg) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TVB2} & Riboflavin (Vitamin B2) (mg) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TNIAC} & Niacin (mg) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TVB6} & Vitamin B6 (mg) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TFOLA} & Total folate (mcg) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TFA} & Folic acid (mcg) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TFF} & Food folate (mcg) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TFDFE} & Folate, DFE (mcg) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TCHL} & Total choline (mg) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TVB12} & Vitamin B12 (mcg) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TB12A} & Added vitamin B12 (mcg) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TVC} & Vitamin C (mg) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TVD} & Vitamin D (D2 + D3) (mcg) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TVK} & Vitamin K (mcg) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TCALC} & Calcium (mg) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TPHOS} & Phosphorus (mg) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TMAGN} & Magnesium (mg) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TIRON} & Iron (mg) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TZINC} & Zinc (mg) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TCOPP} & Copper (mg) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TSODI} & Sodium (mg) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TPOTA} & Potassium (mg) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TSELE} & Selenium (mcg) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TCAFF} & Caffeine (mg) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TTHEO} & Theobromine (mg) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TALCO} & Alcohol (gm) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TMOIS} & Moisture (gm) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TS040} & SFA 4:0 (Butanoic) (gm) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TS060} & SFA 6:0 (Hexanoic) (gm) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TS080} & SFA 8:0 (Octanoic) (gm) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TS100} & SFA 10:0 (Decanoic) (gm) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TS120} & SFA 12:0 (Dodecanoic) (gm) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TS140} & SFA 14:0 (Tetradecanoic) (gm) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TS160} & SFA 16:0 (Hexadecanoic) (gm) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TS180} & SFA 18:0 (Octadecanoic) (gm) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TM161} & MFA 16:1 (Hexadecenoic) (gm) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TM181} & MFA 18:1 (Octadecenoic) (gm) & \texttt{DR1TOT\textunderscore J} \\
  \end{tabular}
  \label{tab:variable_list1}
\end{table}

\begin{table}[ht]
  \centering
  \begin{tabular}{p{1.45cm}p{4cm}p{1.55cm}|p{1.45cm}p{4cm}p{1.55cm}}
    \toprule
    Name & Label & Source & Name & Label & Source \\
    \midrule
    \texttt{DR1TM201} & MFA 20:1 (Eicosenoic) (gm) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TM221} & MFA 22:1 (Docosenoic) (gm) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TP182} & PFA 18:2 (Octadecadienoic) (gm) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TP183} & PFA 18:3 (Octadecatrienoic) (gm) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TP184} & PFA 18:4 (Octadecatetraenoic) (gm) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TP204} & PFA 20:4 (Eicosatetraenoic) (gm) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TP205} & PFA 20:5 (Eicosapentaenoic) (gm) & \texttt{DR1TOT\textunderscore J} & \texttt{DR1TP225} & PFA 22:5 (Docosapentaenoic) (gm) & \texttt{DR1TOT\textunderscore J} \\
    \texttt{DR1TP226} & PFA 22:6 (Docosahexaenoic) (gm) & \texttt{DR1TOT\textunderscore J} & \texttt{DIQ010} & Doctor told you have diabetes & \texttt{DIQ\textunderscore J} \\
    \texttt{DIQ050} & Taking insulin now & \texttt{DIQ\textunderscore J} & \texttt{HIQ011} & Covered by health insurance & \texttt{HIQ\textunderscore J} \\
    \bottomrule
  \end{tabular}
  \label{tab:variable_list2}
\end{table}

## Visualization

In this section, we present two scatterplot matrices that provide a comprehensive visual exploration of the dataset. The first matrix focuses on demographic information, offering insights into the relationships and distributions among key demographic variables. The second matrix encompasses macronutrient intakes. These visualizations aim to reveal potential patterns, correlations, and trends within the dataset.

```{r visualization1, include=FALSE}
scatter.matrix1 = ggpairs(FULLDATA[c("BPXLEVEL", "RIAGENDR", "RIDAGEYR", "RIDRETH3")], 
        aes(color = BPXLEVEL)) + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size = 8),
        axis.text.y = element_text(size = 7))
pdf(file = "scattermatrix1.pdf", width = 8, height = 8)
print(scatter.matrix1)
dev.off()
```
```{r visualization1-embed, echo = FALSE, out.width = "250px", fig.align = "center", fig.cap = "Scatterplot Matrix of BPXLEVEL Against Some Demographic Information\\label{fig:scattermatrix1}"}
knitr::include_graphics("scattermatrix1.pdf")
```

In Figure \ref{fig:scattermatrix1}, we utilized a color-coded scheme to represent different blood pressure levels: red for normal, green for elevated, and blue for hypertension. By examining the relationship between blood pressure levels (`BPXLEVEL`) and gender (`RIAGENDR`), noteworthy patterns emerge. The plot reveals a higher prevalence of elevated blood pressure and hypertension among male participants (coded as 1) compared to their female counterparts (coded as 2).

Further exploration of blood pressure levels against age (`RIAGEYR`) exposes intriguing insights. The distributions indicate a skewed pattern, with individuals younger than 20 predominantly exhibiting normal blood pressure levels. However, a concerning trend is observed among those around 60 years old, who are more likely to have hypertension. Age emerges as a potential influential factor for predicting blood pressure levels in future models.

Analyzing blood pressure levels against race (`RIDRETH3`) unveils distinct prevalence rates. Non-Hispanic White individuals (coded as 3) demonstrate the highest incidence of hypertension, followed by Non-Hispanic Black (coded as 4), Mexican American (coded as 1), and Non-Hispanic Asian (coded as 6) individuals. Categories 2 and 7, representing other Hispanic and other races (including multi-racial), exhibit the lowest hypertension cases.

```{r visualization2, include=FALSE}
scatter.matrix2 = ggpairs(FULLDATA[c("BPXLEVEL", "DR1TCARB", "DR1TPROT", "DR1TFIBE",
                                     "DR1TTFAT", "DR1TCHOL")], 
        aes(color = BPXLEVEL)) + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size = 5),
        axis.text.y = element_text(size = 5))
pdf(file = "scattermatrix2.pdf", width = 10, height = 10)
print(scatter.matrix2)
dev.off()
```
```{r visualization2-embed, echo = FALSE, out.width = "250px", fig.align = "center", fig.cap = "Scatterplot Matrix of BPXLEVEL Against Marco Nutrient Intakes\\label{fig:scattermatrix2}"}
knitr::include_graphics("scattermatrix2.pdf")
```

Figure \ref{fig:scattermatrix2} presents a scatterplot matrix investigating the potential impact of macronutrient intake, including carbohydrates, proteins, fats, and cholesterol [@usda_2022_macronutrients], on blood pressure levels. Histograms of macronutrient distributions across all three blood pressure levels reveal right-skewed patterns, suggesting no single macronutrient significantly influences blood pressure.

Notably, the analysis highlights substantial correlations among macronutrient variables. The highest correlation is observed between protein intake (`DR1TPROT`) and dietary fat intake (`DR1TTFAT`), reaching 0.729. Additional pairs, such as protein intake (`DR1TPROT`) and cholesterol intake (`DR1TCHOL`) with a correlation of 0.684, indicate potential multicollinearity among predictor variables. This observation prompts caution when employing certain parametric modeling methods, such as multinomial regression, which may be sensitive to multicollinearity issues.

These findings lay the groundwork for a nuanced understanding of the dataset and underscore the importance of considering demographic and nutritional factors in predicting blood pressure levels. Subsequent sections will delve deeper into statistical analyses and modeling techniques to derive actionable insights from the presented visualizations.

# Methods

## Data Preparation (One-hot encoding and train test sets spliting)

Categorical predictors often require transformation into numerical format for compatibility with many machine learning algorithms. One-hot encoding is employed to convert categorical variables, `BPACSZ` (4 levels), `BPXPTY` (2 levels), `RIAGENDR` (2 levels), `RIDRETH3` (6 levels), `DIQ010` (4 levels), `DIQ050` (3 levels), and `HIQ011` (4 levels). This technique ensures that the categorical nature of the variables is preserved in the analysis. By applying one-hot encoding to these categorical predictors, our data frame have 91 columns of predictors in total.

In addition to one-hot encoding, we randomly selected 80% of our data (4,900 observations) as the
training set without replacement and the rest 20% of the data (1,225 observations) as the testing set. In this
way, we can evaluate our models objectively.

## Multinomial Regression with Lasso Regularization and 10-fold Cross Validation as a Baseline

Multinomial regression with lasso regularization is chosen as the baseline algorithm for its interpretability and efficacy in handling diverse data characteristics. The model's prediction for a data point $x_i$ is mathematically expressed as:

$$
\hat{y}_i = \text{argmax}_j \left( \frac{e^{\beta_{0j} + \beta_{1j}x_{1i} + \beta_{2j}x_{2i} + \ldots + \beta_{mj}x_{mi}}}{\sum_{k=0}^{N-1} e^{\beta_{0k} + \beta_{1k}x_{1i} + \beta_{2k}x_{2i} + \ldots + \beta_{mk}x_{mi}}} \right)
$$

Here, $Y$ represents the blood pressure level (`BPXLEVEL`) variable, $X$ is the feature matrix with number of columns $m = 91$, $\beta$ denotes the coefficients, and $N = 3$ is the number of classes. The addition of lasso regularization to the objective function introduces a penalty term:

$$
\text{argmax}_{\beta} \sum_{i=0}^{N - 1} \log P(Y = y_i | X) - \lambda \sum_{j=0}^{N - 1} |\beta_j|
$$

where $\lambda$ is the regularization parameter. With the lasso penalty, we will be able to identify which predictors are the most influential in predicting blood pressure types. Moreover, it can help prevent overfitting and control model complexity, potentially deal with the multicollinearity issue

The model's performance is systematically assessed through test accuracy and test ROC-AUC across the 10-fold cross-validation. This thorough evaluation ensures the reliability of our multinomial regression baseline in predicting blood pressure categories.

## eXtreme Gradient Boosting Model with 10-fold Cross Validation

eXtreme Gradient Boosting (XGBoost), stands out as a powerful ensemble learning method, widely recognized for its superior predictive capabilities. The algorithm systematically constructs a collection of weak learners, in the form of decision trees, and amalgamates their predictions to enhance accuracy and generalize well to unseen data.

In this analysis, XGBoost is strategically employed to surpass the baseline set by multinomial regression with lasso regularization. To control model complexity and improve robustness, XGBoost incorporates regularization terms.

The mathematical formulation of the XGBoost algorithm is as follows:

Given a training dataset $\{(x_i, y_i)\}_{i=1}^n$, where $x_i$ represents the predictors of the ith observation and $y_i$ is the corresponding label, XGBoost aims to learn an additive model $F(x)$ of the form:

$$F(x) = \sum_{m=1}^M \gamma_m h_m (x),$$

where
\begin{itemize}
    \item $M$ is the number of weak learners (trees) in the ensemble.
    \item $h_m(x)$ is the mth weak learner (tree).
    \item $\gamma_m$ is the weight (shrinkage) applied to the output of the mth weak learner.
\end{itemize}

The XGBoost algorithm minimizes the following objective function, which comprises a differentiable convex loss function $L(y_i, F_{m-1}(x_i) + \gamma h_m (x_i))$ and a regularization term $\Omega(F_m)$:

$$\mathcal{L}(\{(x_i, y_i)\}_{i=1}^n) = \sum_{i=1}^n L(y_i, F_{m-1}(x_i) + \gamma_m h_m(x_i)) + \Omega(F_m),$$

where:
\begin{itemize}
    \item $L$ is the loss function measuring the difference between predicted values and true labels.
    \item $F_{m-1}$ is the additive model up to the (m-1)th iteration.
    \item $\gamma_m$ is the optimal weight for the mth weak learner.
    \item $h_m(x_i)$ is the prediction of the mth weak learner for the ith observation.
    \item $\Omega(F_m)$ is a regularization term controlling the model's complexity, typically penalizing the number of leaves in the trees and the magnitude of the weights.
\end{itemize}

The objective function is optimized in a stage-wise manner. At each stage, a new weak learner is added to the ensemble by solving for $\gamma_m$ and $h_m (x_i)$, updating the model accordingly:

$$F_m (x) = F_{m-1} (x) + \gamma_m h_m (x)$$

The optimization involves finding values for $\gamma_m$ and the weak learner's parameters that minimize the objective function. This is commonly achieved using gradient boosting, iteratively fitting a weak learner to the negative gradient of the objective function.

The performance of the XGBoost model will be assessed based on test accuracy and the test ROC-AUC. Employing 10-fold cross validation ensures robust estimation of these metrics across different data subsets, enhancing the model's reliability and generalization capabilities.

## Feature Selection with XGBoost Feature Importance

One distinctive feature of XGBoost is its ability to provide valuable insights into feature importance. This is achieved through the computation of importance scores assigned to each predictor, utilizing Gain as the metric, which represents the improvement in accuracy attributed to a specific feature across the model's trees [@chen_2016_xgboost].

This analysis leverages the XGBoost-derived importance scores to discern the most influential predictors. The incorporation of these scores aims to enhance the predictive performance of the model, resulting in improved test accuracy and test ROC-AUC scores. The higher the importance score assigned to a feature, the more impactful it is considered in the overall predictive capacity of the model.

# Results

## Multinomial Regression Model with Lasso Regularization

The Multinomial regression model was trained with various lasso regularization strengths, spanning a range from low to high values, using a 10-fold cross-validation strategy. The model's multinomial deviance was documented for each regularization strength.

```{r sofmax-cv, warning=FALSE, message=FALSE, cache=TRUE, include=FALSE}
set.seed(88)
# Softmax classification with 10 fold CV
cvfit = cv.glmnet(train.X.mm, as.matrix(train.Y), family = "multinomial", 
                  type.multinomial = "grouped", parallel = TRUE,
                  trace.it = TRUE)
pdf(file = "softmax_cv_plot.pdf", width = 5.33, height = 4)
print(plot(cvfit))
dev.off()
```
```{r softmax-cv-plot, fig.align='center', out.width='150px', fig.cap='10-fold Cross Validation Findinig the Best Lasso Regularization Parameter ($\\lambda$)\\label{fig:sofmax-cv}', echo=FALSE, fig.pos = "ht"}
knitr::include_graphics("softmax_cv_plot.pdf")
```

Figure \ref{fig:sofmax-cv} helps identify the optimal $\lambda$ for the lasso penalty term of the multinomial regression. We aim to get the $\lambda$ which minimizes the multinomial deviance. By looking at the figure, we got the smallest multinomial deviance when $\lambda$ = `r cvfit$lambda.min` and only `r cvfit$index[1]` predictors were selected by the lasso penalty, which is indicated by the vertical dash line on the left.

```{r coefficient-table, include=FALSE}
coef.df = as.data.frame(cbind(as.matrix(coef(cvfit, s = cvfit$lambda.min)$`0`), 
                              as.matrix(coef(cvfit, s = cvfit$lambda.min)$`1`),
                              as.matrix(coef(cvfit, s = cvfit$lambda.min)$`2`)))
colnames(coef.df) = c("Coefficient (low - 0)", "Coefficient (elevated - 1)", "Coefficient (high - 2)")
coef.df = coef.df %>% filter_all(all_vars(. != 0))
```

\begin{table}[ht]
  \caption{Coefficients of 43 Selected Predictors and Intercept of the Multinomial Regression Model with Lasso Regularization}
  \centering
  \begin{tabular}{p{1.5cm}p{1.7cm}p{1.7cm}p{1.7cm}|p{1.5cm}p{1.7cm}p{1.7cm}p{1.7cm}}
    \toprule
    Predictor & Coefficient \newline(Y = 0) & Coefficient \newline(Y = 1) & Coefficient \newline(Y = 2) & Predictor & Coefficient \newline(Y = 0) & Coefficient \newline(Y = 1) & Coefficient \newline(Y = 2) \\
    \midrule
    Intercept & 2.2521848 & -0.6889857 & -1.5631991 & \texttt{BPACSZ2} & 0.7865830 & -0.4102865 & -0.3762966 \\
    \texttt{BPACSZ3} & 0.1767172 & -0.1796646 & 0.0029475 & \texttt{BPACSZ5} & -0.4005157 & 0.1731479 & 0.2273678 \\
    \texttt{BPXPLS} & -0.0032227 & -0.0000035 & 0.0032262 & \texttt{BPXPTY2} & 0.1264090 & 0.0395048 & -0.1659138 \\
    \texttt{RIAGENDR2} & 0.2307534 & -0.1249100 & -0.1058434 & \texttt{RIDAGEYR} & -0.0368332 & 0.0039225 & 0.0329107 \\
    \texttt{RIDRETH32} & 0.0241270 & -0.0054927 & -0.0186343 & \texttt{RIDRETH33} & 0.1693800 & -0.0519517 & -0.1174283 \\
    \texttt{RIDRETH34} & -0.1006645 & -0.0656122 & 0.1662767 & \texttt{RIDRETH36} & -0.0342307 & -0.0324576 & 0.0666883 \\
    \texttt{DMDHHSZA} & 0.0403702 & -0.0185643 & -0.0218059 & \texttt{DMDHHSZB} & 0.0348970 & -0.0227166 & -0.0121805 \\
    \texttt{DMDHHSZE} & 0.0189516 & 0.0594003 & -0.0783519 & \texttt{DR1TNUMF} & 0.0064755 & 0.0015578 & -0.0080333 \\
    \texttt{DR1TSUGR} & -0.0000212 & -0.0001419 & 0.0001631 & \texttt{DR1TCHOL} & -0.0000276 & 0.0000052 & 0.0000225 \\
    \texttt{DR1TATOC} & 0.0000225 & 0.0000252 & -0.0000478 & \texttt{DR1TACAR} & -0.0000038 & 0.0000000 & 0.0000038 \\
    \texttt{DR1TLYCO} & 0.0000002 & -0.0000001 & -0.0000001 & \texttt{DR1TNIAC} & -0.0000822 & -0.0000403 & 0.0001225 \\
    \texttt{DR1TVB6} & -0.0030378 & -0.0047378 & 0.0077756 & \texttt{DR1TVB12} & -0.0018620 & 0.0012349 & 0.0006271 \\
    \texttt{DR1TVC} & -0.0000052 & 0.0000135 & -0.0000082 & \texttt{DR1TVD} & 0.0074843 & -0.0026119 & -0.0048724 \\
    \texttt{DR1TZINC} & -0.0004469 & 0.0008037 & -0.0003569 & \texttt{DR1TPOTA} & 0.0000234 & -0.0000043 & -0.0000191 \\
    \texttt{DR1TALCO} & -0.0014905 & 0.0015632 & -0.0000727 & \texttt{DR1TMOIS} & -0.0000438 & 0.0000117 & 0.0000322 \\
    \texttt{DR1TS060} & 0.1333655 & -0.0404757 & -0.0928898 & \texttt{DR1TS080} & 0.0251928 & -0.1485159 & 0.1233232 \\
    \texttt{DR1TS120} & -0.0010639 & -0.0023929 & 0.0034568 & \texttt{DR1TS180} & -0.0027631 & -0.0002637 & 0.0030268 \\
    \texttt{DR1TM161} & -0.0173735 & 0.0056004 & 0.0117731 & \texttt{DR1TM201} & -0.0351797 & 0.0157894 & 0.0193903 \\
    \texttt{DR1TM221} & -0.0781895 & 0.0435119 & 0.0346776 & \texttt{DR1TP184} & -0.3719599 & -0.1414094 & 0.5133693 \\
    \texttt{DR1TP204} & -0.3937469 & 0.2857126 & 0.1080342 & \texttt{DIQ0102} & 0.0059360 & -0.0490002 & 0.0430642 \\
    \texttt{DIQ0502} & -0.0978138 & 0.0281761 & 0.0696377 & \texttt{HIQ0112} & -0.0250448 & -0.0130697 & 0.0381145 \\
    \texttt{HIQ0117} & 0.5921226 & -0.3291087 & -0.2630140 & \texttt{HIQ0119} & -0.0376808 & 0.0126395 & 0.0250413 \\
    \bottomrule
  \end{tabular}
  \label{tab:multinom_coefficients}
\end{table}

Table \ref{tab:multinom_coefficients} shows all coefficients of our 43 predictors and the intercept of the multinomial regression model. The multinomial regression model can be represented as conditional probabilities with respect to each blood pressure level as follows.

$$
\begin{aligned}
\Pr(Y = 0 | X) &= \frac{\exp(2.2521848 + 0.1767172x_{\texttt{BPACSZ3}} - 0.0032227x_{\texttt{BPXPLS}} + \cdots -0.0376808x_{\texttt{HIQ0119}})}{\sum_{i = 0}^{N-1} \exp(\hat{\beta}_{0i} + \hat{\beta}_{1i}x_{\texttt{BPACSZ3}} + \hat{\beta}_{2i}x_{\texttt{BPXPLS}} + \cdots + \hat{\beta}_{43i}x_{\texttt{HIQ0119}})}, \\
\Pr(Y = 1 | X) &= \frac{\exp(-0.6889857 - 0.1796646x_{\texttt{BPACSZ3}} - 0.0000035x_{\texttt{BPXPLS}} + \cdots + 0.0126395x_{\texttt{HIQ0119}})}{\sum_{i = 0}^{N-1} \exp(\hat{\beta}_{0i} + \hat{\beta}_{1i}x_{\texttt{BPACSZ3}} + \hat{\beta}_{2i}x_{\texttt{BPXPLS}} + \cdots + \hat{\beta}_{43i}x_{\texttt{HIQ0119}})}, \\
\Pr(Y = 2 | X) &= \frac{\exp(-1.5631991 + 0.0029475x_{\texttt{BPACSZ3}} +  0.0032262x_{\texttt{BPXPLS}} + \cdots + 0.0250413x_{\texttt{HIQ0119}})}{\sum_{i = 0}^{N-1} \exp(\hat{\beta}_{0i} + \hat{\beta}_{1i}x_{\texttt{BPACSZ3}} + \hat{\beta}_{2i}x_{\texttt{BPXPLS}} + \cdots + \hat{\beta}_{43i}x_{\texttt{HIQ0119}})},
\end{aligned}
$$

where the sum of the three conditional probabilities should be 1 for every input $X \in \mathbb{R}^{1\times 43}$ and

$$
\begin{aligned}
&\sum_{i = 0}^{N-1} \exp(\hat{\beta}_{0i} + \hat{\beta}_{1i}x_{\texttt{BPACSZ3}} + \hat{\beta}_{2i}x_{\texttt{BPXPLS}} + \cdots + \hat{\beta}_{43i}x_{\texttt{HIQ0119}}) \\
=& \exp(2.2521848 + 0.1767172x_{\texttt{BPACSZ3}} - 0.0032227x_{\texttt{BPXPLS}} + \cdots -0.0376808x_{\texttt{HIQ0119}}) + \\
 & \exp(-0.6889857 - 0.1796646x_{\texttt{BPACSZ3}} - 0.0000035x_{\texttt{BPXPLS}} + \cdots + 0.0126395x_{\texttt{HIQ0119}}) + \\
 & \exp(-1.5631991 + 0.0029475x_{\texttt{BPACSZ3}} +  0.0032262x_{\texttt{BPXPLS}} + \cdots + 0.0250413x_{\texttt{HIQ0119}}).
\end{aligned}
$$

```{r softmax-result, include=FALSE, warning=FALSE, message=FALSE}
cvfit$lambda.min
cvfit$index
softmax.predictions = predict(cvfit, newx = test.X.mm, 
                              s = cvfit$lambda.min, type = "class")
mean(softmax.predictions == test.Y)
# ROC - AUC
softmax.rocr = multiclass.roc(as.numeric(test.Y) - 1, 
                              as.numeric(softmax.predictions) - 1)
softmax.rocr$auc
confusionMatrix(as.factor(softmax.predictions), test.Y)
```

Upon inputting the test set data to the multinomial regression model, we achieved a test accuracy of 68% and the test ROC-AUC was calculated as 0.6791, demonstrating the model's robust discriminative ability across multiple classes. These results establish a baseline for our future modeling efforts, showcasing the effectiveness of the multinomial regression approach in capturing and understanding underlying patterns within the dataset.

## XGBoost Model

The XGBoost model was trained using eXtreme Gradient Boosting with exact tree method, a powerful ensemble learning method. The following hyperparameters were modified and utilized in the model:

```{r xgb-model, warning=FALSE, message=FALSE, cache=TRUE, include=FALSE}
# Gradient boosting with 10 fold CV
xgboost.train = function(col) {
  set.seed(88)
  xgb.data = xgb.DMatrix(data = as.matrix(train.X.dummy[col]), 
                         label = recode(train.X.dummy$BPXLEVEL, 
                                        '0'=0, '1'=1, '2'=2))
  xgb.test.X = data.matrix(test.X.dummy[col])
  hyperparameters = list(
    eta = 0.02,
    subsample = 0.75,
    col_subsample = 0.8,
    max_depth = 10
  )
     
  params = list(
    eta = hyperparameters$eta,
    subsample = hyperparameters$subsample,
    colsample_bytree = hyperparameters$col_subsample,  
    max_depth = hyperparameters$max_depth,
    tree_method = "exact",
    objective = "multi:softmax",
    num_class = 3
  )
     
  cv.xgb = xgb.cv(
    params = params,
    data = xgb.data,
    nfold = 10,
    metrics = "mlogloss",
    verbose = 1,
    nrounds = 100
  )
    
  eval.log = as.data.frame(cv.xgb$evaluation_log)
  min.merror = min(eval.log[, 4])
  min.merror.index = which.min(eval.log[, 4])
  xgb.model = xgboost(params = params, 
                      data = xgb.data, 
                      nrounds = min.merror.index, 
                      verbose = 0)
  xgb.predictions = predict(xgb.model, xgb.test.X)
  return(list(xgb.model, xgb.predictions, eval.log))
}
xgboost.all = xgboost.train(colnames(train.X.dummy)[2:92])
mean(xgboost.all[[2]] == test.Y)
multiclass.roc(as.numeric(test.Y) - 1, 
               as.numeric(xgboost.all[[2]]) - 1)$auc
confusionMatrix(as.factor(xgboost.all[[2]]), test.Y)
importance.matrix = xgb.importance(colnames(train.X.dummy)[2:92], model = xgboost.all[[1]])
importance.matrix
```

- Learning Rate (eta): 0.02

- Subsample: 0.75

- Column Subsample: 0.8

- Maximum Depth: 10

- Number of Trees (Rounds): 100

Figure \ref{fig:xgb-train} visually portrays the evolution of training and validation losses over a span of 100 rounds. The training loss, denoted by the blue line, manifests a persistent and incremental reduction throughout the training iterations. Analogously, the validation loss, depicted by the orange line, likewise manifests a sustained and systematic decrease, signifying enhanced generalization performance of the model on previously unseen data.

```{r xgb-loss-plot, include=FALSE}
xgb.loss.plot = ggplot(data = xgboost.all[[3]]) + 
  geom_line(mapping = aes(x = iter, y = train_mlogloss_mean, color = "Train")) +
  geom_line(mapping = aes(x = iter, y = test_mlogloss_mean, color = "Validation")) + 
  scale_color_manual(values = c("blue", "orange"), labels = c("Train", "Validation")) +
  labs(color = "Loss Type") + 
  xlab("Number of Rounds") + 
  ylab("Cross-entropy Loss")
pdf(file = "xgb_loss.pdf", width = 6, height = 2.5)
print(xgb.loss.plot)
dev.off()
```
```{r xgb-loss-plot-embed, fig.align='center', out.width='200px', fig.cap='Training and Validation Losses of XGBoost Model with 10-fold Cross Validation over 100 Rounds of Training\\label{fig:xgb-train}', echo=FALSE, fig.pos = "ht"}
knitr::include_graphics("xgb_loss.pdf")
```

With these hyperparameters, we trained the model with 10-fold cross validation to get a test accuracy of `r round(mean(xgboost.all[[2]] == test.Y) * 100, 5)`% and a test ROC-AUC of 0.6898. Both of these metrics are a little bit higher than that of the multinomial regression model, which indicates that the XGBoost model is slightly better than multimonial regression model with lasso regularization when classifying the blood pressure levels.

## XGBoost Model with Selected Predictors

Figure \ref{fig:importance-plot} shows the Gain scores of the predictors used in the XGBoost model. A higher bar (higher Gain score) represents more important the predictor is. Notably, key predictors such as `RIDAGEYR` (age), `DR1TMOIS` (moisture intake), `BPXPLS` (pulse rate in a minute), `DR1TSUGR` (total sugar intake), and `DR1TCAFF` (caffeine intake) emerged as significant contributors to the predictive power of the model.

```{r importance-plots1, include=FALSE}
importance.plot = xgb.ggplot.importance(importance.matrix, measure = 'Gain') + 
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size = 8),
        axis.text.y = element_text(size = 10)) + 
  geom_bar(aes(fill = "blue"), stat = "identity") +
  ylab("Gain") + coord_cartesian() + ggtitle("")
pdf(file = "importance_plot.pdf", width = 10, height = 3.5)
print(importance.plot)
dev.off()
```
```{r importance-plots2, echo = FALSE, , fig.align='center', out.width='350px', fig.cap='Bar Plots of Gain Score of Each Feature in the XGBoost Model\\label{fig:importance-plot}'}
knitr::include_graphics("importance_plot.pdf")
```

In our pursuit of refining the model and unveiling the most impactful predictors, we executed a meticulous feature selection process. We initiated this process by systematically eliminating predictors, starting with the least important (the one with the lowest Gain score), and subsequently assessed the impact on both test accuracy and test ROC-AUC. This methodical stepwise elimination allowed us to pinpoint a subset of predictors that consistently upheld optimal predictive performance. During this process, we keep using the same hyperparameters we used in the original XGBoost model with 10-fold cross validation at each step.

The results of this feature selection journey revealed a compelling trade-off between the number of predictors and predictive accuracy. Significantly, in figure \ref{fig:xgb-select}, the model showcased a remarkable test accuracy of 68.81633% and a test ROC-AUC of 0.6914502 even with just the top 46 most important predictors. This underscores the efficiency of the selected predictors in encapsulating crucial information for the accurate prediction of health outcomes. As indicated by the red dash lines in figure \ref{fig:xgb-select}, the model achieved the highest test accuracy of 68.89796% and the highest test ROC-AUC of 0.6928988 with the top 64 most important predictors.

```{r xgb-feature-select, warning=FALSE, message=FALSE, cache=TRUE, include=FALSE}
threshold = sort(importance.matrix$Gain)
n.feature = seq(length(importance.matrix$Gain))
threshold.importance = c()
test.accuracy = c()
test.auc = c()
for (i in 1:length(threshold)) {
  selected.features = importance.matrix$Feature[which(importance.matrix$Gain >= threshold[i])]
  xgboost.all = xgboost.train(selected.features)
  accuracy = mean(xgboost.all[[2]] == test.Y)
  auc = multiclass.roc(as.numeric(test.Y) - 1, 
                       as.numeric(xgboost.all[[2]]) - 1)$auc
  threshold.importance[i] = threshold[i]
  test.accuracy[i] = accuracy
  test.auc[i] = auc
}
```

```{r xgb-plot, warning=FALSE, include=FALSE}
df.xgb = data.frame(n.feature = sort(n.feature, decreasing = TRUE),
                    threshold.importance = threshold.importance,
                    test.accuracy = test.accuracy,
                    test.auc = test.auc)
xgb.accuracy = ggplot(df.xgb, aes(x = n.feature, y = test.accuracy)) +
  geom_line() +
  geom_vline(xintercept = df.xgb$n.feature[which.max(df.xgb$test.accuracy)], 
             linetype = "dashed", color = "red") +
  geom_hline(yintercept = mean(df.xgb$test.accuracy), linetype = "dashed", color = "blue") +
  geom_text(aes(label = sprintf("Avg Accuracy: %.3f", mean(df.xgb$test.accuracy))),
            x = max(df.xgb$threshold.importance), y = mean(df.xgb$test.accuracy), 
            vjust = 1, hjust = -1.5, color = "blue") +
  geom_text(aes(label = sprintf("Number of predictors: %d", df.xgb$n.feature[which.max(df.xgb$test.accuracy)])),
            x = df.xgb$n.feature[which.max(df.xgb$test.accuracy)], 
            y = max(df.xgb$test.accuracy), 
            vjust = 20, hjust = 1, color = "red") +
  labs(x = "Number of Predictors",
       y = "Accuracy")
xgb.auc = ggplot(df.xgb, aes(x = n.feature, y = test.auc)) +
  geom_line() +
  geom_vline(xintercept = df.xgb$n.feature[which.max(df.xgb$test.auc)], 
             linetype = "dashed", color = "red") +
  geom_hline(yintercept = mean(df.xgb$test.auc), linetype = "dashed", color = "blue") +
  geom_text(aes(label = sprintf("Avg AUC: %.3f", mean(df.xgb$test.auc))),
            x = max(df.xgb$threshold.importance), y = mean(df.xgb$test.auc), 
            vjust = 1, hjust = -2.5, color = "blue") +
  geom_text(aes(label = sprintf("Number of predictors: %d", df.xgb$n.feature[which.max(df.xgb$test.auc)])),
            x = df.xgb$n.feature[which.max(df.xgb$test.auc)], 
            y = max(df.xgb$test.auc), 
            vjust = 20, hjust = 1, color = "red") +
  labs(x = "Number of Predictors",
       y = "AUC")
xgb.plots = (xgb.accuracy | xgb.auc)
pdf(file = "xgb_plots.pdf", width = 10, height = 4)
print(xgb.plots)
dev.off()
```
```{r xgb-plot-embed, echo = FALSE, out.width = "400px", fig.align = "center", fig.cap = "XGBoost Model Test Accuracy and ROC-AUC from 1 Predictor to 91 Predictors\\label{fig:xgb-select}"}
knitr::include_graphics("xgb_plots.pdf")
```

Table \ref{tab:xgb-select-table} provides a comprehensive overview of the top predictors identified by the XGBoost model with test accuracy higher than or equal to 68.4898% and test ROC-AUC higher than or equal to 0.6898, presenting their corresponding threshold Gain scores, test accuracy, and test ROC-AUC values.

\begin{table}[ht]
  \caption{XGBoost Model Top Predictors and Performance Metrics}
  \centering
  \begin{tabular}{p{2.2cm}p{2cm}p{2.2cm}p{2.4cm}}
    \toprule
    Number of \newline Top Predictors & Threshold Gain Score & Test Accuracy & Test ROC-AUC \\
    \midrule
    64 & 0.0061172 & 68.89796\% & 0.6928988 \\
    58 & 0.0072728 & 68.89796\% & 0.6919519 \\
    46 & 0.0091180 & 68.81633\% & 0.6914502 \\
    53 & 0.0084463 & 68.89796\% & 0.6913785 \\
    50 & 0.0088765 & 68.73469\% & 0.6908387 \\
    56 & 0.0076933 & 68.4898\% & 0.6904320 \\
    \bottomrule
  \end{tabular}
  \label{tab:xgb-select-table}
\end{table}

The table is thoughtfully organized, with entries sorted based on descending test ROC-AUC, prioritizing higher ROC-AUC models. In cases of ties, the sorting is further refined by considering descending test accuracy values and, if necessary, the top number of predictors in ascending order.

```{r xgb-select-coef-table, include=FALSE}
xgb.log = df.xgb %>% 
  arrange(desc(test.auc), desc(test.accuracy), n.feature)
xgb.log$test.accuracy = as.numeric(xgb.log$test.accuracy)
xgb.log$test.accuracy = round(xgb.log$test.accuracy * 100, 5)
xgb.log$test.accuracy = paste0(xgb.log$test.accuracy, "%")
colnames(xgb.log) = c("Number of Predictors",
                      "Threshold Gain Score",
                      "Test Accuracy", 
                      "Test ROC-AUC")
head(xgb.log, 6)
```

## Model Selection

Table \ref{tab:model-summary} provides a summary of the performance of multinomial regression model with lasso penalty, XGBoost model with all predictors, and XGBoost model with selected predictors. As our focus lies on the accuracy and ROC-AUC metrics, and, based on these, the model with the top 64 most important predictors stands out as the preferred choice. This model exhibits higher increases in accuracy, with 0.89796 percentage points and 0.40816 percentage points improvements compared to the multinomial regression model with lasso regularization and the XGBoost model using all predictors, respectively. Moreover, it demonstrates 0.0137988 and 0.0030988 increases in ROC-AUC compared to the multinomial regression model with lasso regularization and the XGBoost model using all predictors, respectively.

\begin{table}[ht]
  \caption{Summary of Three Models’ Test Accuracies and ROC-AUCs}
  \centering
  \begin{tabular}{p{5.5cm}p{2.5cm}p{2.5cm}}
    \toprule
    Model & Test Accuracy & Test ROC-AUC \\
    \midrule
    Multinomial Regression with Lasso & 68\% & 0.67910 \\
    XGBoost Using 91 Predictors & 68.4898\% & 0.6898 \\
    XGBoost Using 64 Predictors & 68.89796\% & 0.6928988 \\
    \bottomrule
  \end{tabular}
  \label{tab:model-summary}
\end{table}

Among these 64 selected predictors, there are 29 predictors selected by both the multinomial model and the XGBoost model. These predictors are `RIDAGEYR`, `DR1TMOIS`, `BPXPLS`, `DR1TSUGR`, `DR1TLYCO`, `RIAGENDR2`, `DR1TCHOL`, `DR1TVC`, `DR1TP204`, `DMDHHSZE`, `BPACSZ5`, `DR1TATOC`, `DR1TM221`, `DR1TM201`, `DR1TACAR`, `DR1TM161`, `DR1TNIAC`, `DR1TVD`, `DR1TNUMF`, `DR1TVB12`, `DR1TVB6`, `DR1TZINC`, `DR1TS120`, `DMDHHSZB`, `DR1TPOTA` `DR1TS180`, `DR1TS060`, `DR1TS080`, and `DR1TALCO`. As both models selected these predictors, indicating the effectiveness of these predictors on predicting blood pressure levels.

```{r predictor-finding, include=FALSE}
intersection = intersect(unname(unlist(importance.matrix[1:64, 1])),
                         rownames(coef.df))
intersection
setdiff(unname(unlist(importance.matrix[1:64, 1])), intersection)
```

Our systematic approach to feature selection not only fine-tuned the model but also provided insightful perspectives on the pivotal factors influencing its predictive power. This enhanced interpretability contributes to a more robust and effective health outcome prediction system.

# Conclusions

Despite the advancements made in developing predictive models, it's crucial to acknowledge certain limitations. One prominent drawback is the challenge of achieving high accuracy and ROC-AUC, particularly in the context of health-related predictions. Accurate blood pressure classification is paramount for providing meaningful health insights, and any inaccuracies in predictions could have significant implications. Notably, discrepancies in predicting health outcomes can impact the reliability of personalized recommendations and interventions, potentially leading to suboptimal health management.

Numerous studies emphasize the importance of accuracy in health-related predictive models. For instance, a study by @sofogianni_2022_cardiovascular highlighted the critical role of accurate predictions in cardiovascular risk assessment models, underscoring the potential consequences of misclassification on patient care. Additionally, research conducted by @grover_2014_an emphasized the need for robust predictive models in chronic disease management, as inaccuracies can compromise the effectiveness of preventive measures and early interventions. These findings underscore the broader concern within the scientific community about the implications of suboptimal accuracy in health-related predictions.

Addressing the aforementioned drawbacks requires a multi-faceted approach. Feature engineering, the process of refining and creating new predictors, could enhance the models' ability to capture intricate patterns in the data, potentially boosting predictive performance. Additionally, acquiring more high-quality data, especially with a focus on diverse demographic groups and health conditions, could contribute to a more comprehensive and representative model. Exploring advanced machine learning techniques, such as deep learning methods like neural networks, holds promise in uncovering complex relationships within the data, potentially elevating predictive accuracy.

\begin{figure}[ht]
  \centering
  \begin{subfigure}{\textwidth}
    \centering
    \includegraphics[width=155px]{web1.png}
    \includegraphics[width=155px]{web2.png}
  \end{subfigure}
  \caption{Website Interface Backend with the Multinomial Regression Model with Lasso Regularization}
  \label{fig:web}
\end{figure}

An exciting application of our predictive models lies in the integration with health apps, such as Apple Health, Samsung Health, and so on. Implementing our models in these platforms could empower individuals to receive personalized daily blood pressure suggestions based on their recorded dietary intakes, known health conditions, and demographic information. To demonstrate, we deployed the multionomial regression model with JavaScript into a website interface at https://bpmodel.ly.gd.edu.kg/, where users can input their demographic, health, and dietary information to get a blood pressure prediction shown in figure \ref{fig:web}. This practical application could serve as a proactive tool for users to manage their health more effectively, offering real-time insights and guidance.

In conclusion, while our predictive models showcase promising results, there is ongoing work to be done in refining their accuracy and applicability. By addressing the identified drawbacks through feature engineering, data enrichment, and the exploration of advanced machine learning techniques, we can move closer to developing highly reliable and impactful predictive models for blood pressure classification. The envisioned integration with health apps presents an exciting avenue for translating our research into actionable insights, fostering proactive health management among individuals.

# Computational Details

The analysis was conducted using R version 4.3.2 for Windows, with the utilization of various R libraries from the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/ to facilitate data manipulation, statistical modeling, and visualization. The following R libraries were employed in this study:

- `caret`: for classification and regression training.
- `dplyr`: for data manipulation and summarization.
- `GGally`: for extension to ggplot2 for correlation plots.
- `ggpubr`: for creating publication-ready plots with ggplot2.
- `glmnet`: for fitting generalized linear models with regularization.
- `grid`: for arranging and combining multiple plots.
- `gridExtra`: for arranging and combining multiple plots.
- `haven`: for reading and writing SPSS, Stata, and SAS files.
- `knitr`: for dynamic report generation in R Markdown.
- `patchwork`: for arranging and combining multiple plots.
- `pROC`: for analyzing ROC curves and assessing model performance.
- `tidyr`: for data tidying and reshaping.
- `xgboost`: for extreme gradient boosting.

The analyses were conducted in the RStudio integrated development environment (IDE) version "Mountain Hydrangea" Release (583b465e, 2023-06-05) for Windows. RStudio can be downloaded at https://posit.co/.

An Intel-compatible 64-bit platform is preferred. At least 2048 MB of RAM is recommended to run the whole script. An operating system of Windows 7 or higher or Mac OS X 10.6 or higher is preferred.

# Reproducibility

Ensuring the reproducibility of this study is of utmost importance. The entire analysis, including data preprocessing, model development, and result generation, is encapsulated in an RMarkdown document. The RMarkdown file, along with the necessary BibTeX and style files, has been made available on GitHub for easy access and replication: https://github.com/lygitdata/bpmodel/.

The RMarkdown file and its relevant files can be downloaded at the following link:

https://bpmodel.ly.gd.edu.kg/manuscript/download.zip

To reproduce the findings and generate the same results presented in this paper, follow these steps:

1. Download the Necessary Files:

  - Navigate to the provided link in your browser.
  - Unzip the downloaded file to a directory of your choice.

2. Open RMarkdown in RStudio:

  - Ensure you have R and RStudio installed on your machine.
  - Open RStudio and navigate to the directory where you unzipped the files.
  - Open the RMarkdown file (`manuscript.Rmd`) in RStudio.

3. Install Required Packages:

  - If not already installed, install the required R packages from the CRAN.

4. Knit the Document:

  - Knit the RMarkdown file to reproduce the analysis. This will execute the code chunks, perform the analysis, and generate the final document.

By following these steps, you can recreate the entire analysis and verify the results presented in this paper. This approach ensures transparency and allows others to validate and build upon the findings of this study.

# Acknowledgement

During the preparation of this work, GPT-4 was used in order to enhance the linguistic quality of some sections. After using this tool, the paper was reviewed and edited the content as needed. The author takes full responsibility for the content of the published article.

\clearpage

# Reference