Divorce Prediction

This project uses machine learning with 2016 American Community Survey Data to predict the probability of divorce.

Environment Set-Up

Variables

First we need to set up the environment

# start clean
rm(list = ls())
# Set working directory
setwd("/Users/cueland/Google Drive/School/670 - Data Science/divorce-prediction/")
# useful libraries
library(plyr)
library(ggplot2)
library(dplyr)
library(reshape2)
library(rpart)
library(rpart.plot)
library(devtools)
library(gridExtra)
library(randomForest)
library(plotROC)

Functions

Set up useful functions.

#Mean F1
meanf1 <- function(actual, predicted){
  #Mean F1 score function
  #actual = a vector of actual labels
  #predicted = predicted labels
  
  classes <- unique(actual)
  results <- data.frame()
  for(k in classes){
    results <- rbind(results, 
                     data.frame(class.name = k,
                                weight = sum(actual == k)/length(actual),
                                precision = sum(predicted == k &
                                                  actual == k) / 
                                  sum(predicted == k),
                                recall = sum(predicted == k &
                                               actual == k) /
                                  sum(actual == k)))
  }
  results$score <- results$weight * 2 *
    (results$precision * results$recall) /
    (results$precision + results$recall) 
  return(sum(results$score))
}
inter <- function(x,y,z) {
  # Create a sequence of numbers based on the number of values you want
  # in the sequence
  # x,y = beginning and end of sequence
  # z = number of intervals to produce
  f <- c(x,y)
  int <- abs(x - y) / z
  d <- seq(min(f),max(f),int)
  return(d)
}

Gather the Data

Data Dictionary

To interpret the actual data download, we first need to download the data dictionary, which will then be used to organize the data we download from the ACS API.

Data Download

ACS data is downloaded by the month, therefore in order to download all of the data from 2016, we need to create a for loop that collects each month’s data from 2016. First we build the URL for each request, download the files, and then process the data using the data dictionary downloaded above.

Process the Data

In order to use the data for divorce predictions, we need to ensure we only use valid data, so we need to drop observations such as surveys that weren’t fully completed, or observations for people under the age of 18.

In addition, we need to transform the data, factorizing the information that is provided in categories, such as marital status.

# load the data from before
pt <- proc.time()[3]
load("/Users/cueland/Google Drive/School/670 - Data Science/divorce-prediction/2016.RData")
proc.time()[3] - pt

elapsed 
 270.59 

# convert all column names to lowercase
names(l) <- tolower(names(l))
# ------------------------------------ Cull Data -----------------------------|
# Select which variables are useful
code <- c("hrhhid", "hrhhid2", "hwhhwtln", "pulineno", "hrmonth", "hryear4", "hurespli", "hufinal", "huspnish", "hetenure", "hehousut", "hefaminc", "hrnumhou", "hubus", "gereg", "gediv", "gtcbsast", "gtmetsta", "gtindvpc", "gtcbsasz", "prtage", "prtfage", "pemaritl", "pesex", "peafever", "peafnow", "peeduca", "ptdtrace", "prdthsp", "pehspnon", "prcitshp", "prinusyr", "pemlr", "puwk", "pubus1", "pubus2ot", "puretot", "pudis", "peret1", "pudis1", "puabsot", "pulay", "pemjot", "pemjnum", "pehrftpt", "pehruslt", "pehrwant", "pehrrsn1", "pehrrsn2", "pehrrsn3", "puhroff1", "puhrot1", "puhrot2", "pehractt", "pulk", "pedwwnto", "prcivlf", "prdisc", "prftlf", "prwksch", "prwkstat", "prmjind1", "prmjocc1", "penlfact", "prchld", "prnmchld", "pedisear", "pediseye", "pedisrem", "pedisphy", "pedisdrs", "pedisout", "pepdemp1")
# Subset to the desired variables
data <- l[code]
# convert the data to all numeric values
for (n in names(data)) {
  data[[n]] <- as.numeric(data[[n]])
}
# Get rid of any observations that don't have a line number
data <- data[data$pulineno > 0, ]
# Keep only complete surveys
data <- data[data$hufinal %in% c(1, 201), ]
# Keep only people 18 years and older
data <- data[data$prtage > 17, ]
# reverse sort by month to only keep the latest observation from 2016
data <- data[order(-as.numeric(data$hrmonth)), ]
data <- data[!duplicated(cbind(data$hrhhid, data$pulineno,
                               data$hrmonth, data$hufinal)),]
# Convert all negative values to NAs
for (n in names(data)) {
  data[[n]][data[[n]] < 0] <- NA
}
# Recode & Factorize marital status
data$mstatus[data$pemaritl %in% c(1, 2, 3)] <- "married"
data$mstatus[data$pemaritl %in% c(4, 5)] <- "divorced"
data$mstatus[data$pemaritl == 6] <- "single"
data$pemaritl <- factor(data$mstatus)
data$mstatus <- NULL
# show distrubution of marriage status
table(data$pemaritl, useNA = "ifany")

divorced  married   single 
  149597   699837   299189 

# Save the numeric data frame before we factorize the rest of the factor vars
data.numeric <- data
# Variables to Factorize (all variables that are not numeric as listed below)
facs <- names(data)[!(names(data) %in% c("hrhhid", "hryear4", "hrnumhou",
                                         "prtage", "pemjnum", "pehruslt",
                                         "puhrot2", "pehractt", "prnmchld"))]
# Loop through each variable to factorize
for (n in facs) {
  data[[n]] <- factor(data[[n]])
}
head(data)

save(data, file = "/Users/cueland/Google Drive/School/670 - Data Science/divorce-prediction/data.RData", compress = TRUE)

Run Decision Tree

First we divide the data randomly into 2 separate parts, train and test, with an approximate 70/30 split. From there, we train the decision tree model using a selected set of variables from the original set. This set of variables was chosen from multiple tests with different combinations, as well as culling certain variables that are explicitly related to marriage status.

set.seed(12345)
train <- runif(nrow(data))
train <- train > 0.7
# Create separate train and test sets
data.train <- data[train, ]
head(data.train[1:6])

data.test <- data[!train, ]
head(data.test[1:6])

Optimization

Because we ran the model with CP = 0 (the most complex), we can subsequently simplify, or “prune”, the decision tree model in order to avoid over-fitting the model to the train data. We use two approaches for this.

XERROR

First, we examine the XERROR of each level of complexity and prune to the level that offers the lowers XERROR value.

data.train$predict <- predict(fit, data.train, type= "class")
meanf1(data.train$predict, data.train$pemaritl)

[1] 0.8713695

# Test the predictions for the test subset for CP = 0
data.test$predict <- predict(fit, data.test, type= "class")
meanf1(data.test$predict, data.test$pemaritl)

[1] 0.7917867

# ------------------------------ min XERROR (mincp) ---------------------------|
# determine the CP that yields the minimum XERROR value 
mincp <- fit$cptable[fit$cptable[, 4] == min(fit$cptable[, 4]), 1]
if (length(mincp) > 1) {
  fit$cptable[fit$cptable[, 4] == min(fit$cptable[, 4]), ]
  mincp <- mincp[1]
}
# prune to the CP that yielded the min XERROR value, mincp
fit.mincp <- prune(fit, cp = mincp)
# Test the predictions for the test subset for CP = mincp
data.test$predict.mincp <- predict(fit.mincp, data.test, type= "class")
meanf1(data.test$predict.mincp, data.test$pemaritl)

[1] 0.8119978

Iteration

Second, we use an iterative approach to determine which level of pruning yields the highest mean F1, which is our measure of best fit. This is done with a series of two loops, the first to find a rough estimate of the best CP value, and the second to further refine the CP value. This yields a mean F1 score for the test subset of 0.812.

# create data frame with results to see how optimization worked
f1s.refined <- data.frame(CP = rev(inter(next.it[1],next.it[2], 10)),
                          mean_F1 = rev(f1s))
f1s.refined

plot(f1s.refined)  # plot optimization

# assign best CP value
cp.best <- f1s.refined[,1][f1s.refined[,2] == max(f1s.refined[,2])]
# find optimal fit
fit.best <- prune(fit, cp = cp.best)
useful.vars.best <- data.frame(name = attr(fit.best$variable.importance,
                                           "names"),
                          importance = fit.best$variable.importance,
                          row.names = NULL)
useful.vars.best

# Show further simplified decision tree graphic
fit.simple <- prune(fit, cp = 0.005)
rpart.plot(fit.simple, shadow.col = "grey", nn = TRUE)

# Test the predictions for the test subset for CP = mincp
data.test$predict.best <- predict(fit.best, data.test, type= "class")
meanf1(data.test$predict.best, data.test$pemaritl)

[1] 0.8119281