Health Data Science Project

Factors influencing obesity among adults in the US

Problem definition

In this project we are trying to identify the various determinants contributing to the prevalence of obesity among adults in the United States, including social, economic, cultural, environmental, and behavioral factors

Abstract and Introduction

We will use the BRFSS (Behavioral Risk Factor Surveillance System) dataset. It is a cross-sectional telephone survey that state health departments conduct monthly over landline telephones and cellular telephones with a standardised questionnaire and technical and methodologic assistance from the CDC (Centers for Disease Control and Prevention). BRFSS is used to collect prevalence data among U.S. residents regarding their risk behaviours and preventive health practices that can affect their health status. We decided on the topic of the factors that influence obesity among adults in the U.S. Studying the factors influencing obesity among adults in the U.S. is crucial. Obesity is a significant public health issue in the U.S., contributing to various chronic diseases like heart disease, diabetes, and certain cancers. Understanding the factors contributing to obesity can inform public health policies and interventions aimed at prevention and management. Obesity-related healthcare costs are substantial, impacting both individuals and the healthcare system as a whole.

WHY IS IT IMPORTANT TO KNOW ABOUT OBESITY? By identifying the factors influencing obesity, healthcare providers can develop more effective strategies for prevention and treatment, potentially reducing healthcare expenditures. Obesity can have profound social and economic consequences, including reduced quality of life, discrimination, and decreased productivity. Exploring the factors behind obesity can help address underlying societal issues such as access to healthy foods, socioeconomic disparities, and cultural influences.

Researching the factors influencing obesity fosters innovation in healthcare and public health. It encourages the development of new technologies, interventions, and policies aimed at combating obesity and promoting healthier lifestyles. Overall, investigating the factors influencing obesity among adults in the U.S. is essential for improving public health outcomes, reducing healthcare costs, addressing social and economic disparities, promoting individual well-being, and driving innovation in healthcare and public health practices.This is a fundamental public health question and can provide insights into preventive strategies and policymaking.

How do we measure obesity in simple terms? Obesity is typically defined by the Body Mass Index (BMI) with a threshold value (e.g., BMI ≥ 30 as obese). Here’s a detailed breakdown of how you might approach analyzing this question using different statistical models, specifying potential variables and offering some analytical tips.

Uploading and preparing the data set

We will upload all the necessary libraries we will use in our analysis here.

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, union

library(readxl)
library(gplots)

## 
## Attaching package: 'gplots'

## The following object is masked from 'package:stats':
## 
##     lowess

library(ggplot2)
library(maps)
library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

library(caret)

## Loading required package: lattice

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

library(boot)

## 
## Attaching package: 'boot'

## The following object is masked from 'package:lattice':
## 
##     melanoma

## The following object is masked from 'package:car':
## 
##     logit

library(geepack)
library(tidyr)
library(nnet)
library(mlogit)

## Loading required package: dfidx

## 
## Attaching package: 'dfidx'

## The following object is masked from 'package:MASS':
## 
##     select

## The following object is masked from 'package:stats':
## 
##     filter

library(RColorBrewer)
library(gridExtra)

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

library(knitr)
library(kableExtra)

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

# Specify the path to your csv file
file_path <- "Demo+Health data.xlsx"

Loading the date

# Read the Excel file
data <- read_excel(file_path)
head(data)

## # A tibble: 6 × 15
##   `_state` genhlth poorhlth cvdcrhd4 diabete4 marital educa employ1 income3
##      <dbl>   <dbl>    <dbl>    <dbl>    <dbl>   <dbl> <dbl>   <dbl>   <dbl>
## 1        1       2       NA        2        1       1     6       7      99
## 2        1       1       NA        2        3       3     4       2       5
## 3        1       2        2        2        3       1     6       7      10
## 4        1       1       NA        2        3       1     4       7      77
## 5        1       4       88        2        3       1     5       5       5
## 6        1       5       88        2        1       1     4       7      99
## # ℹ 6 more variables: smokday2 <dbl>, `_totinda` <dbl>, `_sex` <dbl>,
## #   `_ageg5yr` <dbl>, `_bmi5` <dbl>, `_bmi5cat` <dbl>

dim(data)

## [1] 445132     15

We see that we have 445 132 observations, which is good to work on a big data project. So let us carry on. HOPE YOU WILL ENJOY OUR RIDE :) !

Data Cleaning

# Get the column names
column_names <- names(data)

# Get the data types of each column
data_types <- sapply(data, class)

# Get the number of non-null values in each column
non_null_counts <- sapply(data, function(x) sum(!is.na(x)))

# Create a summary data frame
summary_data <- data.frame(column_names, data_types, non_null_counts)

# Rename the columns of the summary data frame
colnames(summary_data) <- c("Column Name", "Data Type", "Non-Null Count")

# Print the summary data frame
print(summary_data)

##          Column Name Data Type Non-Null Count
## _state        _state   numeric         445132
## genhlth      genhlth   numeric         445129
## poorhlth    poorhlth   numeric         255746
## cvdcrhd4    cvdcrhd4   numeric         445130
## diabete4    diabete4   numeric         445129
## marital      marital   numeric         445124
## educa          educa   numeric         445127
## employ1      employ1   numeric         438936
## income3      income3   numeric         432200
## smokday2    smokday2   numeric         164053
## _totinda    _totinda   numeric         445132
## _sex            _sex   numeric         445132
## _ageg5yr    _ageg5yr   numeric         445132
## _bmi5          _bmi5   numeric         396326
## _bmi5cat    _bmi5cat   numeric         396326

#we see that we have the same number of observations
#non null data, so we do not have missing values.

# Get the number of missing values in each column
missing_counts <- colSums(is.na(data))

# Sort the missing value counts in descending order
missing_counts <- sort(missing_counts, decreasing = TRUE)

# Filter out columns with no missing values
missing_counts <- missing_counts[missing_counts > 0]

# Print the result
print(missing_counts)

## smokday2 poorhlth    _bmi5 _bmi5cat  income3  employ1  marital    educa 
##   281079   189386    48806    48806    12932     6196        8        5 
##  genhlth diabete4 cvdcrhd4 
##        3        3        2

# Create a new column for index
data$index <- seq_len(nrow(data))

# Print the dataframe with the new index
print(data)

## # A tibble: 445,132 × 16
##    `_state` genhlth poorhlth cvdcrhd4 diabete4 marital educa employ1 income3
##       <dbl>   <dbl>    <dbl>    <dbl>    <dbl>   <dbl> <dbl>   <dbl>   <dbl>
##  1        1       2       NA        2        1       1     6       7      99
##  2        1       1       NA        2        3       3     4       2       5
##  3        1       2        2        2        3       1     6       7      10
##  4        1       1       NA        2        3       1     4       7      77
##  5        1       4       88        2        3       1     5       5       5
##  6        1       5       88        2        1       1     4       7      99
##  7        1       2       NA        2        3       2     6       7       8
##  8        1       3       NA        2        3       3     4       7       7
##  9        1       3       NA        2        3       1     6       1       7
## 10        1       3        2        2        1       1     6       7       7
## # ℹ 445,122 more rows
## # ℹ 7 more variables: smokday2 <dbl>, `_totinda` <dbl>, `_sex` <dbl>,
## #   `_ageg5yr` <dbl>, `_bmi5` <dbl>, `_bmi5cat` <dbl>, index <int>

# Remove rows with null values
data <- na.omit(data)

# Print the dataframe after removing null values
print(data)

## # A tibble: 94,694 × 16
##    `_state` genhlth poorhlth cvdcrhd4 diabete4 marital educa employ1 income3
##       <dbl>   <dbl>    <dbl>    <dbl>    <dbl>   <dbl> <dbl>   <dbl>   <dbl>
##  1        1       4       15        2        3       1     4       8       7
##  2        1       4       88        2        1       3     3       7       6
##  3        1       4       88        1        1       2     4       8       2
##  4        1       3       88        1        3       1     3       4      77
##  5        1       5       14        1        1       2     5       8      77
##  6        1       5       30        2        1       1     6       5      99
##  7        1       3       15        2        3       3     5       7       2
##  8        1       3       88        2        3       1     4       7       7
##  9        1       4       88        2        1       1     4       7       1
## 10        1       4        4        2        1       1     4       7      99
## # ℹ 94,684 more rows
## # ℹ 7 more variables: smokday2 <dbl>, `_totinda` <dbl>, `_sex` <dbl>,
## #   `_ageg5yr` <dbl>, `_bmi5` <dbl>, `_bmi5cat` <dbl>, index <int>

#New dimension after deleting the missing values
dim(data)

## [1] 94694    16

# Re_index
data$index <- paste0("row", seq_len(nrow(data)))

We will remove the values for the questions when the respondents let it blank and they are coded under a number, in this case the omi function can not delete them so we had to do it manually for each variabel following the dataset’s codebook.

# Define the columns and values to remove
columns_to_remove <- c('genhlth', 'poorhlth', 'cvdcrhd4', 'diabete4', 'marital', 
                       'educa', 'employ1', 'income3', 'smokday2', '_totinda', '_sex', '_ageg5yr')

#List with values to remove
values_to_remove <- list(
  genhlth = c(9, 7),  
  poorhlth = c(88, 77, 99),  
  cvdcrhd4 = c(7, 9),        
  diabete4 = c(7, 9),        
  marital = c(9),            
  educa = c(9),              
  employ1 = c(9),            
  income3 = c(77, 99),       
  smokday2 = c(7, 9),       
  "_totinda" = c(9),        
  "_ageg5yr"= c(14)        
)

# Function to filter out rows with specific values in each column
filter_values <- function(data, column_name, values) {
  data %>%
    filter(!(!!rlang::sym(column_name) %in% values))
}

# Apply the function to each column
filtered_data <- data
for (column_name in columns_to_remove) {
  filtered_data <- filter_values(filtered_data, column_name, values_to_remove[[column_name]])
}

# Print the first few rows of the filtered data frame
print(head(filtered_data))

## # A tibble: 6 × 16
##   `_state` genhlth poorhlth cvdcrhd4 diabete4 marital educa employ1 income3
##      <dbl>   <dbl>    <dbl>    <dbl>    <dbl>   <dbl> <dbl>   <dbl>   <dbl>
## 1        1       4       15        2        3       1     4       8       7
## 2        1       3       15        2        3       3     5       7       2
## 3        1       4       15        1        3       2     4       7       2
## 4        1       3        1        2        1       1     4       7       3
## 5        1       3        1        2        3       3     4       7       2
## 6        1       4        1        1        3       3     6       7       3
## # ℹ 7 more variables: smokday2 <dbl>, `_totinda` <dbl>, `_sex` <dbl>,
## #   `_ageg5yr` <dbl>, `_bmi5` <dbl>, `_bmi5cat` <dbl>, index <chr>

#The dimension after deleting the missing values with a specific number
dim(filtered_data)

## [1] 39551    16

We will change the names of the columns, so it will be easier for us to understand and work with.

names(filtered_data)

##  [1] "_state"   "genhlth"  "poorhlth" "cvdcrhd4" "diabete4" "marital" 
##  [7] "educa"    "employ1"  "income3"  "smokday2" "_totinda" "_sex"    
## [13] "_ageg5yr" "_bmi5"    "_bmi5cat" "index"

# Change column names
names(filtered_data) <- c("state", "genhlth", "poorhlth", "hrtdisease", "diabetic", "marital", "edu",   "employstatus", "income", "smokeday", "totinda", "sex", "age5cat", "bmi5", "bmi", "index", "obesity_status")

## Warning: The `value` argument of `names<-()` must have the same length as `x` as of
## tibble 3.0.0.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Print new column names
print(names(filtered_data)) 

##  [1] "state"        "genhlth"      "poorhlth"     "hrtdisease"   "diabetic"    
##  [6] "marital"      "edu"          "employstatus" "income"       "smokeday"    
## [11] "totinda"      "sex"          "age5cat"      "bmi5"         "bmi"         
## [16] "index"

summary(filtered_data)

##      state          genhlth         poorhlth       hrtdisease   
##  Min.   : 1.00   Min.   :1.000   Min.   : 1.00   Min.   :1.000  
##  1st Qu.:18.00   1st Qu.:3.000   1st Qu.: 3.00   1st Qu.:2.000  
##  Median :29.00   Median :3.000   Median :10.00   Median :2.000  
##  Mean   :30.67   Mean   :3.277   Mean   :12.67   Mean   :1.885  
##  3rd Qu.:45.00   3rd Qu.:4.000   3rd Qu.:20.00   3rd Qu.:2.000  
##  Max.   :78.00   Max.   :5.000   Max.   :30.00   Max.   :2.000  
##     diabetic        marital           edu         employstatus  
##  Min.   :1.000   Min.   :1.000   Min.   :1.000   Min.   :1.000  
##  1st Qu.:3.000   1st Qu.:1.000   1st Qu.:4.000   1st Qu.:1.000  
##  Median :3.000   Median :2.000   Median :5.000   Median :5.000  
##  Mean   :2.638   Mean   :2.535   Mean   :4.852   Mean   :4.417  
##  3rd Qu.:3.000   3rd Qu.:4.000   3rd Qu.:6.000   3rd Qu.:7.000  
##  Max.   :4.000   Max.   :6.000   Max.   :6.000   Max.   :8.000  
##      income         smokeday        totinda           sex       
##  Min.   : 1.00   Min.   :1.000   Min.   :1.000   Min.   :1.000  
##  1st Qu.: 4.00   1st Qu.:1.000   1st Qu.:1.000   1st Qu.:1.000  
##  Median : 6.00   Median :3.000   Median :1.000   Median :2.000  
##  Mean   : 5.85   Mean   :2.371   Mean   :1.361   Mean   :1.533  
##  3rd Qu.: 8.00   3rd Qu.:3.000   3rd Qu.:2.000   3rd Qu.:2.000  
##  Max.   :11.00   Max.   :3.000   Max.   :2.000   Max.   :2.000  
##     age5cat            bmi5           bmi           index          
##  Min.   : 1.000   Min.   :1205   Min.   :1.000   Length:39551      
##  1st Qu.: 5.000   1st Qu.:2439   1st Qu.:2.000   Class :character  
##  Median : 8.000   Median :2831   Median :3.000   Mode  :character  
##  Mean   : 7.548   Mean   :2959   Mean   :3.097                     
##  3rd Qu.:10.000   3rd Qu.:3330   3rd Qu.:4.000                     
##  Max.   :13.000   Max.   :9566   Max.   :4.000

We will make less categories for the the specific columns and mapping the states as they were coded as numbers from 1 to 78, we mapped them to their corsponding names.

#Mapping the states
state_mapping <- c(
  `1` = "Alabama", `2` = "Alaska", `4` = "Arizona", `5` = "Arkansas", 
  `6` = "California", `8` = "Colorado", `9` = "Connecticut", `10` = "Delaware", 
  `11` = "District of Columbia", `12` = "Florida", `13` = "Georgia", `15` = "Hawaii", 
  `16` = "Idaho", `17` = "Illinois", `18` = "Indiana", `19` = "Iowa", 
  `20` = "Kansas", `21` = "Kentucky", `22` = "Louisiana", `23` = "Maine", 
  `24` = "Maryland", `25` = "Massachusetts", `26` = "Michigan", `27` = "Minnesota", 
  `28` = "Mississippi", `29` = "Missouri", `30` = "Montana", `31` = "Nebraska", 
  `32` = "Nevada", `33` = "New Hampshire", `34` = "New Jersey", `35` = "New Mexico", 
  `36` = "New York", `37` = "North Carolina", `38` = "North Dakota", `39` = "Ohio", 
  `40` = "Oklahoma", `41` = "Oregon", `42` = "Pennsylvania", `44` = "Rhode Island", 
  `45` = "South Carolina", `46` = "South Dakota", `47` = "Tennessee", `48` = "Texas", 
  `49` = "Utah", `50` = "Vermont", `51` = "Virginia", `53` = "Washington", 
  `54` = "West Virginia", `55` = "Wisconsin", `56` = "Wyoming", `66` = "Guam", 
  `72` = "Puerto Rico", `78` = "Virgin Islands"
)

edu_mapping <- c(
  "Never attended school or only kindergarten" = "Elementary or less",
  "Grades 1 through 8 (Elementary)" = "Elementary or less",
  "Grades 9 through 11 (Some high school)" = "High school or equivalent",
  "Grade 12 or GED (High school graduate)" = "High school or equivalent",
  "College 1 year to 3 years (Some college or technical school)" = "College or more",
  "College 4 years or more (College graduate)" = "College or more"
)

# Convert the edu column from numeric codes to descriptive names using the mapping
edu_levels <- c(
  "Never attended school or only kindergarten",
  "Grades 1 through 8 (Elementary)",
  "Grades 9 through 11 (Some high school)",
  "Grade 12 or GED (High school graduate)",
  "College 1 year to 3 years (Some college or technical school)",
  "College 4 years or more (College graduate)"
)

# Ensure age5cat is correctly labeled
age5cat_labels <- c(
  "Age 18 to 24", "Age 25 to 29", "Age 30 to 34", "Age 35 to 39", 
  "Age 40 to 44", "Age 45 to 49", "Age 50 to 54", "Age 55 to 59", 
  "Age 60 to 64", "Age 65 to 69", "Age 70 to 74", "Age 75 to 79", 
  "Age 80 or older"
)

# Convert age5cat to a factor with the correct levels and labels
filtered_data$age5cat <- factor(filtered_data$age5cat, levels = 1:13, labels = age5cat_labels)

# Create the age5cat mapping to fewer categories
age5cat_mapping <- c(
  "Age 18 to 24" = "Age 18 to 29",
  "Age 25 to 29" = "Age 18 to 29",
  "Age 30 to 34" = "Age 30 to 39",
  "Age 35 to 39" = "Age 30 to 39",
  "Age 40 to 44" = "Age 40 to 54",
  "Age 45 to 49" = "Age 40 to 54",
  "Age 50 to 54" = "Age 40 to 54",
  "Age 55 to 59" = "Age 55 to 69",
  "Age 60 to 64" = "Age 55 to 69",
  "Age 65 to 69" = "Age 55 to 69",
  "Age 70 to 74" = "Age 70 and older",
  "Age 75 to 79" = "Age 70 and older",
  "Age 80 or older" = "Age 70 and older"
)

filtered_data <- filtered_data %>% mutate(obesity_status = factor(ifelse(bmi == 4, 1, 0), levels = c(0, 1), labels = c("Not Obese", "Obese")))


# Assuming filtered_data is your dataframe and income is the current income factor
# Define new income categories
filtered_data$income_group <- cut(as.numeric(filtered_data$income), 
                                  breaks = c(0, 4, 7, 11), 
                                  labels = c("Low Income", "Middle Income", "High Income"),
                                  include.lowest = TRUE)

# Create a new variable 'diabetic_status' based on the existing 'diabetic' variable
filtered_data$diabetic_status <- ifelse(filtered_data$diabetic %in% c(1, 2), 1, 
                                        ifelse(filtered_data$diabetic %in% c(3, 4), 2, NA))

# Map the age5cat column to fewer categories
filtered_data$age5cat <- as.factor(age5cat_mapping[as.character(filtered_data$age5cat)])

# Ensure edu is a factor with the correct levels and labels
filtered_data$edu <- factor(filtered_data$edu, levels = 1:6, labels = edu_levels)

# Map the edu column to fewer categories
filtered_data$edu <- as.factor(edu_mapping[as.character(filtered_data$edu)])

# Map state codes to state names
filtered_data$state <- as.factor(state_mapping[as.character(filtered_data$state)])

filtered_data$sex <- factor(filtered_data$sex, levels = c(1, 2), labels = c("Male", "Female"))

# Define the marital status labels as a factor with appropriate levels
filtered_data$marital <- factor(filtered_data$marital, 
                                levels = 1:6, 
                                labels = c("Married", "Divorced", "Widowed", "Separated", "Never married", "Member of an unmarried couple"))

# Define the employment status labels as a factor with appropriate levels
filtered_data$employstatus <- factor(filtered_data$employstatus, 
                                     levels = 1:8, 
                                     labels = c("Employed for wages", "Self-employed", 
                                                "Out of work for 1 year or more", "Out of work for less than 1 year", 
                                                "A homemaker", "A student", "Retired", "Unable to work"))

# Define the heart disease status labels
filtered_data$hrtdisease <- factor(filtered_data$hrtdisease, levels = c(1, 2), labels = c("Yes", "No"))


# Display the structure of the filtered_data data frame to confirm changes
str(filtered_data)

## tibble [39,551 × 19] (S3: tbl_df/tbl/data.frame)
##  $ state          : Factor w/ 54 levels "Alabama","Alaska",..: 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "names")= chr [1:39551] "1" "1" "1" "1" ...
##  $ genhlth        : num [1:39551] 4 3 4 3 3 4 4 5 2 5 ...
##  $ poorhlth       : num [1:39551] 15 15 15 1 1 1 30 30 5 30 ...
##  $ hrtdisease     : Factor w/ 2 levels "Yes","No": 2 2 1 2 2 1 1 1 2 2 ...
##  $ diabetic       : num [1:39551] 3 3 3 1 3 3 1 1 3 3 ...
##  $ marital        : Factor w/ 6 levels "Married","Divorced",..: 1 3 2 1 3 3 1 1 1 1 ...
##  $ edu            : Factor w/ 3 levels "College or more",..: 3 1 3 3 3 1 3 1 3 1 ...
##   ..- attr(*, "names")= chr [1:39551] "Grade 12 or GED (High school graduate)" "College 1 year to 3 years (Some college or technical school)" "Grade 12 or GED (High school graduate)" "Grade 12 or GED (High school graduate)" ...
##  $ employstatus   : Factor w/ 8 levels "Employed for wages",..: 8 7 7 7 7 7 7 8 7 2 ...
##  $ income         : num [1:39551] 7 2 2 3 2 3 5 7 5 8 ...
##  $ smokeday       : num [1:39551] 3 1 1 3 3 1 3 1 3 3 ...
##  $ totinda        : num [1:39551] 2 2 2 1 1 1 2 2 2 2 ...
##  $ sex            : Factor w/ 2 levels "Male","Female": 1 2 2 1 2 1 1 2 2 1 ...
##  $ age5cat        : Factor w/ 5 levels "Age 18 to 29",..: 4 4 4 4 5 5 5 3 5 5 ...
##   ..- attr(*, "names")= chr [1:39551] "Age 55 to 59" "Age 60 to 64" "Age 65 to 69" "Age 65 to 69" ...
##  $ bmi5           : num [1:39551] 2859 4561 3112 3626 3219 ...
##  $ bmi            : num [1:39551] 3 4 4 4 4 2 4 4 3 3 ...
##  $ index          : chr [1:39551] "row1" "row7" "row11" "row17" ...
##  $ obesity_status : Factor w/ 2 levels "Not Obese","Obese": 1 2 2 2 2 1 2 2 1 1 ...
##  $ income_group   : Factor w/ 3 levels "Low Income","Middle Income",..: 2 1 1 1 1 1 2 2 2 3 ...
##  $ diabetic_status: num [1:39551] 2 2 2 1 2 2 1 1 2 2 ...
##  - attr(*, "na.action")= 'omit' Named int [1:350438] 1 2 3 4 5 6 7 8 9 10 ...
##   ..- attr(*, "names")= chr [1:350438] "1" "2" "3" "4" ...

We will convert the column to factors (mainly they are numeric).

# Convert specified columns to factors in the filtered_data data frame
filtered_data$genhlth <- as.factor(filtered_data$genhlth)
filtered_data$poorhlth <- as.factor(filtered_data$poorhlth)
filtered_data$hrtdisease <- as.factor(filtered_data$hrtdisease)
filtered_data$diabetic <- as.factor(filtered_data$diabetic)
filtered_data$marital <- as.factor(filtered_data$marital)
filtered_data$employstatus <- as.factor(filtered_data$employstatus)
filtered_data$income <- as.factor(filtered_data$income)
filtered_data$smokeday <- as.factor(filtered_data$smokeday)
filtered_data$totinda <- as.factor(filtered_data$totinda)
filtered_data$sex <- as.factor(filtered_data$sex)
filtered_data$bmi5 <- as.factor(filtered_data$bmi5)
filtered_data$obesity_status <- factor(filtered_data$obesity_status)
filtered_data$diabetic_status <- as.factor(filtered_data$diabetic_status)

# Display the structure of the filtered_data data frame to confirm changes
str(filtered_data)

## tibble [39,551 × 19] (S3: tbl_df/tbl/data.frame)
##  $ state          : Factor w/ 54 levels "Alabama","Alaska",..: 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "names")= chr [1:39551] "1" "1" "1" "1" ...
##  $ genhlth        : Factor w/ 5 levels "1","2","3","4",..: 4 3 4 3 3 4 4 5 2 5 ...
##  $ poorhlth       : Factor w/ 30 levels "1","2","3","4",..: 15 15 15 1 1 1 30 30 5 30 ...
##  $ hrtdisease     : Factor w/ 2 levels "Yes","No": 2 2 1 2 2 1 1 1 2 2 ...
##  $ diabetic       : Factor w/ 4 levels "1","2","3","4": 3 3 3 1 3 3 1 1 3 3 ...
##  $ marital        : Factor w/ 6 levels "Married","Divorced",..: 1 3 2 1 3 3 1 1 1 1 ...
##  $ edu            : Factor w/ 3 levels "College or more",..: 3 1 3 3 3 1 3 1 3 1 ...
##   ..- attr(*, "names")= chr [1:39551] "Grade 12 or GED (High school graduate)" "College 1 year to 3 years (Some college or technical school)" "Grade 12 or GED (High school graduate)" "Grade 12 or GED (High school graduate)" ...
##  $ employstatus   : Factor w/ 8 levels "Employed for wages",..: 8 7 7 7 7 7 7 8 7 2 ...
##  $ income         : Factor w/ 11 levels "1","2","3","4",..: 7 2 2 3 2 3 5 7 5 8 ...
##  $ smokeday       : Factor w/ 3 levels "1","2","3": 3 1 1 3 3 1 3 1 3 3 ...
##  $ totinda        : Factor w/ 2 levels "1","2": 2 2 2 1 1 1 2 2 2 2 ...
##  $ sex            : Factor w/ 2 levels "Male","Female": 1 2 2 1 2 1 1 2 2 1 ...
##  $ age5cat        : Factor w/ 5 levels "Age 18 to 29",..: 4 4 4 4 5 5 5 3 5 5 ...
##   ..- attr(*, "names")= chr [1:39551] "Age 55 to 59" "Age 60 to 64" "Age 65 to 69" "Age 65 to 69" ...
##  $ bmi5           : Factor w/ 2512 levels "1205","1216",..: 948 2019 1131 1490 1205 383 1291 1459 903 888 ...
##  $ bmi            : num [1:39551] 3 4 4 4 4 2 4 4 3 3 ...
##  $ index          : chr [1:39551] "row1" "row7" "row11" "row17" ...
##  $ obesity_status : Factor w/ 2 levels "Not Obese","Obese": 1 2 2 2 2 1 2 2 1 1 ...
##  $ income_group   : Factor w/ 3 levels "Low Income","Middle Income",..: 2 1 1 1 1 1 2 2 2 3 ...
##  $ diabetic_status: Factor w/ 2 levels "1","2": 2 2 2 1 2 2 1 1 2 2 ...
##  - attr(*, "na.action")= 'omit' Named int [1:350438] 1 2 3 4 5 6 7 8 9 10 ...
##   ..- attr(*, "names")= chr [1:350438] "1" "2" "3" "4" ...

EDA

First we wanted to see the distribution of our variables

# Plot bar chart for gender distribution
ggplot(filtered_data, aes(x = sex, fill = sex)) +
  geom_bar() +
  geom_text(stat = 'count', aes(label = ..count..), 
            vjust = -0.5, size = 5, color = "black", fontface = "bold") +  
  scale_fill_manual(values = c("Male" = "blue4", "Female" = "red4")) +  
  labs(title = "Gender Distribution in Our Dataset", x = "Gender", y = "Count") +
  theme_minimal() +  
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),  
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10),
    legend.position = "none"  # Remove the legend
  ) +
  ylim(0, max(table(filtered_data$sex)) * 1.1)

## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

The bar chart shows that the dataset consists of 18,479 males (46.7%) and 21,072 females (53.3%), indicating a higher representation of females.

# Plot the bar chart
ggplot(filtered_data, aes(x = marital)) +
  geom_bar(fill = "blue4") +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -0.5, size = 5, color = "black", fontface = "bold") +
  labs(title = "Distribution of Marital Status", x = "Marital Status", y = "Count") +
  ylim(0, max(table(filtered_data$marital)) * 1.1) +  
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold", margin = margin(b = 10)),  
    axis.title = element_text(size = 14),  
    axis.text.x = element_text(size = 12, angle = 45, hjust = 1, vjust = 1),  
    axis.text.y = element_text(size = 12),  
    panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "grey80"),  
    panel.grid.minor = element_blank(),  
    legend.position = "none"  
  )

## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

The bar chart shows the distribution of marital status in our dataset, we have 16,238 individuals are married (38.4%), 8,368 are divorced (19.8%), 4,142 are widowed (9.8%), 1,399 are separated (3.3%), 7,152 have never married (16.9%), and 2,252 are members of an unmarried couple (5.3%). The majority of individuals in the dataset are married, followed by those who are divorced and never married.

# Plot the bar chart
ggplot(filtered_data, aes(x = edu)) +
  geom_bar(fill = "blue4") +
  geom_text(stat = 'count', aes(label = ..count..), 
            vjust = -0.5, size = 5, color = "black", fontface = "bold") +
  labs(title = "Distribution of Education Status", x = "Education Level", y = "Count") +
  ylim(0, max(table(filtered_data$edu)) * 1.1) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold", margin = margin(b = 10)),  
    axis.title = element_text(size = 14),  
    axis.text.x = element_text(size = 12, angle = 45, hjust = 1, vjust = 1),  
    axis.text.y = element_text(size = 12),  
    panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "grey80"),  
    panel.grid.minor = element_blank(),  
    legend.position = "none"  
  )

The bar chart shows the distribution of education status in the dataset we have 25,386 individuals have a college education or more (64.8%), 693 have an elementary education or less (1.8%), and 13,472 have a high school education or equivalent (33.4%). The majority of individuals in the dataset have a college education or higher.

# Plot the bar chart
ggplot(filtered_data, aes(x = employstatus)) +
  geom_bar(fill = "blue4") +
  geom_text(stat = 'count', aes(label = ..count..), 
            vjust = -0.5, size = 4, color = "black", fontface = "bold") +
  labs(title = "Distribution of Employment Status", x = "Employment Status", y = "Count") +
  ylim(0, max(table(filtered_data$employstatus)) * 1.1) +  
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold", margin = margin(b = 10)),  
    axis.title = element_text(size = 14),  
    axis.text.x = element_text(size = 10, angle = 45, hjust = 1, vjust = 1),  
    axis.text.y = element_text(size = 12),  
    panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "grey80"),  
    panel.grid.minor = element_blank(),  
    legend.position = "none"  
  )

The bar chart shows the distribution of employment status in the dataset as here we have 13,734 individuals are employed for wages (34.9%), 2,855 are self-employed (7.3%), 1,542 have been out of work for 1 year or more (3.9%), 1,236 have been out of work for less than 1 year (3.1%), 1,330 are homemakers (3.4%), 451 are students (1.1%), 10,888 are retired (27.6%), and 7,515 are unable to work (19.1%). The majority of individuals in the dataset are employed for wages, followed by those who are retired.

# Plot the bar chart
ggplot(filtered_data, aes(x = income_group)) +
  geom_bar(fill = "blue4") +
  geom_text(stat = 'count', aes(label = ..count..), 
            vjust = -0.5, size = 4, color = "black", fontface = "bold") +
  labs(title = "Distribution of Income Status", x = "Income Status", y = "Frequency") +
  ylim(0, max(table(filtered_data$income_group)) * 1.1) +  
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold", margin = margin(b = 10)),  
    axis.title = element_text(size = 14),  
    axis.text.x = element_text(size = 12, angle = 45, hjust = 1, vjust = 1),  
    axis.text.y = element_text(size = 12),  
    panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "grey80"),  
    panel.grid.minor = element_blank(),  
    legend.position = "none"  
  )

The bar chart shows the distribution of income status in the dataset: 11,349 individuals have low income (29.0%), 17,526 have middle income (44.7%), and 10,676 have high income (27.2%). The majority of individuals in the dataset have middle income, followed by low income and high income.

# Plot the pie chart with improved aesthetics
ggplot(filtered_data, aes(x = "", fill = hrtdisease)) +
  geom_bar(width = 1, stat = "count", color = "black") +
  coord_polar("y") +
  geom_text(stat = 'count', aes(label = paste0(round(..count.. / sum(..count..) * 100, 1), "%")), 
            position = position_stack(vjust = 0.5), size = 5, color = "white", fontface = "bold") +
  scale_fill_manual(values = c("Yes" = "red4", "No" = "blue4")) +
  labs(title = "Distribution of Coronary Heart Disease (CHD) or Myocardial Infarction (MI)") +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold", margin = margin(t = 20, b = 10)),  
    axis.title = element_blank(),  
    axis.text = element_blank(),  
    axis.ticks = element_blank(),
    panel.grid = element_blank(),  
    legend.title = element_blank(),  
    legend.position = "bottom"
  )

The pie chart shows the distribution of coronary heart disease (CHD) or myocardial infarction (MI) status in the dataset: 11.5% of individuals have CHD or MI, while 88.5% do not have these conditions. This indicates that the vast majority of individuals in the dataset do not have CHD or MI.

# Plot the pie chart with improved aesthetics
ggplot(filtered_data, aes(x = "", fill = factor(totinda))) +
  geom_bar(width = 1, stat = "count", color = "black") +
  coord_polar("y") +
  geom_text(stat = 'count', aes(label = paste0(round(..count.. / sum(..count..) * 100, 1), "%")), 
            position = position_stack(vjust = 0.5), size = 5, color = "white", fontface = "bold") +
  scale_fill_manual(
    values = c("1" = "blue4", "2" = "red4"), 
    labels = c("1" = "Had physical activity or exercise", 
               "2" = "No physical activity or exercise")
  ) +
  labs(title = "Distribution of Physical Activity or Exercise During the Past 30 Days") +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold", margin = margin(t = 20, b = 10)),  
    axis.title = element_blank(),  
    axis.text = element_blank(),  
    axis.ticks = element_blank(),
    panel.grid = element_blank(),  
    legend.title = element_blank(),  
    legend.position = "bottom"
  )

The pie chart shows the distribution of physical activity or exercise during the past 30 days: 63.9% of individuals had physical activity or exercise, while 36.1% did not have any physical activity or exercise. This indicates that a majority of individuals in the dataset engaged in physical activity or exercise during the past month.

# Plot the pie chart with improved aesthetics
ggplot(filtered_data, aes(x = "", fill = factor(diabetic_status))) +
  geom_bar(width = 1, stat = "count", color = "black") +
  coord_polar("y") +
  geom_text(stat = 'count', aes(label = paste0(round(..count.. / sum(..count..) * 100, 1), "%")), 
            position = position_stack(vjust = 0.5), size = 5, color = "white", fontface = "bold") +
  scale_fill_manual(
    values = c("2" = "blue4", "1" = "red4"), 
    labels = c("1" = "Yes", "2" = "No")
  ) +
  labs(title = "Distribution of Diabetes") +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold", margin = margin(t = 20, b = 10)),  
    axis.title = element_blank(),  
    axis.text = element_blank(),  
    axis.ticks = element_blank(),
    panel.grid = element_blank(),  
    legend.title = element_blank(),  
    legend.position = "bottom"
  )

The pie chart shows the distribution of diabetes status in the dataset: 20.1% of individuals have diabetes, while 79.9% do not have diabetes. This indicates that a significant majority of individuals in the dataset do not have diabetes.

# Plot the bar chart
ggplot(filtered_data, aes(x = factor(genhlth, levels = 1:5, labels = c("Excellent", "Very good", "Good", "Fair", "Poor")))) +
  geom_bar(fill = "blue4") +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -0.5, size = 4, color = "black", fontface = "bold") +
  labs(title = "Distribution of General Health Status", x = "Health Status", y = "Frequency") +
  ylim(0, max(table(filtered_data$genhlth)) * 1.1) +  
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold", margin = margin(b = 10)),  
    axis.title = element_text(size = 14),  
    axis.text.x = element_text(size = 12, angle = 45, hjust = 1, vjust = 1),  
    axis.text.y = element_text(size = 12),  
    panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "grey80"),  
    panel.grid.minor = element_blank(),  
    legend.position = "none"  
  )

The bar chart shows the distribution of general health status in the dataset: 1,927 individuals reported their health as excellent (4.9%), 7,924 as very good (20.1%), 12,899 as good (32.8%), 10,872 as fair (27.6%), and 5,929 as poor (15.1%). The most common health status reported is good, followed by fair.

# Plot the bar chart
ggplot(filtered_data, aes(x = factor(smokeday, levels = 1:3, labels = c("Every day", "Some days", "Not at all")))) +
  geom_bar(fill = "blue4") +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -0.5, size = 4, color = "black", fontface = "bold") +
  labs(title = "Distribution of Smoking Status", x = "Smoking Status", y = "Frequency") +
  ylim(0, max(table(filtered_data$smokeday)) * 1.1) +  
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold", margin = margin(b = 10)),  
    axis.title = element_text(size = 14),  
    axis.text.x = element_text(size = 12, angle = 45, hjust = 1, vjust = 1),  
    axis.text.y = element_text(size = 12),  
    panel.grid.major = element_line(size = 0.5, linetype = 'solid', colour = "grey80"),  
    panel.grid.minor = element_blank(),  
    legend.position = "none"  
  )

The bar chart shows the distribution of smoking status in the dataset: 10,408 individuals smoke every day (21.4%), 4,062 smoke on some days (8.3%), and 25,081 do not smoke at all (70.3%). The majority of individuals in the dataset do not smoke.

Obesity rate over the states

In this part of the project, we wanted to include a map with the obesity rates for each respective state, so we understand the distribution of the obesity in the United States of America for creating a general view.

# Summarize the obesity rate by state
obesity_by_state <- filtered_data %>%
  group_by(state) %>%
  summarize(
    total_individuals = n(),
    obese_individuals = sum(obesity_status == "Obese"),
    obesity_rate = mean(obesity_status == "Obese")
  ) %>%
  arrange(desc(obesity_rate))

# Print the summarized obesity rates per state
print(obesity_by_state)

## # A tibble: 54 × 4
##    state         total_individuals obese_individuals obesity_rate
##    <fct>                     <int>             <int>        <dbl>
##  1 North Dakota                308               149        0.484
##  2 Wisconsin                  1084               509        0.470
##  3 Iowa                        805               377        0.468
##  4 Kentucky                    435               201        0.462
##  5 Ohio                       1676               762        0.455
##  6 South Dakota                603               273        0.453
##  7 West Virginia               630               278        0.441
##  8 Wyoming                     364               160        0.440
##  9 Delaware                    335               147        0.439
## 10 Oklahoma                    582               254        0.436
## # ℹ 44 more rows

# Load US map data
us_map <- map_data("state")

# Merge map data with obesity data
us_map$region <- tolower(us_map$region)
obesity_by_state$state <- tolower(obesity_by_state$state)
map_data <- merge(us_map, obesity_by_state, by.x = "region", by.y = "state", all.x = TRUE)

# Plot the map
ggplot(map_data, aes(x = long, y = lat, group = group, fill = obesity_rate)) +
  geom_polygon(color = "white", size = 0.2) +
  scale_fill_continuous(name = "Obesity Rate", low = "lightblue2", high = "darkred") +
  coord_fixed(1.3) + # Maintain aspect ratio
  theme_minimal() +
  labs(
    title = "Obesity Rates by State",
    x = NULL,
    y = NULL
  ) +
  theme(
    plot.title = element_text(hjust = 0.5, size = 20, face = "bold"),
    legend.position = "bottom",
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 10),
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank(),
    panel.border = element_blank()
  ) +
  guides(fill = guide_colorbar(barwidth = 10, barheight = 1))

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

The map provided above shows the obesity rates by state in the United States and some other parts of countries bordered. It is represented with varying shades of red representing different levels of obesity rates.

GENERAL TREND: Most states have high obesity rates, as indicated by the predominant dark red shading across the map.

States with the darkest red color have the highest obesity rates, close to 0.45. These states are mainly located in the South and parts of the Midwest. Specifically, North Dakota, Wisconsin and Iowa are the states with the highest rate of obesity.

Contingency Tables

Now, we will make some contingency tables. This type of tables show the distribution of individuals observed for the combinations of the categories of the 2 variables or more.

Joint Probability

Let us start with The joint probability

# Joint Probability table for obesity_status and sex
joint_prob_sex <- prop.table(table(filtered_data$obesity_status, filtered_data$sex))
kable(joint_prob_sex, digits = 3, caption = "Joint Probability of Obesity Status and Sex") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

Joint Probability of Obesity Status and Sex

	Male	Female
Not Obese	0.292	0.307
Obese	0.175	0.225

# Joint Probability table for obesity_status and age5cat
joint_prob_age <- prop.table(table(filtered_data$obesity_status, filtered_data$age5cat))
kable(joint_prob_age, digits = 3, caption = "Joint Probability of Obesity Status and Age Category") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

Joint Probability of Obesity Status and Age Category

	Age 18 to 29	Age 30 to 39	Age 40 to 54	Age 55 to 69	Age 70 and older
Not Obese	0.048	0.085	0.140	0.190	0.136
Obese	0.021	0.055	0.114	0.142	0.069

# Joint Probability table for obesity_status and income_group
joint_prob_income <- prop.table(table(filtered_data$obesity_status, filtered_data$income_group))
kable(joint_prob_income, digits = 3, caption = "Joint Probability of Obesity Status and Income Group") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

Joint Probability of Obesity Status and Income Group

	Low Income	Middle Income	High Income
Not Obese	0.167	0.264	0.168
Obese	0.120	0.179	0.102

In the Joint Probability, the probability is calculated of two events occurring simultaneously. It is calculated by dividing the count of occurrences of both events by the total number of observations.

Sex Differences in Obesity: There seems to be a difference in obesity rates between males and females, with females having a higher probability of being obese. (0.225 > 0.175). These probabilities show associations, but they don’t imply causation. For example, being female doesn’t cause obesity, nor does being male prevent it. Other factors not included in this table could contribute to obesity rates.

WE HAVE 5 categories of the age: Generally, the probability of being obese tends to increase with age, with the highest probability among individuals aged 55 to 69 (0.142). Age 70 and older is low rate.

MIDDLE INCOME has the highest obesity: For the probability of being obese, the middle-income group has the highest probability (0.179), followed by the low-income group (0.120) and the high-income group (0.102).

Conditional Probability

# Conditional Probability of obesity_status given sex
conditional_prob_obesity_sex <- prop.table(table(filtered_data$obesity_status, filtered_data$sex), 2)
kable(conditional_prob_obesity_sex, digits = 3, caption = "Conditional Probability of Obesity Status given Sex") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

Conditional Probability of Obesity Status given Sex

	Male	Female
Not Obese	0.626	0.577
Obese	0.374	0.423

# Conditional Probability of obesity_status given age5cat
conditional_prob_obesity_age <- prop.table(table(filtered_data$obesity_status, filtered_data$age5cat), 2)
kable(conditional_prob_obesity_age, digits = 3, caption = "Conditional Probability of Obesity Status given Age Category") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

Conditional Probability of Obesity Status given Age Category

	Age 18 to 29	Age 30 to 39	Age 40 to 54	Age 55 to 69	Age 70 and older
Not Obese	0.699	0.607	0.551	0.572	0.665
Obese	0.301	0.393	0.449	0.428	0.335

# Conditional Probability of obesity_status given income_group
conditional_prob_obesity_income <- prop.table(table(filtered_data$obesity_status, filtered_data$income_group), 2)
kable(conditional_prob_obesity_income, digits = 3, caption = "Conditional Probability of Obesity Status given Income Group") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

Conditional Probability of Obesity Status given Income Group

	Low Income	Middle Income	High Income
Not Obese	0.582	0.596	0.624
Obese	0.418	0.404	0.376

In here we wanted to see the probability of being obese occurring given that another event has already occurred which is Conditional Probability. It is calculated by dividing the joint probability of the two events by the probability of the conditioning event.

Sex Differences in Obesity: The joint probabilities indicated that females have a higher probability of being obese (42.3%) compared to males (37.4%).

Age Differences in Obesity: The conditional probabilities show an increasing trend of obesity with age, peaking at 44.9% for the 40 to 54 age group, and then decreasing to 33.5% for those aged 70 and older.

Income Group Differences in Obesity: The conditional probabilities reveal that the low-income group has the highest probability of being obese (41.8%), followed by the middle-income group (40.4%), and the high-income group (37.6%).

The Models

Multinomial Logistic Regression

library(VGAM)

## Loading required package: stats4

## Loading required package: splines

## 
## Attaching package: 'VGAM'

## The following objects are masked from 'package:mlogit':
## 
##     has.intercept, lrtest

## The following objects are masked from 'package:boot':
## 
##     logit, simplex

## The following object is masked from 'package:caret':
## 
##     predictors

## The following object is masked from 'package:car':
## 
##     logit

# Fit a multinomial logistic regression model
multinom_model <- vglm(bmi ~ hrtdisease + diabetic_status + age5cat + sex + edu + income_group + totinda + smokeday + genhlth, family = multinomial(refLevel = 4), data = filtered_data)

# Summary of the multinomial logistic regression model
summary(multinom_model)

## 
## Call:
## vglm(formula = bmi ~ hrtdisease + diabetic_status + age5cat + 
##     sex + edu + income_group + totinda + smokeday + genhlth, 
##     family = multinomial(refLevel = 4), data = filtered_data)
## 
## Coefficients: 
##                                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1                  -2.4048836  0.2738714  -8.781  < 2e-16 ***
## (Intercept):2                   0.3469836  0.1053895   3.292 0.000993 ***
## (Intercept):3                   0.1628343  0.1023726   1.591 0.111698    
## hrtdiseaseNo:1                  0.1890897  0.1182309   1.599 0.109748    
## hrtdiseaseNo:2                  0.2441503  0.0464568   5.255 1.48e-07 ***
## hrtdiseaseNo:3                  0.0841481  0.0400669   2.100 0.035712 *  
## diabetic_status2:1              1.7364858  0.1311493  13.241  < 2e-16 ***
## diabetic_status2:2              1.2419499  0.0396102  31.354  < 2e-16 ***
## diabetic_status2:3              0.6141595  0.0314411  19.534  < 2e-16 ***
## age5catAge 30 to 39:1          -1.2300101  0.1436864  -8.560  < 2e-16 ***
## age5catAge 30 to 39:2          -0.5111406  0.0587848  -8.695  < 2e-16 ***
## age5catAge 30 to 39:3          -0.0548685  0.0619160  -0.886 0.375522    
## age5catAge 40 to 54:1          -1.5911949  0.1359982 -11.700  < 2e-16 ***
## age5catAge 40 to 54:2          -0.6783984  0.0556025 -12.201  < 2e-16 ***
## age5catAge 40 to 54:3           0.0329873  0.0580859   0.568 0.570100    
## age5catAge 55 to 69:1          -0.9617011  0.1221174  -7.875 3.40e-15 ***
## age5catAge 55 to 69:2          -0.3719874  0.0547527  -6.794 1.09e-11 ***
## age5catAge 55 to 69:3           0.3002255  0.0575674   5.215 1.84e-07 ***
## age5catAge 70 and older:1      -0.0007404  0.1309952  -0.006 0.995490    
## age5catAge 70 and older:2       0.3218763  0.0592592   5.432 5.58e-08 ***
## age5catAge 70 and older:3       0.7882089  0.0615636  12.803  < 2e-16 ***
## sexFemale:1                     0.3233321  0.0754030   4.288 1.80e-05 ***
## sexFemale:2                    -0.0091983  0.0270983  -0.339 0.734278    
## sexFemale:3                    -0.3652744  0.0250112 -14.604  < 2e-16 ***
## eduElementary or less:1         0.6573484  0.1949345   3.372 0.000746 ***
## eduElementary or less:2        -0.1387187  0.1050661  -1.320 0.186735    
## eduElementary or less:3        -0.1111001  0.0970766  -1.144 0.252434    
## eduHigh school or equivalent:1  0.2405442  0.0773762   3.109 0.001879 ** 
## eduHigh school or equivalent:2 -0.0325848  0.0297005  -1.097 0.272592    
## eduHigh school or equivalent:3  0.0153035  0.0275344   0.556 0.578351    
## income_groupMiddle Income:1    -0.3693724  0.0819891  -4.505 6.63e-06 ***
## income_groupMiddle Income:2    -0.1812428  0.0329446  -5.501 3.77e-08 ***
## income_groupMiddle Income:3    -0.0172501  0.0311541  -0.554 0.579782    
## income_groupHigh Income:1      -0.6049170  0.1191839  -5.075 3.86e-07 ***
## income_groupHigh Income:2      -0.2528280  0.0400332  -6.315 2.69e-10 ***
## income_groupHigh Income:3       0.0601733  0.0370133   1.626 0.104009    
## totinda2:1                     -0.1692019  0.0787452  -2.149 0.031656 *  
## totinda2:2                     -0.4930243  0.0301621 -16.346  < 2e-16 ***
## totinda2:3                     -0.3726183  0.0274927 -13.553  < 2e-16 ***
## smokeday2:1                     0.0606895  0.1082555   0.561 0.575061    
## smokeday2:2                     0.0224077  0.0482167   0.465 0.642125    
## smokeday2:3                     0.0778567  0.0466524   1.669 0.095143 .  
## smokeday3:1                    -1.0055988  0.0840500 -11.964  < 2e-16 ***
## smokeday3:2                    -0.5320842  0.0322583 -16.494  < 2e-16 ***
## smokeday3:3                    -0.3145317  0.0304770 -10.320  < 2e-16 ***
## genhlth2:1                     -0.6544326  0.1925174  -3.399 0.000675 ***
## genhlth2:2                     -0.5631199  0.0688547  -8.178 2.88e-16 ***
## genhlth2:3                     -0.2429530  0.0697862  -3.481 0.000499 ***
## genhlth3:1                     -1.1204722  0.1840688  -6.087 1.15e-09 ***
## genhlth3:2                     -1.2546673  0.0669237 -18.748  < 2e-16 ***
## genhlth3:3                     -0.7275233  0.0674643 -10.784  < 2e-16 ***
## genhlth4:1                     -0.8042856  0.1844020  -4.362 1.29e-05 ***
## genhlth4:2                     -1.3016736  0.0691119 -18.834  < 2e-16 ***
## genhlth4:3                     -0.9270128  0.0695222 -13.334  < 2e-16 ***
## genhlth5:1                     -0.1831667  0.1914507  -0.957 0.338703    
## genhlth5:2                     -0.9798165  0.0745895 -13.136  < 2e-16 ***
## genhlth5:3                     -0.8701905  0.0748015 -11.633  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: log(mu[,1]/mu[,4]), log(mu[,2]/mu[,4]), 
## log(mu[,3]/mu[,4])
## 
## Residual deviance: 87081.85 on 118596 degrees of freedom
## 
## Log-likelihood: -43540.93 on 118596 degrees of freedom
## 
## Number of Fisher scoring iterations: 7 
## 
## Warning: Hauck-Donner effect detected in the following estimate(s):
## '(Intercept):1'
## 
## 
## Reference group is level  4  of the response

# Load necessary library
library(VGAM)

# Extract coefficients
coef_multinom <- coef(multinom_model)

# Extract standard errors
std_err <- sqrt(diag(vcov(multinom_model)))

# Convert coefficients to a matrix with 3 columns (for 3 categories: Underweight, Normal Weight, Overweight)
coef_matrix <- matrix(coef_multinom, ncol = 3, byrow = TRUE)
std_err_matrix <- matrix(std_err, ncol = 3, byrow = TRUE)

# Name the columns appropriately
colnames(coef_matrix) <- c("Underweight", "Normal Weight", "Overweight")
colnames(std_err_matrix) <- c("Underweight", "Normal Weight", "Overweight")

# Convert to data frames for easier manipulation
coef_df <- as.data.frame(coef_matrix)
std_err_df <- as.data.frame(std_err_matrix)

# Extract predictor names from the row names of the coefficients matrix
predictor_names <- rownames(coef(multinom_model))

# Assign predictor names to the data frames
rownames(coef_df) <- predictor_names
rownames(std_err_df) <- predictor_names

# Calculate odds ratios and confidence intervals
odds_ratios <- exp(coef_df)
conf_int_lower <- exp(coef_df - 1.96 * std_err_df)
conf_int_upper <- exp(coef_df + 1.96 * std_err_df)

# Convert to long format for plotting
library(reshape2)

## 
## Attaching package: 'reshape2'

## The following object is masked from 'package:tidyr':
## 
##     smiths

# Adding a predictor column for melting
odds_ratios$Predictor <- rownames(odds_ratios)
conf_int_lower$Predictor <- rownames(conf_int_lower)
conf_int_upper$Predictor <- rownames(conf_int_upper)

# Melting the data frames
odds_ratios_long <- melt(odds_ratios, id.vars = "Predictor", variable.name = "BMI_Category", value.name = "Odds_Ratio")
conf_int_lower_long <- melt(conf_int_lower, id.vars = "Predictor", variable.name = "BMI_Category", value.name = "CI_Lower")
conf_int_upper_long <- melt(conf_int_upper, id.vars = "Predictor", variable.name = "BMI_Category", value.name = "CI_Upper")

# Merge the data frames
conf_int_df <- merge(odds_ratios_long, conf_int_lower_long, by = c("Predictor", "BMI_Category"))
conf_int_df <- merge(conf_int_df, conf_int_upper_long, by = c("Predictor", "BMI_Category"))

library(ggplot2)

# Plot
ggplot(conf_int_df, aes(x = Odds_Ratio, y = Predictor, color = BMI_Category)) +
  geom_point(size = 2) +
  geom_errorbarh(aes(xmin = CI_Lower, xmax = CI_Upper), height = 0.2) +
  theme_minimal() +
  labs(title = "Forest Plot of Odds Ratios with Confidence Intervals", x = "Odds Ratio", y = "Predictor") +
  scale_x_continuous(trans = 'log10') +
  facet_wrap(~ BMI_Category, scales = "free_y")

# Load necessary libraries
library(VGAM)
library(reshape2)
library(ggplot2)

# Extract coefficients and standard errors
coef_multinom <- coef(multinom_model)
std_err <- sqrt(diag(vcov(multinom_model)))

# Reshape coefficients and standard errors into matrices
coef_matrix <- matrix(coef_multinom, ncol = 3, byrow = TRUE)
std_err_matrix <- matrix(std_err, ncol = 3, byrow = TRUE)

# Name the columns appropriately
colnames(coef_matrix) <- c("Underweight", "Normal Weight", "Overweight")
colnames(std_err_matrix) <- c("Underweight", "Normal Weight", "Overweight")

# Convert to data frames for easier manipulation
coef_df <- as.data.frame(coef_matrix)
std_err_df <- as.data.frame(std_err_matrix)

# Define predictor names
predictor_names <- c("(Intercept)", "Heart Disease", "Diabetes Status", "Age 30 to 39", "Age 40 to 54", 
                     "Age 55 to 69", "Age 70 and older", "Sex", "Elementary Education", "High School Education", 
                     "Middle Income", "High Income", "Physical Activity", "Smokeday2", "Smokeday3", 
                     "General Health 2", "General Health 3", "General Health 4", "General Health 5")

# Assign predictor names to the data frames
rownames(coef_df) <- predictor_names
rownames(std_err_df) <- predictor_names

# Calculate odds ratios and confidence intervals
odds_ratios <- exp(coef_df)
conf_int_lower <- exp(coef_df - 1.96 * std_err_df)
conf_int_upper <- exp(coef_df + 1.96 * std_err_df)

# Adding a predictor column for melting
odds_ratios$Predictor <- rownames(odds_ratios)
conf_int_lower$Predictor <- rownames(conf_int_lower)
conf_int_upper$Predictor <- rownames(conf_int_upper)

# Melting the data frames
odds_ratios_long <- melt(odds_ratios, id.vars = "Predictor", variable.name = "BMI_Category", value.name = "Odds_Ratio")
conf_int_lower_long <- melt(conf_int_lower, id.vars = "Predictor", variable.name = "BMI_Category", value.name = "CI_Lower")
conf_int_upper_long <- melt(conf_int_upper, id.vars = "Predictor", variable.name = "BMI_Category", value.name = "CI_Upper")

# Merge the data frames
conf_int_df <- merge(odds_ratios_long, conf_int_lower_long, by = c("Predictor", "BMI_Category"))
conf_int_df <- merge(conf_int_df, conf_int_upper_long, by = c("Predictor", "BMI_Category"))

# Function to create individual plots for each BMI category with different colors
create_forest_plot <- function(df, category, color) {
  ggplot(df[df$BMI_Category == category, ], aes(x = Odds_Ratio, y = Predictor)) +
    geom_point(size = 2, color = color) +
    geom_errorbarh(aes(xmin = CI_Lower, xmax = CI_Upper), height = 0.2, color = color) +
    theme_minimal() +
    labs(title = paste("Forest Plot of Odds Ratios with Confidence Intervals\n", category), 
         x = "Odds Ratio", y = "Predictor") +
    scale_x_continuous(trans = 'log10') +
    theme(axis.text.y = element_text(size = 10, angle = 0, hjust = 1)) +
    theme(legend.position = "none")
}

# Create individual plots with different colors
plot_underweight <- create_forest_plot(conf_int_df, "Underweight", "red")
plot_normal_weight <- create_forest_plot(conf_int_df, "Normal Weight", "darkgreen")
plot_overweight <- create_forest_plot(conf_int_df, "Overweight", "blue2")

# Display individual plots
print(plot_underweight)

print(plot_normal_weight)

print(plot_overweight)

Ordinal Logistic Regression

# Fit Model 1 (Age, Gender, Marital, Education)
fit1 <- vglm(bmi ~ age5cat + sex + marital + edu, family = cumulative(parallel = TRUE), data = filtered_data)

# Fit Model 2 (Add health variables)
fit2 <- vglm(bmi ~ age5cat + sex + marital + edu + genhlth + hrtdisease + diabetic_status + smokeday, family = cumulative(parallel = TRUE), data = filtered_data)

# Fit Model 3 (Add Employment Status and income)
fit3 <- vglm(bmi ~ age5cat + sex + marital + edu + genhlth + hrtdisease + diabetic_status + smokeday + income_group + employstatus, family = cumulative(parallel = TRUE), data = filtered_data)

# Calculate AIC and BIC for all three models
aic_fit1 <- AIC(fit1)
aic_fit2 <- AIC(fit2)
aic_fit3 <- AIC(fit3)

bic_fit1 <- BIC(fit1)
bic_fit2 <- BIC(fit2)
bic_fit3 <- BIC(fit3)

# Print AIC and BIC values for comparison
cat("Model 1 AIC:", aic_fit1, "\n")

## Model 1 AIC: 91360.78

cat("Model 2 AIC:", aic_fit2, "\n")

## Model 2 AIC: 88514.37

cat("Model 3 AIC:", aic_fit3, "\n")

## Model 3 AIC: 88376.06

cat("Model 1 BIC:", bic_fit1, "\n")

## Model 1 BIC: 91489.56

cat("Model 2 BIC:", bic_fit2, "\n")

## Model 2 BIC: 88711.83

cat("Model 3 BIC:", bic_fit3, "\n")

## Model 3 BIC: 88650.79

# Determine the best model based on AIC and BIC
if (aic_fit1 < aic_fit2 && aic_fit1 < aic_fit3 && bic_fit1 < bic_fit2 && bic_fit1 < bic_fit3) {
    cat("Model 1 is better based on AIC and BIC.\n")
    summary(fit1)
} else if (aic_fit2 < aic_fit1 && aic_fit2 < aic_fit3 && bic_fit2 < bic_fit1 && bic_fit2 < bic_fit3) {
    cat("Model 2 is better based on AIC and BIC.\n")
    summary(fit2)
} else {
    cat("Model 3 is better based on AIC and BIC.\n")
    summary(fit3)
}

## Model 3 is better based on AIC and BIC.

## 
## Call:
## vglm(formula = bmi ~ age5cat + sex + marital + edu + genhlth + 
##     hrtdisease + diabetic_status + smokeday + income_group + 
##     employstatus, family = cumulative(parallel = TRUE), data = filtered_data)
## 
## Coefficients: 
##                                              Estimate Std. Error z value
## (Intercept):1                                -3.69759    0.08574 -43.126
## (Intercept):2                                -0.72634    0.07975  -9.107
## (Intercept):3                                 0.71596    0.07981   8.971
## age5catAge 30 to 39                          -0.46114    0.04446 -10.373
## age5catAge 40 to 54                          -0.59425    0.04355 -13.647
## age5catAge 55 to 69                          -0.36814    0.04540  -8.109
## age5catAge 70 and older                       0.05793    0.05312   1.090
## sexFemale                                    -0.03967    0.01962  -2.022
## maritalDivorced                               0.12163    0.02721   4.471
## maritalWidowed                                0.16147    0.03507   4.604
## maritalSeparated                              0.12488    0.05356   2.332
## maritalNever married                          0.18672    0.03005   6.213
## maritalMember of an unmarried couple          0.19952    0.04348   4.589
## eduElementary or less                        -0.02087    0.07386  -0.283
## eduHigh school or equivalent                 -0.01607    0.02122  -0.757
## genhlth2                                     -0.36019    0.04707  -7.652
## genhlth3                                     -0.85422    0.04581 -18.646
## genhlth4                                     -0.94138    0.04757 -19.791
## genhlth5                                     -0.71036    0.05184 -13.704
## hrtdiseaseNo                                  0.15178    0.03198   4.746
## diabetic_status2                              0.93148    0.02606  35.745
## smokeday2                                     0.04171    0.03452   1.208
## smokeday3                                    -0.40383    0.02314 -17.450
## income_groupMiddle Income                    -0.08259    0.02559  -3.227
## income_groupHigh Income                      -0.04537    0.03275  -1.385
## employstatusSelf-employed                     0.36847    0.03865   9.535
## employstatusOut of work for 1 year or more    0.25824    0.05226   4.941
## employstatusOut of work for less than 1 year  0.13620    0.05636   2.417
## employstatusA homemaker                       0.22110    0.05511   4.012
## employstatusA student                         0.31627    0.09035   3.501
## employstatusRetired                           0.20287    0.03402   5.963
## employstatusUnable to work                    0.03510    0.03362   1.044
##                                              Pr(>|z|)    
## (Intercept):1                                 < 2e-16 ***
## (Intercept):2                                 < 2e-16 ***
## (Intercept):3                                 < 2e-16 ***
## age5catAge 30 to 39                           < 2e-16 ***
## age5catAge 40 to 54                           < 2e-16 ***
## age5catAge 55 to 69                          5.11e-16 ***
## age5catAge 70 and older                      0.275551    
## sexFemale                                    0.043223 *  
## maritalDivorced                              7.80e-06 ***
## maritalWidowed                               4.15e-06 ***
## maritalSeparated                             0.019727 *  
## maritalNever married                         5.21e-10 ***
## maritalMember of an unmarried couple         4.46e-06 ***
## eduElementary or less                        0.777548    
## eduHigh school or equivalent                 0.448844    
## genhlth2                                     1.98e-14 ***
## genhlth3                                      < 2e-16 ***
## genhlth4                                      < 2e-16 ***
## genhlth5                                      < 2e-16 ***
## hrtdiseaseNo                                 2.08e-06 ***
## diabetic_status2                              < 2e-16 ***
## smokeday2                                    0.226869    
## smokeday3                                     < 2e-16 ***
## income_groupMiddle Income                    0.001249 ** 
## income_groupHigh Income                      0.165965    
## employstatusSelf-employed                     < 2e-16 ***
## employstatusOut of work for 1 year or more   7.75e-07 ***
## employstatusOut of work for less than 1 year 0.015662 *  
## employstatusA homemaker                      6.03e-05 ***
## employstatusA student                        0.000464 ***
## employstatusRetired                          2.48e-09 ***
## employstatusUnable to work                   0.296433    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: logitlink(P[Y<=1]), logitlink(P[Y<=2]), 
## logitlink(P[Y<=3])
## 
## Residual deviance: 88312.06 on 118621 degrees of freedom
## 
## Log-likelihood: -44156.03 on 118621 degrees of freedom
## 
## Number of Fisher scoring iterations: 5 
## 
## No Hauck-Donner effect found in any of the estimates
## 
## 
## Exponentiated coefficients:
##                          age5catAge 30 to 39 
##                                    0.6305625 
##                          age5catAge 40 to 54 
##                                    0.5519764 
##                          age5catAge 55 to 69 
##                                    0.6920195 
##                      age5catAge 70 and older 
##                                    1.0596363 
##                                    sexFemale 
##                                    0.9611102 
##                              maritalDivorced 
##                                    1.1293322 
##                               maritalWidowed 
##                                    1.1752334 
##                             maritalSeparated 
##                                    1.1330172 
##                         maritalNever married 
##                                    1.2052865 
##         maritalMember of an unmarried couple 
##                                    1.2208138 
##                        eduElementary or less 
##                                    0.9793491 
##                 eduHigh school or equivalent 
##                                    0.9840568 
##                                     genhlth2 
##                                    0.6975454 
##                                     genhlth3 
##                                    0.4256160 
##                                     genhlth4 
##                                    0.3900893 
##                                     genhlth5 
##                                    0.4914685 
##                                 hrtdiseaseNo 
##                                    1.1639074 
##                             diabetic_status2 
##                                    2.5382754 
##                                    smokeday2 
##                                    1.0425929 
##                                    smokeday3 
##                                    0.6677603 
##                    income_groupMiddle Income 
##                                    0.9207249 
##                      income_groupHigh Income 
##                                    0.9556413 
##                    employstatusSelf-employed 
##                                    1.4455167 
##   employstatusOut of work for 1 year or more 
##                                    1.2946513 
## employstatusOut of work for less than 1 year 
##                                    1.1459120 
##                      employstatusA homemaker 
##                                    1.2474479 
##                        employstatusA student 
##                                    1.3719977 
##                          employstatusRetired 
##                                    1.2249103 
##                   employstatusUnable to work 
##                                    1.0357211

# Extract intercepts from the best model (fit3)
intercepts <- coef(fit3)[grepl("Intercept", names(coef(fit3)))]
intercept_names <- names(intercepts)

# Create a data frame for plotting
intercept_data <- data.frame(
  Category = factor(intercept_names, levels = intercept_names),
  Estimate = intercepts
)

# Plot the intercepts
library(ggplot2)

ggplot(intercept_data, aes(x = Category, y = Estimate)) +
  geom_point(color = "blue", size = 3) +
  geom_line(group = 1, color = "blue", linetype = "dashed") +
  labs(title = "Intercepts for BMI Categories",
       x = "BMI Category Thresholds",
       y = "Intercept Estimate") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Load necessary libraries
library(VGAM)
library(ggplot2)
library(dplyr)

# Extract coefficients
coef_values <- coef(fit3)
intercepts <- coef_values[grepl("Intercept", names(coef_values))]
betas <- coef_values[!grepl("Intercept", names(coef_values))]
beta_names <- names(betas)

# Create a data frame for plotting betas with different intercepts
plot_data <- data.frame(
  Coefficient = rep(beta_names, 3),
  Estimate = rep(betas, 3),
  Category = rep(c("Underweight", "Normal_weight", "Overweight"), each = length(betas)),
  Intercept = rep(intercepts, each = length(betas))
)

# Calculate the y-values for the lines with different intercepts
plot_data$Fitted_Value <- plot_data$Estimate + plot_data$Intercept

# Plot
ggplot(plot_data, aes(x = Coefficient, y = Fitted_Value, color = Category, group = Category)) +
  geom_point(size = 3) +
  geom_line() +
  labs(title = "Betas for BMI Categories with Different Intercepts",
       x = "Predictor",
       y = "Fitted Value") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))