**r - read and write**, and kindly contributed to R-bloggers)

The following post details how to make a churn model in R. It was part of an interview process for which a take home assignment was one of the stages. The company stated this should take 2hrs, which is entirely unrealistic. To minimise the time cost, my analysis is very succinct and short on the exploratory analysis and amount of models compared. The churn model got me to the final stage, however little in the way of feedback was offered. There is considerable debate in the tech industry as to whether take home exams are a fair assessment or even a reasonable ask on an individual. My assessment is that this is a lazy way to interview and a very high cost on the interviewee. If you really want to work for a particular company, then sure you might be prepared to do what it takes. I doubt I will do another take home assignment.

The analysis was done in Rmarkdown and the output copied into this blog post following these helpful pointers. If you would like to play with the data and full unabridged code please visit this github repo.

## Aim

- Please create a model that predicts which businesses are likely to churn at the start of 2015 based on the
`vertical`

and`incorporation_date`

.- We have an inkling that a dropoff in the number of mandates added might be an advance indicator of someone churning. Please can you assess whether this might be true, and if you think it is useful, incorporate it into your model from part (1).

## The Data

```
# Load packages
suppressPackageStartupMessages({
library(data.table) # Fast I/O
library(dplyr) # Data munging
library(tidyr) # Data munging
library(lubridate) # Makes dates easy
library(plotly) # Interactive charts
library(magrittr) # pipe operators
library(caret) # Handy ML functions
library(rpart) # Decision Trees
library(rpart.plot) # Pretty tree plots
library(ROCR) # ML evaluation
library(e1071) # Misc stat fns
library(randomForest) # rf
})
set.seed(42)
# Read data and drop row number column
df <- fread("monthly_data_(2)_(2).csv", drop = 1)
# Have a glimpse of the data
glimpse(df)
```

```
## Observations: 902
## Variables: 27
## $ company_id
``` 4, 5, 6, 7, 8, 10, 11, 12, 13, 14...
## $ 2014-01-01_payments 8, 0, 2, 3, 0, 0, 0, 0, 0, 0, 4, ...
## $ 2014-02-01_payments 4, 0, 2, 3, 0, 0, 0, 2, 0, 11, 1,...
## $ 2014-03-01_payments 7, 39, 1, 1, 6, 0, 0, 0, 8, 0, 1,...
## $ 2014-04-01_payments 7, 0, 3, 7, 50, 1, 0, 0, 2, 0, 0,...
## $ 2014-05-01_payments 1, 54, 1, 4, 119, 0, 1, 0, 0, 0, ...
## $ 2014-06-01_payments 2, 0, 2, 1, 151, 3, 0, 0, 0, 0, 2...
## $ 2014-07-01_payments 2, 0, 2, 7, 182, 0, 0, 0, 3, 0, 0...
## $ 2014-08-01_payments 4, 22, 1, 2, 167, 0, 0, 0, 2, 0, ...
## $ 2014-09-01_payments 3, 0, 1, 5, 180, 0, 0, 0, 0, 9, 5...
## $ 2014-10-01_payments 5, 0, 2, 8, 157, 1, 0, 0, 0, 2, 1...
## $ 2014-11-01_payments 5, 0, 1, 2, 105, 0, 0, 0, 0, 0, 0...
## $ 2014-12-01_payments 9, 0, 3, 8, 57, 0, 0, 0, 0, 0, 0,...
## $ 2014-01-01_mandates 0, 4, 0, 0, 0, 4, 4, 0, 0, 22, 0,...
## $ 2014-02-01_mandates 0, 31, 1, 0, 2, 3, 8, 0, 0, 20, 1...
## $ 2014-03-01_mandates 0, 24, 0, 0, 0, 5, 19, 0, 0, 11, ...
## $ 2014-04-01_mandates 53, 18, 0, 1, 0, 0, 0, 0, 1, 11, ...
## $ 2014-05-01_mandates 0, 8, 0, 0, 0, 0, 0, 0, 1, 15, 0,...
## $ 2014-06-01_mandates 0, 7, 0, 0, 0, 0, 0, 0, 0, 13, 5,...
## $ 2014-07-01_mandates 0, 0, 0, 0, 0, 0, 0, 0, 0, 13, 16...
## $ 2014-08-01_mandates 0, 0, 0, 0, 0, 0, 0, 0, 0, 14, 8,...
## $ 2014-09-01_mandates 0, 0, 0, 0, 0, 0, 0, 0, 0, 16, 2,...
## $ 2014-10-01_mandates 0, 0, 0, 0, 1, 0, 0, 0, 0, 16, 1,...
## $ 2014-11-01_mandates 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, ...
## $ 2014-12-01_mandates 0, 0, 0, 3, 0, 0, 0, 0, 0, 11, 0,...
## $ vertical "gym/fitness", "freelance d...
## $ incorporation_date "2003-09-25", "2008-10-22", ...

## Data Munging

Some data munging needs to occur for our binary classifiers to make use of the data within. This includes handling dates.

```
# Reshape data and create new columns
df %<>%
gather(key = date, value = quantity, starts_with("20")) %>%
separate(date, c("date","paymentMandate"), "_") %>%
spread(paymentMandate, quantity) %>%
mutate(incorporation_date = as.Date(incorporation_date),
date = as.Date(date),
incorporation_time = round(as.numeric(difftime(as.Date("2014-12-01"),
as.Date(incorporation_date),
unit="weeks")) / 52.25,
digits = 1)) %>%
arrange(date)
# What does the new data look like?
glimpse(df)
```

```
## Observations: 10,824
## Variables: 7
## $ company_id
``` 1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 13, 14,...
## $ vertical "gym/fitness", "gym/fitness", "freelance de...
## $ incorporation_date 2013-05-30, 2003-09-25, 2008-10-22, 2005-0...
## $ date 2014-01-01, 2014-01-01, 2014-01-01, 2014-0...
## $ mandates 1, 0, 0, 0, 4, 0, 0, 0, 4, 4, 0, 0, 22, 0, ...
## $ payments 0, 1, 6, 8, 0, 2, 3, 0, 0, 0, 0, 0, 0, 4, 4...
## $ incorporation_time 1.5, 11.2, 6.1, 9.4, 1.2, 9.1, 4.0, 2.9, 1....

### What is Churn

For the purposes of this project I have defined ‘churn’ as zero mandates and payments for the first instance of that after the last mandate/payment made.

```
# Create binary 'churn' column
df$churn <- 0
# For use in for loop - upper bound of data
max.date <- max(df$date)
# Mark all companies as churned in the month immediately their last activity
for (company in unique(df$company_id)) {
# Subset data to company
df.sub <- subset(df, df$company_id == company)
# Index of last positive mandate OR payment
last.pos.idx <- tail(which(df.sub$mandates != 0 | df.sub$payments != 0), 1)
# Get date of last activity
last.activity.date <- df.sub$date[last.pos.idx]
# If less than max.date of dataset mark churn ELSE do nothing i.e. positive at end of period
if (last.activity.date < max.date) {
# Get churn date (last positive month plus 1mth)
churn.date <- last.activity.date %m+% months(1)
# Mark month of churn as 1
df[df$date == churn.date & df$company_id == company, ]$churn <- 1
}
}
# Multiple rows per company, filter for last month or churn month values...
# Get churners
df %>% filter(churn == 1) -> churners
# Get max date row of remainers (non-churners)
df %>%
filter(churn == 0 & !(company_id %in% churners$company_id) & date == max(date)) -> remainers
# Combine and variables coded ready for modelling
churners %>%
rbind(remainers) %>%
mutate(vertical = as.factor(vertical),
churn = as.factor(churn)) -> model.df
```

### Churns over the year

```
# Plot churners
model.df %>%
filter(churn == 1) %>%
group_by(date) %>%
summarise(n = n()) %>%
plot_ly( x = ~date, y = ~n, type = 'scatter', mode = 'lines')
```

Could do more on exploratory but this is the easiest to cut back for this report.

### Balance of the data

The proportion of churn is 32.04%. Representing an imbalanced dataset. Accuracy is an inappropriate measure (I could get 67.96% accuracy predicting no businesses leave), so I will focus on recall and accuracy.

```
# Loyal vs Churn
table(model.df$churn)
##
## 0 1
## 613 289
```

## Model

Survival models and binary classifiers are common approaches to ‘Churn’ models. I will approach the model using the latter, though if I had more time I would investigate other classifiers and a survival model. I will limit the range of models to a logistic, a decision tree and an ensemble.

### Split Data

```
# 80/20 train test split
index <- createDataPartition(model.df$churn, p = 0.8, list = FALSE)
train.df <- model.df[index, ]
test.df <- model.df[-index, ]
# Check balance of the training split
table(train.df$churn)
##
## 0 1
## 491 232
```

```
# Check balance of the test split
table(test.df$churn)
##
## 0 1
## 122 57
```

### Logistic

```
# Run model
Logistic.model <- glm(churn ~ incorporation_time + vertical,
data = train.df,
family = binomial(link = 'logit'))
# Predict
log.pred <- predict(Logistic.model, newdata = test.df, type = 'response')
# Convert probs to binary
log.pred <- as.factor(ifelse(log.pred > 0.5, 1, 0))
# Evaluation Metrics
log.result <- confusionMatrix(data = log.pred, test.df$churn)
log.precision <- log.result$byClass['Pos Pred Value']
log.recall <- log.result$byClass['Sensitivity']
log.F1 <- log.result$byClass['F1']
```

### Decision Tree

```
# Train model
tree.model <- rpart(churn ~ incorporation_time + vertical,
data = train.df,
method = "class",
control = rpart.control(xval = 10))
# Plot
rpart.plot(tree.model)
```

```
# Evaluation metrics
tree.pred <- predict(tree.model, newdata = test.df, type = "class")
tree.result <- confusionMatrix(data = tree.pred, test.df$churn)
tree.precision <- tree.result$byClass['Pos Pred Value']
tree.recall <- tree.result$byClass['Sensitivity']
tree.F1 <- tree.result$byClass['F1']
```

### Random Forest (Ensemble)

```
# Train model
forest.model <- randomForest(churn ~ incorporation_time + vertical,
data = train.df,
ntree=200,
type="classification")
# See error reduction with number of trees ( not much gained beyond ~25 trees)
plot(forest.model)
```

```
# Look at the variable Importance from the random forest
varImpPlot(forest.model, sort = T, main="Variable Importance")
```

```
# Evaluation metrics
forest.pred <- predict(forest.model, newdata = test.df, type = "class")
forest.result <- confusionMatrix(data = forest.pred, test.df$churn)
forest.precision <- forest.result$byClass['Pos Pred Value']
forest.recall <- forest.result$byClass['Sensitivity']
forest.F1 <- forest.result$byClass['F1']
```

### Evaluation Metrics

```
log.precision;
## Pos Pred Value
## 0.8333333
tree.precision;
## Pos Pred Value
## 0.8088235
forest.precision;
## Pos Pred Value
## 0.8134328
```

```
log.recall;
## Sensitivity
## 0.9016393
tree.recall;
## Sensitivity
## 0.9016393
forest.recall;
## Sensitivity
## 0.8934426
```

```
log.F1;
## F1
## 0.8661417
tree.F1;
## F1
## 0.8527132
forest.F1;
## F1
## 0.8515625
```

Surprisingly, the logistic regression model performs the best, with the top precision score and equal recall score with that of the decision tree. With more time, I’d see if tweaks to the decision tree and random forest models could change this. Its also possible over/undersampling could help. From here on I will use the logistic regression model.

## Examine the incorporation of time information

Per the second aim, examine the inclusion of time information. This requires more data munging. Taking code used from above, I will extend to include the derivation of a very basic time period variable.

```
# Create binary 'leading_indicator' column
model.df$leading_indicator <- 0
# Min date for which a leading_indicator can be calculated (lower limit of data)
min.date <- min(df$date)
# If month before 'churn' (churn-1), is below the level of mandates of the month 2 months prior (churn-2) then
# make leading_indicator == 1
for (company in churners$company_id) {
# Subset data to company
df.sub <- subset(df, df$company_id == company)
# Get month prior to churn
month.prior <- df.sub$date[df.sub$churn == 1] %m-% months(1)
# Get two months prior to churn
two.month.prior <- df.sub$date[df.sub$churn == 1] %m-% months(2)
# If two months prior is within dataset date range and level of mandates is greater than 0
if ((two.month.prior > min.date) && (df.sub$mandates[df.sub$date == two.month.prior] > 0)) {
# Compare number of mandates 1 month prior to 2 months prior, if less, mark 'leading_indicator' as '1'
if (df.sub$mandates[df.sub$date == month.prior] < df.sub$mandates[df.sub$date == two.month.prior]) {
model.df[model.df$company_id == company, ]$leading_indicator <- 1
}
}
}
```

### Of the churners, how many have a leading indicator?

```
model.df %>%
filter(churn == 1) %>%
group_by(leading_indicator) %>%
summarise(n=n()) %>%
mutate(percent = round(n / sum(n) * 100, 1))
```

```
## # A tibble: 2 x 3
## leading_indicator n percent
##
```
## 1 0 264 91.3
## 2 1 25 8.7

### Re-train model and evaluate

```
# Re-split data so 'leading_indicator' is in columns (the index remains the same)
train.df <- model.df[index, ]
test.df <- model.df[-index, ]
```

### Logistic

```
# Run model
lead.logistic.model <- glm(churn ~ incorporation_time + vertical + leading_indicator,
data = train.df,
family = binomial(link = 'logit'))
# Predict
log.pred <- predict(lead.logistic.model, newdata = test.df, type = 'response')
# Convert probs to binary
log.pred <- as.factor(ifelse(log.pred > 0.5, 1, 0))
# Evaluation Metrics
lead.log.result <- confusionMatrix(data = log.pred, test.df$churn)
lead.log.precision <- log.result$byClass['Pos Pred Value']
lead.log.recall <- log.result$byClass['Sensitivity']
lead.log.F1 <- log.result$byClass['F1']
```

### Evaluation Metrics

Compare the precision and recall of the logistic model with and without the `lead_indicator`

.

```
log.precision
## Pos Pred Value
## 0.8333333
```

```
lead.log.precision
## Pos Pred Value
## 0.8333333
```

```
log.recall
## Sensitivity
## 0.9016393
```

```
lead.log.recall
## Sensitivity
## 0.9016393
```

```
# Have a look at the model coefficients and p-values
summary(lead.logistic.model)
## Call:
## glm(formula = churn ~ incorporation_time + vertical + leading_indicator,
## family = binomial(link = "logit"), data = train.df)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4429 -0.8252 -0.4957 0.9938 2.5564
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.74066 0.20007 3.702 0.000214 ***
## incorporation_time -0.26140 0.03045 -8.584 < 2e-16 ***
## verticalfreelance developer 0.12614 0.21737 0.580 0.561714
## verticalgym/fitness -0.91099 0.22342 -4.077 4.55e-05 ***
## leading_indicator 18.06219 482.97671 0.037 0.970168
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 907.42 on 722 degrees of freedom
## Residual deviance: 736.18 on 718 degrees of freedom
## AIC: 746.18
##
## Number of Fisher Scoring iterations: 15
```

Generally, it is best to minimise the `AIC`

. The `Logistic.model`

model AIC is 800.5 compared with the `lead.logistic.model`

incorporating the `lead_indicator`

is 746.18. The `lead_indicator`

has not changed the predictive power on this dataset, but since the `AIC`

favours a better model fit whilst penalising for additional predictors, the model to choose is the `lead.logistic.model`

.

#### Re-train model and save

```
# Re-train on all data
final.model <- glm(churn ~ incorporation_time + vertical + leading_indicator,
data = model.df,
family = binomial(link = 'logit'))
save(final.model, file = "model.rda")
```

## Recommendations for further investigation/Comments

- The
`lead_indicator`

was a quick and dirty test, I’d spend more time looking at a better construction (e.g. moving avg) - If the cost of rentention activities vs losing a customer was known then an the optimal trade-off in terms of business cost could be found.
- Incorporating additional data, e.g. CRM data showing interactions. Potentially assessing sales/account staff.
- If it were a production model prompting staff to do retention calls, evaluate the impact of such calls through modelling e.g. a/b testing
- Test model performance with a change to balancing the classes e.g. under/oversampling, boostrap samples. This may explain the relative underperformance of tree based models in this exercise.
- Try other binary classifier models.

**leave a comment**for the author, please follow the link and comment on their blog:

**r - read and write**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...