The four pillars:

  • Descriptive Analytics
  • Diagnostic Analytics
  • Predictive Analytics
  • Prescriptive Analytics

Prescriptive Analytics | What should we do?

Prescriptive analytics seeks to provide a data-driven course of action, that is, recommendations that are likely to result in the best outcome as determined from the data. These recommendations may take the form of decision rules, prioritized lists of actions, or specific prescriptions for how to allocate resources or make trade-offs.

Identifying anomalies

What days had an unusually high number of downtimes?

Let’s make use of our prediction models to identify days with unusual downtimes.

tb_subset = tb %>% 
  filter(day >= as.Date("2020/11/01")) %>% 
  mutate(uid = paste0(plant, "_", machine_id))

tb_formatted_data = tb_subset %>% 
  group_by(day, plant) %>% 
  summarize(tot_dt = sum(dt), dt_instances = n()) %>% 
  ungroup() %>% 
  group_by(plant) %>% 
  complete(day = seq(min(tb_subset$day), max(tb_subset$day), by = 'day')) %>% 
  mutate(wk = epiweek(day)) %>% 
  mutate(yr = year(day)) %>% 
  mutate(wid = paste0(yr, "_", wk)) %>% 
  replace_na(list(tot_dt = 0, dt_instances = 0)) %>%
  mutate(wday = weekdays(day))

tb_reg_data = tb_formatted_data %>% 
  dummy_cols(select_columns = c("plant", "wday")) %>% 
  select(dt_instances, tot_dt, day, wid, plant, wday, 9:18) %>% 
  mutate(l1_dt_instances = lag(dt_instances, 1)) %>% 
  mutate(l7_dt_instances = lag(dt_instances, 7)) %>% 
  drop_na() %>% 
  filter(dt_instances > 5)

cat_nms = names(tb_reg_data)[7:16]

cat_vars = expand_grid(plant = cat_nms[1:3], wday = cat_nms[4:10]) %>% 
  mutate(vars = paste0(plant, ":", wday)) %>% 
  pull(vars) %>% 
  paste0(collapse = " + ")

tot_inst_formula = paste0(
  "dt_instances ~ ", cat_vars, "+ l1_dt_instances + l7_dt_instances - 1"
) %>% as.formula()

mod_inst = lm(tot_inst_formula, tb_reg_data)

tb_inst_results = tb_reg_data %>% 
  mutate(residuals = mod_inst$residuals)

tb_inst_results %>% 
  mutate(formatted_date = paste0(wday, ", ", format.Date(day))) %>% 
  plot_ly(
    color = ~plant, x = ~residuals, text = ~formatted_date, type = "box", 
    boxpoints = "all", quartilemethod = "inclusive",
    pointpos = 0, jitter = 0.5, alpha = 0.5
  )

What days had unusually high total downtimes?

tb_reg_data2 = tb_formatted_data %>% 
  dummy_cols(select_columns = c("plant", "wday")) %>% 
  select(dt_instances, tot_dt, day, wid, plant, wday, 9:18) %>%  
  mutate(l1_tot_dt = lag(tot_dt, 1)) %>% 
  mutate(l7_tot_dt = lag(tot_dt, 7)) %>% 
  drop_na() %>% 
  filter(dt_instances > 5)

tot_dt_formula = paste0(
  "tot_dt ~ ", cat_vars, "+ l1_tot_dt + l7_tot_dt - 1"
) %>% as.formula()

mod_tot_dt = lm(tot_dt_formula, tb_reg_data2)

tb_tot_dt_results = tb_reg_data2 %>% 
  mutate(residuals = mod_tot_dt$residuals)

tb_tot_dt_results %>% 
  mutate(formatted_date = paste0(wday, ", ", format.Date(day))) %>% 
  plot_ly(
    color = ~plant, x = ~residuals, text = ~formatted_date, type = "box", 
    boxpoints = "all", quartilemethod = "inclusive",
    pointpos = 0, jitter = 0.5, alpha = 0.5
  )
# Grab all outliers above the mean
plants = tb_inst_results$plant %>% unique()
tb_date_outliers = lapply(plants, function(p) {
  inst_vec = tb_inst_results %>% 
    filter(plant == p) %>% 
    pull(residuals)
  
  totd_vec = tb_tot_dt_results %>% 
    filter(plant == p) %>% 
    pull(residuals)
  
  inst_out = boxplot.stats(inst_vec)$out[boxplot.stats(inst_vec)$out > mean(inst_vec)]
  totd_out = boxplot.stats(totd_vec)$out[boxplot.stats(totd_vec)$out > mean(totd_vec)]
  
  tmp1 = tb_inst_results %>% 
    filter(plant == p) %>% 
    filter(residuals %in% inst_out) %>% 
    select(day, wday, plant, dt_instances:tot_dt)
  
  tmp2 = tb_tot_dt_results %>% 
    filter(plant == p) %>% 
    filter(residuals %in% totd_out) %>% 
    select(day, wday, plant, dt_instances:tot_dt)
  
  bind_rows(tmp1, tmp2) %>% 
    mutate(day = format.Date(day))
}) %>% bind_rows() %>% 
  distinct()

tb_date_outliers %>% 
  DT::datatable()

What machines had unusually high downtime instances & downtimes?

We’ll use our behavior models from the previous section to determine outliers.

First let’s look at downtime instances.

unique_states = tb_subset$reason_code_id %>% unique()
unique_machines = tb_subset %>% 
  group_by(uid) %>% 
  summarize(n = n()) %>% 
  filter(n > 100) %>% 
  pull(uid)

tm_list = lapply(unique_machines, function(u_id) {
  state_seq = tb_subset %>% 
    filter(uid == u_id) %>% 
    pull(reason_code_id)
  
  mc_fit = markovchainFit(data = state_seq, possibleStates = unique_states)
  
  return(mc_fit$estimate@transitionMatrix)
}) %>% setNames(unique_machines)

# Calculate the distance matrix
tb_machine_dist = expand_grid(row = unique_machines, col = unique_machines) %>% 
  rowwise() %>% 
  mutate(d = norm(tm_list[[row]] - tm_list[[col]], 'F'))

# Create the distance matrix
M = tb_machine_dist %>%
  pivot_wider(names_from = col, values_from = d) %>%
  select(-row) %>%
  as.matrix() %>%
  magrittr::set_rownames(colnames(.))


# Get the hierarchical tree
M_dist = as.dist(M)

M_tree = hclust(M_dist)

M_dend = as.dendrogram(M_tree)

K = 4

clusters    = cutree(M_dend, k=K, order_clusters_as_data = FALSE)
tb_clusters = tibble(uid = names(clusters), cluster = clusters)

uid_order = tb_clusters$uid

tb_per_machine = tb_subset %>%
  group_by(uid, plant, category, description) %>%
  summarize(dt_instances = n(), avg_dt = mean(dt), med_dt = median(dt), tot_dt = sum(dt)) %>%
  filter(dt_instances > 100) %>%
  ungroup() %>%
  left_join(tb_clusters, by = 'uid') %>%
  mutate(cluster = factor(cluster))

tb_per_machine %>%
  mutate(detail_description = paste0(description, " (", category, ")")) %>% 
  mutate(cluster = fct_reorder(cluster, tot_dt)) %>%
  plot_ly(
    color = ~cluster, x = ~dt_instances, type = "box", 
    boxpoints = "all", text = ~detail_description, quartilemethod = "inclusive",
    pointpos = 0, jitter = 0.5, alpha = 0.5
  ) %>% 
  layout(
    xaxis = list(title = "downtime instances"),
    yaxis = list(title = "behavior cluster"),
    title = "downtime instances per machine"
  )

Next, let’s look at total downtime.

tb_per_machine %>%
  mutate(detail_description = paste0(description, " (", category, ")")) %>% 
  mutate(cluster = fct_reorder(cluster, tot_dt)) %>%
  mutate(tot_dt = tot_dt/(60*24)) %>% 
  plot_ly(
    color = ~cluster, x = ~tot_dt, type = "box", 
    boxpoints = "all", text = ~detail_description, quartilemethod = "inclusive",
    pointpos = 0, jitter = 0.5, alpha = 0.5
  ) %>% 
  layout(
    xaxis = list(title = "total downtime (days)"),
    yaxis = list(title = "behavior cluster"),
    title = "total downtime per machine"
  )

Let’s put this information together.

# Get all anomalies (outliers above the mean, cluster 1)
tb_anomalies = lapply(1:4, function(c) {
  tmp = tb_per_machine %>% 
    filter(cluster == c)
  
  out_inst = boxplot.stats(tmp$dt_instances)$out
  out_inst = out_inst[out_inst > mean(tmp$dt_instances)]
  
  out_totd = boxplot.stats(tmp$tot_dt)$out
  out_totd = out_totd[out_totd > mean(tmp$tot_dt)]
  
  if(c == 1) {
    tmp
  } else {
    tmp %>% 
      filter(dt_instances %in% out_inst | tot_dt %in% out_totd)
  }
}) %>% bind_rows() %>% 
  mutate(anomaly = T)

plants = c("all", tb_subset$plant %>% unique())

tb_tmp = lapply(plants, function(p) {
  
  tmp = tb_subset %>% 
    mutate(uid = paste0(plant, "_", machine_id)) %>% 
    group_by(uid, category, plant) %>% 
    summarize(tot_dt = sum(dt)/(60*24), dt_inst = n()) %>% 
    mutate(group = p)
  
  # To help in redraw of visualization
  dummy_data = tibble(category = tb$category %>% unique()) %>% 
    mutate(tot_dt = 0) %>% 
    mutate(dt_inst = 0) %>% 
    mutate(plant = p)
  
  if (p != "all") {
    tmp = tmp %>% filter(plant == p)
  }
  
  tmp %>% 
    select(-plant) %>% 
    rename(plant = group) %>% 
    bind_rows(dummy_data) %>% 
    return()
}) %>% bind_rows() %>% 
  left_join(tb_anomalies %>% select(uid, anomaly), by = "uid") %>% 
  replace_na(list(anomaly = F)) %>% 
  left_join(tb_per_machine %>% select(uid, description, cluster), by = "uid") %>% 
  mutate(
    detail_description = paste0(
      description, " (", category, ")(cluster = ", cluster, ")"
      )
    )

tb_tmp %>% 
  plot_ly(x = ~tot_dt, y = ~dt_inst, text = ~detail_description) %>% 
  add_markers(data = tb_tmp, frame = ~plant, color = ~anomaly, alpha = 0.7) %>% 
  layout(
    xaxis = list(
      title = "total downtime (days)"
      ),
    yaxis = list(
      title = "total downtime instances"
      ),
    title = "identifying anomalous behavior in machines"
  ) %>% 
  animation_button(visible = F) %>% 
  animation_opts(redraw = F, transition = 0) %>% 
  animation_slider(currentvalue = list(font = list(size = 0, color = "white")))
tb_anomalies %>% 
  select(-anomaly) %>% 
  DT::datatable()