# ML clustering analysis # load packages and prep data---- pacman::p_load(here, rio, here, janitor, tidyverse, mclust, vcd, broom, gtsummary, patchwork, sjPlot, cluster, ggdist, ggridges, tidyquant) ml_prep <- import(here("ML_analysis", "Vitality","ml_preprocessed.dta")) # variables in order to make it easier to select variables for clustering analysis and drop NAs ml_data <- ml_prep |> select(-c(num_timepoints)) |> select(tot1:tot24, contains("s_"), everything()) |> drop_na() # analyses done on dataset including n=86 individuals with no missing data # check variance apply(ml_data[,2:25], 2, var) # 1. K-means clustering algorithm ---- ### selecting k ---- # z.perc.between <- numeric(1) # for (i in 1:5) { # z <- kmeans(ml_data[,2:25], centers = i, nstart = 1) # z.perc.between[i] <- 100*z$betweenss/z$totss # } # z.perc.between[1:5] # plot(1:5, z.perc.between[1:5], type = "l", ylab = "Percentage of variability between clusters",xlab = "k") ### run kmeans algorithm ---- set.seed(7) kmeans <- kmeans(ml_data[,1:24], centers = 3, nstart = 10) kmeans # add cluster id to master dataset ml_data$cluster_id_kmeans <- kmeans$cluster kmeans[["centers"]] kmeans$centers[2,] #overall adherence mean(kmeans$centers[1,1:24]) # 64% mean(kmeans$centers[2,1:24]) # 24% mean(kmeans$centers[3,1:24]) # 88% ml_data |> filter(cluster_id_kmeans == 1) |> select(tot1:tot24) |> colMeans(1:24) |> mean() ml_data |> filter(cluster_id_kmeans == 2) |> select(tot1:tot24) |> colMeans(1:24) |> mean() ml_data |> filter(cluster_id_kmeans == 3) |> select(tot1:tot24) |> colMeans(1:24) |> mean() ml_data |> filter(cluster_id_kmeans == 1) |> select(age) |> pull() |> mean() ml_data |> filter(cluster_id_kmeans == 2) |> select(age) |> pull() |> mean() ml_data |> filter(cluster_id_kmeans == 3) |> select(age) |> pull() |> mean() ml_data |> filter(cluster_id_kmeans == 1) |> select(sesquint) |> pull() |> median() ml_data |> filter(cluster_id_kmeans == 2) |> select(sesquint) |> pull() |> median() ml_data |> filter(cluster_id_kmeans == 3) |> select(sesquint) |> pull() |> median() ml_data |> filter(cluster_id_kmeans == 1) |> select(male) |> pull() |> mean() ml_data |> filter(cluster_id_kmeans == 2) |> select(male) |> pull() |> mean() ml_data |> filter(cluster_id_kmeans == 3) |> select(male) |> pull() |> mean() ml_data |> filter(cluster_id_kmeans == 1) |> select(orphan) |> pull() |> mean() ml_data |> filter(cluster_id_kmeans == 2) |> select(orphan) |> pull() |> mean() ml_data |> filter(cluster_id_kmeans == 3) |> select(orphan) |> pull() |> mean() # K means # cluster 1 mean age 15.9 ses 2 male 46% orphan 58% adherence 64% # cluster 2 mean age 18.0 ses 4 male 64% orphan 79% adherence 24% # cluster 3 mean age 14.7 ses 2 male 44% orphan 50% adherence 88% # 2. Hierarchical clustering algorithm ---- library(cluster) ml_ag <- agnes(ml_data[,1:24], diss = FALSE, metric = "euclidean", method = "ward") plot(ml_ag,main = "Average Link, 26 Variables", labels = ml_data$study_id, which.plot = 2) ## dendrogram rect.hclust(ml_ag, k = 3, border = "blue") groups.3 <- cutree(ml_ag,3) # store the results groups.3 ml_data$cluster_id_agg <- groups.3 # use aggregate to get summaries by cluster id - quicker than method using tidyverse aggregate(ml_data,list(ml_data$cluster_id_kmeans),mean) # list needed to specify subgroups aggregate(ml_data,list(ml_data$cluster_id_agg),mean) # list needed to specify subgroups aggregate(ml_data,list(ml_data$cluster_id_agg),median) # list needed to specify subgroups ml_data |> filter(cluster_id_agg==1) |> select(tot1:tot24) |> colMeans(1:24) |> mean() ml_data |> filter(cluster_id_agg==2) |> select(tot1:tot24) |> colMeans(1:24) |> mean() ml_data |> filter(cluster_id_agg==3) |> select(tot1:tot24) |> colMeans(1:24) |> mean() # Hierarchical clustering # cluster 1 mean age 14.8 ses 2 male 48% orphan 50% adherence 88% # cluster 2 mean age 15.7 ses 2 male 38% orphan 58% adherence 67% # cluster 3 mean age 17.6 ses 4 male 63% orphan 75% adherence 28% # 3. Gaussian mixture modelling using mclust package ---- m1 <- Mclust(ml_data[,1:24]) summary(m1) ml_data$cluster_id_gmm <- m1$classification m1[["classification"]] # adherence mean(m1$data[m1$classification==1,1:24]) # 64% mean(m1$data[m1$classification==2,1:24]) # 88% mean(m1$data[m1$classification==3,1:24]) # 24% aggregate(ml_data,list(ml_data$cluster_id_gmm),median) # list needed to specify subgroups aggregate(ml_data,list(ml_data$cluster_id_gmm),mean) # list needed to specify subgroups # GMM # cluster 1 mean age 15.9 ses 2 male 44% orphan 58% adherence 64% # cluster 2 mean age 14.7 ses 2 male 46% orphan 50% adherence 88% # cluster 3 mean age 18.0 ses 4 male 64% orphan 78% adherence 28% # 4. K-medoids clustering ---- kmedoids <- pam(ml_data[,2:25], 3) summary(kmedoids) # add cluster id to master dataset ml_data$cluster_id_kmedoids <- kmedoids$clustering ml_data |> filter(cluster_id_kmedoids == 1) |> select(tot1:tot24) |> colMeans(1:24) |> mean() # 85% ml_data |> filter(cluster_id_kmedoids == 2) |> select(tot1:tot24) |> colMeans(1:24) |> mean() # 25% ml_data |> filter(cluster_id_kmedoids == 3) |> select(tot1:tot24) |> colMeans(1:24) |> mean() # 65% ml_data |> filter(cluster_id_kmedoids == 1) |> select(age) |> pull() |> mean() # 14.8 ml_data |> filter(cluster_id_kmedoids == 2) |> select(age) |> pull() |> mean() # 18 ml_data |> filter(cluster_id_kmedoids == 3) |> select(age) |> pull() |> mean() # 16.1 ml_data |> filter(cluster_id_kmedoids == 1) |> select(sesquint) |> pull() |> median() # 2 ml_data |> filter(cluster_id_kmedoids == 2) |> select(sesquint) |> pull() |> median() # 4 ml_data |> filter(cluster_id_kmedoids == 3) |> select(sesquint) |> pull() |> median() # 2 ml_data |> filter(cluster_id_kmedoids == 1) |> select(male) |> pull() |> mean() # 45% ml_data |> filter(cluster_id_kmedoids == 2) |> select(male) |> pull() |> mean() # 60% ml_data |> filter(cluster_id_kmedoids == 3) |> select(male) |> pull() |> mean() # 47% ml_data |> filter(cluster_id_kmedoids == 1) |> select(orphan) |> pull() |> mean() # 52% ml_data |> filter(cluster_id_kmedoids == 2) |> select(orphan) |> pull() |> mean() # 73% ml_data |> filter(cluster_id_kmedoids == 3) |> select(orphan) |> pull() |> mean() # 60% # K-medoids # cluster 1 mean age 14.8 ses 2 male 45% orphan 52% adherence 85% # cluster 2 mean age 18 ses 4 male 60% orphan 73% adherence 25% # cluster 3 mean age 16.1 ses 2 male 47% orphan 60% adherence 65% # active Vit D arm ---- tapply(ml_data$base25OHD, ml_data$arm, summary) # $`1` # Min. 1st Qu. Median Mean 3rd Qu. Max. # 37.90 55.40 63.20 64.77 71.10 120.80 # # $`2` # Min. 1st Qu. Median Mean 3rd Qu. Max. # 42.00 57.60 65.30 67.71 77.60 107.30 tapply(ml_data$wk4825OHD, ml_data$arm, summary) # $`1` # Min. 1st Qu. Median Mean 3rd Qu. Max. # 51.30 74.40 94.70 91.43 107.40 138.40 # # $`2` # Min. 1st Qu. Median Mean 3rd Qu. Max. # 39.20 59.00 65.40 66.55 73.40 96.00 ml_data_vitD <- ml_data |> filter(arm == 1) |> mutate(adherence = factor( case_when( cluster_id_gmm == 2 ~ "High", cluster_id_gmm == 1 ~ "Moderate", cluster_id_gmm == 3 ~ "Low" )) |> fct_relevel("Low", "Moderate", "High") ) # change the reference group ml_data_vitD$cluster_id_gmm <- factor(ml_data_vitD$cluster_id_gmm) ml_data_vitD$cluster_id_gmm <- relevel(ml_data_vitD$cluster_id_gmm, ref = "3") m1 <- lm(wk4825OHD ~ factor(cluster_id_gmm), data = ml_data_vitD) summary(m1) tapply(ml_data_vitD$wk4825OHD, ml_data_vitD$cluster_id_gmm, summary) # raincloud col <- thematic::okabe_ito(4) colors <- c("#F0E442", "#56B4E9","#D55E00") vitD <- ml_data_vitD |> ggplot(aes(x = adherence, y = wk4825OHD, fill = adherence)) + # add half-violin from ggdist package ggdist::stat_halfeye( ## custom bandwidth adjust = 0.5, ## move geom to the right justification = -.2, ## remove slab interval .width = 0, point_colour = NA, scale = 0.5 ) + geom_boxplot( width = .12, ## remove outliers outlier.color = NA, alpha = 0.5 ) + stat_summary( geom = "text", fun = "median", aes(label = round(..y.., 2)), fontface = "bold", vjust = 2.35 ) + #Adjust theme scale_fill_tq() + theme_tq(18) + labs( title = "", subtitle = "25(OH)D by adherence cluster", y = "25(OH)D", x = "", fill = "Adherence cluster" ) + coord_flip() + #theme(legend.position="none") theme(axis.text.y = element_blank(), #remove x axis labels axis.ticks.y = element_blank()) +#remove x axis ticks scale_colour_manual(values = colors) + theme(legend.position = "right") + guides(fill = guide_legend(reverse = T)) vitD