"
#' date: "`r format(Sys.time(), '%B %Y')`"
#' tags: [DSPA, SOCR, MIDAS, Big Data, Predictive Analytics]
#' output:
#' html_document:
#' theme: spacelab
#' highlight: tango
#' includes:
#' before_body: SOCR_header.html
#' toc: true
#' number_sections: true
#' toc_depth: 2
#' toc_float:
#' collapsed: false
#' smooth_scroll: true
#' code_folding: show
#' self_contained: yes
#' ---
#'
#' # Motivation
#'
#' [HTTP cookies](https://en.wikipedia.org/wiki/HTTP_cookie) are used to track web-surfing and Internet traffic. We often notice that promotions (ads) on websites tend to match our needs, reveal our prior browsing history, or may reflect our interests. That is not an accident. Nowadays, [recommendation systems](https://en.wikipedia.org/wiki/Recommender_system) are highly based on machine learning methods that can learn the behavior, e.g., purchasing patterns, of individual consumers. In this chapter, we will uncover some of the mystery behind recommendation systems like market basket analysis. Specifically, we will (1) discuss association rules and their support and confidence, (2) present the *Apriori algorithm* for association rule learning, and (3) cover step-by-step a set of case-studies, including a toy example, Head and Neck Cancer Medications, Grocery purchases, and survival of Titanic passengers.
#'
#' # Association Rules
#'
#' Association rules are the result of process analytics (e.g., market basket analysis) that specify patterns of relationships among items. One specific example would be:
#' $$\{charcoal, \, lighter, \, chicken\, wings\}\rightarrow\{barbecue\, sauce\}$$
#' In words, charcoal, lighter and chicken wings imply barbecue sauce. The curly brackets indicate that we have a set of items and the arrow suggest a direction of the association. Items in a set are called *elements.* When an item-set like $\{charcoal, \, lighter, \, chicken\, wings, \, barbecue\, sauce\}$ appears in our dataset with some regularity, we can mine and discover pattern of association with other item-sets.
#'
#' Association rules are commonly used for unsupervised discovery of knowledge rather than prediction of pre-specified outcomes. In biomedical research, association rules are widely used to:
#'
#' * searching for interesting or frequently occurring patterns of DNA,
#' * searching for protein sequences in an analysis of cancer data,
#' * finding patterns of medical claims that occur in combination with credit card or insurance fraud.
#'
#' # The Apriori algorithm for association rule learning
#'
#' Association rules are mostly applied to *transactional data*, which is usually records of trade, exchange, or arrangement, e.g., medical records. These datasets are typically very large in number of transactions and features, e.g., [electronic health record (EHR)](https://en.wikipedia.org/wiki/Electronic_health_record). In such data archives, the number and complexity of transactions is high, which complicates efforts to extract patterns, conduct market analysis, or predict basket purchases.
#'
#' **Apriori association rules** help untangle such difficult modeling problems. If we have a simple prior belief about the properties of frequent elements, we may be able to efficiently reduce the number of features or combinations that we need to look at.
#'
#' The Apriori algorithm is based on a simple `apriori` belief that *all subsets of a frequent item-set must also be frequent*. This is known as the **Apriori property**. In the last example, the full set $\{charcoal, \, lighter, \, chicken\, wings, \, barbecue\, sauce\}$ is *frequent* if and only if itself and all of its subsets, including single elements, pairs, and triples, occur frequently. Naturally, the apriori rule is designed for finding patterns in large datasets where patterns that appear frequently are considered "interesting", "valuable", or "important".
#'
#' # Rule **support** and **confidence**
#'
#' We can measure rule's importance by computing its **support** and **confidence** metrics. The support and confidence represent two criteria useful in deciding whether a pattern is "valuable". By setting thresholds for these two criteria, we can easily limit the number of interesting rules or item-sets reported.
#'
#' For item-sets $X$ and $Y$, the `support` of an item-set measures how (relatively) frequently it appears in the data:
#' $$support(X)=\frac{count(X)}{N},$$
#' where *N* is the total number of transactions in the database and *count(X)* is the number of observations (transactions) containing the item-set *X*.
#'
#' In a set-theoretic sense, the union of item-sets is an item-set itself. In other words, if $Z={X,Y}={X}\cup{Y}$, then $$support(Z)=support(X,Y).$$
#'
#' For a given rule $X \rightarrow Y$, the `rule's confidence` measures the relative accuracy of the rule:
#' $$confidence(X \rightarrow Y)=\frac{support(X, Y)}{support(X)}.$$
#'
#' The `confidence` measures the joint occurrence of *X* and *Y* over the *X* domain. If whenever *X* appears *Y* tends to also be present, then we will have a high $confidence(X\rightarrow Y)$.
#'
#' Note that the ranges of the support and the confidence are $0 \leq support,\ confidence \leq 1$.
#'
#'
#' One intuitive example of a strong association rule is $\{peanut\, butter\}\rightarrow\{bread\}$, because it has high *support* as well as high *confidence* in grocery store transactions. Shoppers tend to purchase bread when they get peanut butter. These items tend to appear in the same baskets, which yields high confidence for the rule $\{peanut\, butter\}\rightarrow\{bread\}$. Of course, there may be cases where $\{peanut\, butter\}\rightarrow\{celery\}$. These recommendation systems may not be perfect and can fail in unexpected ways, see the ["How Target Figured Out A Teen Girl Was Pregnant Before Her Father Did?" article](https://www.forbes.com/sites/kashmirhill/2012/02/16/how-target-figured-out-a-teen-girl-was-pregnant-before-her-father-did).
#'
#' # Building a set of rules with the Apriori principle
#'
#' Remember that the number of arrangements of $n$ elements taken $k$ at a time, i.e., the [number of combinations](https://en.wikipedia.org/wiki/Combination), increases exponentially with the size of the item inventory ($n$), see [this SOCR activity](http://wiki.stat.ucla.edu/socr/index.php/AP_Statistics_Curriculum_2007_Prob_Count). This is precisely why if a restaurant only uses 10 ingredients, there are ${10\choose 5}=253$ possible menu items for customers to order. clearly the complexity of the number of "baskets" rapidly increases with the inventory of the available items or ingredients. For instance,
#'
#' + $n=100$ (ingredients) and $k=50$ (menus of 50 ingredients), yields ${100\choose 50}=100891344545564193334812497256\gt 10^{29}$ possible orders, and
#' + $n=100$ (ingredients) and $k=5$ (menus of only 5 ingredients), yields ${100\choose 5}=75287520\gt 7M$ possible orders.
#'
#' To avoid this complexity, we will introduce a two-step process of building few, simple, and informative sets of rules:
#'
#' * **Step 1**: Filter all item-sets with a minimum *support* threshold. This is accomplished iteratively by increasing the size of the item-sets. In the first iteration, we compute the support of singletons, 1-item-sets. Next iteration, we compute the support of pairs of items, then, triples of items, etc. Item-sets passing iteration *i* could be considered as candidates for the next iteration, *i+1*. If *{A}*, *{B}*, *{C}* are all frequent singletons, but *D* is not frequent in the first singleton-selection round, then in the second iteration we only consider the support of these pairs *{A, B}*, *{A,C}*, *{B,C}*, ignoring all pairs including *D*. This substantially reduces the cardinality of the potential item-sets and ensures the feasibility of the algorithm. At the third iteration, if *{A,C}*, and *{B,C}* are frequently occurring, but *{A, B}* is not, then the algorithm may terminate, as the support of *{A,B,C}* is trivial (does not pass the support threshold), given that *{A, B}* was not frequent enough.
#'
#' * **Step 2**: Using the item-sets selected in step 1, generate new rules with *confidence* larger than a predefined minimum confidence threshold. The candidate item-sets that passed step 1 would include all frequent item-sets. For the highly-supported item-set *{A, C}*, we would compute the confidence measures for $\{A\}\rightarrow\{C\}$ as well as $\{C\}\rightarrow\{A\}$ and compare these against the minimum confidence threshold. The *surviving rules are the ones with confidence levels exceeding that minimum threshold*.
#'
#' # A toy example
#' Assume that a large supermarket tracks sales data by stock-keeping unit (SKU) for each item, i.e., each item, such as "butter" or "bread", is identified by an SKU number. The supermarket has a database of transactions where each transaction is a set of SKUs that were bought together.
#'
#' Suppose the database of transactions consist of following item-sets, each representing a purchasing order:
#'
#'
require(knitr)
item_table = as.data.frame(t(c("{1,2,3,4}","{1,2,4}","{1,2}","{2,3,4}","{2,3}","{3,4}","{2,4}")))
colnames(item_table) <- c("choice1","choice2","choice3","choice4","choice5","choice6","choice7")
kable(item_table, caption = "Item table")
#'
#'
#' We will use *Apriori* to determine the frequent item-sets of this database. To do so, we will say that an item-set is frequent if it appears in at least $3$ transactions of the database, i.e., the value $3$ is the support threshold.
#'
#' The first step of Apriori is to count up the number of occurrences, i.e., the support, of each member item separately. By scanning the database for the first time, we obtain get:
#'
#'
item_table = as.data.frame(t(c(3,6,4,5)))
colnames(item_table) <- c("items: {1}","{2}","{3}","{4}")
rownames(item_table) <- "(N=7)*support"
kable(item_table,caption = "Size 1 Support")
#'
#'
#' All the singletons, item-sets of size 1, have a support of at least 3, so they are all frequent. The next step is to generate a list of all pairs of frequent items.
#'
#' For example, regarding the pair $\{1,2\}$: the first table of Example 2 shows items 1 and 2 appearing together in three of the item-sets; therefore, we say that the support of the item $\{1,2\}$ is $3$.
#'
#'
item_table = as.data.frame(t(c(3,1,2,3,4,3)))
colnames(item_table) <- c("{1,2}","{1,3}","{1,4}","{2,3}","{2,4}","{3,4}")
rownames(item_table) <- "N*support"
kable(item_table,caption = "Size 2 Support")
#'
#'
#' The pairs $\{1,2\}$, $\{2,3\}$, $\{2,4\}$, and $\{3,4\}$ all meet or exceed the minimum support of $3$, so they are *frequent*. The pairs $\{1,3\}$ and
#' $\{1,4\}$ are not and any larger set which contains $\{1,3\}$ or $\{1,4\}$ cannot be frequent. In this way, we can prune sets: we will now look for frequent triples in the database, but we can already exclude all the triples that contain one of these two pairs:
#'
#'
item_table = as.data.frame(t(c(2)))
colnames(item_table) <- c("{2,3,4}")
rownames(item_table) <- "N*support"
kable(item_table,caption = "Size 3 Support")
#'
#'
#' In the example, there are no frequent triplets -- the support of the item-set $\{2,3,4\}$ is below the minimal threshold, and the other triplets were excluded because they were super sets of pairs that were already below the threshold. We have thus determined the frequent sets of items in the database, and illustrated how some items were not counted because some of their subsets were already known to be below the threshold.
#'
#' # Case Study 1: Head and Neck Cancer Medications
#'
#' ## Step 1 - collecting data
#'
#' To demonstrate the *Apriori* algorithm in a real biomedical case-study, we will use a transactional healthcare data representing [a subset of the Head and Neck Cancer Medication data](https://umich.instructure.com/files/1678540/download?download_frd=1), which it is available in [our case-studies collection](https://umich.instructure.com/courses/38100/files/folder/data) as `10_medication_descriptions.csv`. It consists of inpatient medications for head and neck cancer patients.
#'
#' The data is a wide format, see [Chapter 1](https://www.socr.umich.edu/people/dinov/courses/DSPA_notes/01_Foundation.html), where each row represents a patient. During the study period, each patient had records for a maximum of 5 encounters. *NA* represents no medication administration records in this specific time point for the specific patient. This dataset contains a total of 528 patients.
#'
#'
#subsetting patients with 2-5 encounters
in_medication<-read.csv("https://umich.instructure.com/files/1678540/download?download_frd=1", header=T)
library(plyr)
a<-count(inmedwpid$PID)
b<-a[a$freq %in% c(2:5), ]
hn_med1<-inmedwpid[inmedwpid$PID%in% b$x, ]
table(hn_med1$PID)
# To save the results locally
# write.csv(hn_med1, "D:\\folder\\hn_med1.csv")
#'
#'
#'
#make a csv file that is read.transactions() friendly.
# hn_med1<-read.csv("D:\\folder\\hn_med1.csv")
hn_med1<-hn_med1[, -1]
hn_med1<-hn_med1[, c(3, 9)]
med.l<-count(hn_med1, c("PID", "MEDICATION_DESC"))
m_count<-count(med.l$PID)
enc<-c()
for (i in 1:length(m_count$x)){
for (j in 1:m_count$freq[i]){
enc=c(enc, j)
}
}
med.l<-cbind(med.l, enc)
med.l$MEDICATION_DESC<-gsub("[[:punct:]]", " ", med.l$MEDICATION_DESC)
med.l$MEDICATION_DESC<-casefold(med.l$MEDICATION_DESC, upper=F)
med.w<-reshape(med.l, timevar="enc", idvar =c("PID", "freq"), direction="wide")
med.w<-med.w[, -c(1, 2)]
write.csv(med.w, " med.w.csv")
#'
#'
#' ## Step 2 - exploring and preparing the data
#'
#' Different from our data imports in the previous chapters, transactional data need to be ingested in R using the `read.transactions()` function. This function will store data as a matrix with each row representing a basket (or transaction example) and each column representing items in the transaction.
#'
#' Let's load the dataset and delete the irrelevant *index* column. Using the `write.csv(R data, "path")` function, we can output our R data file into a local CSV file. To avoid generating another index column in the output CSV file, we can use the `row.names=F` option.
#'
#'
med<-read.csv("https://umich.instructure.com/files/1678540/download?download_frd=1", stringsAsFactors = FALSE)
med<-med[, -1]
write.csv(med, "medication.csv", row.names=F)
library(knitr)
kable(med[1:5, ])
#'
#'
#' Now we can use `read.transactions()` in the `arules` package to read the CSV file we just outputted.
#'
#'
# install.packages("arules")
library(arules)
med<-read.transactions("medication.csv", sep = ",", skip = 1, rm.duplicates=TRUE)
summary(med)
#'
#'
#' Here we use the option `rm.duplicates=T` because we may have similar medication administration records for two different patients. The option `skip=1` means we skip the heading line in the CSV file. Now we get a transactional data with unique rows.
#'
#' The summary of a transactional data contains rich information. The first block of information tells us that we have 528 rows and 88 different medicines in this matrix. Using the density number we can calculate how many non *NA* medication records are in the data. In total, we have $528\times 88=46,464$ positions in the matrix. Thus, there are $46,464\times 0.0209=971$ medicines prescribed during the study period.
#'
#' The second block lists the most frequent medicines and their frequencies in the matrix. For example `fentanyl injection uh` appeared 211 times, which represents a proportion of $211/528=0.4$ of the (treatment) transactions. Since [fentanyl](https://en.wikipedia.org/wiki/Fentanyl) is frequently used to help prevent pain after surgery or other medical treatments, we can see that many of these patients may have undergone some significant medical procedures that require post-operative pain management.
#'
#' The last block shows statistics about the size of the transaction. 248 patients had only one medicine in the study period, while 12 of them had 5 medication records one for each time point. On average, the patients are having 1.8 different medicines.
#'
#' ### Visualizing item support - item frequency plots
#'
#' The summary might may still be fairly abstract, so we can visualize the information.
#'
#'
inspect(med[1:5,])
#'
#'
#' The `inspect()` call shows the transactional dataset. We can see that the medication records of each patient are nicely formatted as item-sets.
#'
#' We can further analyze the frequent terms using `itemFrequency()`. This will show all item frequencies alphabetically ordered from the first five outputs.
#'
#'
library(plotly)
itemFrequency(med[, 1:5])
# itemFrequencyPlot(med, topN=20)
# convert the data to a matrix
mat <- as.matrix(med@data)
# View matrix as image
# image(mat)
# capture the row/column names
rowNames <- med@itemInfo$labels
colNames <- paste0("S", c(1:dim(mat)[2]))
rownames(mat) <- rowNames
colnames(mat) <- colNames
# convert matrix to DF for processing, order rows based on their average (across subjects/cases/columns), back to matrix for display
df <- as.data.frame(1*mat)
df$avg <- rowMeans(df)
dfOrdered <- df[order(df$avg, decreasing = T), ]
matOrdered <- as.matrix(dfOrdered)
# track the ordered row names
rowNames <- rownames(dfOrdered)
colNames <- colnames(dfOrdered)
# 2D top 20 terms bar plot
# To order the meds based on "avg", instead of alphabetically (mind the "-" sign to order Large to small!)
plot_ly(x = reorder(rowNames[c(1:20)], -dfOrdered[1:20, "avg"]), y=dfOrdered[1:20, "avg"], name="Top 20 Meds", type="bar") %>%
layout(title='Frequency of Medications (Top 20) based on averaging across Cases',
xaxis = list(title="Term"),
yaxis = list(title="Relative Frequency"))
# 3D surface plot
plot_ly(x = colNames, y = rowNames, z = 2*matOrdered, type = "surface") %>%
layout(title='Term (X) by Sample (Y) Frequency (Z) Plot',
xaxis = list(title="Term"),
yaxis = list(title="Sample ID")) %>% hide_colorbar()
#'
#'
#' We can only display the top 20 medicines that are most frequently present in this dataset. Consistent with the prior `summary()` output, `fentanyl` is still the most frequent item. You can also try to plot the items with a threshold for support. Instead of `topN=20`, just use the option `support=0.1`, which will give you all the items have a support greater or equal to $0.1$.
#'
#' Let's generalize this process and define a new function, `itemFrequencyPlotly(transactionObject, numTopItemps = 10)`, that can generate frequency plots for any transaction object.
#'
#'
# define a generic plot_ly ItemFrequency plotting function
itemFrequencyPlotly <- function(transactionObject, numTopItemps = 10, name="") {
name <- ifelse(name=="",
paste0('Frequency of Items (Top ', numTopItemps,
' ) based on averaging across Cases'),
paste0('Frequency of Items (Top ', numTopItemps,
' ) based on averaging across Cases (Data=', name, ')'))
mat <- as.matrix(transactionObject@data)
# View matrix as image
# image(mat)
# capture the row/column names
rowNames <- transactionObject@itemInfo$labels
colNames <- paste0("S", c(1:dim(mat)[2]))
rownames(mat) <- rowNames
colnames(mat) <- colNames
# convert matrix to DF for processing, order rows based on their average (across subjects/cases/columns), back to matrix for display
df <- as.data.frame(1*mat)
df$avg <- rowMeans(df)
dfOrdered <- df[order(df$avg, decreasing = T), ]
matOrdered <- as.matrix(dfOrdered)
# track the ordered row names
rowNames <- rownames(dfOrdered)
colNames <- colnames(dfOrdered)
# 2D top "20"numTopItemps" terms bar plot
# To order the meds based on "avg", instead of alphabetically (mind the "-" sign to order Large to small!)
plot_ly(x = reorder(rowNames[c(1:numTopItemps)], -dfOrdered[1:numTopItemps, "avg"]),
y=dfOrdered[1:numTopItemps, "avg"], name=paste0("Top ", numTopItemps, " Meds"), type="bar") %>%
layout(title=name,
xaxis = list(title="Terms"),
yaxis = list(title="Relative Frequency"))
}
#'
#'
#' ### Visualizing transaction data - plotting the sparse matrix
#'
#' The sparse matrix will show what medications were prescribed for each patient. Below we only show the top-20 medications for the first 15 cases.
#'
#'
# image(med[1:5, ])
plot_ly(x=reorder(rowNames[c(1:20)], -dfOrdered[1:20, "avg"]), y=colNames[1:15],
z=2*matOrdered[1:15, 1:20], type="heatmap") %>%
layout(title='Heatmap - Top-20 Medications for the first 15 Cases') %>% hide_colorbar()
#'
#'
#' This images has 15 rows (translations), as we only specified the first 15 patients) and 20 columns (top-20 meds). Although the picture may be a little hard to interpret, it gives a sense of what kind of medicine is prescribed for each patient in the study.
#'
#' Let's see an expanded graph including the top 30 medications for a random roster of 50 patients.
#'
#'
subset_int <- sample(ncol(matOrdered), 50, replace = F)
# image(med[subset_int, ])
plot_ly(x=rowNames[1:30], y=colNames[subset_int],
z=2*t(matOrdered[1:30, subset_int]), type="heatmap") %>%
layout(title='Heatmap - Bottom-30 Medications for a Random set of 50 Patients') %>% hide_colorbar()
#'
#'
#' It shows us clearly that some medications are more popular than others. Now, let's fit the *Apriori* model.
#'
#' ## Step 3 - training a model on the data
#'
#' With the data in place, we can build the *association rules* using the `arules::apriori()` function.
#'
#' `myrules <- apriori(data=mydata, parameter=list(support=0.1, confidence=0.8, minlen=1))`
#'
#' * *data*: a sparse matrix created by `read.transacations()`.
#' * *support*: minimum threshold for support.
#' * *confidence*: minimum threshold for confidence.
#' * *minlen*: minimum required rule items (in our case, medications).
#'
#' Setting up the threshold could be hard. You don't want it to be too high so that you get no rules or rules that everyone knows. You don't want to set it too low either, to avoid too many rules present. Let's see what we get under the default setting `support=0.1, confidence=0.8`:
#'
#'
apriori(med)
#'
#' Not surprisingly, we have 0 rules. The default setting is too high. In practice, we might need some time to fine-tune these thresholds, which may require certain familiarity with the underlying process or clinical phenomenon.
#'
#' In this case study, we set `support=0.01` and `confidence=0.25`. This requires rules that have appeared in at least 10% of the head and neck cancer patients in the study. Also, the rules have to have least 25% accuracy. Moreover, `minlen=2` would be a very helpful option because it removes all rules that have fewer than two items.
#'
#'
med_rule<-apriori(med, parameter=list(support=0.01, confidence=0.25, minlen=2))
med_rule
#'
#'
#' The result suggest we have a new `rules` object consisting of 29 rules.
#'
#' ## Step 4 - evaluating model performance
#'
#' First, we can obtain the overall summary of this set of rules.
#'
#'
summary(med_rule)
#'
#'
#' We have 13 rules that contain two items; 12 rules containing 3 items, and the remaining 4 rules contain 4 items.
#'
#' The `lift` column shows how much more likely one medicine is to be prescribed to a patient given another medicine is prescribed. It is obtained by the following formula:
#' $$lift(X\rightarrow Y)=\frac{confidence(X\rightarrow Y)}{support(Y)}$$
#' Note that $lift(X\rightarrow Y)$ is the same as $lift(Y\rightarrow X)$. The range of $lift$ is $[0,\infty)$ and higher $lift$ is better. We don't need to worry about the support, since we already set a threshold that the support must exceed.
#'
#' Using the `arulesViz` package we can visualize the confidence and support scatter plots for all the rules.
#'
#'
# install.packages("arulesViz")
library(arulesViz)
# plot(sort(med_rule))
sortedRule <- sort(med_rule)
x1 <- sortedRule@quality$support
y1 <- sortedRule@quality$confidence
z1 <- sortedRule@quality$lift
col1 <- sortedRule@quality$count
ruleNames <- paste0("Rule", c(1:length(sortedRule@quality$support)))
plot_ly(x = ~x1, y = ~y1, z = ~z1, color = ~z1, name=ruleNames) %>%
add_markers() %>%
layout(title=paste0("Arule Support-Confidence-Lift Plot (for all ", length(sortedRule@quality$support), " rules)"),
scene = list(xaxis = list(title = 'Support'),
yaxis = list(title = 'Confidence'),
zaxis = list(title = 'Lift'))) %>%
hide_colorbar()
#'
#'
#' Again, we can utilize the `inspect()` function to see exactly what are these rules.
#'
#'
inspect(med_rule[1:3])
#'
#'
#' Here, `lhs` and `rhs` refer to "left hand side" and "right hand side" of the rule, respectively. `lhs` is the given condition and `rhs` is the predicted result. Using the first row as an example: If a head-and-neck patient has been prescribed acetaminophen (pain reliever and fever reducer), it is likely that the patient is also prescribed *cefazolin* (antibiotic prescribed for treatment of resistant bacterial infections); bacterial infections are associated with fevers and some cancers.
#'
#' ## Step 5 - sorting the set of association rules
#'
#' Sorting the resulting association rules corresponding to high **lift** values will help us identify the most useful rules.
#'
#'
inspect(sort(med_rule, by="lift")[1:3])
#'
#'
#' These rules may need to be interpreted by clinicians and experts in the specific context of the study. For instance, the first row, *{fentanyl, heparin, hydrocodone acetaminophen}* implies *{cefazolin}*. *Fentanyl* and hydrocodone acetaminophen are both pain relievers that may be prescribed after surgery to relieve moderate to severe pain based on an narcotic opioid pain reliever (*hydrocodone*) and a non-opioid pain reliever (acetaminophen). *Heparin* is usually used before surgery to reduce the risk of blood clots. This rule may suggest patients that have undergone surgical treatments and are likely that they will need *cefazolin* to prevent post-surgical bacterial infection. *Cefazolin* is an antibiotic used for the treatment of bacterial infections and also to prevent group *B streptococcal disease* around the time of delivery and before general surgery.
#'
#' ## Step 6 - taking subsets of association rules
#'
#' If we are more interested in investigating associations that are linked to a specific medicine, we can narrow the rules down by making subsets. Let us try investigating rules related to [fentanyl](https://en.wikipedia.org/wiki/Fentanyl), since it appears to be the most frequently prescribed medicine. [Fentanyl](https://en.wikipedia.org/wiki/Fentanyl) is used in the management of chronic cancer pain.
#'
#'
fi_rules<-subset(med_rule, items %in% "fentanyl injection uh")
inspect(fi_rules)
#'
#'
#' Earlier, we saw that $\%in\%$ is a simple intuitive interface to `match()` and is used as a binary operator returning a logical vector indicating if there is a match (T) or not (F) for its left operand within the right table object. In total, there are 14 rules related to this item. Let's plot them.
#'
#' ## Graphical depiction of association rules
#'
#'
# ?arulesViz::plot()
# plot(sort(fi_rules, by="lift"), method="grouped", control=list(type="items"),
# main = "Grouped Matrix for the 14 Fentanyl-associated Rules")
# plot(sort(fi_rules, by="lift"), method = "graph", engine = "plotly")
# plot(sort(fi_rules, by="lift"), method = "graph", engine = "igraph")
# plot(sort(fi_rules, by="lift"), method = "graph", engine = "paracoord")
plot(sort(fi_rules, by="lift"), method = "graph", engine = "htmlwidget",
control=list(main = list(title="Grouped Matrix for the 14 Fentanyl-associated Rules")))
subrules2 <- sample(subset(fi_rules, lift > 2), 5)
plot(sort(subrules2, by="lift"), method="grouped", control=list(type="items"), engine = "htmlwidget")
plot(sort(fi_rules, by="lift"), method="grouped", k = 7, control=list(type="items"), engine = "htmlwidget")
m <- rules2matrix(sort(fi_rules, by="lift")[1:10], measure = "lift")
# m
plot(fi_rules[1:10], method = "matrix", engine = "htmlwidget")
## Grouped matrix
# create a matrix with LHSs grouped in k = 10 groups
m <- rules2groupedMatrix(sort(fi_rules, by="lift")[1:10], k = 10)
# m$m
# number of rules per group
# table(m$clustering_rules)
# get rules for group 1
inspect(fi_rules[m$clustering_rules == 1])
# the corresponding plot
plot(fi_rules, method = "grouped matrix", k = 7, engine = "htmlwidget")
# plot(fi_rules, method = "grouped matrix", k = 7)
#### For interactive runs try the RShinyApp:
# ruleExplorer(fi_rules)
#'
#'
#'
plot(fi_rules, method="graph", measure = "support", engine="htmlwidget", # nodeCol=rainbow(14),
shading = "lift", control = list(verbose = TRUE))
#'
#'
#' ## Saving association rules to a file or data frame
#'
#' We can save these rules into a CSV file using `write()`. It is similar with the function `write.csv()` that we have mentioned in the beginning of this case study.
#'
#'
write(med_rule, file = "medrule.csv", sep=",", row.names=F)
#'
#'
#' Sometimes it is more convenient to convert the rules into a data frame.
#'
#'
med_df <- as(med_rule, "data.frame")
str(med_df)
#'
#'
#' As we can see, the rules are converted into a factor vector. Finally, remember that matrices and data-frames can be converted to `arules` *transactions* format using `arulesObject <- as(input.df, "transactions")`.
#'
#' # Practice Problems
#'
#' ## Groceries
#'
#' In this practice problem, we will investigate the associations of frequently purchased groceries using the *grocery* dataset in the R base. Firstly, let's load the data.
#'
#'
data("Groceries")
summary(Groceries)
#'
#'
#' We will try to find out the top 5 frequent grocery items and plot them.
#'
#'
# itemFrequencyPlot(Groceries, topN=5)
itemFrequencyPlotly(Groceries, 10, "grocieries")
#'
#'
#' Then, try to use `support = 0.006, confidence = 0.25, minlen = 2` to set up the grocery association rules. Sort the top 3 rules with highest lift.
#'
#'
groceryrules <- apriori(Groceries, parameter = list(support =
0.006, confidence = 0.25, minlen = 2))
groceryrules
inspect(sort(groceryrules, by = "lift")[1:3])
#'
#'
#' The number of rules ($463$) appears excessive. We can try stringer parameters. In practice, it's more possible to observe underlying rules if you set a higher confidence. Here we set the $confidence=0.6$.
#'
#'
groceryrules <- apriori(Groceries, parameter = list(support = 0.006, confidence = 0.6, minlen = 2))
groceryrules
inspect(sort(groceryrules, by = "lift")[1:3])
#'
#'
#' We observe mainly rules between dairy products. It makes sense that customers pick up milk when they walk down the dairy products isle. Experiment further with various parameter settings and try to interpret the results in the context of this grocery case-study.
#'
#' ## Titanic Passengers
#'
#' Next we'll use the [Titanic Passengers Dataset](https://umich.instructure.com/courses/38100/files/folder/data/16_TitanicPassengerSurvivalDataset). Let's start by loading the data first.
#'
#'
dat <- read.csv("https://umich.instructure.com/files/9372716/download?download_frd=1")
# Choose only the key data features
dat <- dat[ , c("pclass", "survived", "sex", "age", "fare", "cabin")]
# Factorize categorical features
dat$pclass <- as.factor(dat$pclass)
dat$survived <- factor(dat$survived, levels = c("0", "1"))
dat$sex <- as.factor(dat$sex)
# Convert the Cabin number to a character A-F, Z=missing cabin ID
dat$cabin <- substring(dat$cabin,1,1)
for (i in 1:length(dat$cabin))
if ((dat$cabin[i]=="")) dat$cabin[i] <- "Z"
dat$cabin <- as.factor(dat$cabin)
# Convert the Ticket Fair from numeric to categorical label
f <- as.character(dat$fare)
for (i in 1:length(dat$fare)) {
if (is.na(dat$fare[i])) f[i] <- as.character("low")
else if (dat$fare[i]<50) f[i] <- as.character("low")
else if (50<=dat$fare[i] && dat$fare[i]<100) f[i] <- as.character("medium")
else if (100<=dat$fare[i] && dat$fare[i]<200) f[i] <- as.character("high")
else f[i] <- as.character("extreme") # if (200<=dat$fare[i])
}
dat$fare <- as.factor(f)
table(as.factor(dat$fare))
# Convert Age from numeric to categorical (Decade-of-life) label
f <- as.character(dat$age)
for (i in 1:length(dat$age)) {
if (is.na(dat$age[i])) f[i] <- as.character("1")
else {
a = 1 + dat$age[i] %/% 10 # integer division by 10 (per decade of life)
f[i] <- as.character(a)
}
}
dat$age <- as.factor(f)
table(as.factor(dat$age))
str(dat)
#'
#'
#' Next, we can mine the association rules.
#'
#'
library(arules)
rules <- apriori(dat, parameter = list(minlen=5, supp=0.02, conf=0.8))
inspect(rules[1:20])
#'
#'
#' We can focus on the binary *survival* outcome (Survived=1/0).
#'
#'
# examine the rules with containing "Survived=1" in the RHS
rules <- apriori(dat,
parameter = list(minlen=3, supp=0.02, conf=0.7),
appearance = list(rhs="survived=1", default="lhs"), control = list(verbose=F))
rules.sorted <- sort(rules, by="lift")
inspect(head(rules.sorted, 30))
#'
#'
#' Prune any redundant association rules. For instance, some rules may provide no extra knowledge. Rules that are highly related with prior rules may be redundant. Pruning reduces the number of rules from *27 to 18*.
#'
#'
# search for redundant rules
rules_lift <- sort(rules, by = 'lift')
rules_pruned <- rules_lift[!is.redundant(rules_lift, measure="lift")]
inspect(head(rules_pruned, 30))
#'
#'
#' The package `arulesViz` supplies rule-visualization routines using scatter, bubble, and parallel coordinates plots. The visualization of the frequent itemsets using parallel coordinates allows all items to be placed on the vertical axis with their position determined by their group and by the frequency of the item in descending order, vertical ranking based the support of the 1-itemset containing the item, [see this paper](https://link.springer.com/chapter/10.1007/978-3-030-04921-8_12). To render the maximal frequent itemsets, all subsets are implicitly drawn as sub-segments of one polyline. All itemsets are visualized as a parallel coordinate polylines. The *support* can be mapped to line color or width. The visualization effectively . Itemsets that share common parts may have overlapping polylines.
#'
#'
#'
library(arulesViz)
plot(rules)
# scatter plot of association rules
plot(rules_pruned, method="graph", control=list(verbose = FALSE), engine="htmlwidget")
# display association rules
# plot(rules_pruned, method="paracoord", control=list(reorder=TRUE, verbose = FALSE), engine="htmlwidget")
# interactive association network graph
plot(rules_pruned, method="graph", measure = "support", engine="htmlwidget", nodeCol=rainbow(14),
shading = "lift", control = list(verbose = FALSE))
# <- rules2matrix(sort(rules_pruned, by="lift")[1:10], measure = "lift")
# m
# plot(rules_pruned, method = "matrix", engine = "htmlwidget")
## Grouped matrix
# create a matrix with LHSs grouped in k = 10 groups
# m <- rules2groupedMatrix(sort(rules_pruned, by="lift")[1:10], k = 10)
# m$m
# number of rules per group
# table(m$clustering_rules)
# get rules for group 1
inspect(rules_pruned[m$clustering_rules == 1])
# the corresponding plot
plot(rules_pruned, method = "grouped matrix", k = 7, engine = "htmlwidget")
#'
#'
#' Consider plotting *all rules* (not only the pruned list of association rules) and examining rules predicting passenger death (use `"survived=0"`).
#'
#'
#' # Summary
#'
#' * The Apriori algorithm for association rule learning is only suitable for large transactional data. For some small datasets, it might not be very helpful.
#' * It is useful for discovering associations, mostly in early phases of an exploratory study.
#' * Some rules can be built due to chance and may need further verification.
#' * [See also Chapter 19 (Text Mining and NLP)](https://www.socr.umich.edu/people/dinov/courses/DSPA_notes/19_NLP_TextMining.html).
#'
#' Try to replicate these results with [other data from the list of our Case-Studies](https://umich.instructure.com/courses/38100/files/).
#'
#'
#'