Association Rules and Anomaly Detection

These are two different unsupervised tasks. Association rules examine which items occur together; anomaly detection assigns high scores to observations that look unusual under a model. Neither directly proves disease or fraud.

1. Association rule mining

Suppose a fictional dataset has binary indicators for services A, B, and C on each visit. A rule A → B says: among rows with A, how often was B also present? It does not mean A caused B.

services <- data.frame(
  A = c(1,1,1,1,0,0,1,0,1,0),
  B = c(1,1,1,0,1,0,0,0,1,1),
  C = c(0,1,0,1,1,1,0,1,1,0)
)
support_A <- mean(services$A == 1)
support_B <- mean(services$B == 1)
support_AB <- mean(services$A == 1 & services$B == 1)
confidence_A_to_B <- support_AB / support_A
lift_A_to_B <- confidence_A_to_B / support_B
c(support_AB = support_AB,
  confidence = confidence_A_to_B,
  lift = lift_A_to_B)

support_AB confidence       lift 
 0.4000000  0.6666667  1.1111111 

Definitions:

• Support(A and B): proportion of all transactions containing A and B.
• Confidence(A → B): proportion of A transactions that also contain B.
• Lift(A → B): confidence divided by the overall prevalence of B. A value greater than one means B occurs more often after observing A than its marginal frequency, but confounding and selection can explain the association.

These formulas require denominators greater than zero. For sparse clinical events, support may be tiny and rules unstable.

2. Apriori and FP-growth

Apriori finds frequent item sets by extending frequent smaller sets while pruning candidates. FP-growth compresses transaction patterns into a tree representation and often avoids Apriori’s explicit candidate enumeration. Both require thresholds to manage the enormous number of possible rules.

Example: Apriori works in passes. First count single items and drop the ones below the minimum support; then build pairs only from the survivors. With the ten fictional visits above and a minimum support of 0.4 (4 visits), the passes look like this in base R.

colSums(services)                 # pass 1: visits containing each single item

A B C 
6 6 6 

crossprod(as.matrix(services))    # pass 2: visits containing each pair (diagonal = single items)

  A B C
A 6 4 3
B 4 6 3
C 3 3 6

Each single service appears in 6 visits, so all three pass 4. For pairs, A and B occur together in 4 visits, which meets the minimum, so the pair passes, while A with C and B with C occur together in only 3, so those two pairs are pruned and no set containing them is tried later. The saving is small here and large when there are thousands of items. FP-growth reaches the same frequent sets by a different route and gives the same rules.

In R, the arules package runs Apriori (and Eclat, a related method). The next five sections use it on the synthetic EHR data. An optional companion, arulesViz, draws rules as scatter plots and graphs.

3. Turn the data into transactions

Apriori needs yes/no items, not numbers. Numeric columns must be cut into labeled bins first, using cut(). Each bin then becomes one item, such as age=60+. Choose bin edges that a reader can understand; here are simple ones.

library(arules)
ehr <- readr::read_csv("data/synthetic_ehr.csv", show_col_types = FALSE)
bins <- data.frame(
  age = cut(ehr$age, c(0, 40, 60, Inf), labels = c("<40", "40-59", "60+"),
            right = FALSE),
  bmi = cut(ehr$bmi, c(0, 25, 30, Inf), labels = c("<25", "25-29", "30+"),
            right = FALSE),
  sbp = cut(ehr$sbp, c(0, 120, 140, Inf), labels = c("<120", "120-139", "140+"),
            right = FALSE),
  family_history = factor(ehr$family_history, levels = 0:1, labels = c("No", "Yes")),
  diabetes = factor(ehr$diabetes, levels = c("No", "Yes"))
)
tx <- as(bins, "transactions")
tx

transactions in sparse format with
 600 transactions (rows) and
 13 items (columns)

itemLabels(tx)

 [1] "age=<40"            "age=40-59"          "age=60+"           
 [4] "bmi=<25"            "bmi=25-29"          "bmi=30+"           
 [7] "sbp=<120"           "sbp=120-139"        "sbp=140+"          
[10] "family_history=No"  "family_history=Yes" "diabetes=No"       
[13] "diabetes=Yes"      

as(bins, "transactions") turns each factor level into an item named column=level, so 600 rows become 600 transactions over 13 items. A missing value (NA) just produces no item for that column in that row. Each row is one visit, which is fine for this exercise but see the pitfalls at the end of this chapter.

4. Mine rules with apriori() and read them

apriori() needs a minimum support (supp), a minimum confidence (conf) and, to keep the output small, a longest rule length (maxlen). Here a rule must cover at least 5% of the visits and be right at least half the time, with at most three items in total.

rules <- apriori(tx,
                 parameter = list(supp = 0.05, conf = 0.5, minlen = 2, maxlen = 3),
                 control = list(verbose = FALSE))
summary(rules)

set of 142 rules

rule length distribution (lhs + rhs):sizes
  2   3 
 34 108 

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.000   3.000   3.000   2.761   3.000   3.000 

summary of quality measures:
    support          confidence        coverage            lift       
 Min.   :0.05167   Min.   :0.5000   Min.   :0.05333   Min.   :0.6394  
 1st Qu.:0.07333   1st Qu.:0.6912   1st Qu.:0.09667   1st Qu.:1.0016  
 Median :0.12333   Median :0.7364   Median :0.15667   Median :1.0700  
 Mean   :0.15577   Mean   :0.7697   Mean   :0.20899   Mean   :1.2536  
 3rd Qu.:0.21417   3rd Qu.:0.9129   3rd Qu.:0.29167   3rd Qu.:1.4000  
 Max.   :0.64833   Max.   :1.0000   Max.   :0.91000   Max.   :2.6272  
     count       
 Min.   : 31.00  
 1st Qu.: 44.00  
 Median : 74.00  
 Mean   : 93.46  
 3rd Qu.:128.50  
 Max.   :389.00  

mining info:
 data ntransactions support confidence
   tx           600    0.05        0.5
                                                                                                                   call
 apriori(data = tx, parameter = list(supp = 0.05, conf = 0.5, minlen = 2, maxlen = 3), control = list(verbose = FALSE))

summary(rules) reports how many rules were found, how many items they have, and the range of support, confidence and lift. Sort by lift to see the strongest associations first:

inspect(head(sort(rules, by = "lift"), 5))

    lhs                              rhs        support    confidence
[1] {sbp=140+, family_history=No} => {age=60+}  0.07000000 0.9545455 
[2] {sbp=140+}                    => {age=60+}  0.09166667 0.9482759 
[3] {bmi=30+, sbp=140+}           => {age=60+}  0.06500000 0.9285714 
[4] {sbp=140+, diabetes=No}       => {age=60+}  0.05333333 0.9142857 
[5] {age=<40, bmi=25-29}          => {sbp=<120} 0.06833333 0.8913043 
    coverage   lift     count
[1] 0.07333333 2.627189 42   
[2] 0.09666667 2.609934 55   
[3] 0.07000000 2.555701 39   
[4] 0.05833333 2.516383 32   
[5] 0.07666667 2.209845 41   

Each row reads “if the left side is present, the right side tends to be present”. For the first rule, 7% of all visits contain all three items (support), and 95% of the visits with sbp=140+ and family_history=No are also age 60 or older (confidence). lift is about 2.6: age 60 or older is 2.6 times as common in these visits as in all visits. count is the number of rows behind the rule. The left side has only 44 rows (coverage times 600), which is small, so treat the exact numbers loosely.

5. Read support, confidence and lift together

No one of the three is enough. Sorting by confidence instead of lift gives a different list:

inspect(head(sort(rules, by = "confidence"), 5))

    lhs                       rhs           support    confidence coverage  
[1] {age=40-59, bmi=25-29} => {diabetes=No} 0.06166667 1          0.06166667
[2] {bmi=25-29, sbp=<120}  => {diabetes=No} 0.09833333 1          0.09833333
[3] {age=40-59, bmi=<25}   => {diabetes=No} 0.08333333 1          0.08333333
[4] {age=<40, bmi=<25}     => {diabetes=No} 0.08000000 1          0.08000000
[5] {bmi=<25, sbp=<120}    => {diabetes=No} 0.13500000 1          0.13500000
    lift     count
[1] 1.098901 37   
[2] 1.098901 59   
[3] 1.098901 50   
[4] 1.098901 48   
[5] 1.098901 81   

base_rate <- mean(bins$diabetes == "No")
round(base_rate, 3)

[1] 0.91

These rules have confidence 1 (the right side is always present when the left side is), yet the lift is only about 1.1. The reason is the base rate: 91% of all visits have diabetes=No, so a rule that predicts diabetes=No is almost always right without knowing anything. Lift divides confidence by that base rate, so it exposes this.

A lift near 1 (roughly 0.9 to 1.1) means no association, whatever the confidence says. A high lift with a small count can be chance, so read count too.

6. Rules for one outcome with appearance

Often one item is the question, for example diabetes=Yes. Use appearance to keep only rules with that item on the right and everything else on the left. Lower support is needed because only 9% of visits have the outcome.

rules_dm <- apriori(tx,
                    parameter = list(supp = 0.02, conf = 0.2, minlen = 2, maxlen = 3),
                    appearance = list(rhs = "diabetes=Yes", default = "lhs"),
                    control = list(verbose = FALSE))
length(rules_dm)

[1] 6

inspect(head(sort(rules_dm, by = "lift"), 5))

    lhs                              rhs            support    confidence
[1] {bmi=30+, sbp=140+}           => {diabetes=Yes} 0.03500000 0.5000000 
[2] {age=60+, sbp=140+}           => {diabetes=Yes} 0.03833333 0.4181818 
[3] {sbp=140+}                    => {diabetes=Yes} 0.03833333 0.3965517 
[4] {sbp=140+, family_history=No} => {diabetes=Yes} 0.02666667 0.3636364 
[5] {age=60+, bmi=30+}            => {diabetes=Yes} 0.05166667 0.2897196 
    coverage   lift     count
[1] 0.07000000 5.555556 21   
[2] 0.09166667 4.646465 23   
[3] 0.09666667 4.406130 23   
[4] 0.07333333 4.040404 16   
[5] 0.17833333 3.219107 31   

Here the top rule has a lift above 5: diabetes=Yes is more than 5 times as common among visits with bmi=30+ and sbp=140+ as among all visits. It still describes only 21 rows, and it is a pattern in fake data, not a risk score.

7. When apriori() returns no rules

An empty result usually means the thresholds are too strict for the data. Compare the same question at a strict and a relaxed setting.

strict <- apriori(tx,
                  parameter = list(supp = 0.3, conf = 0.95, minlen = 2, maxlen = 3),
                  appearance = list(rhs = "diabetes=Yes", default = "lhs"),
                  control = list(verbose = FALSE))
length(strict)

[1] 0

length(rules_dm)

[1] 6

sum(ehr$diabetes == "Yes")

[1] 54

The strict run finds nothing. Only 54 of the 600 visits (9%) have diabetes=Yes, and every rule here has it on the right, so no rule can have support above 0.09 and none can meet supp = 0.3. When this happens, lower supp first (it sets how rare an item set may be), then conf. Also check that the bins are not too fine, which spreads rows across many tiny items, and that the right-hand item exists in itemLabels(tx). Do not lower the thresholds until rules appear and then report them as findings: very low support gives many rules based on a handful of rows.

8. Anomaly and novelty detection

Anomaly detection asks: Which points are unusually different from a defined reference?

It is not a binary classifier unless verified anomaly labels exist. If a machine flags a rare systolic BP value, that record may represent a true unusual patient, a measurement problem, or a unit error. Do not automatically discard it.

Example: look at the five highest systolic pressures in the synthetic data and decide what a flag would mean.

sort(ehr$sbp, decreasing = TRUE)[1:5]

[1] 154 152 151 151 150

The highest values are 154, 152, 151, 151 and 150 mmHg. A detector would rank these first, but all are plausible readings, so the sensible action is to leave them in. A slip such as 15.4 or 1540 would look very different, and even then the score cannot tell you whether it is a typing error or a real patient. Only a look at the record, or the person who entered it, can decide.

9. Standard scores and robust scores

A simple z-score measures how far a value is from the average in standard-deviation units. A robust approach uses medians and median absolute deviation (MAD), often less sensitive to extreme points.

ehr <- readr::read_csv("data/synthetic_ehr.csv", show_col_types = FALSE)
values <- ehr$sbp
z <- as.numeric(scale(values))
med <- median(values, na.rm = TRUE)
madv <- mad(values, constant = 1, na.rm = TRUE)
robust_z <- (values - med) / madv
summary(z)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max.      NAs 
-2.71227 -0.70967 -0.01647  0.00000  0.67674  2.44827       20 

summary(robust_z)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max.      NAs 
-3.88889 -1.00000  0.00000  0.02375  1.00000  3.55556       20 

sum(abs(z) > 3, na.rm = TRUE)

[1] 0

Because constant = 1 leaves the MAD unscaled, robust_z is about 1.48 times larger than a MAD score scaled to match the standard deviation, so compare the two scores by rank rather than raw size. A rule like |z| > 3 is a statistical screening heuristic, not an automatic clinical cutoff. A high score simply means distant from the center (mean for z, median for robust_z) of this dataset.

10. Local and model-based anomaly methods

Optional R packages: dbscan (LOF), isotree (isolation forest), kernlab (one-class SVM), robustbase (robust covariance). The added unsupervised script demonstrates safe, limited anomaly scores without needing specialized packages.

Example: base R has mahalanobis(), which (with the ordinary cov() used here, the classical version) scores a whole row at once and accounts for correlation between the columns. It returns squared distances, and a common cut-off is a high quantile of the chi-square distribution with one degree of freedom per column.

three <- na.omit(ehr[, c("age", "bmi", "sbp")])
d2 <- mahalanobis(three, colMeans(three), cov(three))
round(sort(d2, decreasing = TRUE)[1:3], 1)

[1] 8.0 7.8 7.3

three[order(-d2)[1:3], ]

# A tibble: 3 x 3
    age   bmi   sbp
  <dbl> <dbl> <dbl>
1    18  18.8    87
2    81  43     137
3    22  42.8   106

sum(d2 > qchisq(0.999, df = 3))

[1] 0

The three largest squared distances are 8.0, 7.8 and 7.3, and the rows behind them are unusual combinations (for example age 18 with a BMI of 18.8 and an SBP of 87). None passes the cut-off (qchisq(0.999, 3) is about 16.3, and the last line counts 0 rows above it), so this method flags nothing in these fake data. That is a reasonable outcome: a method that always finds anomalies is not detecting anything. The robust version swaps cov() for a robust covariance estimate (the robustbase package, listed above, offers one); it is an option, not used here.

11. Evaluation when anomalies are rare or unlabeled

If there are verified anomaly labels, evaluate precision/recall at inspection capacity and assess false alarm workload. Without ground truth, use simulated corruption, domain expert audits, stability studies, and prospective review rather than claiming high “sensitivity.” A feature indicating the clinician already corrected a record would leak the answer.

Example: simulated corruption gives you labels you control. Take the recorded sbp values, divide 10 randomly chosen ones by 10 (a decimal slip), and check how many a robust score catches.

set.seed(27)
sbp_all <- na.omit(ehr$sbp)
slipped <- sample(length(sbp_all), 10)
sbp_bad <- sbp_all
sbp_bad[slipped] <- sbp_bad[slipped] / 10
score <- abs(sbp_bad - median(sbp_bad)) / mad(sbp_bad)     # MAD scaled to match an SD
flag <- score > 3.5
c(caught = sum(flag[slipped]), false_alarms = sum(flag[-slipped]), total_flagged = sum(flag))

       caught  false_alarms total_flagged 
           10             0            10 

All 10 slipped values are flagged and none of the 570 untouched values is, so on this test the score has perfect recall and no false alarms. That is an easy test: a value ten times too small is far from the rest. Real errors are subtler, so repeat the exercise with smaller changes before trusting a threshold, and treat the numbers as a check on the method, not as sensitivity in real data.

12. Association mining and EHR pitfalls

• One patient with ten visits could generate ten correlated transactions.
• Frequent services may reflect workflows rather than disease mechanisms.
• Small groups can create misleadingly high lift.
• Searching thousands of rules creates multiple-testing and adaptive selection issues.
• Rare combinations may re-identify individuals even if names were removed.
• Selecting a model or threshold based on test outcomes destroys independence.

Example: the first pitfall in numbers. Here are five fictional visits from three patients. Patient p1 comes three times and always gets service A.

visits <- data.frame(patient = c("p1", "p1", "p1", "p2", "p3"),
                     A = c(1, 1, 1, 0, 0))
mean(visits$A == 1)                                  # support counted over visits

[1] 0.6

mean(tapply(visits$A, visits$patient, max) == 1)     # support counted over patients

[1] 0.3333333

length(unique(ehr$patient_id)); nrow(ehr)               # the course data

[1] 600

[1] 600

Counted over visits, A has support 0.6; counted over patients it has 0.33, because one patient supplies three of the three A visits. In the synthetic data every patient_id appears once (600 patients, 600 rows), so this problem does not arise here, but it would in a real visit table. Decide whether a transaction is a visit or a patient before mining.

c(support_A = support_A, support_B = support_B, support_AB = support_AB)

 support_A  support_B support_AB 
       0.6        0.6        0.4 

support_AB / support_A

[1] 0.6666667

sum(services$A == 1 & services$B == 1) / sum(services$A == 1)

[1] 0.6666667

rule_stats <- function(d) {
  s_A <- mean(d$A == 1); s_B <- mean(d$B == 1); s_AB <- mean(d$A == 1 & d$B == 1)
  round(c(support_B = s_B, confidence = s_AB / s_A, lift = s_AB / s_A / s_B), 3)
}
changed <- services
changed$B[5] <- 0
rbind(original = rule_stats(services), changed = rule_stats(changed))

         support_B confidence  lift
original       0.6      0.667 1.111
changed        0.5      0.667 1.333

cor(z, robust_z, method = "spearman", use = "complete.obs")

[1] 1

c(z_over_3 = sum(abs(z) > 3, na.rm = TRUE),
  robust_over_3 = sum(abs(robust_z) > 3, na.rm = TRUE),
  robust_scaled_over_3 = sum(abs(robust_z / 1.4826) > 3, na.rm = TRUE))

            z_over_3        robust_over_3 robust_scaled_over_3 
                   0                   13                    0 

bins2 <- bins
bins2$bmi <- cut(ehr$bmi, c(0, 30, Inf), labels = c("<30", "30+"), right = FALSE)
tx2 <- as(bins2, "transactions")
par_a <- list(supp = 0.05, conf = 0.5, minlen = 2, maxlen = 3)
c(three_bins = length(rules),
  two_bins = length(apriori(tx2, parameter = par_a, control = list(verbose = FALSE))))

three_bins   two_bins 
       142        135 

rules_age <- apriori(tx,
                     parameter = list(supp = 0.05, conf = 0.3, minlen = 2, maxlen = 3),
                     appearance = list(rhs = "age=60+", default = "lhs"),
                     control = list(verbose = FALSE))
inspect(head(sort(rules_age, by = "lift"), 3))

    lhs                              rhs       support    confidence coverage  
[1] {sbp=140+, family_history=No} => {age=60+} 0.07000000 0.9545455  0.07333333
[2] {sbp=140+}                    => {age=60+} 0.09166667 0.9482759  0.09666667
[3] {bmi=30+, sbp=140+}           => {age=60+} 0.06500000 0.9285714  0.07000000
    lift     count
[1] 2.627189 42   
[2] 2.609934 55   
[3] 2.555701 39   

mean(bins$age == "60+")

[1] 0.3633333

inspect(head(sort(rules, by = "confidence"), 3))

    lhs                       rhs           support    confidence coverage  
[1] {age=40-59, bmi=25-29} => {diabetes=No} 0.06166667 1          0.06166667
[2] {bmi=25-29, sbp=<120}  => {diabetes=No} 0.09833333 1          0.09833333
[3] {age=40-59, bmi=<25}   => {diabetes=No} 0.08333333 1          0.08333333
    lift     count
[1] 1.098901 37   
[2] 1.098901 59   
[3] 1.098901 50   

round(mean(bins$diabetes == "No"), 3)

[1] 0.91

strict_rules <- apriori(tx, parameter = list(supp = 0.4, conf = 0.9, minlen = 2, maxlen = 3),
                        control = list(verbose = FALSE))
length(strict_rules)

[1] 2

inspect(strict_rules)

    lhs                    rhs           support   confidence coverage 
[1] {sbp=120-139}       => {diabetes=No} 0.4250000 0.9107143  0.4666667
[2] {family_history=No} => {diabetes=No} 0.6483333 0.9067599  0.7150000
    lift      count
[1] 1.0007849 255  
[2] 0.9964395 389