#' ---	
#' title: "Data Science and Predictive Analytics (UMich HS650)"	
#' subtitle: "<h2><u>Apriori Association Rules Learning</u></h2>"	
#' author: "<h3>SOCR/MIDAS (Ivo Dinov)</h3>"	
#' 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	
#'       after_body: SOCR_footer_tracker.html	
#'     toc: true	
#'     number_sections: true	
#'     toc_depth: 2	
#'     toc_float:	
#'       collapsed: false	
#'       smooth_scroll: true	
#' ---	
#' 	
#' # Motivation	
#' 	
#' [HTTP cookies](https://en.wikipedia.org/wiki/HTTP_cookie) are used to track web-surfing the 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) 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, and Grocery purchases.	
#' 	
#' # 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. Those curly brackets indicate that we have a set. 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 discover the above pattern.	
#' 	
#' Association rules are commonly used for unsupervised discovery of knowledge rather than prediction of 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 fraudulent credit card or insurance use.	
#' 	
#' # The Apriori algorithm for association rule learning	
#' 	
#' Association rules are mostly applied to transactional data, which is usually transactional records like medical records. These datasets are typically very large in number of transactions and features. This will add lots of possible orders and patterns when we try to do a market basket analysis, which makes data mining a very hard task. 	
#' 	
#' With the **Apriori** rule, this problem is easily solved. If we have a simple prior (belief about the properties of frequent elements), we can efficiently reduce the number of features or combinations that we need to look at.	
#' 	
#' The Apriori algorithm has a simple `apriori` belief that *all subsets of a frequent item-set must also be frequent*. This is known as the **Apriori property**. For the last example $\{charcoal, \, lighter, \, chicken\, wings, \, barbecue\, sauce\}$, this full set can be frequent if and only if itself and all its subsets of single elements, pairs and triples occur frequently. We can see that this algorithm is designed for finding patterns in large datasets. If a pattern happens frequently, it is considered "interesting".	
#' 	
#' # Measuring rule importance by using **support** and **confidence**	
#' 	
#' Support and confidence are the two criteria to help us decide whether a pattern is "interesting". 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 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*. Of course, the union of item-sets is an item-set itself, i.e., if $Z={X,Y}$, then $$support(Z)=support(X,Y).$$	
#' 	
#' For a 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)}$$	
#' 	
#' This measures the joint occurrence of *X* and *Y* over the *X* domain. If whenever *X* appears *Y* tends to be present too, we will have a high $confidence(X\rightarrow Y)$. The ranges of the support and confidence are $0 \leq support,\ confidence \leq 1$.	
#' 	
#' 	
#' $\{peanut\, butter\}\rightarrow\{bread\}$ would be an example of a strong rule 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\}$.	
#' 	
#' # Building a set of rules with the Apriori principle	
#' 	
#' To build a set of rules, we need to go through two steps:	
#' 	
#' * **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, etc. Item-sets passing iteration *i* could be considered as candidates for the next iteration, *i+1*. If *{A}*, *{B}*, *{C}* are all frequent, 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("item1","item2","item3","item4")	
rownames(item_table) <- "support"	
kable(item_table,caption = "Size 1 Support")	
#' 	
#' 	
#' All the 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) <- "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) <- "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](http://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/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 an example and each column representing a feature.	
#' 	
#' Let's load the dataset and delete the irrelevant *index* column. With 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)	
#' 	
#' 	
#' 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 that is $211/528=40%$ of the (treatment) transactions. Since [fentanyl](https://en.wikipedia.org/wiki/Fentanyl) is frequently used to help prevent pain after surgery or other medical procedure, we can see that many of these patients were going through some painful medical procedures.	
#' 	
#' 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 still be fairly abstract, let's visualize the data.	
#' 	
#' 	
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.	
#' 	
#' 	
itemFrequency(med[, 1:5])	
itemFrequencyPlot(med, topN=20)	
#' 	
#' 	
#' The above graph is showing us 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$. 	
#' 	
#' ### Visualizing transaction data - plotting the sparse matrix	
#' 	
#' The sparse matrix will show what mediations were prescribed for each patient.	
#' 	
#' 	
image(med[1:5, ])	
#' 	
#' 	
#' This images has 5 rows (we only requested the first 5 patients) and 88 columns (88 different medicines). 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 100 randomly chosen patients.	
#' 	
#' 	
subset_int <- sample(nrow(med), 100, replace = F)  	
image(med[subset_int, ])	
#' 	
#' 	
#' 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 `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.1` 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 contains two items; 12 rules contains 3 items and the rest 4 rules contains 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 support since we already set a threshold that the support will exceed.	
#' 	
#' Using hte `arugleViz` package we can visualize the confidence and support scatter plots for all the rules.	
#' 	
#' 	
# install.packages("arulesViz")	
library(arulesViz)	
plot(sort(med_rule))	
#' 	
#' 	
#' 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 that resist bacterial infections); bacterial infections are associated with fevers and some cancers.	
#' 	
#' ## Step 5 - improving model performance	
#' 	
#' ### Sorting the set of association rules	
#' 	
#' Sorting the resulting association rules corresponding to high **lift** values will help us select 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. *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.	
#' 	
#' ### 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, since it appears to be the most frequently prescribed medicine.	
#' 	
#' 	
fi_rules<-subset(med_rule, items %in% "fentanyl injection uh")	
inspect(fi_rules)	
#' 	
#' 	
#' `%in%` means "belongs to" in R language. There are 14 rules related to this item. Let's plot them.	
#' 	
#' 	
plot(sort(fi_rules, by="lift"), method="grouped", control=list(type="items"), main = "Grouped Matrix for the 14 Fentanyl-associated Rules")	
#' 	
#' 	
#' ### 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.	
#' 	
#' # 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)	
#' 	
#' 	
#' 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.	
#' 	
#' # 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 verifications.	
#' * [See also Chapter 19 (Text Mining and NLP)](http://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/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/).