Multiple Regression
Prasanna Sadashivaiah
Overview:
The dataset used for this study is the National Health and Nutrition Examination Survey (NHANES) adult sample. The dataset contains approximately 500 cases (participants) and over 70 variables that describe demographic characteristics, health outcomes, lifestyle behaviors, and physical measurements. Each row represents an individual adult participant, and variables include age, gender, income, physical activity, blood pressure, diabetes status, and Body Mass Index (BMI).
You’ll load it from a URL, clean it, fit simple and multiple linear regression models, interpret the outputs, check assumptions, and diagnose model fit using residuals, plots, and metrics like R² and RMSE.
Step 1: Examine Dataset
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.1 ✔ readr 2.2.0
## ✔ ggplot2 4.0.3 ✔ stringr 1.6.0
## ✔ lubridate 1.9.5 ✔ tibble 3.3.1
## ✔ purrr 1.2.2 ✔ tidyr 1.3.2## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsnhanes <- read.csv("C:/PMS/Data_Science/Data101/Final_Project/nhanes.samp.adult.500.csv",
stringsAsFactors = FALSE)
dim(nhanes)## [1] 500 76str(nhanes)## 'data.frame': 500 obs. of 76 variables:
## $ ID : int 63106 67820 57178 68693 69465 61505 67341 69498 57323 56025 ...
## $ SurveyYr : chr "2011_12" "2011_12" "2009_10" "2011_12" ...
## $ Gender : chr "male" "female" "male" "male" ...
## $ Age : int 50 47 46 28 50 39 74 31 21 80 ...
## $ AgeDecade : chr " 50-59" " 40-49" " 40-49" " 20-29" ...
## $ AgeMonths : int NA NA 561 NA NA 471 NA NA 260 NA ...
## $ Race1 : chr "White" "Black" "White" "White" ...
## $ Race3 : chr "White" "Black" NA "White" ...
## $ Education : chr "9 - 11th Grade" "College Grad" "Some College" "Some College" ...
## $ MaritalStatus : chr "Divorced" "Separated" "Married" "NeverMarried" ...
## $ HHIncome : chr "10000-14999" "35000-44999" "more 99999" "more 99999" ...
## $ HHIncomeMid : int 12500 40000 100000 100000 40000 NA 50000 100000 50000 40000 ...
## $ Poverty : num 0.95 1.74 4.99 4.14 2.16 NA 3.3 5 0.92 3.51 ...
## $ HomeRooms : int 7 6 6 8 10 4 4 4 4 6 ...
## $ HomeOwn : chr "Own" "Rent" "Own" "Own" ...
## $ Work : chr "NotWorking" "Working" "Working" "Working" ...
## $ Weight : num 82.8 79.9 73.7 80.9 70.5 ...
## $ Length : logi NA NA NA NA NA NA ...
## $ HeadCirc : logi NA NA NA NA NA NA ...
## $ Height : num 172 165 170 177 162 ...
## $ BMI : num 27.9 29.4 25.4 25.8 26.9 ...
## $ BMICatUnder20yrs: logi NA NA NA NA NA NA ...
## $ BMI_WHO : chr "25.0_to_29.9" "25.0_to_29.9" "25.0_to_29.9" "25.0_to_29.9" ...
## $ Pulse : int 58 70 74 58 76 58 68 60 84 82 ...
## $ BPSysAve : int 125 121 120 132 152 148 122 114 108 168 ...
## $ BPDiaAve : int 86 68 74 74 103 88 65 76 46 52 ...
## $ BPSys1 : int 122 124 120 134 144 150 120 114 108 170 ...
## $ BPDia1 : int 88 66 70 72 104 94 68 82 44 76 ...
## $ BPSys2 : int 124 120 118 130 150 148 118 114 108 168 ...
## $ BPDia2 : int 86 66 74 72 106 92 62 76 46 52 ...
## $ BPSys3 : int 126 122 122 134 154 148 126 114 108 168 ...
## $ BPDia3 : int 86 70 74 76 100 84 68 76 46 52 ...
## $ Testosterone : num 525.37 5.98 NA 653.19 8.17 ...
## $ DirectChol : num 1.29 1.22 1.4 1.84 2.43 NA 1.14 1.29 1.34 1.6 ...
## $ TotChol : num 5.07 3.7 6.03 4.55 5.92 NA 3.21 5.02 6.15 5.02 ...
## $ UrineVol1 : int 244 65 105 51 30 NA 39 120 81 47 ...
## $ UrineFlow1 : num 1.683 0.442 0.682 0.464 1.304 ...
## $ UrineVol2 : int NA NA NA NA 114 NA NA NA NA NA ...
## $ UrineFlow2 : num NA NA NA NA 1.12 ...
## $ Diabetes : chr "No" "No" "No" "No" ...
## $ DiabetesAge : int NA NA NA NA NA NA NA NA NA NA ...
## $ HealthGen : chr "Fair" "Vgood" "Vgood" "Excellent" ...
## $ DaysPhysHlthBad : int 5 2 0 0 0 NA 30 0 0 0 ...
## $ DaysMentHlthBad : int 30 5 0 0 0 NA 0 10 30 0 ...
## $ LittleInterest : chr "Most" "None" "None" "None" ...
## $ Depressed : chr "Several" "Several" "None" "None" ...
## $ nPregnancies : int NA 2 NA NA 3 NA NA NA NA 1 ...
## $ nBabies : int NA 2 NA NA 3 NA NA NA NA 1 ...
## $ Age1stBaby : int NA 21 NA NA 27 NA NA NA NA NA ...
## $ SleepHrsNight : int 4 6 5 7 6 6 8 7 5 7 ...
## $ SleepTrouble : chr "Yes" "No" "No" "No" ...
## $ PhysActive : chr "No" "Yes" "Yes" "Yes" ...
## $ PhysActiveDays : int 3 NA 3 NA NA 4 3 NA NA NA ...
## $ TVHrsDay : chr "2_hr" "3_hr" NA "0_to_1_hr" ...
## $ CompHrsDay : chr "0_hrs" "2_hr" NA "4_hr" ...
## $ TVHrsDayChild : logi NA NA NA NA NA NA ...
## $ CompHrsDayChild : logi NA NA NA NA NA NA ...
## $ Alcohol12PlusYr : chr "Yes" "Yes" "Yes" "No" ...
## $ AlcoholDay : int 1 2 NA NA 1 NA 1 2 2 1 ...
## $ AlcoholYear : int 24 3 0 0 364 NA 12 104 2 1 ...
## $ SmokeNow : chr "No" NA NA NA ...
## $ Smoke100 : chr "Yes" "No" "No" "No" ...
## $ Smoke100n : chr "Smoker" "Non-Smoker" "Non-Smoker" "Non-Smoker" ...
## $ SmokeAge : int 18 NA NA NA NA NA NA NA NA 21 ...
## $ Marijuana : chr "No" "Yes" "Yes" "No" ...
## $ AgeFirstMarij : int NA 19 14 NA NA NA NA 20 19 NA ...
## $ RegularMarij : chr "No" "Yes" "Yes" "No" ...
## $ AgeRegMarij : int NA 20 16 NA NA NA NA NA NA NA ...
## $ HardDrugs : chr "No" "No" "Yes" "No" ...
## $ SexEver : chr "Yes" "Yes" "Yes" "No" ...
## $ SexAge : int 16 17 14 NA 17 NA NA 19 13 NA ...
## $ SexNumPartnLife : int 26 10 50 0 4 NA NA 3 15 NA ...
## $ SexNumPartYear : int 2 2 1 0 1 NA NA 1 2 NA ...
## $ SameSex : chr "No" "No" "No" "No" ...
## $ SexOrientation : chr "Heterosexual" "Heterosexual" "Heterosexual" "Heterosexual" ...
## $ PregnantNow : chr NA NA NA NA ...summary(nhanes$BMI)## Min. 1st Qu. Median Mean 3rd Qu. Max. NAs
## 17.80 24.67 28.20 29.39 32.93 68.63 3Step 2: Select Relevant Variables
nhanes_clean <- nhanes %>%
select(BMI, Age, HHIncomeMid, PhysActive)#Step 3: Remove Missing Values
nhanes_clean <- nhanes_clean %>%
filter(!is.na(BMI),
!is.na(Age),
!is.na(HHIncomeMid),
!is.na(PhysActive))Step 4: Create New Variable
nhanes_clean <- nhanes_clean %>%
mutate(
Active = ifelse(PhysActive == "Yes", 1, 0)
)Step 5: Run Multiple Regression
model <- lm(BMI ~ Age + HHIncomeMid + Active,
data = nhanes_clean)
summary(model)##
## Call:
## lm(formula = BMI ~ Age + HHIncomeMid + Active, data = nhanes_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.838 -4.662 -1.115 3.829 37.517
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.128e+01 1.214e+00 25.772 < 2e-16 ***
## Age 4.620e-03 1.925e-02 0.240 0.81042
## HHIncomeMid -1.845e-05 1.006e-05 -1.834 0.06724 .
## Active -2.051e+00 6.594e-01 -3.110 0.00199 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.813 on 451 degrees of freedom
## Multiple R-squared: 0.03536, Adjusted R-squared: 0.02894
## F-statistic: 5.511 on 3 and 451 DF, p-value: 0.001006#Step 6: Select Variables and Clean Data
nhanes_clean <- nhanes %>%
select(BMI, Age, HHIncomeMid, PhysActive) %>%
filter(
!is.na(BMI),
!is.na(Age),
!is.na(HHIncomeMid),
!is.na(PhysActive)
)#Step 7: Create New Variable
nhanes_clean <- nhanes_clean %>%
mutate(
Active = ifelse(PhysActive == "Yes", 1, 0)
)Step 8: Run Multiple Regression Model
model <- lm(BMI ~ Age + HHIncomeMid + Active,
data = nhanes_clean)
summary(model)##
## Call:
## lm(formula = BMI ~ Age + HHIncomeMid + Active, data = nhanes_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.838 -4.662 -1.115 3.829 37.517
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.128e+01 1.214e+00 25.772 < 2e-16 ***
## Age 4.620e-03 1.925e-02 0.240 0.81042
## HHIncomeMid -1.845e-05 1.006e-05 -1.834 0.06724 .
## Active -2.051e+00 6.594e-01 -3.110 0.00199 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.813 on 451 degrees of freedom
## Multiple R-squared: 0.03536, Adjusted R-squared: 0.02894
## F-statistic: 5.511 on 3 and 451 DF, p-value: 0.001006Step 8: table of coefficients
summary(model)##
## Call:
## lm(formula = BMI ~ Age + HHIncomeMid + Active, data = nhanes_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.838 -4.662 -1.115 3.829 37.517
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.128e+01 1.214e+00 25.772 < 2e-16 ***
## Age 4.620e-03 1.925e-02 0.240 0.81042
## HHIncomeMid -1.845e-05 1.006e-05 -1.834 0.06724 .
## Active -2.051e+00 6.594e-01 -3.110 0.00199 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.813 on 451 degrees of freedom
## Multiple R-squared: 0.03536, Adjusted R-squared: 0.02894
## F-statistic: 5.511 on 3 and 451 DF, p-value: 0.001006Step 9: Check Regression Assumption
plot(model, which = 1)
Step 10: Normal Q-Q Plot (Normality)
plot(model, which = 2)
Scale-Location Plot (Homoscedasticity)
plot(model, which = 3)
Residuals vs Leverage Plot
plot(model, which = 5)
Multicollinearity
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:purrr':
##
## some## The following object is masked from 'package:dplyr':
##
## recodemodel <- lm(BMI ~ Age + HHIncomeMid + Active,
data = nhanes_clean)
vif(model)## Age HHIncomeMid Active
## 1.022579 1.042409 1.065456