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

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.

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
##               <dbl> <int>   <dbl>
## 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 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")