161.777 Practical Data Mining
Project 2025
Esha Kiran Kumar Kalyanam
# Add any other packages you need to load here.
library(tidyverse)
library(arules)
# Read in the data
ice.train <- read_csv("https://www.massey.ac.nz/~jcmarsha/161777/assessment/data/ice-train.csv")
ice.test <- read_csv("https://www.massey.ac.nz/~jcmarsha/161777/assessment/data/ice-test.csv")
campy.train <- read_csv("https://www.massey.ac.nz/~jcmarsha/161777/assessment/data/campy-train.csv") |>
mutate(Source = factor(Source))
campy.test <- read_csv("https://www.massey.ac.nz/~jcmarsha/161777/assessment/data/campy-test.csv")
wv_survey <- read_csv("https://www.massey.ac.nz/~jcmarsha/161777/assessment/data/wv_survey2014.csv")
flags <- read.transactions("https://www.massey.ac.nz/~jcmarsha/161777/assessment/data/flags.txt", cols=1)Exercise 1: Predicting ice thickness
# Putting all the libraries that i can use for this exercise
library(tidymodels)## ── Attaching packages ────────────────────────────────────── tidymodels 1.3.0 ──## ✔ broom 1.0.8 ✔ rsample 1.3.0
## ✔ dials 1.4.0 ✔ tune 1.3.0
## ✔ infer 1.0.8 ✔ workflows 1.2.0
## ✔ modeldata 1.4.0 ✔ workflowsets 1.1.1
## ✔ parsnip 1.3.2 ✔ yardstick 1.3.2
## ✔ recipes 1.3.1## ── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
## ✖ scales::discard() masks purrr::discard()
## ✖ recipes::discretize() masks arules::discretize()
## ✖ Matrix::expand() masks tidyr::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ recipes::fixed() masks stringr::fixed()
## ✖ dplyr::lag() masks stats::lag()
## ✖ Matrix::pack() masks tidyr::pack()
## ✖ arules::recode() masks dplyr::recode()
## ✖ yardstick::spec() masks readr::spec()
## ✖ recipes::step() masks stats::step()
## ✖ Matrix::unpack() masks tidyr::unpack()
## ✖ recipes::update() masks Matrix::update(), stats::update()#library(janitor)
#library(GGally)
#library(skimr)Histogram of the thickness
# Distribution of target
ggplot(ice.train, aes(x = Thickness)) + geom_histogram(bins = 30)
This graph shows us that the target variable is skewed (the distribution
is basically not normal) which mean that linear regression models might
not work that well because it usually expects a normally distributed
residuals, and works best when the target variable is roughly
symmetric.
So we should work with a model like Random Forest which doesnt assume normality.
- Thickness Trend over the Years
ggplot(ice.train, aes(Year, Thickness)) +
geom_smooth(method = "loess", color = "red") +
labs(title = "Thickness Trend Over Years")## `geom_smooth()` using formula = 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 2007## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 1.01## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 9.3734e-28## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 1## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 2007## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 1.01## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 9.3734e-28## Warning in predLoess(object$y, object$x, newx = if (is.null(newdata)) object$x
## else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 1
A curved trend in the red line (LOESS fit), suggesting non-linear
relationships over time. The gray shading shows uncertainty increasing,
possibly due to fewer data points or higher variance in those years.
Again, suggests that a flexible, non-linear model like Random Forest
will perform better.
- Max Latitude vs Thickness (hexbin plot)
ggplot(ice.train, aes(MaxLat, Thickness)) +
geom_hex(bins = 30) +
scale_fill_viridis_c() +
labs(title = "Thickness by Latitude")
A light center and dark edges — that is, a concentrated density of
moderate Thickness around a specific latitude, and more spread/extremes
elsewhere.
Non-linear spatial pattern: Linear regression might fail to capture this structure, whereas Random Forest is good at capturing these kinds of conditional interactions and localized behavior.
sapply(ice.train, function(x) sum(is.na(x)))## Year MinDay MaxDay MinLat MaxLat MinLon MaxLon Length
## 0 0 0 0 0 0 0 0
## Nsamps Thickness
## 0 0No missingness shows us that we dont need any imputation for this exercise
1.1 Ice Thickness Methodology
WRITE UP YOUR METHODOLOGY HERE
# Recipe: Normalize numeric predictors
ice_recipe <- recipe(Thickness ~ ., data = ice.train) %>%
step_normalize(all_predictors())
# Model spec: Random Forest
rf_spec <- rand_forest(mtry = tune(), min_n = tune(), trees = 500) %>%
set_engine("ranger") %>%
set_mode("regression")
# Workflow
rf_workflow <- workflow() %>%
add_recipe(ice_recipe) %>%
add_model(rf_spec)
# Cross-validation setup
cv_folds <- vfold_cv(ice.train, v = 5)
# Tuning grid
rf_grid <- grid_regular(
mtry(range = c(2, 8)),
min_n(range = c(2, 10)),
levels = 4
)
# Tune the model
rf_tune <- tune_grid(
rf_workflow,
resamples = cv_folds,
grid = rf_grid,
metrics = metric_set(rmse)
)
# Select best parameters
best_rf <- rf_tune %>%
select_best(metric = "rmse")
# Finalize workflow with best params
final_rf_workflow <- finalize_workflow(rf_workflow, best_rf)
# Fit on full training data
final_rf_fit <- fit(final_rf_workflow, data = ice.train)
# Predict on test data
ice.test.with.predictions <- ice.test %>%
bind_cols(predict(final_rf_fit, new_data = ice.test))
# Write the result to a CSV (change studentid to your actual ID)
write_csv(ice.test.with.predictions, "ex1_24020251.csv")# Recipe: Normalize numeric predictors
ice_recipe <- recipe(Thickness ~ ., data = ice.train) %>%
step_normalize(all_predictors())
# Random Forest model spec with fixed parameters (no tuning)
rf_spec_simple <- rand_forest(mtry = 5, min_n = 5, trees = 500) %>%
set_engine("ranger") %>%
set_mode("regression")
# Workflow
rf_workflow_simple <- workflow() %>%
add_recipe(ice_recipe) %>%
add_model(rf_spec_simple)
# Fit on full training data
final_rf_fit_simple <- fit(rf_workflow_simple, data = ice.train)
# Predict on test data
ice.test.with.predictions <- ice.test %>%
bind_cols(predict(final_rf_fit_simple, new_data = ice.test))
# Write result to CSV (change studentid to your actual ID)
write_csv(ice.test.with.predictions, "ex1_24020251_simple.csv")We have coded using random forest both with tunning and a simple model to compare, the tunning definitely took more time (~approx 20 minutes) Now we will compare the 2 models and see which one performs better
set.seed(123)
data_split <- initial_split(ice.train, prop = 0.8)
train_data <- training(data_split)
valid_data <- testing(data_split)Testing for Untuned model
# testing on untuned model
# Refit simple model on train_data
final_rf_fit_simple <- fit(rf_workflow_simple, data = train_data)
# Predict and evaluate on validation data
simple_preds <- predict(final_rf_fit_simple, new_data = valid_data) %>%
bind_cols(valid_data)
metrics(simple_preds, truth = Thickness, estimate = .pred)## # A tibble: 3 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 rmse standard 0.349
## 2 rsq standard 0.873
## 3 mae standard 0.232# testing on tuned model
# Use same folds for tuning but only from train_data
cv_folds <- vfold_cv(train_data, v = 5)
# (tune_rf... same as before)
# After tuning:
final_rf_fit <- fit(final_rf_workflow, data = train_data)
# Predict and evaluate on validation data
tuned_preds <- predict(final_rf_fit, new_data = valid_data) %>%
bind_cols(valid_data)
metrics(tuned_preds, truth = Thickness, estimate = .pred)## # A tibble: 3 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 rmse standard 0.347
## 2 rsq standard 0.873
## 3 mae standard 0.230Here we need Mae to reach 0, and the smaller value of rmse performs better, and larger value of rsq performs better, this case tunned model performs slightly better, To note this isnt the test where we have tested it on tho.
# Get metrics from both
untuned_metrics <- metrics(simple_preds, truth = Thickness, estimate = .pred) %>%
mutate(Model = "Untuned RF")
tuned_metrics <- metrics(tuned_preds, truth = Thickness, estimate = .pred) %>%
mutate(Model = "Tuned RF")
# Combine both
all_metrics <- bind_rows(tuned_metrics, untuned_metrics)
# Plot
ggplot(all_metrics, aes(x = .metric, y = .estimate, fill = Model)) +
geom_bar(stat = "identity", position = position_dodge(), width = 0.6) +
scale_fill_manual(values = c("Tuned RF" = "#4CAF50", "Untuned RF" = "#FFC107")) +
labs(title = "Tuned vs Untuned Random Forest Model Performance",
x = "Metric",
y = "Score",
fill = "Model") +
theme_minimal() +
theme(text = element_text(size = 12))
## Exercise 2: Predicting the source of campylobacter
# Check structure
glimpse(campy.train)## Rows: 400
## Columns: 51
## $ Source <fct> Sheep, Beef, Beef, Poultry, Beef, Poultry, Poultry, Poultr…
## $ CAMP0003.V1 <dbl> -185.92723, 77.99245, 82.51236, 81.49695, 81.79943, 81.496…
## $ CAMP0003.V2 <dbl> 2.3053030, -5.2710045, 6.4827317, 5.1527925, -11.8247562, …
## $ CAMP0012.V1 <dbl> 36.71589, -15.06629, -13.12199, -15.06629, -15.06629, -15.…
## $ CAMP0012.V2 <dbl> -4.406882, -4.317577, 9.310441, -4.317577, -4.317577, -4.3…
## $ CAMP0074.V1 <dbl> 127.28044, -27.22901, -29.92261, -36.57515, -28.32487, -28…
## $ CAMP0074.V2 <dbl> 3.484098, -8.794580, -9.946344, -3.173427, -8.516313, -8.5…
## $ CAMP0133.V1 <dbl> 192.70503, -49.42997, -59.27090, -53.74389, -56.27665, -56…
## $ CAMP0133.V2 <dbl> 7.7104440, 7.8912927, 4.4158963, 17.0999832, 19.2778281, 7…
## $ CAMP0201.V1 <dbl> -265.6023, 161.5823, 161.3653, 164.5170, 164.5170, 163.520…
## $ CAMP0201.V2 <dbl> 19.0953706, -7.1245339, 4.9481621, -11.3762227, -11.376222…
## $ CAMP0226.V1 <dbl> -218.7569, 110.4446, 106.0126, 103.1607, 105.4473, 104.045…
## $ CAMP0226.V2 <dbl> -11.560149, -10.886575, -5.472559, -3.802787, -2.101317, -…
## $ CAMP0236.V1 <dbl> -53.23589, 20.28444, 21.42368, 21.22623, 17.68375, 21.2262…
## $ CAMP0236.V2 <dbl> 10.98413393, 4.96094073, -8.60420238, 4.10837595, -5.20575…
## $ CAMP0289.V1 <dbl> -40.936758, 9.498759, 8.538700, 8.907599, 289.749193, 9.19…
## $ CAMP0289.V2 <dbl> -383.75218, -260.29853, -261.19871, -261.00437, 579.24317,…
## $ CAMP0609.V1 <dbl> 266.49397, -97.36202, -101.19318, -98.29937, -99.36610, -9…
## $ CAMP0609.V2 <dbl> -5.2663075, 2.1039444, 6.7873064, 5.1851643, 4.7805647, 5.…
## $ CAMP0879.V1 <dbl> -110.76238, 49.89587, 52.54514, 54.70570, 55.53442, 55.535…
## $ CAMP0879.V2 <dbl> -0.6630628, 9.4071670, 22.3415935, 22.0737632, 23.1426925,…
## $ CAMP0976.V1 <dbl> 15.540253, -4.270418, -5.510377, -5.673055, -4.843925, -5.…
## $ CAMP0976.V2 <dbl> -0.2775330, -3.6683760, -3.9919041, -4.4073620, -3.1113960…
## $ CAMP1015.V1 <dbl> -4.847205, 1.750955, 2.437035, 1.750955, 1.750955, 1.75095…
## $ CAMP1015.V2 <dbl> -0.1836036, -2.1642475, -2.9888490, -2.1642475, -2.1642475…
## $ CAMP1076.V1 <dbl> 242.50513, -82.19588, -97.21802, -93.17179, -93.17179, -96…
## $ CAMP1076.V2 <dbl> 3.8424376, -7.5214780, -25.3930349, 25.8398837, 25.8398837…
## $ CAMP1080.V1 <dbl> -445.4994, 182.8193, 192.6265, 192.6265, 193.6830, 193.819…
## $ CAMP1080.V2 <dbl> 13.9994478, -13.6545235, -57.8575436, -57.8575436, -56.571…
## $ CAMP1082.V1 <dbl> 140.06821, -37.53879, -36.27916, -34.29013, -35.34534, -36…
## $ CAMP1082.V2 <dbl> 1.801343, -7.144345, 3.154446, -10.798173, 4.048094, 2.467…
## $ CAMP1162.V1 <dbl> 114.56253, -28.19458, -31.90352, -30.18738, -29.43875, -24…
## $ CAMP1162.V2 <dbl> 2.1005307, -5.4344872, -0.5582965, 0.9994499, 4.3748606, -…
## $ CAMP1176.V1 <dbl> 153.76031, -57.79709, -56.87719, -58.44376, -55.03340, -58…
## $ CAMP1176.V2 <dbl> -2.3040141, 0.3310222, -12.0987121, 13.4442918, -10.280233…
## $ CAMP1178.V1 <dbl> -213.86438, 88.46626, 92.26185, 90.00291, 90.80221, 80.342…
## $ CAMP1178.V2 <dbl> 24.198082, 35.019066, 9.164906, 8.643049, 8.179239, 12.356…
## $ CAMP1179.V1 <dbl> 133.49628, -75.64334, -73.87632, -73.87632, -73.87632, -64…
## $ CAMP1179.V2 <dbl> 0.8648811, 0.1344619, -0.1911970, -0.1911970, -0.1911970, …
## $ CAMP1213.V1 <dbl> 15.58629, 18.32932, 20.55584, 24.20177, 23.24912, 22.02397…
## $ CAMP1213.V2 <dbl> 13.20111, 20.47641, 13.26037, -78.64188, -78.36459, -78.37…
## $ CAMP1251.V1 <dbl> -151.5062075, -150.6651212, 9.4286645, 14.7790781, 7.96852…
## $ CAMP1251.V2 <dbl> -257.57188, -256.76081, -154.36687, -132.52259, -131.79333…
## $ CAMP1300.V1 <dbl> 28.20891, -11.80734, -13.71199, -12.76947, -12.76947, -12.…
## $ CAMP1300.V2 <dbl> 0.03383027, 0.45099922, -2.65296220, -0.24293368, -0.24293…
## $ CAMP1473.V1 <dbl> -92.81715, 17.75843, 16.67712, 19.52353, 18.48182, 19.1287…
## $ CAMP1473.V2 <dbl> 1.2781163, 2.3199203, 1.9282542, -0.8038236, -0.7653928, -…
## $ CAMP1538.V1 <dbl> 24.38298, -12.71836, -15.11377, -12.38350, -12.38350, -12.…
## $ CAMP1538.V2 <dbl> -0.04545584, 2.23272134, 3.35528459, 3.00514028, 3.0051402…
## $ CAMP1613.V1 <dbl> -46.93894, 19.28144, 16.24487, 19.28139, 19.28139, 19.2813…
## $ CAMP1613.V2 <dbl> 8.2435334, 0.8618736, 0.6620816, 1.5784402, 1.5784402, 1.5…# Count of each class
campy.train %>% count(Source)## # A tibble: 3 × 2
## Source n
## <fct> <int>
## 1 Beef 128
## 2 Poultry 139
## 3 Sheep 133# Check for missing values
sum(is.na(campy.train))## [1] 0# Summary stats
summary(campy.train)## Source CAMP0003.V1 CAMP0003.V2 CAMP0012.V1
## Beef :128 Min. :-204.46 Min. :-16.0597 Min. :-17.4899
## Poultry:139 1st Qu.: 74.46 1st Qu.:-11.0379 1st Qu.:-15.0961
## Sheep :133 Median : 77.99 Median : 0.9881 Median :-15.0663
## Mean : 14.13 Mean : -1.7141 Mean : -0.8513
## 3rd Qu.: 81.50 3rd Qu.: 5.0889 3rd Qu.:-12.3948
## Max. : 85.33 Max. : 28.2928 Max. : 51.7786
## CAMP0012.V2 CAMP0074.V1 CAMP0074.V2 CAMP0133.V1
## Min. :-10.377 Min. :-88.348 Min. :-18.264 Min. :-59.271
## 1st Qu.: -4.318 1st Qu.:-29.106 1st Qu.: -8.795 1st Qu.:-53.930
## Median : -3.795 Median :-28.325 Median : -8.516 Median :-52.647
## Mean : -0.953 Mean : 9.105 Mean : -4.691 Mean : 7.547
## 3rd Qu.: 1.733 3rd Qu.:-25.595 3rd Qu.: 3.450 3rd Qu.:-47.989
## Max. : 9.469 Max. :150.906 Max. : 47.311 Max. :195.332
## CAMP0133.V2 CAMP0201.V1 CAMP0201.V2 CAMP0226.V1
## Min. :-25.9698 Min. :-268.41 Min. :-12.5693 Min. :-267.72
## 1st Qu.: -0.6180 1st Qu.: 155.68 1st Qu.:-11.3762 1st Qu.: 99.79
## Median : -0.3057 Median : 163.59 Median : -4.9698 Median : 104.05
## Mean : 3.4674 Mean : 57.57 Mean : -5.0215 Mean : 26.35
## 3rd Qu.: 17.1000 3rd Qu.: 164.52 3rd Qu.: -0.1926 3rd Qu.: 106.01
## Max. : 20.8938 Max. : 168.43 Max. : 19.0954 Max. : 162.03
## CAMP0226.V2 CAMP0236.V1 CAMP0236.V2 CAMP0289.V1
## Min. :-23.572 Min. :-77.2692 Min. :-17.678 Min. :-561.117
## 1st Qu.: -9.077 1st Qu.: 17.6838 1st Qu.:-12.091 1st Qu.: -44.026
## Median : -4.369 Median : 19.2390 Median : -0.129 Median : 8.539
## Mean : -5.133 Mean : 0.3673 Mean : -1.505 Mean : -4.590
## 3rd Qu.: -3.264 3rd Qu.: 21.1276 3rd Qu.: 4.108 3rd Qu.: 9.499
## Max. :168.274 Max. : 37.0678 Max. : 18.413 Max. : 592.067
## CAMP0289.V2 CAMP0609.V1 CAMP0609.V2 CAMP0879.V1
## Min. :-384.42 Min. :-108.248 Min. :-23.013 Min. :-115.31
## 1st Qu.:-261.77 1st Qu.:-100.266 1st Qu.: -5.118 1st Qu.: 48.95
## Median :-261.00 Median : -98.299 Median : 3.444 Median : 54.71
## Mean :-172.16 Mean : -8.253 Mean : 1.353 Mean : 13.62
## 3rd Qu.: -82.38 3rd Qu.: -90.562 3rd Qu.: 5.539 3rd Qu.: 55.08
## Max. : 779.03 Max. : 269.499 Max. : 15.856 Max. : 63.91
## CAMP0879.V2 CAMP0976.V1 CAMP0976.V2 CAMP1015.V1
## Min. :-31.729 Min. :-9.1648 Min. :-5.1697 Min. :-5.9622
## 1st Qu.: -1.424 1st Qu.:-5.9109 1st Qu.:-4.4074 1st Qu.: 1.0101
## Median : 22.074 Median :-4.8439 Median :-3.1114 Median : 1.7510
## Mean : 9.022 Mean :-0.3095 Mean :-2.0925 Mean : 0.7256
## 3rd Qu.: 22.342 3rd Qu.: 2.0623 3rd Qu.: 0.5721 3rd Qu.: 1.7510
## Max. : 24.595 Max. :15.5887 Max. : 8.1437 Max. : 2.4370
## CAMP1015.V2 CAMP1076.V1 CAMP1076.V2 CAMP1080.V1
## Min. :-2.9888 Min. :-97.218 Min. :-28.842 Min. :-453.22
## 1st Qu.:-2.1642 1st Qu.:-93.369 1st Qu.: -5.382 1st Qu.: 169.96
## Median :-2.1642 Median :-84.672 Median : 6.326 Median : 179.29
## Mean :-1.3428 Mean : -7.145 Mean : 6.665 Mean : 28.43
## 3rd Qu.:-0.1808 3rd Qu.:-69.149 3rd Qu.: 23.282 3rd Qu.: 192.63
## Max. : 2.8211 Max. :245.375 Max. : 25.840 Max. : 194.58
## CAMP1080.V2 CAMP1082.V1 CAMP1082.V2 CAMP1162.V1
## Min. :-59.08 Min. :-40.671 Min. :-13.420 Min. :-54.667
## 1st Qu.:-57.86 1st Qu.:-36.279 1st Qu.: -8.009 1st Qu.:-29.537
## Median :-12.23 Median :-35.130 Median : 1.438 Median :-28.466
## Mean :-17.35 Mean : 7.849 Mean : -1.431 Mean : 5.457
## 3rd Qu.: 6.60 3rd Qu.:-28.946 3rd Qu.: 3.154 3rd Qu.:-23.170
## Max. : 38.44 Max. :141.063 Max. : 12.675 Max. :116.532
## CAMP1162.V2 CAMP1176.V1 CAMP1176.V2 CAMP1178.V1
## Min. :-29.0788 Min. :-68.919 Min. :-24.8112 Min. :-266.269
## 1st Qu.: -0.2941 1st Qu.:-58.444 1st Qu.:-11.4941 1st Qu.:-126.235
## Median : 2.1026 Median :-55.971 Median : -2.9233 Median : 86.134
## Mean : 0.3327 Mean : -5.393 Mean : -0.3658 Mean : -8.906
## 3rd Qu.: 4.3242 3rd Qu.:-49.526 3rd Qu.: 13.4443 3rd Qu.: 90.378
## Max. : 22.6198 Max. :154.828 Max. : 18.8680 Max. : 109.235
## CAMP1178.V2 CAMP1179.V1 CAMP1179.V2 CAMP1213.V1
## Min. :-256.045 Min. :-86.57 Min. :-51.6123 Min. :-200.73
## 1st Qu.: 8.179 1st Qu.:-71.82 1st Qu.:-11.6794 1st Qu.: 19.31
## Median : 12.977 Median :-64.05 Median : -0.1912 Median : 22.02
## Mean : 12.426 Mean :-17.92 Mean : -1.0622 Mean : 14.75
## 3rd Qu.: 41.763 3rd Qu.:-59.31 3rd Qu.: 3.6769 3rd Qu.: 23.25
## Max. : 91.072 Max. :143.60 Max. : 34.2917 Max. : 54.62
## CAMP1213.V2 CAMP1251.V1 CAMP1251.V2 CAMP1300.V1
## Min. :-125.38 Min. :-1155.36 Min. :-257.57 Min. :-13.738
## 1st Qu.: -78.36 1st Qu.: -151.51 1st Qu.:-136.41 1st Qu.:-12.769
## Median : 13.26 Median : 13.96 Median :-130.51 Median :-11.807
## Mean : -17.90 Mean : -54.57 Mean : -76.09 Mean : -2.347
## 3rd Qu.: 21.35 3rd Qu.: 19.65 3rd Qu.: -26.32 3rd Qu.:-10.781
## Max. : 77.69 Max. : 435.94 Max. : 977.97 Max. : 29.204
## CAMP1300.V2 CAMP1473.V1 CAMP1473.V2 CAMP1538.V1
## Min. :-2.65296 Min. :-93.907 Min. :-47.7272 Min. :-17.999
## 1st Qu.:-0.24293 1st Qu.: 13.903 1st Qu.: -0.8415 1st Qu.:-12.718
## Median :-0.24293 Median : 18.063 Median : -0.7654 Median :-12.383
## Mean :-0.36992 Mean : -7.906 Mean : 0.2930 Mean : -4.263
## 3rd Qu.: 0.03383 3rd Qu.: 19.524 3rd Qu.: 1.5623 3rd Qu.:-11.597
## Max. : 2.09736 Max. : 24.204 Max. : 48.6252 Max. : 25.319
## CAMP1538.V2 CAMP1613.V1 CAMP1613.V2
## Min. :-8.8594 Min. :-66.6940 Min. :-36.7314
## 1st Qu.:-5.4720 1st Qu.: 16.2449 1st Qu.: 0.6712
## Median : 2.2327 Median : 19.2814 Median : 1.5784
## Mean :-0.6817 Mean : -0.2004 Mean : -0.5082
## 3rd Qu.: 3.0051 3rd Qu.: 20.0791 3rd Qu.: 1.5784
## Max. : 5.3694 Max. : 22.6006 Max. : 8.3241library(GGally)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2ggpairs(campy.train, columns = 1:5, aes(color = Source))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
2.1 Campylobacter Methodology
Data Splitting for cross-validation
set.seed(123)
campy_split <- initial_split(campy.train, prop = 0.8, strata = Source)
campy_train_split <- training(campy_split)
campy_valid_split <- testing(campy_split)Preprocessing Recipe
campy_recipe <- recipe(Source ~ ., data = campy_train_split) %>%
step_normalize(all_numeric_predictors())Trying different Models 1. Random Forest
rf_spec <- rand_forest(mtry = tune(), min_n = tune(), trees = 500) %>%
set_mode("classification") %>%
set_engine("ranger")
rf_workflow <- workflow() %>%
add_recipe(campy_recipe) %>%
add_model(rf_spec)
rf_res <- tune_grid(
rf_workflow,
resamples = vfold_cv(campy_train_split, v = 5, strata = Source),
grid = 10,
metrics = metric_set(accuracy)
)## i Creating pre-processing data to finalize unknown parameter: mtrybest_rf <- rf_res %>%
select_best(metric = "accuracy")
final_rf <- finalize_workflow(rf_workflow, best_rf)print("Best Random Forest configuration:")## [1] "Best Random Forest configuration:"print(best_rf)## # A tibble: 1 × 3
## mtry min_n .config
## <int> <int> <chr>
## 1 1 14 Preprocessor1_Model01# Showing all performances metrics from the tuning
rf_metrics <- collect_metrics(rf_res)
print("Random Forest Cross-Validation Results:")## [1] "Random Forest Cross-Validation Results:"print(rf_metrics)## # A tibble: 10 × 8
## mtry min_n .metric .estimator mean n std_err .config
## <int> <int> <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 1 14 accuracy multiclass 0.602 5 0.0120 Preprocessor1_Model01
## 2 6 31 accuracy multiclass 0.583 5 0.0230 Preprocessor1_Model02
## 3 11 2 accuracy multiclass 0.589 5 0.0218 Preprocessor1_Model03
## 4 17 18 accuracy multiclass 0.577 5 0.0222 Preprocessor1_Model04
## 5 22 35 accuracy multiclass 0.567 5 0.0227 Preprocessor1_Model05
## 6 28 6 accuracy multiclass 0.567 5 0.0194 Preprocessor1_Model06
## 7 33 23 accuracy multiclass 0.583 5 0.0196 Preprocessor1_Model07
## 8 39 40 accuracy multiclass 0.570 5 0.0201 Preprocessor1_Model08
## 9 44 10 accuracy multiclass 0.586 5 0.0216 Preprocessor1_Model09
## 10 50 27 accuracy multiclass 0.583 5 0.0216 Preprocessor1_Model10- SVM
svm_spec <- svm_rbf(cost = tune(), rbf_sigma = tune()) %>%
set_mode("classification") %>%
set_engine("kernlab")
svm_workflow <- workflow() %>%
add_recipe(campy_recipe) %>%
add_model(svm_spec)
svm_res <- tune_grid(
svm_workflow,
resamples = vfold_cv(campy_train_split, v = 5, strata = Source),
grid = 10,
metrics = metric_set(accuracy)
)
best_svm <- svm_res %>%
select_best(metric = "accuracy")
final_svm <- finalize_workflow(svm_workflow, best_svm)
print("Best SVM configuration:")## [1] "Best SVM configuration:"print(best_svm)## # A tibble: 1 × 3
## cost rbf_sigma .config
## <dbl> <dbl> <chr>
## 1 32 0.000464 Preprocessor1_Model10# Showing all performances metrics from the tuning
svm_metrics <- collect_metrics(svm_res)
print("SVM Cross-Validation Results:")## [1] "SVM Cross-Validation Results:"print(svm_metrics)## # A tibble: 10 × 8
## cost rbf_sigma .metric .estimator mean n std_err .config
## <dbl> <dbl> <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 0.000977 0.000000215 accuracy multiclass 0.348 5 0.00106 Preprocessor…
## 2 0.00310 0.00599 accuracy multiclass 0.348 5 0.00106 Preprocessor…
## 3 0.00984 0.0000000001 accuracy multiclass 0.348 5 0.00106 Preprocessor…
## 4 0.0312 0.00000278 accuracy multiclass 0.348 5 0.00106 Preprocessor…
## 5 0.0992 0.0774 accuracy multiclass 0.551 5 0.0155 Preprocessor…
## 6 0.315 0.00000000129 accuracy multiclass 0.348 5 0.00106 Preprocessor…
## 7 1 0.0000359 accuracy multiclass 0.348 5 0.00106 Preprocessor…
## 8 3.17 1 accuracy multiclass 0.533 5 0.0203 Preprocessor…
## 9 10.1 0.0000000167 accuracy multiclass 0.348 5 0.00106 Preprocessor…
## 10 32 0.000464 accuracy multiclass 0.605 5 0.0300 Preprocessor…best_rf_acc <- best_rf %>% inner_join(rf_metrics, by = c("mtry", "min_n"))
best_svm_acc <- best_svm %>% inner_join(svm_metrics, by = c("cost", "rbf_sigma"))
print("Top Random Forest Accuracy:")## [1] "Top Random Forest Accuracy:"print(best_rf_acc)## # A tibble: 1 × 9
## mtry min_n .config.x .metric .estimator mean n std_err .config.y
## <int> <int> <chr> <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 1 14 Preprocessor1_Mo… accura… multiclass 0.602 5 0.0120 Preproce…print("Top SVM Accuracy:")## [1] "Top SVM Accuracy:"print(best_svm_acc)## # A tibble: 1 × 9
## cost rbf_sigma .config.x .metric .estimator mean n std_err .config.y
## <dbl> <dbl> <chr> <chr> <chr> <dbl> <int> <dbl> <chr>
## 1 32 0.000464 Preprocessor… accura… multiclass 0.605 5 0.0300 Preproce…Fit final Model on full train set
final_model <- last_fit(final_rf, campy_split)
collect_metrics(final_model)## # A tibble: 3 × 4
## .metric .estimator .estimate .config
## <chr> <chr> <dbl> <chr>
## 1 accuracy multiclass 0.679 Preprocessor1_Model1
## 2 roc_auc hand_till 0.857 Preprocessor1_Model1
## 3 brier_class multiclass 0.204 Preprocessor1_Model1Predict on Test Set and written in the csv file
final_rf_fit <- fit(final_rf, data = campy.train)
campy.test.with.classifications <- predict(final_rf_fit, new_data = campy.test) %>%
bind_cols(campy.test)
# Ensure column name is .pred_class
head(campy.test.with.classifications$.pred_class)## [1] Sheep Sheep Sheep Sheep Sheep Sheep
## Levels: Beef Poultry Sheep# Save to CSV
write_csv(campy.test.with.classifications, "ex2_24020251.csv")Exercise 3: Clustering the New Zealand World Values Survey
3.1 Initial exploration of the dataset
# View structure and glimpse the data
glimpse(wv_survey)## Rows: 841
## Columns: 10
## $ Female <dbl> 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, …
## $ Age <dbl> 69, 29, NA, 22, 27, 61, 25, 56, 54, 35, 72, 37, 76,…
## $ StateofHealth <dbl> 2, 2, 2, 1, 2, 1, 4, 1, 2, 1, 1, 1, 2, 2, 1, 1, 1, …
## $ MaritalStatus <dbl> 1, 2, 1, 1, 6, 1, 2, 1, 1, 6, 1, 1, 5, 1, 6, 1, 2, …
## $ HighestEducation <dbl> 5, 3, NA, 4, 4, 2, 4, 5, 5, 4, 5, 5, 2, 2, 3, 5, 2,…
## $ LifeSatisfaction <dbl> 8, 3, 8, 9, 8, 10, 3, 7, 10, 7, 9, 10, 5, 10, 9, 10…
## $ IncomeDecile <dbl> 8, 8, 8, 3, 4, 3, 5, 10, 3, 6, 5, 10, 1, 4, 3, 8, 2…
## $ SocialClass <dbl> 4, 4, 3, 3, 4, 4, 4, 3, 3, NA, 2, 2, 4, 4, 5, NA, 3…
## $ BelieveinGod <dbl> NA, 1, NA, 1, 1, NA, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, …
## $ FeelingofHappiness <dbl> 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 1, 1, 2, …# Summary statistics for numeric variables
summary(wv_survey)## Female Age StateofHealth MaritalStatus
## Min. :0.0000 Min. :18.00 Min. :1.000 Min. :1.000
## 1st Qu.:0.0000 1st Qu.:39.00 1st Qu.:1.000 1st Qu.:1.000
## Median :1.0000 Median :51.50 Median :2.000 Median :1.000
## Mean :0.5772 Mean :51.44 Mean :1.854 Mean :2.271
## 3rd Qu.:1.0000 3rd Qu.:64.00 3rd Qu.:2.000 3rd Qu.:3.000
## Max. :1.0000 Max. :90.00 Max. :4.000 Max. :6.000
## NA's :6 NA's :13 NA's :13 NA's :29
## HighestEducation LifeSatisfaction IncomeDecile SocialClass
## Min. :1.000 Min. : 1.000 Min. : 1.000 Min. :1.00
## 1st Qu.:2.000 1st Qu.: 7.000 1st Qu.: 3.000 1st Qu.:2.00
## Median :4.000 Median : 8.000 Median : 6.000 Median :3.00
## Mean :3.502 Mean : 7.649 Mean : 5.722 Mean :3.12
## 3rd Qu.:5.000 3rd Qu.: 9.000 3rd Qu.: 8.000 3rd Qu.:4.00
## Max. :5.000 Max. :10.000 Max. :10.000 Max. :5.00
## NA's :34 NA's :12 NA's :71 NA's :99
## BelieveinGod FeelingofHappiness
## Min. :0.0000 Min. :1.000
## 1st Qu.:0.0000 1st Qu.:1.000
## Median :1.0000 Median :2.000
## Mean :0.7063 Mean :1.713
## 3rd Qu.:1.0000 3rd Qu.:2.000
## Max. :1.0000 Max. :4.000
## NA's :160 NA's :22# Distribution plots for each variable
# Numeric variables: Age, LifeSatisfaction, IncomeDecile
numeric_vars <- c("Age", "LifeSatisfaction", "IncomeDecile")
# Categorical variables: Female, StateofHealth, MaritalStatus, HighestEducation, SocialClass, BelieveinGod, FeelingofHappiness
categorical_vars <- c("Female", "StateofHealth", "MaritalStatus", "HighestEducation", "SocialClass", "BelieveinGod", "FeelingofHappiness")
# Histograms for numeric variables
wv_survey %>%
select(all_of(numeric_vars)) %>%
pivot_longer(everything()) %>%
ggplot(aes(x = value)) +
geom_histogram(bins = 30, fill = "skyblue", color = "black") +
facet_wrap(~ name, scales = "free_x") +
labs(title = "Histograms of Numeric Variables")## Warning: Removed 96 rows containing non-finite outside the scale range
## (`stat_bin()`).
# Bar plots for categorical variables
for (var in categorical_vars) {
print(
wv_survey %>%
filter(!is.na(.data[[var]])) %>%
count(.data[[var]]) %>%
ggplot(aes(x = factor(.data[[var]]), y = n)) +
geom_bar(stat = "identity", fill = "coral") +
labs(title = paste("Distribution of", var),
x = var,
y = "Count")
)
}
# You can also do frequency tables:
for (var in c(categorical_vars, numeric_vars)) {
print(paste("Summary for", var))
print(table(wv_survey[[var]], useNA = "ifany"))
print(summary(wv_survey[[var]]))
}## [1] "Summary for Female"
##
## 0 1 <NA>
## 353 482 6
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0000 0.0000 1.0000 0.5772 1.0000 1.0000 6
## [1] "Summary for StateofHealth"
##
## 1 2 3 4 <NA>
## 296 380 129 23 13
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1.000 1.000 2.000 1.854 2.000 4.000 13
## [1] "Summary for MaritalStatus"
##
## 1 2 3 4 5 6 <NA>
## 477 113 34 20 49 119 29
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1.000 1.000 1.000 2.271 3.000 6.000 29
## [1] "Summary for HighestEducation"
##
## 1 2 3 4 5 <NA>
## 25 209 150 182 241 34
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1.000 2.000 4.000 3.502 5.000 5.000 34
## [1] "Summary for SocialClass"
##
## 1 2 3 4 5 <NA>
## 1 216 240 263 22 99
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1.00 2.00 3.00 3.12 4.00 5.00 99
## [1] "Summary for BelieveinGod"
##
## 0 1 <NA>
## 200 481 160
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0000 0.0000 1.0000 0.7063 1.0000 1.0000 160
## [1] "Summary for FeelingofHappiness"
##
## 1 2 3 4 <NA>
## 283 494 36 6 22
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1.000 1.000 2.000 1.713 2.000 4.000 22
## [1] "Summary for Age"
##
## 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
## 6 11 7 4 9 8 7 7 9 12 11 9 9 11 6 11
## 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
## 9 15 7 19 16 15 15 13 14 18 17 16 12 21 18 22
## 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
## 14 16 23 18 16 17 20 15 13 11 13 20 17 12 13 14
## 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
## 18 7 14 8 16 16 8 13 9 13 14 9 8 5 6 4
## 82 83 84 85 86 87 88 89 90 <NA>
## 4 3 2 3 4 2 1 3 2 13
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 18.00 39.00 51.50 51.44 64.00 90.00 13
## [1] "Summary for LifeSatisfaction"
##
## 1 2 3 4 5 6 7 8 9 10 <NA>
## 9 8 20 24 50 72 126 234 136 150 12
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1.000 7.000 8.000 7.649 9.000 10.000 12
## [1] "Summary for IncomeDecile"
##
## 1 2 3 4 5 6 7 8 9 10 <NA>
## 54 79 71 87 72 71 86 82 91 77 71
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1.000 3.000 6.000 5.722 8.000 10.000 71Some groups may be more or less willing to participate in the World Values Survey, leading to bias. For example, older adults or those with lower education may be less likely to respond, skewing results toward younger or more educated respondents. Sensitive questions like income or belief in God may have higher non-response rates. This can bias the sample away from representing the full population’s true distribution of values.
For missings values
missing_counts <- sapply(wv_survey, function(x) sum(is.na(x)))
missing_counts_sorted <- sort(missing_counts, decreasing = TRUE)
missing_counts_sorted[1:3] # top 3 variables with most missing## BelieveinGod SocialClass IncomeDecile
## 160 99 713.2 Normalised dataset
library(caret)## Loading required package: lattice##
## Attaching package: 'caret'## The following objects are masked from 'package:yardstick':
##
## precision, recall, sensitivity, specificity## The following object is masked from 'package:purrr':
##
## lift# Select variables for clustering
vars_to_normalize <- c("Age", "Female", "MaritalStatus", "StateofHealth", "HighestEducation", "LifeSatisfaction", "IncomeDecile", "SocialClass", "BelieveinGod", "FeelingofHappiness")
wv_survey_norm <- wv_survey %>%
select(all_of(vars_to_normalize)) %>%
na.omit() %>%
mutate(across(everything(), scale))Variables have different units and ranges — e.g., Age can be from 18 to 90+, LifeSatisfaction from 1 to 10, Female is binary (0 or 1), and others are categorical codes. Without normalization, variables with larger scales (like Age) would dominate distance calculations in clustering. Normalization puts all variables on the same scale so no single variable dominates the clustering outcome.
3.3 Dendrogram
library(cluster)
library(factoextra)## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa# Prepare data
cluster_vars <- c("Age", "Female", "MaritalStatus", "StateofHealth", "HighestEducation", "LifeSatisfaction")
wv_cluster <- wv_survey %>%
select(all_of(cluster_vars)) %>%
na.omit() %>%
mutate(across(everything(), scale))
# Distance matrix using Euclidean distance
dist_matrix <- dist(wv_cluster)
# Hierarchical clustering with complete linkage
hc_complete <- hclust(dist_matrix, method = "complete")
# Plot dendrogram
fviz_dend(hc_complete, k = 4, # to highlight 4 clusters, just an example
cex = 0.5, # label size
main = "Dendrogram of WVS data (Complete linkage)")## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
## Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
fviz_dend(hc_complete, k = 4,
cex = 0.5,
main = "Dendrogram of WVS data (Complete linkage)",
rect = TRUE)
### 3.4 k-means
set.seed(123)
wss <- numeric(10)
avg_sil <- numeric(10)
for (k in 1:10) {
km <- kmeans(wv_cluster, centers = k, nstart = 25)
wss[k] <- km$tot.withinss
# silhouette width only for k > 1
if (k > 1) {
ss <- silhouette(km$cluster, dist(wv_cluster))
avg_sil[k] <- mean(ss[, 3])
} else {
avg_sil[k] <- NA
}
}
# Plot total within sum of squares (Elbow method)
plot(1:10, wss, type = "b", xlab = "Number of clusters K", ylab = "Total within-cluster sum of squares", main = "Elbow Method")
# Plot average silhouette width
plot(2:10, avg_sil[2:10], type = "b", xlab = "Number of clusters K", ylab = "Average silhouette width", main = "Silhouette Method")
# Choose k based on elbow and silhouettekm_final <- kmeans(wv_cluster, centers = 4, nstart = 25)
fviz_silhouette(silhouette(km_final$cluster, dist(wv_cluster)))## cluster size ave.sil.width
## 1 1 136 0.07
## 2 2 228 0.19
## 3 3 112 0.19
## 4 4 282 0.29
### 3.5 Visualisation of heirarchical clustering
fviz_cluster(list(data = wv_cluster, cluster = cutree(hc_complete, k = 4)),
geom = "point",
ellipse.type = "convex",
main = "Clusters from hierarchical clustering (k=4)")
fviz_cluster(list(data = wv_cluster, cluster = cutree(hc_complete, k = 5)),
geom = "point",
ellipse.type = "convex",
main = "Clusters from hierarchical clustering (k=5)")
fviz_cluster(list(data = wv_cluster, cluster = cutree(hc_complete, k = 3)),
geom = "point",
ellipse.type = "convex",
main = "Clusters from hierarchical clustering (k=3)")
fviz_cluster(list(data = wv_cluster, cluster = cutree(hc_complete, k = 2)),
geom = "point",
ellipse.type = "convex",
main = "Clusters from hierarchical clustering (k=2)")
### 3.6 Visualisation of k-means clusters
# Add cluster membership to data
wv_cluster$cluster <- factor(km_final$cluster)
# Gather data for boxplot
wv_long <- wv_cluster %>%
pivot_longer(cols = -cluster, names_to = "variable", values_to = "value")
ggplot(wv_long, aes(x = cluster, y = value, fill = cluster)) +
geom_boxplot() +
facet_wrap(~variable, scales = "free") +
theme_minimal() +
labs(title = "Boxplots of variables by k-means clusters")
Cluster 1 (Red): Older, mostly married individuals in poorer health and with lower life satisfaction. Moderately educated.
Cluster 2 (Green): Younger, highly satisfied with life, less educated, more likely single or cohabiting, and in decent health.
Cluster 3 (Cyan): Middle-aged, highly educated, moderate health, moderate satisfaction. Mixed marital status.
Cluster 4 (Purple): Well-educated individuals with high life satisfaction, good health, and a balanced age group.
Exercise 4
4.1 Frequency of colour use
# Frequency plot
itemFrequencyPlot(flags, topN = 10, type = "absolute", main = "Flag Colour Frequencies")
# Find top 3 colours
item_freq <- itemFrequency(flags, type = "absolute")
sort(item_freq, decreasing = TRUE)[1:3]## Red White Blue
## 165 151 954.2 Number of rules
rules <- apriori(flags, parameter = list(support = 0.023, confidence = 0.5))## Apriori
##
## Parameter specification:
## confidence minval smax arem aval originalSupport maxtime support minlen
## 0.5 0.1 1 none FALSE TRUE 5 0.023 1
## maxlen target ext
## 10 rules TRUE
##
## Algorithmic control:
## filter tree heap memopt load sort verbose
## 0.1 TRUE TRUE FALSE TRUE 2 TRUE
##
## Absolute minimum support count: 4
##
## set item appearances ...[0 item(s)] done [0.00s].
## set transactions ...[12 item(s), 212 transaction(s)] done [0.00s].
## sorting and recoding items ... [9 item(s)] done [0.00s].
## creating transaction tree ... done [0.00s].
## checking subsets of size 1 2 3 4 5 done [0.00s].
## writing ... [77 rule(s)] done [0.00s].
## creating S4 object ... done [0.00s].length(rules)## [1] 774.3 Highest support rules
# Sort rules by support
rules_sorted_support <- sort(rules, by = "support", decreasing = TRUE)
# Get top 5
inspect(rules_sorted_support[1:5])## lhs rhs support confidence coverage lift count
## [1] {} => {Red} 0.7783019 0.7783019 1.0000000 1.000000 165
## [2] {} => {White} 0.7122642 0.7122642 1.0000000 1.000000 151
## [3] {White} => {Red} 0.5566038 0.7814570 0.7122642 1.004054 118
## [4] {Red} => {White} 0.5566038 0.7151515 0.7783019 1.004054 118
## [5] {Blue} => {White} 0.3349057 0.7473684 0.4481132 1.049285 714.4 Lift greater than 2
high_lift_rules <- subset(rules, subset = lift > 2)
length(high_lift_rules)## [1] 2inspect(high_lift_rules)## lhs rhs support confidence coverage
## [1] {Red, White, Yellow} => {Black} 0.05660377 0.5454545 0.10377358
## [2] {Green, Red, White, Yellow} => {Black} 0.02830189 0.5000000 0.05660377
## lift count
## [1] 2.569697 12
## [2] 2.355556 6There are 2 rules with lift greater than 2.
These rules suggest that Black appears much more frequently than expected when certain combinations of other colours (e.g., Red, White, Yellow) are present.
For instance, in the rule {Red, White, Yellow} => {Black}, the lift of 2.57 indicates that Black is 2.57 times more likely to appear with those colours than by random chance. This suggests a strong association between Black and those colours in flags.
4.5 Choiseul independence flag
black_rules <- subset(rules, lhs %in% "Black")
inspect(sort(black_rules, by = "confidence", decreasing = TRUE)[1:5])## lhs rhs support confidence coverage
## [1] {Black, Blue, Red} => {White} 0.02358491 1.0000000 0.02358491
## [2] {Black, Green, White, Yellow} => {Red} 0.02830189 1.0000000 0.02830189
## [3] {Black, Green, White} => {Red} 0.08018868 0.9444444 0.08490566
## [4] {Black, White, Yellow} => {Red} 0.05660377 0.9230769 0.06132075
## [5] {Black, Green, Yellow} => {Red} 0.04716981 0.9090909 0.05188679
## lift count
## [1] 1.403974 5
## [2] 1.284848 6
## [3] 1.213468 17
## [4] 1.186014 12
## [5] 1.168044 10itemLabels(flags)## [1] "Aquamarine Blue" "Black" "Blue" "Dark Blue"
## [5] "Gold" "Green" "Light Blue" "Maroon"
## [9] "Orange" "Red" "White" "Yellow"black_rules <- subset(rules, lhs %in% "Black")
inspect(sort(black_rules, by = "confidence", decreasing = TRUE)[1:5])## lhs rhs support confidence coverage
## [1] {Black, Blue, Red} => {White} 0.02358491 1.0000000 0.02358491
## [2] {Black, Green, White, Yellow} => {Red} 0.02830189 1.0000000 0.02830189
## [3] {Black, Green, White} => {Red} 0.08018868 0.9444444 0.08490566
## [4] {Black, White, Yellow} => {Red} 0.05660377 0.9230769 0.06132075
## [5] {Black, Green, Yellow} => {Red} 0.04716981 0.9090909 0.05188679
## lift count
## [1] 1.403974 5
## [2] 1.284848 6
## [3] 1.213468 17
## [4] 1.186014 12
## [5] 1.168044 10#combo_rules <- subset(rules, lhs %ain% c("Black", "Blue", "Red") & length(lhs) == 3)
#inspect(sort(combo_rules, by = "confidence", decreasing = TRUE)[1:5])
# Subset rules where lhs contains ALL of Black, Blue, and Red (but maybe more too)
#combo_rules <- subset(rules, lhs %ain% c("Black", "Blue", "Red"))
#combo_rules <- subset(rules, lhs %in% list(c("Black", "Blue", "Red")))
#inspect(combo_rules)
# Now sort and inspect top 5
#inspect(sort(combo_rules, by = "confidence", decreasing = TRUE)[1:5])