Targeting Housing Subsidy with People-Based Machine Learning

Version 1.0 | First Created Oct 22, 2023 | Updated Nov 1, 2023

Keywords: Logistic Regression, Confusion Matrix, ROC Curve, Sensitivity, Specificity, Cost-Benefit Analysis, McFadden R-Squared

Introduction

People-based machine learning, a field harnessing the power of data-driven insights into human behaviors and preferences, has become pivotal for businesses and organizations. With algorithms trained on human-centric data, businesses get to personalize products, services, and experiences based on individual customer preferences. At any point, the bottom line of any businesses and organizations may be affected, if a client no longer wish to participate, donate, or purchase a good, which makes it extra important to have a predictive model that captures these potential changes. Like any machine learning applications, however, it is difficult to capture the nuances in human behavior with a single model, nor is it possible to have a unbiased sample dataset that truly reflect human behavior. This means that continuous model adaptation and recognition of these complexities are essential for effective decision-making based on these models.

The goal of this report is to fit a logistic regression model to predict whether homeowners in Emil City will enter and take a home repair tax credit. To contextualize, the city is considering a more proactive approach for targeting homeowners who would quality for this tax credit program and the Department of Housing and Community Development (HCD) have been trying to proactively reach out to eligible homeowners ever year. It is predicted that houses that transacted after taking the credit, will sell with a $10,000 premium, on average. Homes surrounding the repaired home see an aggregate premium of $56,000, on average. The thorny task for the HCD is to target homeowners with the highest likelihood of accepting the credit. In this case, they could maximize the benefits of the aggregate premium of houses while reduce the marketing resources cost. As of now, typically not all eligible homeowners they reach out and enter the program ultimately takes the credit.

That said, we will need to train a better model to ensure that the HCD will not be reaching out to eligible homeowners at random. Logistic regression is used here due to its interpretability and suitability for binary classification problems (in this case, our dependent variable indicate whether the a homeowner will take the credit or not). We will first start with data exploration to determine variables most explains homeowner’s behavior. Then, we split the data into training and testing set, fit a regression model, and constructed confusion matrix at 50% threshold to see the extent to which our model can correctly predict homeowner entering the credit program (true positive), or correctly predict homeowner not taking the credit (true negative). This is followed by a cross validation test to compute ROC, sensitivity, and specificity. We conclude our analysis with a cost/benefit analysis and determine the optimal threshold of delineating positive and negative prediction that would maximize the revenue.

The GitHub repository for this report is here.

Variables Examination

The first part of this report focuses on data exploration to see how each factor in our dataset relates to homeowner’s decision to accept home repair tax credit. Variables are split into categorical variables and numeric variables respectively, for which we conducted chi-square and anova tests to examine if there’s a significant association between dependent variable and independent variables.

Categorical Variables

The figure below shows the difference in the number of homeowners who accepted and did not accept tax credit based on a different categorical features. We observed the following. There seems to be a significant difference in the number of people accepting credit who were contacted via cellular, who had a university degree, work as admin or technician, were married, were previously contacted in May, were contacted recently, were Philadelphia residents, and had previously accepted the credit. On the other hand, there doesn’t seem to be a significant difference in the number of people accepting credit depending on the day of the week they were contacted, their mortgage status, and presence of tax lien.

cat_var <- housing %>% 
  select(-c(X, age, y_numeric, previous, unemploy_rate, cons.price.idx, cons.conf.idx, inflation_rate, spent_on_repairs, campaign))
cat_var %>%
    gather(Variable, value, -y) %>%
    count(Variable, value, y) %>%
      ggplot(., aes(value, n, fill = y)) +   
        geom_bar(position = "dodge", stat="identity") +
        facet_wrap(~Variable, scales="free", ncol = 4, 
                   labeller= labeller(Variable = c(
                     `contact` = "Ways of Contact",
                     `day_of_week` = "Day of the Week",
                     `education` = "Educational Attainment",
                     `pdays` = "Days After Contact",
                     `job` = "Occupation Indicator",
                     `marital` = "Marital Status",
                     `month` = "Month of Last Contact",
                     `mortgage` = "Mortgage",
                     `poutcome` = "Previous Campaign Outcome",
                     `taxbill_in_phl` = "Philadelphia Residents",
                     `taxLien` = "Lien"))) +
        scale_fill_manual(values = palette2, name = "Credit") +
        labs(x="Credit", y="Count",
             title = "Feature Associations with the Likelihood of Credit Acceptance",
             subtitle = "Categorical features") +
        theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  theme(plot.subtitle = element_text(size = 12,face = "italic"),
        plot.title = element_text(size = 18, face = "bold"), 
        axis.text.x=element_text(size=8),
        axis.text.y=element_text(size=8), 
        axis.title=element_text(size=9), 
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, linewidth =0.8))

Using pure visual evidence is not enough to explain the significance of our categorical variables in predicting whether homeowners will take the credit. To support the above visualization, we conducted chi-square test iteratively for all categorical variables. The chi-square test determines if there is a significant association between categorical variables. It compares the observed frequencies of different categories in a contingency table with the frequencies that would be expected if the variables were independent. The table below shows that Ways of Contact, Edcucation, Marital, Occupation, Month of Last Contact, Previoust Campaign Outcome, Lien Presence and Days After Last Contact are significant predictors in determining homeowner’s likelihood of accepting the credit.

chi_var <- c("contact", "day_of_week", "education", "job", "marital", "month", "mortgage", "poutcome", "taxbill_in_phl", "taxLien", "pdays")

Label = c("contact" = "Ways of Contact",
          "day_of_week" = "Day of the Week",
          "education" = "Educational Attainment",
          "job" = "Occupation Indicator",
          "marital" = "Marital Status",
          "month" = "Month of Last Contact",
          "mortgage" = "Mortgage",
          "poutcome" = "Previous Campaign Outcome",
          "taxbill_in_phl" = "Philadelphia Residents",
          "taxLien" = "Lien", 
          "pdays" = "Days After Contact")

chi_square_results <- data.frame(
  Df = integer(),
  X_Squared = numeric(),
  P_Value = numeric(),
  stringsAsFactors = FALSE
)

for (var in chi_var) {
  contingency_table <- table(cat_var$y, cat_var[[var]])
  chi_square <- chisq.test(contingency_table)
  chi_square_results <- rbind(chi_square_results, data.frame(
    Df = chi_square$parameter,
    X_Squared = chi_square$statistic,
    P_Value = chi_square$p.value
  ))
}

rownames(chi_square_results) <- Label


chi_square_results %>% 
 kable() %>% 
 kable_styling( bootstrap_options = c("striped", "hover", "condensed")) %>% 
    footnote(general_title = "\n", general = "Table 1")

	Df	X_Squared	P_Value
Ways of Contact	1	76.8462374	0.0000000
Day of the Week	4	0.5121761	0.9723049
Educational Attainment	7	22.2918199	0.0022621
Occupation Indicator	11	69.9786737	0.0000000
Marital Status	3	10.2859188	0.0162857
Month of Last Contact	9	299.7881167	0.0000000
Mortgage	2	0.6273842	0.7307440
Previous Campaign Outcome	2	454.4865243	0.0000000
Philadelphia Residents	1	0.7595710	0.3834628
Lien	2	24.1856637	0.0000056
Days After Contact	1	448.2212094	0.0000000
Table 1

Numeric Variables

The figure below shows the difference in the number of homeowners who accepted and did not accept tax credit based on a different numeric features. There seems to be a significant difference in the number of people accepting credit based on their age, number of contacts for this campaign and before this campaign, as well as inflation and unemployment rate of that time.

numeric_var <- housing %>% 
  select(age, y, previous, unemploy_rate, cons.price.idx, cons.conf.idx, inflation_rate, spent_on_repairs, campaign)


numeric_var %>%
  gather(Variable, value, -y) %>%
    ggplot(aes(y, value, fill=y)) + 
      geom_bar(position = "dodge", stat = "summary", fun = "mean") + 
      facet_wrap(~Variable, scales = "free", ncol = 4, labeller= labeller(Variable = c(
                     `age` = "Age",
                     `previous` = "Contact before Campaign",
                     `unemploy_rate` = "Unemployment Rate",
                     `cons.price.idx` = "Consumer Price Index",
                     `cons.conf.idx` = "Consumer Confidence Index",
                     `campaign` = "Contacts for Campaign",
                     `inflation_rate` = "Inflation Rate",
                     `spent_on_repairs` = "Amount Spent on Repairs"))) +
      scale_fill_manual(values = palette2) +
      labs(x="Credit", y="Value", 
           title = "Feature Associations with the Likelihood of Credit Acceptance",
           subtitle = "Continous features") +
      theme(legend.position = "none") +
    theme(plot.subtitle = element_text(size = 12,face = "italic"),
        plot.title = element_text(size = 18, face = "bold"), 
        axis.text.x=element_text(size=10),
        axis.text.y=element_text(size=10), 
        axis.title=element_text(size=9), 
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, linewidth =0.8))

Because these numeric features all have different scales and are difficult to compare across, we made another density plot, which conveys the same information, but shows distribution of a continuous variable. Surprisingly, we see that the distribution of Consumer Confidence Index, Consumer Price Index, and Amount Spent on Repairs all significantly differ between homeowners who accepted and did not accept the tax credit.

numeric_var %>%
    gather(Variable, value, -y) %>%
    ggplot() + 
    geom_density(aes(value, color=y), fill = "transparent") + 
   facet_wrap(~Variable, scales = "free", ncol = 4, labeller= labeller(Variable = c(
                     `age` = "Age",
                     `previous` = "Contact before Campaign",
                     `unemploy_rate` = "Unemployment Rate",
                     `cons.price.idx` = "Consumer Price Index",
                     `cons.conf.idx` = "Consumer Confidence Index",
                     `inflation_rate` = "Inflation Rate",
                     `campaign` = "Contacts for Campaign",
                     `spent_on_repairs` = "Amount Spent on Repairs")))+
  scale_color_manual(values = palette1) +
    labs(title = "Feature Distributions Based on Credit Acceptance",
         subtitle = "Continous features") +
  theme(plot.subtitle = element_text(size = 12,face = "italic"),
        plot.title = element_text(size = 18, face = "bold"), 
        axis.text.x=element_text(size=9),
        axis.text.y=element_text(size=9), 
        axis.title=element_text(size=9), 
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, linewidth =0.8))

We further conducted ANOVA test for all numeric features to see if there’s a significant difference in mean age, for example, across people who take and did not take the credit. The result is in concur with our expectations: all numeric variables are significant predictors.

anova_var <- c("age", "previous", "unemploy_rate", "cons.price.idx", "cons.conf.idx", "inflation_rate", "spent_on_repairs", "campaign")

new_names <- c("age" = "Age",
                "previous" = "Contact before Campaign",
                "unemploy_rate" = "Unemployment Rate",
                "cons.price.idx" = "Consumer Price Index",
                "cons.conf.idx" = "Consumer Confidence Index",
                "inflation_rate" = "Inflation Rate",
                "spent_on_repairs" = "Amount Spent on Repairs",
               "campaign" = "Contacts for Campaign")

anova_results <- data.frame(
                            Df = integer(), 
                            Sum_Sq = numeric(), 
                            Mean_Sq = numeric(), 
                            F_value = numeric(), 
                            P_Value = numeric(), stringsAsFactors = FALSE)

for (var in anova_var) {
  anova_result <- aov(numeric_var[[var]] ~ numeric_var$y)
  summary_data <- summary(anova_result)[[1]][, c("Df", "Sum Sq", "Mean Sq", "F value", "Pr(>F)")][1:1, ]
  anova_results <- rbind(anova_results, summary_data)
}


rownames(anova_results) <- new_names


anova_results %>% 
 kable() %>% 
 kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>% 
    footnote(general_title = "\n", general = "Table 2")

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Age	1	1596.56982	1596.56982	15.06149	0.0001057
Contact before Campaign	1	79.03066	79.03066	288.00244	0.0000000
Unemployment Rate	1	807.05379	807.05379	359.02715	0.0000000
Consumer Price Index	1	13.36299	13.36299	40.19177	0.0000000
Consumer Confidence Index	1	257.19577	257.19577	12.21670	0.0004786
Inflation Rate	1	1103.20812	1103.20812	402.90944	0.0000000
Amount Spent on Repairs	1	2725800.03197	2725800.03197	571.90285	0.0000000
Contacts for Campaign	1	157.25175	157.25175	23.97551	0.0000010
Table 2

Logistic Regression

We then fit a logistic regression containing all variables in our dataset. To do so, we split our data randomly into 65% of training dataset and 35% of testing dataset. Some of our variables, such as Lien and Education, have single occurrence and therefore needs to be put into the training set to avoid new levels appearing in the testing set. This model shows that some indicators are significantly correlated with credit acceptance with p < 0.05. This include age, some categories of job, ways of contact, month last contacted, previous campaign outcome, unemployment rate, consumer price and confidence indices. However, the majority of predictors are not significant.

set.seed(479)
trainIndex <- createDataPartition(y = paste(housing$taxLien, housing$education) , p = .65,
                                  list = FALSE,
                                  times = 1)
housingTrain <- housing[ trainIndex,]
housingTest  <- housing[-trainIndex,]

housingModel <- glm(y_numeric ~ .,
                  data=housingTrain %>% 
                    dplyr::select(age, job, marital, education, taxLien, mortgage, taxbill_in_phl, contact, month, day_of_week, campaign, pdays, poutcome, unemploy_rate, cons.conf.idx, cons.price.idx, inflation_rate, spent_on_repairs, y_numeric), family="binomial" (link="logit"))

housing_sum <- summary(housingModel)
coefficients_table <- as.data.frame(housing_sum$coefficients)

coefficients_table$significance <- ifelse(coefficients_table$`Pr(>|z|)` < 0.001, '***',
                                         ifelse(coefficients_table$`Pr(>|z|)` < 0.01, '**',
                                                ifelse(coefficients_table$`Pr(>|z|)` < 0.05, '*',
                                                       ifelse(coefficients_table$`Pr(>|z|)` < 0.1, '.', ''))))

coefficients_table$p_value <- paste0(round(coefficients_table$`Pr(>|z|)`, digits = 3), coefficients_table$significance)

coefficients_table %>%
  select(-significance, -`Pr(>|z|)`) %>% 
  kable(align = "r") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))  %>%
  footnote(general_title = "\n", general = "Table 3")

	Estimate	Std. Error	z value	p_value
(Intercept)	-335.1496539	132.8764110	-2.5222660	0.012*
age	0.0204678	0.0086973	2.3533692	0.019*
jobblue-collar	-0.5325684	0.2907044	-1.8319927	0.067.
jobentrepreneur	-0.7143944	0.5111034	-1.3977491	0.162
jobhousemaid	-0.3226287	0.5351631	-0.6028605	0.547
jobmanagement	-0.6072725	0.3018896	-2.0115715	0.044*
jobretired	-0.5592430	0.3903821	-1.4325527	0.152
jobself-employed	-0.8229489	0.4298447	-1.9145261	0.056.
jobservices	-0.3355901	0.3051141	-1.0998838	0.271
jobstudent	0.1281199	0.4256531	0.3009960	0.763
jobtechnician	0.0225732	0.2395056	0.0942491	0.925
jobunemployed	0.0109616	0.4385126	0.0249973	0.98
jobunknown	-0.3124032	0.7707176	-0.4053407	0.685
maritalmarried	0.1371049	0.2520088	0.5440481	0.586
maritalsingle	0.2632753	0.2887379	0.9118142	0.362
maritalunknown	-13.8144973	556.2455938	-0.0248352	0.98
educationbasic.6y	0.1448773	0.4419572	0.3278083	0.743
educationbasic.9y	0.1232291	0.3452445	0.3569327	0.721
educationhigh.school	-0.0210035	0.3362686	-0.0624604	0.95
educationilliterate	-14.0457400	1455.3976440	-0.0096508	0.992
educationprofessional.course	-0.0864436	0.3719707	-0.2323935	0.816
educationuniversity.degree	0.2385109	0.3323036	0.7177499	0.473
educationunknown	0.2113254	0.4195492	0.5036964	0.614
taxLienunknown	-0.0509971	0.2272948	-0.2243655	0.822
taxLienyes	-12.2808134	1455.3976455	-0.0084381	0.993
mortgageunknown	0.0597981	0.5034481	0.1187771	0.905
mortgageyes	-0.2120102	0.1481208	-1.4313331	0.152
taxbill_in_phlyes	0.1796557	0.2094041	0.8579381	0.391
contacttelephone	-0.9822124	0.3039599	-3.2313878	0.001**
monthaug	0.4298805	0.4366130	0.9845802	0.325
monthdec	1.8238219	0.7846121	2.3244886	0.02*
monthjul	0.1025432	0.3828448	0.2678453	0.789
monthjun	-0.1581236	0.4690165	-0.3371387	0.736
monthmar	1.5457954	0.5865365	2.6354634	0.008**
monthmay	-0.2394845	0.3070199	-0.7800292	0.435
monthnov	-0.0057614	0.4441257	-0.0129726	0.99
monthoct	0.4016013	0.5628718	0.7134863	0.476
monthsep	0.6429439	0.6474800	0.9929943	0.321
day_of_weekmon	-0.1183577	0.2296430	-0.5153985	0.606
day_of_weekthu	0.0174090	0.2259178	0.0770591	0.939
day_of_weektue	0.0213120	0.2298068	0.0927387	0.926
day_of_weekwed	-0.0203213	0.2397434	-0.0847626	0.932
campaign	-0.0598325	0.0404233	-1.4801468	0.139
pdaysYes	0.6436547	0.6282873	1.0244592	0.306
poutcomenonexistent	0.5081719	0.2280928	2.2279169	0.026*
poutcomesuccess	1.1689464	0.6464320	1.8083054	0.071.
unemploy_rate	-1.2113721	0.4927379	-2.4584514	0.014*
cons.conf.idx	0.0633293	0.0296774	2.1339226	0.033*
cons.price.idx	2.5636069	0.8719343	2.9401378	0.003**
inflation_rate	-0.6576623	0.4702556	-1.3985209	0.162
spent_on_repairs	0.0187400	0.0108730	1.7235321	0.085.
Table 3

The model is being tested on the test set and a mosaic plot is created to visualize the confusion matrix. Unfortunately, the model has low sensitivity rate of only 0.23. Sensitivity is the true-positive rate, which is the proportion of actual positive cases that were correctly identified by the model. Since our model’s accuracy rate is high (0.89), this means that our model missing a significant number of positive cases. These are cases where we correctly predict that homeowners will accept the credit and they actually did.

testProbs <- data.frame(Outcome = as.factor(housingTest$y_numeric),
                        Probs = predict(housingModel, housingTest, type= "response"))

testProbs <- 
  testProbs %>%
  mutate(predOutcome  = as.factor(ifelse(testProbs$Probs > 0.5 , 1, 0)))  

 
mosaicplot(confusionMatrix(testProbs$predOutcome, testProbs$Outcome, 
                       positive = "1")$table, color=c("#DBC2CF","#9FA2B2"), main = "Mosaic Plot for Original Confusion Matrix",
           xlab = "Prediction", ylab = "Reference")

housingmatrix <- as.data.frame(as.table(confusionMatrix(testProbs$predOutcome, testProbs$Outcome, 
                       positive = "1")$table))
housingmatrix %>% 
  kable() %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))  %>%
  footnote(general_title = "\n", general = "Table 4")

Prediction	Reference	Freq
0	0	1251
1	0	22
0	1	122
1	1	37
Table 4

Improving Model

However, just because a variable is significant at predicting credit acceptance does not mean that it is a relevant or meaningful predictor. To improve this model, not only do we need to perform feature engineering to increase sensitivity while maintain overall accuracy, but also do we need to exclude irrelevant features. To fit a new model:

• Firstly, we collapsed the categories in education and employment. For education, considering that having a university degree serves a benchmark in determining credit acceptance, we recode it into two categories: highschool.below and highschool.above. For employment, we recode it into four major categories: unemployed, white-collar, blue-collar, and self-employed, considering that employment status and employment sector will both affect credit acceptance.

• Secondly, we create interaction features between education and employment considering that different job may require different levels or type of education.

• Thirdly, we computed a new variable called campaign efficiency, which is the ratio of previous campaign successes to the number of contacts made, This will indicate the efficiency of previous campaigns in terms of successful credit acceptance. We also computed a composite feature representing economic conditions based on unemploy_rate, cons.price.idx, cons.conf.idx, and inflation_rate.

• Fourthly, carefully reviewing the dataset leads us to realize that it is imbalance. Usually, imbalanced datasets can lead the model to be biased towards the majority class, causing it to perform poorly on the minority class. In our case, a significant number of attributes we have are cases when homeowner declined tax credit. This makes feature weighting extra important. By assigning appropriate weights to features, we can ensure that the model gives sufficient importance to the minority class, helping it learn the underlying patterns in both classes more effectively.

• Finally, we excluded variable Month we last contacted individual, Day of the week we last contacted individual, Lien and marital, not merely because these variables are not significant predictors, but because the former two features won’t be relevant to credit acceptance, variable lien contains too many unknown values, and there might be colinearity between martial, age, job, and education.

housing <- housing %>% 
  mutate(edu = case_when(
    education == "high.school" | education == "professional.course" | education == "university.degree" ~ "highschool.above", 
    TRUE ~ "highschool.below"),
    employment = case_when(
      job == "self-employed" ~ "self-employed", 
      job == "retired" | job == "unemployed" | job == "student" | job == "unknown" | job == "housemaid" ~ "not-employed",
      job == "blue-collar" ~ "blue-collar", 
      TRUE ~ "white-collar"))
  
housing$job_education_interaction <- paste(housing$edu, housing$employment, sep="_")
housing$campaign_efficiency <- housing$previous / housing$campaign
housing$economic_conditions <- (housing$unemploy_rate + housing$cons.price.idx + housing$cons.conf.idx - housing$inflation_rate) / 4
imbalance_weight <- sum(housing$y_numeric == 0) / sum(housing$y_numeric == 1)

Our improved model shows that most indicators are significantly correlated with credit acceptance with p < 0.05, except for some categories of job-education interaction as well as pdays, which is whether an individual has been contacted before. This is probably because pday is correlated with poutcome, another variable indicating success of previous contact.

set.seed(479)
trainIndex <- createDataPartition(housing$y , p = .65,
                                  list = FALSE,
                                  times = 1)
housingTrain <- housing[ trainIndex,]
housingTest  <- housing[-trainIndex,]

improvedModel <- glm(y_numeric ~ .,
                  data=housingTrain %>% 
                    dplyr::select(y_numeric, age, job_education_interaction, mortgage, pdays, contact, poutcome, economic_conditions, previous, campaign_efficiency),  weights = ifelse(y_numeric == 1, imbalance_weight, 1), 
                  family="binomial" (link="logit"))

improved_sum <- summary(improvedModel)
coefficients_table <- as.data.frame(improved_sum$coefficients)

coefficients_table$significance <- ifelse(coefficients_table$`Pr(>|z|)` < 0.001, '***',
                                         ifelse(coefficients_table$`Pr(>|z|)` < 0.01, '**',
                                                ifelse(coefficients_table$`Pr(>|z|)` < 0.05, '*',
                                                       ifelse(coefficients_table$`Pr(>|z|)` < 0.1, '.', ''))))

coefficients_table$p_value <- paste0(round(coefficients_table$`Pr(>|z|)`, digits = 3), coefficients_table$significance)

coefficients_table %>%
  select(-significance, -`Pr(>|z|)`) %>% 
  kable(align = "r") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))  %>%
  footnote(general_title = "\n", general = "Table 5")

	Estimate	Std. Error	z value	p_value
(Intercept)	-2.4483602	0.4237519	-5.7778149	0***
age	0.0067819	0.0029976	2.2624372	0.024*
job_education_interactionhighschool.above_not-employed	0.5082476	0.2026669	2.5077973	0.012*
job_education_interactionhighschool.above_self-employed	-0.0568971	0.2457110	-0.2315608	0.817
job_education_interactionhighschool.above_white-collar	0.0943712	0.1708588	0.5523346	0.581
job_education_interactionhighschool.below_blue-collar	-0.1894807	0.1833597	-1.0333825	0.301
job_education_interactionhighschool.below_not-employed	0.6054083	0.2088516	2.8987492	0.004**
job_education_interactionhighschool.below_self-employed	-0.8584100	0.4391580	-1.9546725	0.051.
job_education_interactionhighschool.below_white-collar	-0.1024480	0.2011602	-0.5092856	0.611
mortgageunknown	-0.4819300	0.2161021	-2.2301038	0.026*
mortgageyes	-0.1803795	0.0638614	-2.8245484	0.005**
pdaysYes	0.7306544	0.4883692	1.4961108	0.135
contacttelephone	-0.8533536	0.0753342	-11.3275664	0***
poutcomenonexistent	0.8229103	0.2266384	3.6309397	0***
poutcomesuccess	1.3750342	0.5028598	2.7344284	0.006**
economic_conditions	0.1140043	0.0269205	4.2348516	0***
previous	1.1581184	0.2597144	4.4591994	0***
campaign_efficiency	-0.4983248	0.2125183	-2.3448563	0.019*
Table 5

We tested our model and a mosaic plot is created to visualize the confusion matrix. Now, the model has much higher sensitivity rate of 0.61. The accuracy rate dropped to 0.70, but is still acceptable. This means that now, our model is able to capture a significantly more number of positive cases.

testProbs <- data.frame(Outcome = as.factor(housingTest$y_numeric),
                        Probs = predict(improvedModel, housingTest, type= "response"))

testProbs <- 
  testProbs %>%
  mutate(predOutcome  = as.factor(ifelse(testProbs$Probs > 0.5 , 1, 0)))  


mosaicplot(confusionMatrix(testProbs$predOutcome, testProbs$Outcome, 
                       positive = "1")$table, color=c("#DBC2CF","#9FA2B2"), main = "Mosaic Plot for Improved Confusion Matrix",
           xlab = "Prediction", ylab = "Reference")

improvedmatrix <- as.data.frame(as.table(confusionMatrix(testProbs$predOutcome, testProbs$Outcome, 
                       positive = "1")$table))

improvedmatrix %>% 
  kable() %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))  %>%
  footnote(general_title = "\n", general = "Table 6")

Prediction	Reference	Freq
0	0	915
1	0	368
0	1	60
1	1	97
Table 6

We computed the Area Under the Receiver Operating Characteristic Curve (AUC-ROC) metric to further evaluate the performance of our improved model. The AUC value ranges from 0 to 1, where 0.5 means a model with no discriminatory power (it performs as well as random chance), and 1 means a perfect model (it perfectly distinguishes between positive and negative cases).Generally, an AUC between 0.7 and 0.8 is considered acceptable and indicates a fair to good performance of the model. While there is room for improvement, a score of 0.74 suggests that the model is capturing important patterns in the data and making predictions better than random guessing.

ggplot(testProbs, aes(d = as.numeric(Outcome), m = Probs)) +
  geom_roc(n.cuts = 50, labels = FALSE, colour = "#2E4756") +
    labs(title = "ROC Curve for Improved Model") +
  style_roc(theme = theme_grey) +
  geom_abline(slope = 1, intercept = 0, size = 1.5, color = "#DBC2CF") +
  theme(plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"), 
        axis.text.x=element_text(size=8),
        axis.text.y=element_text(size=8), 
        axis.title=element_text(size=9), 
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, linewidth =0.8))

Cross Validation

Next, we performed cross validation test for both of original model and improved model. The trainControl parameter is set to run 100 k-folds and to output predicted probabilities, classProbs, for two classes, those who accepted the credit and those who did not. Additional parameters output AUC (the train function refers to this as ‘ROC’) and confusion metrics for each fold. The three metrics in the cvFit output are for mean AUC, Sensitivity, and Specificity across all 100 folds. Here, the ROC for our original model and improved model are 0.76 and 0.70, the Sensitivity are 0.97 and 0.98, while the Specificity are 0.21 and 0.17 respectively.

ctrl <- trainControl(method = "cv", number = 100, classProbs=TRUE, summaryFunction=twoClassSummary)

cvFit_original <- train(y ~ .,
                  data=housing %>% 
                    dplyr::select(age, job, marital, education, taxLien, mortgage, taxbill_in_phl, contact, month, day_of_week, campaign, pdays, poutcome, unemploy_rate, cons.conf.idx, cons.price.idx, inflation_rate, spent_on_repairs, y), 
                method="glm", family="binomial",
   
                             metric="ROC", trControl = ctrl)
cvFit_improved <- train(y ~ .,
                  data=housing %>% 
                    dplyr::select(y, age, job_education_interaction, mortgage, pdays, contact, poutcome, economic_conditions, previous, campaign_efficiency),   weights = ifelse(y == "Yes", imbalance_weight, 1),
                method="glm", family="binomial",
                metric="ROC", trControl = ctrl)

The figure below plots the distribution of AUC, Sensitivity, and Specificity across the 100 folds. The tighter each distribution is to its mean, the more generalizable the model. Both models generalizes well to sensitivity – the rate it correctly predicts credit acceptance, instead of specificity. It seems our would-be decision-making tool is inconsistent in how it predicts the credit acceptance. This inconsistency will have a direct effect on the marketing and resource allocation process should this algorithm be put into production.

dplyr::select(cvFit_original$resample, -Resample) %>%
  gather(metric, value) %>%
  left_join(gather(cvFit_original$results[2:4], metric, mean)) %>%
  ggplot(aes(value)) + 
    geom_histogram(bins=20, fill = "#DBC2CF") +
    facet_wrap(~metric) +
    geom_vline(aes(xintercept = mean), colour = "#2E4756", linetype = 2, size = 0.6) +
    scale_x_continuous(limits = c(0, 1)) +
    labs(x="Goodness of Fit", y="Count", title="CV Goodness of Fit Metrics Original Model",
         subtitle = "Across-fold mean reprented as dotted lines") + 
   theme(plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"), 
        axis.text.x=element_text(size=8),
        axis.text.y=element_text(size=8), 
        axis.title=element_text(size=9), 
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, linewidth =0.8))

dplyr::select(cvFit_improved$resample, -Resample) %>%
  gather(metric, value) %>%
  left_join(gather(cvFit_improved$results[2:4], metric, mean)) %>%
  ggplot(aes(value)) + 
    geom_histogram(bins=20, fill = "#DBC2CF") +
    facet_wrap(~metric) +
    geom_vline(aes(xintercept = mean), colour = "#2E4756", linetype = 2, size = 0.8) +
    scale_x_continuous(limits = c(0, 1)) +
    labs(x="Goodness of Fit", y="Count", title="CV Goodness of Fit Metrics for Improved Model",
         subtitle = "Across-fold mean reprented as dotted lines") + 
  theme(plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"), 
        axis.text.x=element_text(size=8),
        axis.text.y=element_text(size=8), 
        axis.title=element_text(size=9), 
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, linewidth =0.8))

Cost Benefit Analysis

In this final section, we conducted a cost benefit analysis to help HCD maximize their revenue. If we predict that a household will take the credit, then HCD allocates (on average) $2,850 per homeowner which includes staff and resources to facilitate mailers, phone calls, and information/counseling sessions at the HCD offices. And among those eligible homeowners, 25% of them will take the $5000 credit. This wil lead to a $10,000 premium, on average, for each houses that transacted after taking the credit and an aggregate premium of $56,000, on average, for homes surrounding the repaired home. Based on this, we outlined the cost/benefit for each scenario in the confusion matrix, where Count refers to the number of occurrences of each scenario:

• True Positive: We predicted correctly homeowner would enter credit program; allocated the marketing resources, and 25% ultimately achieved the credit: \(-2850 * Count - 5000 * Count * 0.25 + (10000 + 56000) * Count\)

• True Negative: We predicted correctly homeowner would not take the credit, no marketing resources were allocated, and no credit was allocated: \(0*Count\)

• False Positive: We predicted incorrectly homeowner would take the credit; allocated marketing resources; no credit allocated: \(-2850 * Count\)

• False Negative: We predicted that a homeowner would not take the credit but they did. These are likely homeowners who signed up for reasons unrelated to the marketing campaign. Thus, we ‘0 out’ this category.

We may see that

cost_benefit_table <-
   testProbs %>%
      count(predOutcome, Outcome) %>%
      summarize(True_Negative = sum(n[predOutcome==0 & Outcome==0]),
                True_Positive = sum(n[predOutcome==1 & Outcome==1]),
                False_Negative = sum(n[predOutcome==0 & Outcome==1]),
                False_Positive = sum(n[predOutcome==1 & Outcome==0])) %>%
       gather(Variable, Count) %>%
       mutate(Revenue =
               ifelse(Variable == "True_Negative", Count * 0,
               ifelse(Variable == "True_Positive", 61900 * Count,
               ifelse(Variable == "False_Negative", 0 * Count,
               ifelse(Variable == "False_Positive", (-2850) * Count, 0)))))
   
kable(cost_benefit_table) %>% 
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))  %>%
  footnote(general_title = "\n", general = "Table 7")

Variable	Count	Revenue
True_Negative	915	0
True_Positive	97	6004300
False_Negative	60	0
False_Positive	368	-1048800
Table 7

Because a different confusion matrix, set of errors and cost/benefit calculation exists for each threshold, we need to compute an ideal threshold that returns the greatest cost/benefit. The figure below visualizes how each of the four scenarios in our confusion matrix change along with the threshold. As we increase the threshold, our true positive and false positive rate both increases. This means that our likelihood of getting benefits or seeing costs both increases. The goal is therefore to find an optimal threshold that maximizes benefits and minimizes costs along this line.

all_threshold %>% 
  ggplot() +
  geom_line(aes(x = threshold, y=Count, color = Variable), size = 1.5, linetype = 1) + 
  scale_color_manual(values = palette4, guide=FALSE) + 
  facet_wrap(~Variable) + 
  labs(title = "Confusion Metric Outcomes for Each Threshold") +
  xlab("Threshold") + 
  theme(plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"), 
        axis.text.x=element_text(size=8),
        axis.text.y=element_text(size=8), 
        axis.title=element_text(size=9), 
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, linewidth =0.8))

The iterateThresholds function is based on a while loop. The threshold x, starts at 0.01 and while it is less than 1, a predicted classification, predOutcome, is mutated given x; confusion metrics are calculated; the cost benefit is performed and is appended to a data frame called all_prediction. Finally x is iterated by adding 0.01 and the process continues until x = 1.

iterateThresholds <- function(data) {
  x = .01
  all_prediction <- data.frame()
  while (x <= 1) {
  
  this_prediction <-
      testProbs %>%
      mutate(predOutcome = ifelse(Probs > x, 1, 0)) %>%
      count(predOutcome, Outcome) %>%
      summarize(True_Negative = sum(n[predOutcome==0 & Outcome==0]),
                True_Positive = sum(n[predOutcome==1 & Outcome==1]),
                False_Negative = sum(n[predOutcome==0 & Outcome==1]),
                False_Positive = sum(n[predOutcome==1 & Outcome==0])) %>%
     gather(Variable, Count) %>%
     mutate(Revenue =
               ifelse(Variable == "True_Negative", Count * 0,
               ifelse(Variable == "True_Positive", 61900 * Count,
               ifelse(Variable == "False_Negative", 0 * Count,
               ifelse(Variable == "False_Positive", (-2850) * Count, 0)))),
            Threshold = x)
  
  all_prediction <- rbind(all_prediction, this_prediction)
  x <- x + .01
  }
return(all_prediction)
}

Using this function, we computed our optimal threshold, which is 0.4 in this case.

whichThreshold <- iterateThresholds(testProbs)

whichThreshold_revenue <- 
whichThreshold %>% 
    group_by(Threshold) %>% 
    summarize(Revenue = sum(Revenue))

plot1 <- ggplot(whichThreshold_revenue)+
  geom_line(aes(x = Threshold, y = Revenue), linewidth = 2, color = "#DBC2CF")+
  geom_vline(xintercept =  pull(arrange(whichThreshold_revenue, -Revenue)[1,1]), color = "#9FA2B2", linewidth = 2.5)+
    labs(title = "Model Revenues By Threshold and Count",
         subtitle = "Vertical Line Denotes Optimal Threshold") + 
  theme(plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"), 
        axis.text.x=element_text(size=8),
        axis.text.y=element_text(size=8), 
        axis.title=element_text(size=9), 
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, linewidth =0.8))
  
credict_count <-whichThreshold %>% 
  filter(Variable == "True_Positive") %>% 
  mutate(Credit_Count = Count * 0.25)

plot2 <- credict_count %>% 
  ggplot() +
  geom_line(aes(x = Threshold, y = Credit_Count), linewidth = 2, color = "#DBC2CF") + 
  geom_vline(xintercept =  pull(arrange(whichThreshold_revenue, -Revenue)[1,1]), color = "#9FA2B2", linewidth = 2.5) +
   theme(plot.subtitle = element_text(size = 9,face = "italic"),
        plot.title = element_text(size = 12, face = "bold"), 
        axis.text.x=element_text(size=8),
        axis.text.y=element_text(size=8), 
        axis.title=element_text(size=9), 
        panel.background = element_blank(),
        panel.border = element_rect(colour = "grey", fill=NA, linewidth =0.8))

plot1 / plot2 

Finally, we compare the revenue produced with optimal threshold and a 50% threshold. It shows that a threshold of 40% is much more optimal in that it yields the greatest revenue at $6,403,3900 while a 50% threshold will only yield $4,817,450. In addition, the total count of credits earned is much higher when using the optimal threshold.

final_table<- data.frame(
  Threshold = c("Optimal_Threshold", "50%_Threshold"),
  Credit = c("34.75", "23.75"),
  Revenue = c("6403900", "4817450")
)

final_table %>% 
  kbl() %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))  %>%
  footnote(general_title = "\n", general = "Table 8")

Threshold	Credit	Revenue
Optimal_Threshold	34.75	6403900
50%_Threshold	23.75	4817450
Table 8

Conclusions

In conclusion, we have designed a model to predict the likelihood of Emil City homeowners accepting home repair tax credit with acceptable accuracy and sensitivity, but this does not mean that this model can be put into production at this point. Our primary concerns with the model include the following. Firstly, there’s significant discrepancies in terms of the number of outcomes of our predictor variables. In logistic regression, unbalanced sample sizes can lead the model to be biased towards the majority class. If the minority class is of particular interest (in our case, people who accepted the credit is the minority class), an unbalanced data set can result in poor predictive performance for the minority class, leading to a model with high accuracy because it predicts the majority class well, but performs poorly on the minority class. Despite that we have applied weight on our minority class to alleviate the effect of imbalanced dataset, it will be good to either oversample the minority class, undersample the majority class, or use alternative algorithm approaches to assign different classifications costs to different classes.

In addition, the current data set contains a lot of field categorized as unknown, and sometimes we see a statistically significant differences in mean between our unknown group and other group in a paired wise comparison. This is leading to misinterpretation in our results. In most cases, models are the most accurate when built using previous data on credit acceptance, and having unknown fields in these variables is likely lead to poor performance of our model.

Moreover, human behavior is influenced by a myriad of factors, including emotions, culture, social norms, personal experiences, and cognitive processes, and those dimensions are not fully included in our model. While capturing this complexity of human behavior in a mathematical model is inherently challenging, it is good practice to include variables such as lifestyle choices, emotional state, device usage, and most importantly, results from any feedback surveys into this model.

And finally, we need to think about what is the most “optimal” model for HCD. For organizations whose goal is solely to achieve greater economic benefits, than we have provided an optimized model for them. But for government agencies like HCD who has many unquantifiable bottom lines like equity and social cohesion, this might not be the case. Whether or not to allocated credit and to whom they allcoate credit should not be based entirely on revenues made, but on people who need those resources.