Aging and Depression: Random Forest & Neural Network Analyses (in progress)
Shannon, Choi, Hites
2024-11-08
Predicting Depression in Older Adults in the Southeastern USA: Comparing Random Forest and Linear Regression Models
Introduction
This analysis explores the predictors of depressive symptoms among older adults using machine learning techniques, including Random Forest and Multiple Regression. (Update from mmanuscript)
# Install required packages
#install.packages("haven")
#install.packages("psych")
#install.packages("ggplot2")
#install.packages("reshape2")
#install.packages("corrplot")
#install.packages("randomForest")
#install.packages("nnet")
#install.packages("pROC")
#install.packages("patchwork")
#install.packages("iml")
#install.packages("grid")
#install.packages("gridExtra")
#install.packages("NeuralNetTools")
#install.packages("png")
#install.packages("themis")
#install.packages("smotefamily")
# Load libraries
library(dplyr)
library(smotefamily)
library(themis)
library(haven)
library(psych)
library(ggplot2)
library(reshape2)
library(corrplot)
library(randomForest)
library(nnet)
library(pROC)
library(grid)
library(gridExtra)
library(NeuralNetTools)
library(png)
library(caret)Data Cleaning and Preparation
The data is imported, cleaned, and variables are transformed or categorized for analysis.
#Load the data
data_ESC <- read_sav("C:/Users/nhshannoniv/Downloads/Choi_10-20-24_EastSouthCentral Respondents only_deleteSpouseVar_2012-2020LBQ+randhrs1992_2018v1 (1).sav")
#Start Data Cleaning
data_ESC <- data.frame(
depresv_symp = data_ESC$R14CESD, # Number of depressive symptoms
diff_day_actt = data_ESC$R14ADLA, # Difficulty in daily activities
age506580 = data_ESC$r14Age506580, # Age categorized into 50, 65, and 80 years
self_rt_help = data_ESC$R14SRHpf, # Self-rated health perception
num_comorbidity = data_ESC$r14comorbidityNopsych, # Number of comorbid conditions excluding psychological disorders
cog_status = data_ESC$r13CogStat3, # Cognitive status (3 levels)
female = data_ESC$female, # Gender indicator (1 = female, 0 = male)
marital_stat = data_ESC$r14MS, # Marital status
num_ppl_in_house = data_ESC$H14HHRES, # Number of people living in the household
livealone = data_ESC$h14Livealone, # Indicator for living alone
yrs_of_edu = data_ESC$RAEDYRS, # Years of education completed
employment_type = data_ESC$r14employment, # Employment type (e.g., full-time, part-time)
non_house_wealth_log = data_ESC$LognohsngW, # Log-transformed non-housing wealth
house_wealth_log = data_ESC$LoghousingW, # Log-transformed housing wealth
income_log = data_ESC$Loghhincome, # Log-transformed household income
insurance = data_ESC$HealthInsu # Health insurance coverage
)
#Additional cleaning
data_ESC$age506580[data_ESC$age506580 == 0] <- NA
data_ESC$num_comorbidity[data_ESC$num_comorbidity %in% c(5, 6, 7)] <- 4
# Transform variables
data_ESC$depresv_symp <- as.numeric(data_ESC$depresv_symp)
data_ESC$depresv_symp <- cut(data_ESC$depresv_symp, breaks = c(-Inf, 2, Inf),
labels = c("0-2 Symptoms", "3+ Symptoms"))
# Convert variables to factors
categorical_vars <- c("depresv_symp", "age506580", "cog_status", "female", "marital_stat",
"livealone", "employment_type", "insurance")
data_ESC[categorical_vars] <- lapply(data_ESC[categorical_vars], as.factor)
#changing the reference categories
data_ESC$marital_stat <- relevel(data_ESC$marital_stat, ref = 2) #Single
data_ESC$employment_type <- relevel(data_ESC$employment_type, ref = 1) #Full time
data_ESC$insurance <- relevel(data_ESC$insurance, ref = 2) #Employer
# Handling missing data
# Pairwise deletion is applied here to exclude rows with any missing values. While this simplifies the analysis,
# it may introduce bias if the excluded data is not missing completely at random. Future work might consider
# using imputation techniques to preserve sample size and reduce potential biases.
data_ESC <- na.omit(data_ESC)#Splt the data into Training and Test dataset
# Train-Test Split
set.seed(1989)
sample_index <- sample(1:nrow(data_ESC), size = 0.7 * nrow(data_ESC))
train_data <- data_ESC[sample_index, ]
test_data <- data_ESC[-sample_index, ]Correlation Table
# Compute and visualize correlations for numerical predictors
numeric_data <- train_data[, c("yrs_of_edu", "non_house_wealth_log", "house_wealth_log",
"income_log", "num_comorbidity",
"num_ppl_in_house")]
correlation_matrix <- cor(numeric_data, use = "complete.obs")
corrplot(correlation_matrix, method = "color", col = colorRampPalette(c("blue", "white", "red"))(200),
type = "upper", addCoef.col = "black", tl.col = "black", tl.srt = 45,
title = "Correlation Matrix for Training Data")
correlation_matrix## yrs_of_edu non_house_wealth_log house_wealth_log
## yrs_of_edu 1.0000000 0.29461099 0.21452628
## non_house_wealth_log 0.2946110 1.00000000 0.33065844
## house_wealth_log 0.2145263 0.33065844 1.00000000
## income_log 0.3508956 0.33876979 0.29939443
## num_comorbidity -0.1928991 -0.08071173 -0.06263273
## num_ppl_in_house -0.1309180 -0.06397789 0.01963504
## income_log num_comorbidity num_ppl_in_house
## yrs_of_edu 0.35089557 -0.192899144 -0.130917977
## non_house_wealth_log 0.33876979 -0.080711730 -0.063977887
## house_wealth_log 0.29939443 -0.062632729 0.019635043
## income_log 1.00000000 -0.120470815 -0.022499988
## num_comorbidity -0.12047082 1.000000000 -0.005705693
## num_ppl_in_house -0.02249999 -0.005705693 1.000000000Random Forest Analysis
Classification Table
A Random Forest model is used to predict depressive symptom categories.
# Train Random Forest on training data
rf_model <- randomForest(depresv_symp ~ ., data = train_data, importance = TRUE, ntree = 500)
rf_model##
## Call:
## randomForest(formula = depresv_symp ~ ., data = train_data, importance = TRUE, ntree = 500)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 3
##
## OOB estimate of error rate: 22.79%
## Confusion matrix:
## 0-2 Symptoms 3+ Symptoms class.error
## 0-2 Symptoms 500 46 0.08424908
## 3+ Symptoms 124 76 0.62000000# Evaluate performance on test data
rf_predictions <- predict(rf_model, newdata = test_data)
confusion_matrix_rf <- table(Predicted = rf_predictions, Actual = test_data$depresv_symp)
confusion_matrix_rf## Actual
## Predicted 0-2 Symptoms 3+ Symptoms
## 0-2 Symptoms 232 49
## 3+ Symptoms 14 26Variable Importance Plot
# Display model summary and plot variable importance
varImpPlot(rf_model, main = "Variable Importance in Predicting Depressive Symptoms (depresv_symp)")
# Interpretation
# The variable importance plot highlights which predictors contribute the most to the model's performance. For example, variables with high importance scores like `yrs_of_edu` or `income_log` indicate a stronger relationship with depressive symptoms. This insight helps to prioritize factors for potential interventions or further study.
varImp(rf_model)## 0-2 Symptoms 3+ Symptoms
## diff_day_actt 18.52687980 18.52687980
## age506580 5.89468221 5.89468221
## self_rt_help 19.91811770 19.91811770
## num_comorbidity 6.45406990 6.45406990
## cog_status -0.41670847 -0.41670847
## female 0.28213990 0.28213990
## marital_stat 6.36843836 6.36843836
## num_ppl_in_house 0.03640394 0.03640394
## livealone 1.41967920 1.41967920
## yrs_of_edu -0.81597597 -0.81597597
## employment_type 3.61980679 3.61980679
## non_house_wealth_log 11.37233905 11.37233905
## house_wealth_log 7.78917495 7.78917495
## income_log 3.89733283 3.89733283
## insurance 3.89741615 3.89741615# Compute predicted probabilities for the test data
rf_probabilities <- predict(rf_model, newdata = test_data, type = "prob")[, 2]
# Convert the actual outcome to a numeric binary variable (if necessary)
# Assuming "3+ Symptoms" is the positive class and "0-2 Symptoms" is the negative class
test_data$depresv_symp_numeric <- ifelse(test_data$depresv_symp == "3+ Symptoms", 1, 0)
# Generate the ROC curve and compute AUC
roc_rf <- roc(test_data$depresv_symp_numeric, rf_probabilities)## Setting levels: control = 0, case = 1## Setting direction: controls < cases# Explain ROC and AUC
# ROC (Receiver Operating Characteristic) curves illustrate the diagnostic ability of a binary classifier system as its discrimination threshold is varied.
# AUC (Area Under the Curve) provides a single scalar value to summarize the performance of the model, with higher values indicating better performance.
# Plot the ROC curve
plot.roc(roc_rf, main = "ROC Curve for Random Forest Model",
col = "blue", lwd = 2, print.auc = TRUE, auc.polygon = TRUE,
auc.polygon.col = "lightblue", grid = TRUE)
#Comparing RF to linear regression
# Train a Multiple Linear Regression model
# Convert depresv_symp to a factor if not already
train_data$depresv_symp <- as.factor(train_data$depresv_symp)
test_data$depresv_symp <- as.factor(test_data$depresv_symp)
# Fit the logistic regression model
logistic_model <- glm(depresv_symp ~ ., data = train_data, family = binomial)
# Display model summary
summary(logistic_model)##
## Call:
## glm(formula = depresv_symp ~ ., family = binomial, data = train_data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.525905 1.093076 -3.226 0.00126 **
## diff_day_actt 0.450734 0.101297 4.450 8.60e-06 ***
## age5065802 -0.784446 0.287453 -2.729 0.00635 **
## age5065803 -1.071947 0.384635 -2.787 0.00532 **
## self_rt_help 1.245741 0.219893 5.665 1.47e-08 ***
## num_comorbidity 0.176562 0.091821 1.923 0.05449 .
## cog_status1 0.027516 0.268488 0.102 0.91837
## cog_status2 0.351085 0.486109 0.722 0.47015
## female1 0.009433 0.211809 0.045 0.96448
## marital_stat1 -0.412535 0.395394 -1.043 0.29679
## marital_stat3 0.604520 0.456501 1.324 0.18542
## marital_stat4 0.983545 0.477457 2.060 0.03940 *
## marital_stat5 0.849041 0.509426 1.667 0.09558 .
## num_ppl_in_house 0.038308 0.122042 0.314 0.75360
## livealone1 -0.187347 0.350696 -0.534 0.59319
## yrs_of_edu 0.026245 0.042262 0.621 0.53460
## employment_type2 -0.355501 0.409028 -0.869 0.38477
## employment_type3 0.372871 0.331065 1.126 0.26005
## employment_type4 1.142393 0.429886 2.657 0.00787 **
## non_house_wealth_log -0.070654 0.033891 -2.085 0.03709 *
## house_wealth_log 0.044272 0.044753 0.989 0.32254
## income_log 0.225221 0.161399 1.395 0.16288
## insurance0 0.337294 0.392787 0.859 0.39050
## insurance2 -0.133175 0.375620 -0.355 0.72293
## insurance3 0.095173 0.464556 0.205 0.83767
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 867.38 on 745 degrees of freedom
## Residual deviance: 671.17 on 721 degrees of freedom
## AIC: 721.17
##
## Number of Fisher Scoring iterations: 5# Predict on test data
lm_predictions <- predict(logistic_model, newdata = test_data)
# Evaluate the performance of the model
mse_lm <- mean((test_data$depresv_symp - lm_predictions)^2)## Warning in Ops.factor(test_data$depresv_symp, lm_predictions): '-' not
## meaningful for factorscat("Mean Squared Error (Linear Regression):", mse_lm, "\n")## Mean Squared Error (Linear Regression): NA# Generate the ROC curve and compute AUC for regression (using predicted probabilities)
roc_lm <- roc(test_data$depresv_symp, lm_predictions)## Setting levels: control = 0-2 Symptoms, case = 3+ Symptoms## Setting direction: controls < cases# Plot the ROC curve
plot.roc(roc_lm, main = "ROC Curve for Linear Regression Model",
col = "red", lwd = 2, print.auc = TRUE, auc.polygon = TRUE,
auc.polygon.col = "pink", grid = TRUE)
# Define a threshold for classification
# Assuming a cutoff of 2.5 to separate "0-2 Symptoms" from "3+ Symptoms"
threshold <- 2.5
lm_predictions_class <- ifelse(lm_predictions >= threshold, "3+ Symptoms", "0-2 Symptoms")
# Create a confusion matrix
confusion_matrix_lm <- table(Predicted = lm_predictions_class, Actual = test_data$depresv_symp)
# Display the confusion matrix
cat("Confusion Matrix for Linear Regression Model:\n")## Confusion Matrix for Linear Regression Model:print(confusion_matrix_lm)## Actual
## Predicted 0-2 Symptoms 3+ Symptoms
## 0-2 Symptoms 245 75
## 3+ Symptoms 1 0# Compute overall accuracy
accuracy_lm <- sum(diag(confusion_matrix_lm)) / sum(confusion_matrix_lm)
cat("Accuracy (Linear Regression):", round(accuracy_lm * 100, 2), "%\n")## Accuracy (Linear Regression): 76.32 %