A national veterans organization is trying to develop a model that helps their marketing team predict donors vs non-donors on a mailing campaign. Their main database contains observations on 13 million donors and recent analytics have uncovered then only 5.1% of those individuals actually donate anything. With average donations being $13.00 and the organization mailing costs being $0.68, they wish to maximize their profit margin by removing non-donors from their mailing lists.
The goal of this project will be to fit three classification models on a sampled data set consisting of 3,000 donors and chosen which model performs the best. The highest performing model will then be used to predict future donors from non-donors using a threshold of 32% on its predictions. The final model will categorizing a person with a likelihood higher than 32% as a donor and less than 32% as a non-donor.
Two data files are being used within this project -
Fundraising.csv - The first data set is a CSV file of fundraising data containing a weighted sample from a database used by a national veterans organization consisting of 3000 observations across 21 columns with 50% of the data being donors and the other 50% non-donors. The weighted sample containing a 50/50 split is necessary because it prevents imbalances and bias in the final results of the model. Fitting a sample that had a higher ratio of donor than non-donors would over represent one class and under represent another.
Future_fundraising.csv - The second data set contains data related to future candidates that the campaign wishes to target.
This data was loaded in from a local file directory on my desktop -
the original file is a download from Canvas. I will store the main data
set in an object titled data and the future data set in an
object titled future_data
# Load data from a local working directory
data <- read.csv("C:/Users/justi/OneDrive/Desktop/Grad School/UTSA 3rd semester/STA 6543 - Predictive Modeling/Week 10 - Final Project/fundraising.csv")
future_data <- read_csv("C:/Users/justi/OneDrive/Desktop/Grad School/UTSA 3rd semester/STA 6543 - Predictive Modeling/Week 10 - Final Project/future_fundraising.csv",
show_col_types = FALSE)
The main data contains 3000 rows and 21 columns and the future data contains 120 rows and 20 columns. The main data contains 7 character variables, 10 integer, and 4 numeric variables. Similarly, the future data consists of 6 character variables and 14 numeric.
glimpse(data)
## Rows: 3,000
## Columns: 21
## $ zipconvert2 <chr> "Yes", "No", "No", "No", "No", "No", "No", "Yes", …
## $ zipconvert3 <chr> "No", "No", "No", "Yes", "Yes", "No", "No", "No", …
## $ zipconvert4 <chr> "No", "No", "No", "No", "No", "No", "Yes", "No", "…
## $ zipconvert5 <chr> "No", "Yes", "Yes", "No", "No", "Yes", "No", "No",…
## $ homeowner <chr> "Yes", "No", "Yes", "Yes", "Yes", "Yes", "Yes", "Y…
## $ num_child <int> 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ income <int> 1, 5, 3, 4, 4, 4, 4, 4, 4, 1, 4, 5, 2, 3, 4, 4, 2,…
## $ female <chr> "No", "Yes", "No", "No", "Yes", "Yes", "No", "Yes"…
## $ wealth <int> 7, 8, 4, 8, 8, 8, 5, 8, 8, 5, 5, 8, 8, 5, 6, 9, 7,…
## $ home_value <int> 698, 828, 1471, 547, 482, 857, 505, 1438, 1316, 42…
## $ med_fam_inc <int> 422, 358, 484, 386, 242, 450, 333, 458, 541, 203, …
## $ avg_fam_inc <int> 463, 376, 546, 432, 275, 498, 388, 533, 575, 271, …
## $ pct_lt15k <int> 4, 13, 4, 7, 28, 5, 16, 8, 11, 39, 6, 8, 5, 3, 13,…
## $ num_prom <int> 46, 32, 94, 20, 38, 47, 51, 21, 66, 73, 59, 25, 27…
## $ lifetime_gifts <dbl> 94, 30, 177, 23, 73, 139, 63, 26, 108, 161, 84, 40…
## $ largest_gift <dbl> 12, 10, 10, 11, 10, 20, 15, 16, 12, 6, 5, 10, 20, …
## $ last_gift <dbl> 12, 5, 8, 11, 10, 20, 10, 16, 7, 3, 3, 10, 20, 7, …
## $ months_since_donate <int> 34, 29, 30, 30, 31, 37, 37, 30, 31, 32, 30, 32, 37…
## $ time_lag <int> 6, 7, 3, 6, 3, 3, 8, 6, 1, 7, 12, 2, 7, 1, 10, 3, …
## $ avg_gift <dbl> 9.400000, 4.285714, 7.080000, 7.666667, 7.300000, …
## $ target <chr> "Donor", "Donor", "No Donor", "No Donor", "Donor",…
glimpse(future_data)
## Rows: 120
## Columns: 20
## $ zipconvert2 <chr> "No", "Yes", "No", "No", "No", "Yes", "No", "Yes",…
## $ zipconvert3 <chr> "Yes", "No", "No", "No", "Yes", "No", "Yes", "No",…
## $ zipconvert4 <chr> "No", "No", "No", "Yes", "No", "No", "No", "No", "…
## $ zipconvert5 <chr> "No", "No", "Yes", "No", "No", "No", "No", "No", "…
## $ homeowner <chr> "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "…
## $ num_child <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ income <dbl> 5, 1, 4, 4, 2, 4, 2, 3, 4, 2, 6, 1, 4, 4, 1, 4, 4,…
## $ female <chr> "Yes", "No", "Yes", "No", "Yes", "Yes", "Yes", "No…
## $ wealth <dbl> 9, 7, 1, 8, 7, 8, 1, 8, 3, 5, 8, 8, 8, 4, 8, 8, 7,…
## $ home_value <dbl> 1399, 1355, 835, 1019, 992, 834, 639, 457, 349, 69…
## $ med_fam_inc <dbl> 637, 411, 310, 389, 524, 371, 209, 253, 302, 335, …
## $ avg_fam_inc <dbl> 703, 497, 364, 473, 563, 408, 259, 285, 324, 348, …
## $ pct_lt15k <dbl> 1, 9, 22, 15, 6, 10, 36, 25, 19, 14, 4, 40, 7, 5, …
## $ num_prom <dbl> 74, 77, 70, 21, 63, 35, 72, 68, 55, 59, 19, 47, 31…
## $ lifetime_gifts <dbl> 102.00, 249.00, 126.00, 26.00, 100.00, 92.00, 146.…
## $ largest_gift <dbl> 6.00, 15.00, 6.00, 16.00, 20.00, 37.00, 12.00, 5.0…
## $ last_gift <dbl> 5, 7, 6, 16, 3, 37, 11, 3, 5, 13, 30, 5, 15, 12, 1…
## $ months_since_donate <dbl> 29, 35, 34, 37, 21, 37, 36, 32, 30, 33, 37, 17, 30…
## $ time_lag <dbl> 3, 3, 8, 5, 6, 5, 5, 9, 9, 10, 5, 1, 6, 1, 6, 3, 1…
## $ avg_gift <dbl> 4.857143, 9.576923, 4.344828, 13.000000, 7.692308,…
table(sapply(data, class))
##
## character integer numeric
## 7 10 4
table(sapply(future_data, class))
##
## character numeric
## 6 14
To begin my EDA I scanned for null values within the data set and found that there were none. A univariate and bivariate analysis were used to investigate single variables using histogram plots to evaluate the distribution across individual variables and scatter plots to evaluate dual relationships. Next, I used correlation plots and the plot() function on the holistic data set to scan for additional univariate relationships concluding by selecting variables that seemed to have co linear relationships or unique patterns.
The following code reveals that the future and main dataset contained no null values giving the green light to proceed.
# Scan all columns for NAs
colSums(is.na(data))
## zipconvert2 zipconvert3 zipconvert4 zipconvert5
## 0 0 0 0
## homeowner num_child income female
## 0 0 0 0
## wealth home_value med_fam_inc avg_fam_inc
## 0 0 0 0
## pct_lt15k num_prom lifetime_gifts largest_gift
## 0 0 0 0
## last_gift months_since_donate time_lag avg_gift
## 0 0 0 0
## target
## 0
colSums(is.na(future_data))
## zipconvert2 zipconvert3 zipconvert4 zipconvert5
## 0 0 0 0
## homeowner num_child income female
## 0 0 0 0
## wealth home_value med_fam_inc avg_fam_inc
## 0 0 0 0
## pct_lt15k num_prom lifetime_gifts largest_gift
## 0 0 0 0
## last_gift months_since_donate time_lag avg_gift
## 0 0 0 0
The univariate analysis uncovered many insights on the data’s individual variables as follows…
# Separate the numeric variables from the character
num <- select(data, where(is.numeric))
# plot all numeric variables
num |>
pivot_longer(everything()) |>
ggplot(aes(value)) +
geom_histogram(bins = 30, fill = "darkblue") +
facet_wrap(~name, scales = "free")
# Separate categorical variables
cat <- select(data, where(is.character))
# Plot all categorical variables
cat |>
pivot_longer(everything()) |>
ggplot(aes(x = value)) +
geom_bar(fill = "darkblue") +
facet_wrap(~name, scales = "free") +
theme_bw()
Observations: Several of the variable contain right skewed distributions and outlier categories highlighted below:
Avg_fam_inc and med_fam_inc - These variables contained a right skewed distribution showing a significant portion of the data was gathered on individuals with average income around $31,800 - $51,600 indicating this data was calculated on a lower middle class area.
home_value - This variable also contained a major right skew revealing a mean home value of $114,300.
income - This variable contained a major outlier category in bracket number 4, indicating that individuals in this income class were heavily sampled amongst the data’s observations.
time_lag - This variable contained a large spike in its right skewed pattern indicating that a large majority of individuals in this data set often spend around 7 months between giving their 1st and 2nd gift.
pct_lt15k - This variable contained a strong right skew in the 0-21 count section indicating that a large majority of the donor neighborhoods contain individuals that earn less than $15,000 annually.
wealth - This variable contained a large outlier spike in the 8th wealth class indicating this group was heavily sampled in the data.
Females - This categorical variable showed a YES category that was about 30% larger than the NO category indicating that females are sampled more heavily in this data set.
Target - The response variable that defines whether or not an individual donated is perfectly evenly distributed.
zipconvert5 - This categorical variable shows that contains the highest level of individuals with higher ZIP codes (80000–99999) are more heavily sampled in this data set.
Next, a bivariate analysis was performed to investigate variable relationships. The results are as follows…
# Use a correlation plot to scan for relationships
corrplot(cor(num), method = "color")
# Create multiple uni variate plots
p1 <- data |> ggplot(aes(home_value, med_fam_inc)) + geom_point(color="darkblue") + theme_bw()
p2 <- data |> ggplot(aes(home_value, avg_fam_inc)) + geom_point(color="darkblue") + theme_bw()
p3 <- data |> ggplot(aes(avg_fam_inc, pct_lt15k)) + geom_point(color = "darkblue") + theme_bw()
p4 <- data |> ggplot(aes(num_child, time_lag)) + geom_point(color = "darkblue") + theme_bw()
# Add the plots together
(p1 | p2) / (p3 | p4) +
plot_annotation(title = "Relationships Amongst Predictors")
# Create multiple uni variate plots
p1 <- data |> ggplot(aes(num_child, last_gift)) + geom_point(color="darkblue") + theme_bw()
p2 <- data |> ggplot(aes(home_value, wealth)) + geom_point(color="darkblue") + theme_bw()
p3 <- data |> ggplot(aes(last_gift, pct_lt15k)) + geom_point(color = "darkblue") + theme_bw()
p4 <- data |> ggplot(aes(num_child, time_lag)) + geom_point(color = "darkblue") + theme_bw()
# Add the plots together
(p1 | p2) / (p3 | p4) +
plot_annotation(title = "Relationships Amongst Predictors")
Observations: The correlation plot contained evidence that there was significant correlation between income, home, value, pct_lt15k, numbers of children, and gifts so I started investigating these.
There is a linear relationship between home value and median/average family income showing that the more expensive a home was the more a families median/average income was.
There is a downward trending relationship on pct_lt15k and avg_fam_inc showing that the higher a families average income the less likely they were to live next to someone earning less than $15,000 annually.
Time_lag decreased as number of children decreased indicating that the less children a family has, the more likely they are to give gifts.
The response variable that I will be using to indicate whether a
person is a donor or non-donor is target. Prior to building
the models to begin predictions the data set must be transformed.
To begin the variable transformations, the response variable is converted into a factor variable containing two levels “Non_Donor” and “Donor.” Further, Binary transformations were performed on the categorical variables within the data. For example, instead of the Yes/No condition used in the raw data for zipconvert1-5, I used 1 to indicate yes and 0 to indicate no. This change is critical for running classification models because R must have factor inputs rather than textual. Lastly, the data was scaled by using a standardization method in R to ensure that the results aren’t skewed. To demonstrate the importance of scaling, consider a variable called income on a scale of 100,000 - 500,000 and wealth is on a scale of 1 - 10. Due to the vast difference in evenness between the two scales the model could be thrown off balance, hence why scaling is needed prior to running the model.
# Target is the binary response variable
y_train <- factor(data$target, labels = c("Non_Donor", "Donor"))
# Convert factor variables into dummy variables
x_train <- model.matrix(~ . - 1, data = data |> select(-target)) |> as.data.frame()
x_future <- model.matrix(~ . - 1, data = future_data) |> as.data.frame()
# Ensure columns match between future and training data
x_future <- x_future[, names(x_train)]
# Standardize Predictors
preproc <- preProcess(x_train, method = c("center", "scale"))
# Apply Pre-processing steps to the main data
train_processed <- cbind(
predict(preproc, x_train),
target = y_train
)
# Apply these steps to the future data
future_processed <- predict(preproc, x_future)
K-fold cross validation was used to ensure a clean test and training split was prior to modeling. Since the value of K that was chosen is 10 the cross-validation method is effectively splits the 3000 records into 10 folds consisting of 300 records each. The first fold might be the test set and the other 9 folds might be the train set demonstrated in a pattern as follows: test,train,train,train,train,train,train,train,train,train. In the next run, it might look something like train,test,train,train,train,train,train,train,train,train continuing to iterate through the folds until each fold has been tested on. 10-fold CV then takes the average accuracy rate of the 10 folds and averages them to give a performance score on the model.
# Set the seed to ensure results are reproducible
set.seed(12345)
# Use 10-fold cross validation on the data
fit_control <- trainControl(
method = "cv",
number = 10,
classProbs = TRUE,
summaryFunction = twoClassSummary,
savePredictions = "final"
)
The three models chosen to predict donors vs non-donors were Random Forest, Linear Discriminant Analysis, and Gradient Boosting. The success metric that was chosen was the ROC AUC metric to determine which model was the most effective at making predictions. The ROC AUC (area under the curve) metric is used to determine how good of a job the model does at predicting donors vs non-donors. ROC AUC can be compared to an extremely detailed confusion matrix that tracks how well the model classifies its predictions, instead of evaluating the model at one specific threshold like a confusion matrix, a ROC evaluates it at every possible threshold. Generally, a ROC AUC of .50 is practically guessing whereas a ROC of .90 is doing an excellent job of predicting. Anything above a .5 ROC AUC is better than random guessing.
Prior to fitting these three models several other classification models such as Quadratic Discriminant Analysis, K-nearest neighbors, Logistic Regression, and Naive Bayes were used, all of which underperformed the models shown. To fit these three models the caret package in R is utilized to re sample the test/train data, fit, and evaluate the models. Although these models are among some of the best classification models there are some notable limitations as follows.
Random Forest: Random Forest does a great job creating predictions that are extremely diversified due to how many trees it reliably fits. However, one weakness of RF is that the correct value of M (one of Random Forest’s hyperparameters) must be selected - often referred as ‘mtry’ in this section. If the correct M value is not selected than the performance of the model can decrease (James et al - ISLR, 2021).
Linear Discriminant Analysis: LDA is known for working extremely well on models that have normally distributed predictors often found when the predictors form a bell shaped curve. However, a weakness of LDA is if normal distribution is not the case, than LDA still relies on the normal distribution and linear boundary assumptions resulting in less flexibility in the model. (James et al - ISLR, 2021).
Gradient Boosting: Gradient Boosting models can achieve higher performance than RF due to its unique tuning options. However, a weakness of Gradient Boosting is its hyperparameters are more complex. Much like RF, the performance of the GB model is dependent on how the hyperparameters or tuned.
Description: The first classification model used is Random Forest. With the combined CV resampling work above the fit_control object is called in the trControl argument of the RF model. This Random Forest model will create 300 decision trees where each tree is trained on a different sample of the data. To diversify the trees individual nodes the mtry function is used to tell the model to consider different splits of tree sets at 1-split, 2-splits, and 3-splits. This diversity will allow the model more flexibility to find the best predictors that influence the response variable and showcases one of RF strengths.
# Fit the Random Forest
rf_fit <- train(
target ~ .,
data = train_processed,
method = "rf",
metric = "ROC",
trControl = fit_control, # use the CV object created above
tuneGrid = expand.grid(mtry = c(1, 2, 3)),
ntree = 300
)
# View the model performance
datatable(round(rf_fit$results[, c("mtry",
"ROC",
"Sens",
"Spec")], 3) ,
caption = "Random Forest Performance Metrics"
)
# View the RF importance plot
plot(varImp(rf_fit))
# View the RF ROC AUC plot
plot(rf_fit,
main = "What Value of mtry Gives the Highest ROC?",
xlab = "Random Selected Predictors (mtry)")
# View the model performance
datatable(round(rf_fit$results[, c("mtry",
"ROC",
"Sens",
"Spec")], 3) ,
caption = "Random Forest Performance Metrics",
options = list(dom = "t")
)
Observations: After fitting the random forest the highest value of ROC was achieved on mtry=2 with a ROC value of 0.5686 (shown in the graph below). The importance plot shows that home_value, avg_gift, and med_fam_inc are the highest predictors in whether or not a person will donate. For the target hyper parameter of the Random Forest, several variable combinations were attempted the underperformed using all of the variables together. This demonstrates that just because a variable is low on the importance plot does not mean it is useless in making predictions. Additionally, the number of trees and M value were testing using several different variations which all underperformed the final output that is shown above.
The RF model performance can be evaluated by looking at the table above. The mtry, ROC, Sens, and Spec columns of the results table give a comprehensive breakdown explained as follows: First, the mtry column shows the three variations of trees that were selected to be a part of the model. The ROC values displayed are the average ROC result of the 10-folds - mtry=2 appears to be the highest with a value of 0.5686354. Next, sens and spec values at mtry=2 show that the model correctly identifies about 56.6% of the actual donors and 52.7% of the non-donors.
Model Assessment - see the model comparison section below for the full assessment
Description: For the second classification model a Linear Discriminant Analysis model is fit using the same workflow that used on the Random Forest. A key difference to point out in Random Forest and LDA is that the RF model fit hundreds of decision trees to make its predictions while LDA fits a linear boundary between donors and non-donors offering a much less flexible approach. Another distinction between the two models is their variable chosen shown in the importance plots.
# Fit the LDA
lda_fit <- train(
target ~ .,
data = train_processed,
method = "lda",
metric = "ROC",
trControl = fit_control # use the CV object created above
)
# View the LDA importance plot
plot(varImp(lda_fit))
# View the model performance
datatable(round(lda_fit$results[, c("ROC",
"Sens",
"Spec")], 3) ,
caption = "LDA Performance Metrics",
options = list(dom = "t")
)
Observations: After fitting the LDA model the best ROC value was a ROC value of 0.578. The importance plot shows that months_since_donate, last_gift, and avg_gift are the highest predictors in whether or not a person will donate. Notice that the important variables are different than the RF indicating the differences in the LDA models preference. LDA also does not have an “mtry” argument due to its more rigid model parameters so only one average ROC value is outputted
Model Assessment - see the model comparison section below for the full assessment
Description: For the third classification model, a Gradient Boosting model is fit. Much like Random Forest, Gradient Boosting is a tree-based classification model that uses a non-linear decision boundary. Unlike the Random Forest, Gradient Boosting builds its tree sequentially designing each tree to improve upon the error within the previous trees. This sequential improvement allows Gradient Boosting to learn complex relationships very well but can raise the risk of the model over fitting if it is not tuned properly (GeeksforGeeks, 2024).
# Fit the gradient boosting model
gbm_fit <- train(
target ~ .,
data = train_processed,
method = "gbm",
metric = "ROC",
trControl = fit_control, # use the CV object created above
verbose = FALSE
)
# View the GB importance plot
plot(varImp(gbm_fit))
# View the model performance
datatable(round(gbm_fit$results[, c("ROC",
"Sens",
"Spec",
"n.trees",
"interaction.depth",
"shrinkage")], 4) ,
caption = "GBM Performance Metrics",
options = list(dom = "t")
)
Observations: The Gradient Boosting model performed the best of the three models with a ROC of .5895. Unlike the other models, the importance plot ruled out several of the predictor variables as insignificant in its predictions showing a little more rigidity than RF. Uniquely, its top three predictors were months_since_donate, largest_gift, and avg_gift. Due to Gradient Boosting unique style of fitting sequences of trees, the models output table is more extensive. The highest ROC occurred with 50 trees with an interaction depth of 1 indicating that a less complex combination of tree variables at 50 trees won over the more complex options at 100 and 150 trees.
Model Assessment - see the model comparison section below for the full assessment
# Create objects for each of the models
rf <- rf_fit
lda <- lda_fit
gb <- gbm_fit
# Create a table that shows all three models side-by-side
datatable(
data.frame(
Model = c("Random Forest",
"Linear Discriminant Analysis",
"Gradient Boosting"),
ROC = round(c(
max(rf_fit$results$ROC),
max(lda_fit$results$ROC),
max(gbm_fit$results$ROC)
), 4),
Sens = round(c(
rf_fit$results$Sens[which.max(rf_fit$results$ROC)],
lda_fit$results$Sens[which.max(lda_fit$results$ROC)],
gbm_fit$results$Sens[which.max(gbm_fit$results$ROC)]
), 3),
Spec = round(c(
rf_fit$results$Spec[which.max(rf_fit$results$ROC)],
lda_fit$results$Spec[which.max(lda_fit$results$ROC)],
gbm_fit$results$Spec[which.max(gbm_fit$results$ROC)]
), 3)
),
caption = "Top 3 Performing Models",
options = list(dom = "t")
)
RF Model Analysis: The Random Forest model performed the worst of the three models with a ROC of 0.5686. The Sens and Spec scores show that the model correctly identifies about 56.6% of the actual donors and 52.7% of the non-donors.
LDA Model Analysis: The Linear Discriminant Analysis model was the second best of the three models with a ROC of 0.5783. The Sens and Spec scores show that the model correctly identifies about 57.4% of the actual donors and 54.2% of the non-donors.
GBM Model Analysis: The Linear Discriminant Analysis model was the second best of the three models with a ROC of 0.5859. The Sens and Spec scores show that the model correctly identifies about 53.6% of the actual donors and 58.6% of the non-donors.
Lastly, taking the most effective model predictions are made to predict future donors…
# Future Donor Predictions
preds <- predict(gbm_fit,
newdata = future_processed,
type="prob")
# Apply a Cutoff
cutoff <- 0.32
results <- data.frame(
donor_probability = preds$Donor,
Prediction = ifelse(
preds$Donor >= cutoff,
"Donor",
"Non_Donor"
)
)
# View Results
results |> head(120)
# Write output for the CSV
submission <- data.frame(
value = results$Prediction
)
# Write to a CSV
write.csv(submission, "future_fundraising_predictions.csv", row.names = FALSE)
cat("Predictions exported to future_fundraising_predictions.csv\n")
## Predictions exported to future_fundraising_predictions.csv
Conclusion: With the 32% threshold, if there is a chance that someone might donate that is larger than 32% they will be marked as a donor, if not then they will be marked as a non donor. Evidently, there is already a large number non-donors that the model predicted that prove useful in helping the veterans organization’s marketing team maximize their profit margins.
James, G., Witten, D., Hastie, T., Tibshirani, R., & Taylor, J. (2021). An introduction to statistical learning: With applications in R (2nd ed.). Springer. https://doi.org/10.1007/978-1-0716-1418-1
GeeksforGeeks. (2024, April 9). Gradient boosting vs random forest. https://www.geeksforgeeks.org/machine-learning/gradient-boosting-vs-random-forest/