STA Project 1 - MLR & Bootstrapping

1 Introduction

This project will demonstrate the use of variable selection to determine the final model. Bootstrapping procedures will be used to generate estimated confidence intervals for the coefficients of the regression model identified in a previous assignment. Specifically a model that looks at how different variables like age, hypertension, and diabetes can help predict a patients bmi.

1.1 Description of the Dataset

The data set being used is called Diabetes Prediction. The data set contains over 100,000 patient observations on both their medical and demographic information. This includes a patients diabetic status, hypertension status, heart disease status, as either having or not having the medical ailment. The full in order list of the data sets variables is as follows,

gender: chr

age: num

hypertension: int

heart_disease: int

smoking_history: chr

bmi: num

HbA1c_level: num

blood_glucose_level: int

diabetes: int

2 Materials

2.1 Data Preparation

url = "https://raw.githubusercontent.com/ChloeWinters79/STA321/refs/heads/main/Data/diabetes_prediction_dataset.csv"
bmi.diabetes = read.csv(url, header = TRUE)

HbA1c <- bmi.diabetes$HbA1c_level
Diabetic <- bmi.diabetes$diabetes
BloodGlucose <- bmi.diabetes$blood_glucose_level


prediabetic = (HbA1c > 5.69)  & (BloodGlucose > 99) & (Diabetic < 1 ) #create the variable prediabetic from those with blood glucose and HbA1c levels that fall above the normal range (determined by the American Diabetes Association) for non diabetic patients. 
bmi.diabetes$prediabetic = as.character(prediabetic) #convert from logic to character

The variables HbA1c_level, blood_glucose_level, and diabetes are used to define the variable prediabetic. The variable prediabetic is defined as follows, prediabetic = TRUE if HbA1c_level > 5.69, blood_glucose_level > 99 & diabetes < 1; prediabetic = FALSE otherwise.

American Diabetes Association (ADA) published that a normal HbA1c level is under 5.7, while a prediabetic HbA1c is between 5.7 and 6.5, and anything above 6.5 is considered diabetic. Additionally, the ADA finds that blood glucose (bg) levels under 100 are withing the normal range, bg levels from 100 to 125 are prediabetic, and bg levels 126 or over are diabetic.

These parameters from the ADA were used to create the variable prediabetic. Any patient with HbA1c and blood glucose levels outside of the normal range that are not diagnosed diabetics were flagged at prediabetic patients.

final.data = bmi.diabetes[, -c(1,5)] #remove gender and smoking_history, the variables not used in the initial candidate model
kable(head(final.data))

age	hypertension	heart_disease	bmi	HbA1c_level	blood_glucose_level	diabetes	prediabetic
80	0	1	25.19	6.6	140	0	TRUE
54	0	0	27.32	6.6	80	0	FALSE
28	0	0	27.32	5.7	158	0	TRUE
36	0	0	23.45	5.0	155	0	FALSE
76	1	1	20.14	4.8	155	0	FALSE
20	0	0	27.32	6.6	85	0	FALSE

bmi.log.age = lm((bmi)^-0.5 ~ log(age) + hypertension + heart_disease +  HbA1c_level
       + blood_glucose_level + diabetes + prediabetic, data = final.data)

kable(summary(bmi.log.age)$coef, caption = "Inferential Statistics of Model")

Inferential Statistics of Model

	Estimate	Std. Error	t value	Pr(>|t|)
(Intercept)	0.2421664	0.0005924	408.7822241	0.0000000
log(age)	-0.0133497	0.0000692	-192.8523791	0.0000000
hypertension	-0.0027677	0.0002457	-11.2636679	0.0000000
heart_disease	0.0040241	0.0003292	12.2222941	0.0000000
HbA1c_level	0.0000510	0.0000879	0.5801947	0.5617846
blood_glucose_level	-0.0000001	0.0000019	-0.0617243	0.9507825
diabetes	-0.0079106	0.0003583	-22.0761603	0.0000000
prediabeticTRUE	-0.0000032	0.0001923	-0.0163902	0.9869231

Considering the fact that three variables have high p-values and low t-values, it seems like some backwards selection is necessary before moving forward. In this case we are going to use the t-statistic to decide which value to throw away during each step. The variable selection will stop when all variables have a t-value higher than 1. Since blood_glucose_level has the smallest t-value we will start by removing that variable from the model.

bmi.log.age.6 = lm((bmi)^-0.5 ~ log(age) + hypertension + heart_disease +  HbA1c_level
                 + diabetes + prediabetic, data = final.data)
summary(bmi.log.age.6)

## 
## Call:
## lm(formula = (bmi)^-0.5 ~ log(age) + hypertension + heart_disease + 
##     HbA1c_level + diabetes + prediabetic, data = final.data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.114252 -0.011085  0.000521  0.011184  0.132117 
## 
## Coefficients:
##                   Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)      2.421e-01  4.726e-04  512.317   <2e-16 ***
## log(age)        -1.335e-02  6.922e-05 -192.854   <2e-16 ***
## hypertension    -2.768e-03  2.457e-04  -11.264   <2e-16 ***
## heart_disease    4.024e-03  3.292e-04   12.223   <2e-16 ***
## HbA1c_level      5.265e-05  8.362e-05    0.630    0.529    
## diabetes        -7.923e-03  3.003e-04  -26.382   <2e-16 ***
## prediabeticTRUE -8.424e-06  1.723e-04   -0.049    0.961    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01974 on 99993 degrees of freedom
## Multiple R-squared:  0.3084, Adjusted R-squared:  0.3084 
## F-statistic:  7433 on 6 and 99993 DF,  p-value: < 2.2e-16

There are still variables with a t-statistic less than zero so the backwards selection continues until no more variable with a low t-statistic remain.

bmi.log.age.5 = lm((bmi)^-0.5 ~ log(age) + hypertension + heart_disease +  HbA1c_level
                 + diabetes, data = final.data)
summary(bmi.log.age.5)

## 
## Call:
## lm(formula = (bmi)^-0.5 ~ log(age) + hypertension + heart_disease + 
##     HbA1c_level + diabetes, data = final.data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.114255 -0.011083  0.000519  0.011182  0.132114 
## 
## Coefficients:
##                 Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)    2.422e-01  4.206e-04  575.682   <2e-16 ***
## log(age)      -1.335e-02  6.922e-05 -192.855   <2e-16 ***
## hypertension  -2.768e-03  2.457e-04  -11.264   <2e-16 ***
## heart_disease  4.024e-03  3.292e-04   12.223   <2e-16 ***
## HbA1c_level    5.000e-05  6.364e-05    0.786    0.432    
## diabetes      -7.915e-03  2.533e-04  -31.248   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01974 on 99994 degrees of freedom
## Multiple R-squared:  0.3084, Adjusted R-squared:  0.3084 
## F-statistic:  8920 on 5 and 99994 DF,  p-value: < 2.2e-16

# Output shows 1 remaining variable with a t-statistic under 1 so we continue with backwards selection. 

bmi.log.age.4 = lm((bmi)^-0.5 ~ log(age) + hypertension + heart_disease
                 + diabetes, data = final.data)
summary(bmi.log.age.4)

## 
## Call:
## lm(formula = (bmi)^-0.5 ~ log(age) + hypertension + heart_disease + 
##     diabetes, data = final.data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.114220 -0.011113  0.000545  0.011190  0.132144 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.424e-01  2.422e-04 1000.92   <2e-16 ***
## log(age)      -1.335e-02  6.922e-05 -192.86   <2e-16 ***
## hypertension  -2.767e-03  2.457e-04  -11.26   <2e-16 ***
## heart_disease  4.024e-03  3.292e-04   12.22   <2e-16 ***
## diabetes      -7.838e-03  2.336e-04  -33.55   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01974 on 99995 degrees of freedom
## Multiple R-squared:  0.3084, Adjusted R-squared:  0.3084 
## F-statistic: 1.115e+04 on 4 and 99995 DF,  p-value: < 2.2e-16

Now that all the variables have a t-statistic above 1, we can continue onto determining the final model.

bmi.coef.4 <- summary(bmi.log.age.4)$coef
kable(summary(bmi.log.age.4)$coef, caption = "Inferential Statistics of Final Model")

Inferential Statistics of Final Model

	Estimate	Std. Error	t value	Pr(>|t|)
(Intercept)	0.2424251	0.0002422	1000.92041	0
log(age)	-0.0133499	0.0000692	-192.85792	0
hypertension	-0.0027673	0.0002457	-11.26224	0
heart_disease	0.0040239	0.0003292	12.22173	0
diabetes	-0.0078379	0.0002336	-33.55337	0

The final model is,

\[ \text{bmi} ^ {-0.5} = 0.2424251 - 0.0133499\times \log(\text{age}) - 0.0027673\times \text{hypertension} + 0.0040239\times \text{heart_disease} - 0.0078379\times \text{diabetes} \]

3 Methodology and Analysis

3.1 Bootstrapping

Now that the final model has been determine, bootstrapping can be used to help find the confidence intervals for the coefficients in the final regression model.

B = 1000       # choose the number of bootstrap replicates.
## 
num.p = dim(model.frame(bmi.log.age.4))[2]  # returns number of parameters in the model
smpl.n = dim(model.frame(bmi.log.age.4))[1] # sample size
## zero matrix to store bootstrap coefficients 
coef.mtrx = matrix(rep(0, B*num.p), ncol = num.p)       
## 
for (i in 1:B){
  bootc.id = sample(1:smpl.n, smpl.n, replace = TRUE) # fit final model to the bootstrap sample
  bmi.log.age.boot = lm((bmi)^-0.5 ~ log(age) + hypertension + heart_disease
                 + diabetes, data = final.data[bootc.id,])     
  coef.mtrx[i,] = coef(bmi.log.age.boot)}    # extract coefs from bootstrap regression model    

3.1.1 Histograms

An R function was defined to make the histograms of the bootstrap regression coefficients. This function is also used later on to create the histograms for the residual bootstrap estimated regression coefficients.

boot.hist = function(bmi.coef.4, bt.coef.mtrx, var.id, var.nm){
  x1.1 <- seq(min(bt.coef.mtrx[,var.id]), max(bt.coef.mtrx[,var.id]), length=300 )
  y1.1 <- dnorm(x1.1, mean(bt.coef.mtrx[,var.id]), sd(bt.coef.mtrx[,var.id]))
  highestbar = max(hist(bt.coef.mtrx[,var.id], plot = FALSE)$density) 
  ylimit <- max(c(y1.1,highestbar))
  hist(bt.coef.mtrx[,var.id], probability = TRUE, main = var.nm, xlab="", 
       col = "azure1",ylim=c(0,ylimit), border="lightseagreen")
  lines(x = x1.1, y = y1.1, col = "red3")
  lines(density(bt.coef.mtrx[,var.id], adjust=2), col="blue")}

par(mfrow=c(2,3))  # histograms of bootstrap coefs
boot.hist(bt.coef.mtrx=coef.mtrx, var.id=1, var.nm ="Intercept" )
boot.hist(bt.coef.mtrx=coef.mtrx, var.id=2, var.nm ="Log Age" )
boot.hist(bt.coef.mtrx=coef.mtrx, var.id=3, var.nm ="Hypertension" )
boot.hist(bt.coef.mtrx=coef.mtrx, var.id=4, var.nm ="Heart Disease" )
boot.hist(bt.coef.mtrx=coef.mtrx, var.id=5, var.nm ="Diabetes" )

The histograms were created to depict to normal density curves on top of each of the histograms. The first curve is red and it uses the estimated regression coefficients and their standard error in the output of the regression procedure. The red curve is the basis for the reported p-values. The second curve on the histograms is blue. The blue curve is a non-parametric estimate of the density of the bootstrap sampling distribution. The blue curve is the basis for the reported confidence intervals.

All of the density curves on the histograms seem to be quite close, with some of them even sharing a substantial amount of overlap. None of the curves appear to show any skewness, they are all normally distributed.

3.1.2 Confidence Intervals

num.p = dim(coef.mtrx)[2]
boot.ci = NULL
boot.wd = NULL
for (i in 1:num.p){
  lowerci0.025 = round(quantile(coef.mtrx[, i], 0.025, type = 2), 8)
  upperci0.975 = round(quantile(coef.mtrx [, i], 0.975, type =2), 8)
  boot.wd[i] = upperci0.975 - lowerci0.025
  boot.ci[i] = paste ("[", round(lowerci0.025, 4), " , ", round(upperci0.975,4), "]")}

kable(as.data.frame(cbind(formatC(bmi.coef.4, 4, format = "f"), boot.ci.95=boot.ci)), caption = "Regression Coefficient Matrix With 95% Residiual Bootstrap Confidence Interval")

Regression Coefficient Matrix With 95% Residiual Bootstrap Confidence Interval

	Estimate	Std. Error	t value	Pr(>|t|)	boot.ci.95
(Intercept)	0.2424	0.0002	1000.9204	0.0000	[ 0.2418 , 0.243 ]
log(age)	-0.0133	0.0001	-192.8579	0.0000	[ -0.0135 , -0.0132 ]
hypertension	-0.0028	0.0002	-11.2622	0.0000	[ -0.0033 , -0.0023 ]
heart_disease	0.0040	0.0003	12.2217	0.0000	[ 0.0034 , 0.0046 ]
diabetes	-0.0078	0.0002	-33.5534	0.0000	[ -0.0083 , -0.0074 ]

Looking at the confidence intervals for the variables, we want to confirm that none of them contain the value 0. From the output, it seems that none of the 95% confidence intervals span 0 and they are all consistent. Additionally, all the p-values are smaller than 0.05. This output aligns with what is depicted in the normal density curves on the histograms above.

3.2 Bootstrapping Residuals

Now we are going use the bootstrap residual methods to create bootstrap confidence intervals.

3.2.1 Residuals

hist(sort(bmi.log.age.4$residuals),
     xlab = "Residuals",
     main = "Histogram of Residuals")

The residual plot shows that the distribution of the residuals has no outlines and appears to have a normal distribution.

Now we are going to use the below code to generate the bootstrap confidence intervals of the regression coefficients

final.bmi.model = lm((bmi)^-0.5 ~ log(age) + hypertension + heart_disease
                 + diabetes, data = final.data)

model.residual = final.bmi.model$residuals
num.p = dim(model.matrix(final.bmi.model))[2]
sample.num = dim(model.matrix(final.bmi.model))[1]
boot.matrix = matrix(rep(0,5*B), ncol = num.p)

for (i in 1:B) {
  boot.log.bmi = final.bmi.model$fitted.values +
  sample(final.bmi.model$residuals, sample.num, replace = TRUE)
  final.data$boot.log.bmi = boot.log.bmi
  bootr.model = lm(boot.log.bmi ~ log(age) + hypertension + heart_disease
                 + diabetes, data = final.data)   
  boot.matrix[i, ] = bootr.model$coefficients}

3.2.2 Histograms

Now using the function that was mentioned earlier we are going to make the histograms of the residual bootstrap estimates of the regression coefficients.

boot.res.hist = function(bt.res.mtrx, var.id, var.nm){
  x1.1 <- seq(min(bt.res.mtrx[,var.id]), max(bt.res.mtrx[,var.id]), length=300 )
  y1.1 <- dnorm(x1.1, mean(bt.res.mtrx[,var.id]), sd(bt.res.mtrx[,var.id]))
  highestbar = max(hist(bt.res.mtrx[,var.id], plot = FALSE)$density) 
  ylimit <- max(c(y1.1,highestbar))
  hist(bt.res.mtrx[,var.id], probability = TRUE, main = var.nm, xlab="", 
       col = "azure1",ylim=c(0,ylimit), border="lightseagreen")
  lines(x = x1.1, y = y1.1, col = "red3")
  lines(density(bt.res.mtrx[,var.id], adjust=2), col="blue")}

par(mfrow=c(2,3))  # histograms of bootstrap residuals
boot.res.hist(bt.res.mtrx=boot.matrix, var.id=1, var.nm ="Intercept" )
boot.res.hist(bt.res.mtrx=boot.matrix, var.id=2, var.nm ="Log Age" )
boot.res.hist(bt.res.mtrx=boot.matrix, var.id=3, var.nm ="Hypertension" )
boot.res.hist(bt.res.mtrx=boot.matrix, var.id=4, var.nm ="Heart Disease" )
boot.res.hist(bt.res.mtrx=boot.matrix, var.id=5, var.nm ="Diabetes" )

The curves on all the histograms seem to be very close and have good overlap similar to the earlier histograms. There does seems to be some slight skewness from the center for both the Intercept and Log Age histograms, however the skewness seems to be very slight and nothing to be extremely concerned with. The confidence intervals should yield similar information to the histograms so we can confirm that the skewness is minor when analyzing the confidence intervals. The 95% residual bootstrap confidence intervals are given below,

3.2.3 Confidence Intervals

num.p = dim(boot.matrix)[2]
boot.r.ci = NULL
boot.r.wd = NULL
for (i in 1:num.p){
  lowerci0.025.res = round(quantile(boot.matrix[, i], 0.025, type = 2), 8)
  upperci0.975.res = round(quantile(boot.matrix [, i], 0.975, type =2), 8)
  boot.r.wd[i] = upperci0.975.res - lowerci0.025.res
  boot.r.ci[i] = paste ("[", round(lowerci0.025.res, 4), " , ", round(upperci0.975.res,4), "]")}

kable(as.data.frame(cbind(formatC(bmi.coef.4, 4, format = "f"), boot.r.ci.95=boot.r.ci)), caption = "Regression Coefficient Matric")

Regression Coefficient Matric

	Estimate	Std. Error	t value	Pr(>|t|)	boot.r.ci.95
(Intercept)	0.2424	0.0002	1000.9204	0.0000	[ 0.2419 , 0.2429 ]
log(age)	-0.0133	0.0001	-192.8579	0.0000	[ -0.0135 , -0.0132 ]
hypertension	-0.0028	0.0002	-11.2622	0.0000	[ -0.0033 , -0.0023 ]
heart_disease	0.0040	0.0003	12.2217	0.0000	[ 0.0033 , 0.0047 ]
diabetes	-0.0078	0.0002	-33.5534	0.0000	[ -0.0083 , -0.0074 ]

The residual bootstrap confidence intervals yield the same results as the p-values. This is to be expected considering the significantly large sample size. The large sample size allows for the distribution of the estimated coefficients to have a good approximation of the normal distributions.

4 Results and Conclusions

4.1 Inferential Statistics

kable(as.data.frame(cbind(formatC(bmi.coef.4[, -3], 4, format ="f"), btc.ct.95 = boot.ci, btr.ci.95 = boot.r.ci)),
      caption = "Final Combined Inferential Statistics with p-values and Bootstrap Confidence Intervals")

Final Combined Inferential Statistics with p-values and Bootstrap Confidence Intervals

	Estimate	Std. Error	Pr(>|t|)	btc.ct.95	btr.ci.95
(Intercept)	0.2424	0.0002	0.0000	[ 0.2418 , 0.243 ]	[ 0.2419 , 0.2429 ]
log(age)	-0.0133	0.0001	0.0000	[ -0.0135 , -0.0132 ]	[ -0.0135 , -0.0132 ]
hypertension	-0.0028	0.0002	0.0000	[ -0.0033 , -0.0023 ]	[ -0.0033 , -0.0023 ]
heart_disease	0.0040	0.0003	0.0000	[ 0.0034 , 0.0046 ]	[ 0.0033 , 0.0047 ]
diabetes	-0.0078	0.0002	0.0000	[ -0.0083 , -0.0074 ]	[ -0.0083 , -0.0074 ]

Looking at the above output it is clear that the three methods yield the same results, which is the significance of the individual explanatory variables. This makes sense considering that the final model \[ \text{bmi} = 0.2424251 - 0.0133499\times \log(\text{age}) - 0.0027673\times \text{hypertension} + 0.0040239\times \text{heart_disease} - 0.0078379\times \text{diabetes} \] does not have serious violations of the model assumptions.

kable(round(cbind(boot.wd, boot.r.wd), 4), caption = "Widths of both Bootstrap Confidence Intervals")

Widths of both Bootstrap Confidence Intervals

boot.wd	boot.r.wd
0.0012	0.0010
0.0003	0.0003
0.0010	0.0009
0.0012	0.0014
0.0009	0.0009

The above table shows that the widths of the residual and case bootstrap confidence intervals are either very similar or exactly the same as each other.

kable(bmi.coef.4, caption = "Inferential Statistics of Final Model")

Inferential Statistics of Final Model

	Estimate	Std. Error	t value	Pr(>|t|)
(Intercept)	0.2424251	0.0002422	1000.92041	0
log(age)	-0.0133499	0.0000692	-192.85792	0
hypertension	-0.0027673	0.0002457	-11.26224	0
heart_disease	0.0040239	0.0003292	12.22173	0
diabetes	-0.0078379	0.0002336	-33.55337	0

The regression coefficients explain the correlation between the adjusted bmi and the corresponding variables, with log(age) instead of age. Taking into consideration that both the confidence intervals and the p-values depicted that these variables had a significant ability to determine the response variable that was not 0. We can determine that this model can be used to help predict a patients expected bmi based on their other factors like hypertension, and age.

Additionally, we can explain the estimated regression coefficient of -0.0078379 as follows, since diabetes is a binary variable it can only be expressed as 0 or 1. Due to the transformation on BMI during model building the variable diabetes would be expressed as follows,

\[ Diabetes_1^{-0.5} - Diabetes_0^{-0.5} = -0.0078379 \] Since, 1 raised to any power will always be 1, and 0 raised to -0.5 is undefined it simplifies to the following.

\[ Diabetes_1 = -0.0078379 \] Using this model, patients with diabetes with see a 0.78% decrease in the bmi of the patient. The other binary variables in the model can be interpreted similarly.

5 General Discussion

While the model depicted several significant predictive variables for the response variable in the context of this data set it is important to note that while the variables are significant their impact on the value of bmi itself is rather small. Unfortunately there is the misconception when it comes to the medical field that everyone who has issues like diabetes, hypertension, or heart disease must be overweight and all their problems come from being overweight. However, this model depicts only a small change to the predicted bmi if a patient has diabetes, hypertension, or heart disease or not. While the American Diabetes Association does note that having a higher bmi does put patients at risk for having diabetes, anyone can become diabetic, regardless of their weight. Any misconceptions on a exclusive relationship between bmi and diabetic status are simply that, misconceptions.

Also, it is important to note that this data set is incredibly large with over 100,000 observations. Large data sets make it very easy to get p-values that are close to if not zero. Ideally, it would be better to use a smaller data set when analyzing the variables for the final model to see what the p-values look like when they are not pushed to zero by a large sample size.

However, the model was not determined solely on the p-value, and there was additional testing that took place. Additionally, while ethically using factors like hypertension, heart disease and diabetes should not be the main aspects of trying to determine a patients bmi. Due to the fact that is would most likely further the connotation that only people with those ailments have high bmi’s, the model could technically be used in that way.

6 References

American Diabetes Association. (2023). Understanding Diabetes Diagnosis. Diabetes.org. https://diabetes.org/about-diabetes/diagnosis

American Diabetes Association. Extra Weight, Extra Risk Diabetes.org. https://diabetes.org/health-wellness/weight-management/extra-weight-extra-risk