reporting in word

Reporting Analysis done with R into a Word Document

We will use R to analyse the Bank Marketing dataset and apply a Logistic Regression model on it to predict if the client will subscribe a term deposit (variable y), then deliver results in a Word document using R officer & flextable libraries.

options(warn = -1)
library(officer)
library(flextable)
library(ggplot2)
library(reshape)
library(dplyr)
library(gridExtra)
library(ROCR)

data <- read.csv("bank.csv", sep = ";")
str(data)

## 'data.frame':    4521 obs. of  17 variables:
##  $ age      : int  30 33 35 30 59 35 36 39 41 43 ...
##  $ job      : Factor w/ 12 levels "admin.","blue-collar",..: 11 8 5 5 2 5 7 10 3 8 ...
##  $ marital  : Factor w/ 3 levels "divorced","married",..: 2 2 3 2 2 3 2 2 2 2 ...
##  $ education: Factor w/ 4 levels "primary","secondary",..: 1 2 3 3 2 3 3 2 3 1 ...
##  $ default  : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ balance  : int  1787 4789 1350 1476 0 747 307 147 221 -88 ...
##  $ housing  : Factor w/ 2 levels "no","yes": 1 2 2 2 2 1 2 2 2 2 ...
##  $ loan     : Factor w/ 2 levels "no","yes": 1 2 1 2 1 1 1 1 1 2 ...
##  $ contact  : Factor w/ 3 levels "cellular","telephone",..: 1 1 1 3 3 1 1 1 3 1 ...
##  $ day      : int  19 11 16 3 5 23 14 6 14 17 ...
##  $ month    : Factor w/ 12 levels "apr","aug","dec",..: 11 9 1 7 9 4 9 9 9 1 ...
##  $ duration : int  79 220 185 199 226 141 341 151 57 313 ...
##  $ campaign : int  1 1 1 4 1 2 1 2 2 1 ...
##  $ pdays    : int  -1 339 330 -1 -1 176 330 -1 -1 147 ...
##  $ previous : int  0 4 1 0 0 3 2 0 0 2 ...
##  $ poutcome : Factor w/ 4 levels "failure","other",..: 4 1 1 4 4 1 2 4 4 1 ...
##  $ y        : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...

summary(data)

##       age                 job          marital         education   
##  Min.   :19.00   management :969   divorced: 528   primary  : 678  
##  1st Qu.:33.00   blue-collar:946   married :2797   secondary:2306  
##  Median :39.00   technician :768   single  :1196   tertiary :1350  
##  Mean   :41.17   admin.     :478                   unknown  : 187  
##  3rd Qu.:49.00   services   :417                                   
##  Max.   :87.00   retired    :230                                   
##                  (Other)    :713                                   
##  default       balance      housing     loan           contact    
##  no :4445   Min.   :-3313   no :1962   no :3830   cellular :2896  
##  yes:  76   1st Qu.:   69   yes:2559   yes: 691   telephone: 301  
##             Median :  444                         unknown  :1324  
##             Mean   : 1423                                         
##             3rd Qu.: 1480                                         
##             Max.   :71188                                         
##                                                                   
##       day            month         duration       campaign     
##  Min.   : 1.00   may    :1398   Min.   :   4   Min.   : 1.000  
##  1st Qu.: 9.00   jul    : 706   1st Qu.: 104   1st Qu.: 1.000  
##  Median :16.00   aug    : 633   Median : 185   Median : 2.000  
##  Mean   :15.92   jun    : 531   Mean   : 264   Mean   : 2.794  
##  3rd Qu.:21.00   nov    : 389   3rd Qu.: 329   3rd Qu.: 3.000  
##  Max.   :31.00   apr    : 293   Max.   :3025   Max.   :50.000  
##                  (Other): 571                                  
##      pdays           previous          poutcome      y       
##  Min.   : -1.00   Min.   : 0.0000   failure: 490   no :4000  
##  1st Qu.: -1.00   1st Qu.: 0.0000   other  : 197   yes: 521  
##  Median : -1.00   Median : 0.0000   success: 129             
##  Mean   : 39.77   Mean   : 0.5426   unknown:3705             
##  3rd Qu.: -1.00   3rd Qu.: 0.0000                            
##  Max.   :871.00   Max.   :25.0000                            
## 

docx <- read_docx()
officer::body_add_par(docx, value = "Table of Content", 
                      style = "heading 1")
body_add_toc(docx, level = 2)
body_add_break(docx)

title_case <- function(x) {
  lower_case <- tolower(x)
  without_underscore <- gsub("_", " ", lower_case, perl = TRUE)
  titled <- gsub("(^|[[:space:]])([[:alpha:]])",
                 "\\1\\U\\2",
                 without_underscore,
                 perl = TRUE)
  titled
}

# trim <- function(x) {
#   gsub("^\\s+|\\s+$", "", x)
# }

create_flextable <- function(df, theme = "box") {
  thm <- paste0("flextable::theme_", theme, "(ft)")
  ft <- flextable::flextable(df)
  ns <- title_case(names(df))
  ex <- paste0('"', names(df), '"', "=", '"', ns, '"', collapse = ",")
  expr <- paste0("flextable::set_header_labels(ft,")
  exp <- paste0(expr, ex, ")")
  ft <- eval(parse(text = exp))
  ft <- eval(parse(text = thm))
  ft <- ft %>%
    flextable::bold(part = "header") %>%
    flextable::bold(j = 1) %>%
    flextable::font(fontname = "Arrial Narrow", part = "all") %>%
    flextable::fontsize(size = 10, part = "all") %>%
    flextable::align(align = "left", part = "all") %>%
    flextable::autofit()
  ft
}

body_add_par(docx, "Exploring Data", style = "heading 1")
body_add_par(docx, "Summarizing Data", style = "heading 2")
body_add_par(docx, '- Examining the "age" and "balance" averages of each "Marital Status" group',
             style = "Normal")
age_balance_marital <- data %>% group_by(marital) %>%
  summarise(age_average = mean(age),
            balance_average = mean(balance))
abm <- create_flextable(age_balance_marital, "box")
flextable::body_add_flextable(docx, value = abm, align = "left")

body_add_par(docx, "", style = "Normal")
body_add_par(docx, '- Examining the "age" and "balance" averages of each "Job" who subscribed in this campaign',
             style = "Normal")
age_balance_job <-
  data %>% group_by(job, y) %>% filter(y == "yes") %>%
  summarise(age_average = mean(age),
            balance_average = mean(balance)) %>%
  arrange(desc(balance_average))
abj <- create_flextable(age_balance_job, "zebra")
flextable::body_add_flextable(docx, value = abj, align = "left")

body_add_par(docx, "", style = "Normal")
body_add_par(docx, '- Examining the "duration" of each "Contact Method"',
             style = "Normal")
contact_duration <- data %>% group_by(contact) %>%
  summarise(duration_average_min = mean(duration)) %>%
  mutate(duration_average_min = duration_average_min / 60)
cd <- create_flextable(contact_duration, "vanilla")
flextable::body_add_flextable(docx, value = cd, align = "left")

cursor_end(docx)
body_add_break(docx)
body_add_par(docx, "Visualizing Data", style = "heading 2")
body_add_par(docx, "- Visualizing distribution of numeric columns: ", 
             style = "Normal")

# in case we don't know the exact spelling of the column names we can 
# extract them with this regex
regex <- "(?=.*age)|(?=.*balanc)|(?=.*durat)"
cols <- colnames(data)
nums <- grep(regex, cols, perl = TRUE, value = TRUE)

age_balance_duration <- data[, c(nums, "y")] %>%
  mutate(balance = log(balance, base = 10),
         duration = duration / 60) %>%
  dplyr::rename(log10_balance = balance,
                duration_min = duration) %>%
  melt(id.vars = "y")

ggplot(age_balance_duration, aes(x = value, fill = y)) +
  geom_density(color = "white", alpha = 0.6) +
  facet_wrap( ~ variable, scales = "free") +
  theme_minimal() + scale_fill_brewer(palette = "Set1") +
  theme(legend.position = "top")

filepath <- "."
filename <- paste0(dirname(filepath), "/density.png")
ggsave(filename = filename, height = 4, width = 6)
body_add_img(docx, src = filename, height = 4, width = 6)

body_add_par(docx, "", style = "Normal")
body_add_par(docx, "- Visualizing relationship between numeric columns and y: ", style = "Normal")
ggplot(age_balance_duration, aes(x = value, y = as.numeric(y))) +
  geom_point(position = position_jitter(w = 0.05, h = 0.05),
             alpha = 0.3) +
  geom_smooth(method = "loess", se = FALSE) +
  facet_wrap( ~ variable, scales = "free", ncol = 2) +
  theme_minimal() + ylab("Subscription")

filename <- paste0(dirname(filepath), "/scatter.png")
ggsave(filename = filename, height = 4, width = 6)
body_add_img(docx, src = filename, height = 4, width = 6)

body_add_par(docx, "", style = "Normal")
body_add_par(docx, "- spotting days, months and jobs in which users tend to subscribe", style = "Normal")
data <- data %>%
  mutate(month =
           factor(gsub("(^[[:alpha:]])", "\\U\\1",
                       month, perl = TRUE),
                  month.abb))
gg <- ggplot(data, aes(x = as.numeric(month), color = y)) +
  geom_freqpoly(stat = "count") +
  xlab("month") + 
  scale_x_continuous(breaks = seq(1, 12, 1), labels = month.abb) +
  scale_y_continuous(breaks = seq(0, 1500, 250)) +
  theme_minimal() + scale_color_brewer(palette = "Set1") +
  xlab("last contact month of the year")

gg2 <- ggplot(data, aes(x = job, fill = y)) +
  geom_bar(width = 0.8,
           color = "white",
           alpha = 0.8) +
  theme_minimal() + scale_fill_brewer(palette = "Set1")

gg3 <- ggplot(data, aes(x = day, fill = y)) +
  geom_histogram(bins = 30, col = "white") + theme_minimal() +
  scale_fill_brewer(palette = "Accent") +
  xlab("last contact day of the month")

ggsave("daymonthjob.png", 
       arrangeGrob(gg, gg2, gg3), height = 4, width = 7)
body_add_img(docx, src = "daymonthjob.png", height = 4, width = 7)

body_add_par(docx, "", style = "Normal")
body_add_par(docx, "- Visualizing frequency of categorical variables and their relationship with y: ",
             style = "Normal")

# function (not in)
'%!in%' <- function(x, y)
  ! ('%in%'(x, y))

vars <- colnames(data)
# summary of numeric data
summary(data[, vars[sapply(data[, vars], class) %in% 
                      c("numeric", "integer")]])

##       age           balance           day           duration   
##  Min.   :19.00   Min.   :-3313   Min.   : 1.00   Min.   :   4  
##  1st Qu.:33.00   1st Qu.:   69   1st Qu.: 9.00   1st Qu.: 104  
##  Median :39.00   Median :  444   Median :16.00   Median : 185  
##  Mean   :41.17   Mean   : 1423   Mean   :15.92   Mean   : 264  
##  3rd Qu.:49.00   3rd Qu.: 1480   3rd Qu.:21.00   3rd Qu.: 329  
##  Max.   :87.00   Max.   :71188   Max.   :31.00   Max.   :3025  
##     campaign          pdays           previous      
##  Min.   : 1.000   Min.   : -1.00   Min.   : 0.0000  
##  1st Qu.: 1.000   1st Qu.: -1.00   1st Qu.: 0.0000  
##  Median : 2.000   Median : -1.00   Median : 0.0000  
##  Mean   : 2.794   Mean   : 39.77   Mean   : 0.5426  
##  3rd Qu.: 3.000   3rd Qu.: -1.00   3rd Qu.: 0.0000  
##  Max.   :50.000   Max.   :871.00   Max.   :25.0000

# change pdays variable to factor of intervals and delete it
data$days_after_last_contact <-
  cut(data$pdays,
      breaks = c(-1, 0, 150, seq(300, 900, 300)),
      include.lowest = T)
data$pdays <- NULL

# change campaign and previous to factor of intervals and rename them
data <- data %>% mutate(
  campaign = cut(
    campaign,
    breaks = c(seq(1, 10, 2),
               seq(25, 50, 25)),
    include.lowest = TRUE
  ),
  previous = cut(
    previous,
    breaks = c(seq(0, 6, 2),
               seq(15, 25, 10)),
    include.lowest = TRUE
  )
) %>% dplyr::rename(
    number_contacts_this_campaign = campaign,
    number_contacts_before_this_campaign = previous
  )

# new column names
vars <- colnames(data)
# categorical variables
cat_vars <-
  vars[sapply(data[, vars], class) %in% c("factor", "character")]
# we have already examined month and job
cat_vars <- cat_vars[cat_vars %!in% c("month", "job")]

cat_melted <- data[, cat_vars] %>%
  melt(id.vars = "y")
ggplot(cat_melted, aes(x = value, fill = y)) + 
  geom_bar(width = 0.5, color = "white") +
  theme_minimal() + theme(
    legend.position = c(0.6, 0.1),
    legend.direction = "horizontal",
    legend.key.size = unit(10, "mm"),
    legend.title = element_text("Subscription")
  ) +
  facet_wrap(~ variable, scales = "free", ncol = 3) +
  scale_fill_brewer(palette = "Paired")

ggsave("catvars.png", height = 5, width = 7)
body_add_img(docx, src = "catvars.png", height = 5, width = 7)

body_add_break(docx)
body_add_par(docx, "Modeling Data", style = "heading 1")
body_add_par(docx, "Preparing data for modeling", 
             style = "heading 2")
body_add_par(docx, "- one hot encoding of categorical variables and normalizing numeric ones", 
             style = "Normal")
body_add_par(docx, paste("dimenstions of original data:",
                         paste(dim(data), collapse = " , ")),
             style = "Normal")

data$day <- factor(data$day)
# predictors and outcome
pred_vars <- setdiff(vars, "y")

## to one hot encode factor values and normalize numeric ones if needed
cat <-
  pred_vars[sapply(data[, pred_vars], class) %in% c("factor", "character")]
num <-
  pred_vars[sapply(data[, pred_vars], class) %in% c("numeric", "integer")]

for (i in cat) {
  dict <- unique(data[, i])
  for (key in dict) {
    data[[paste0(i, "_", key)]] <- 1.0 * (data[, i] == key)
  }
}

data[, num] <- apply(data[, num], 2, function(x) {
  (x - min(x)) / (max(x) - min(x))
})

data <- data[, -which(colnames(data) %in% cat)]

body_add_par(docx, paste("dimenstions of new encoded data:",
                         paste(dim(data), collapse = " , ")),
             style = "Normal")

cursor_end(docx)
body_add_par(docx, "Train/Test Splitting", 
             style = "heading 2")
body_add_par(docx, "- splitting data to train and test datasets (80/20)", style = "Normal")

sample_size <- floor(0.8 * nrow(data))
## set the seed to make your splits reproducible
set.seed(13)
train_indices <- sample(seq(nrow(data)), size = sample_size)
train <- data[train_indices,]
test <- data[-train_indices,]

body_add_par(docx, paste("nrows of training dataset:", 
                         nrow(train)), 
             style = "Normal")

body_add_par(docx, paste("nrows of test dataset:", 
                         nrow(test)), 
             style = "Normal")

body_add_par(docx, "Fitting Classification Model", 
             style = "heading 2")
body_add_par(docx, "- Applying Logistic Regression algorithm to data: ", style = "Normal")

glmmodel <- glm(y ~ ., data = train, family = binomial(link = "logit"))

# how well the model fits the data it has seen
train$pred <- predict(glmmodel, newdata = train, type = "response")

# visualizing classification thresheold
body_add_par(docx, "", style = "Normal")
body_add_par(docx, "Spotting Classification Threshold Probability Value: ",
             style = "Normal")

ggplot(train, aes(x = pred, fill = y)) + 
  geom_density(alpha = 0.6, col = "white") +
  scale_fill_brewer(palette = "Set1") + theme_minimal()

ggsave("threshold.png", height = 4, width = 5)
body_add_img(docx, src = "threshold.png", height = 4, width = 5)

body_add_par(docx, "", style = "Normal")
predobj <- prediction(train$pred, train$y)
perf <- performance(predobj, "sens", "spec")
sensitivity <- perf@y.values[[1]]
specificity <- perf@x.values[[1]]
acc <- performance(predobj, "acc")
accuracy <- acc@y.values[[1]]
error.rate <- 1 - accuracy
threshold <- acc@x.values[[1]]
errors <- data.frame(cbind(threshold,
                           cbind(error.rate,
                                 cbind(
                                   accuracy,
                                   cbind(sensitivity, specificity)
                                 ))))

error.data <- melt(errors, id.vars = "threshold")
ggplot(error.data, aes(x = threshold, y = value, 
                       col = variable, linetype = variable)) + 
  geom_line() + theme_minimal() +
  scale_color_brewer(palette = "Set1") + 
  scale_x_continuous(breaks = seq(0, 1, 0.1)) +
  theme(legend.position = "top")

ggsave("threshold2.png", height = 4, width = 5)
body_add_img(docx, src = "threshold2.png", height = 4, width = 5)

body_add_par(docx, "The model seems to classify well at probablity: 0.12",
             style = "Normal")
body_add_par(docx, "", style = "Normal")
body_add_par(docx, "Evaluating the Model", 
             style = "heading 2")
body_add_par(docx, "- Measuring Model Accuracy: ", style = "Normal")

(tab <- table(train$pred >= 0.12, train$y))

##        
##           no  yes
##   FALSE 2690   74
##   TRUE   508  344

(`train accuracy` <- round((tab[1, 1] + tab[2, 2]) / sum(tab), 2))

## [1] 0.84

test$pred <- predict(glmmodel, newdata = test, type = "response")
(tab <- table(test$pred >= 0.12, test$y))

##        
##          no yes
##   FALSE 683  15
##   TRUE  119  88

(`test accuracy` <- round((tab[1, 1] + tab[2, 2]) / sum(tab), 2))

## [1] 0.85

accuracy <- t(cbind(`train accuracy`, `test accuracy`))
accuracy <- as.data.frame(cbind(rownames(accuracy), accuracy))

ns <- names(accuracy)
ft <- flextable(accuracy)
# to exclude headers in flextable
# set_header_labels(ft, V1 = NULL, V2 = NULL)
# names(df) <- paste0("X", 1:ncol(df))

# so we need to know in advance the column names but there's a trick 
# to apply it without knowing them
ex <- paste0('"', ns, '"', "=NULL", collapse = ",")
expr <- paste0('flextable::set_header_labels(ft,')
exp <- paste0(expr, ex, ")")
ft <- eval(parse(text = exp))

ft <- ft %>% bold(j = 1) %>% font(fontname = "Arial") %>%
  fontsize(size = 11) %>% border_remove() %>% autofit() %>%
  align(align = "left", part = "all") %>% 
  color(i = 2, j = 2, color = "darkblue") %>%
  bold(j = 2, i = 2)
flextable::body_add_flextable(docx, value = ft, align = "left")

print(docx, target = "report.docx")