--- title: "VietNamI Dataset Analysis" author: "Analysis" date: "`r Sys.Date()`" output: "html_notebook" --- # Import libraries ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) # Import libraries library(ggplot2) library(tidyverse) library(stringr) library(plotly) library(factoextra) library(cluster) library(mclust) library(kernlab) library(GGally) library(parallel) library(doParallel) library(foreach) library(psych) library(bootstrap) library(boot) ``` \# ------------------------------------------------- \# 1. PREPROCESSING \# ------------------------------------------------- ### 1. Quick study and inspection of the dataset ```{r Quick-study-and-inspection} readLines("VietNamI.csv",2) db = read.csv(file ="VietNamI.csv", header = TRUE, sep = ",", dec=".", stringsAsFactors = TRUE) summary(db) str(db) glimpse(db) ``` ## 1.1 Remove useless rows, NAs and duplicated rows ```{r Remove-rows} db$rownames <- NULL colSums(is.na(db)) sum(duplicated(db)) db <- db[!duplicated(db), ] ``` ## 1.2 Convert categorical variables to factor ```{r Convert-to-factor} db$married <- as.factor(db$married) db$injury <- as.factor(db$injury) db$insurance <- as.factor(db$insurance) db$commune <- as.factor(db$commune) ``` ## 1.3 Age and household total expenditure are in logarithm form, so transform them # into real age and expenses ```{r Transform-log-variables} db$age_real <- round(exp(db$age)) summary(db$age_real) db$hhexp <- exp(db$lnhhexp) summary(db$hhexp) ``` ## 1.4 Look for impossible values and remove them ```{r Check-impossible-values} summary(db$age_real) summary(db$educ) summary(db$illdays) summary(db$actdays) summary(db$pharvis) summary(db$illness) table(db$insurance) table(db$sex) table(db$married) ``` ### 2 Visualize outliers # Pharvis ```{r Plot-Pharvis-Outliers} ggplot(db, aes(x = pharvis)) + geom_boxplot(fill = "lightblue", color = "blue", outlier.color = "red", outlier.shape = 16) + theme_minimal() + labs(title = "Outliers in Pharvis", x = "Pharvis", y = "") ``` # Illdays ```{r Plot-Illdays-Outliers} ggplot(db, aes(x = illdays)) + geom_boxplot(fill = "lightgreen", color = "darkgreen", outlier.color = "red", outlier.shape = 16) + theme_minimal() + labs(title = "Outliers in Illdays", x = "Illdays", y = "") ``` # Actdays ```{r Plot-Actdays-Outliers} ggplot(db, aes(x = actdays)) + geom_boxplot(fill = "lightpink", color = "darkred", outlier.color = "red", outlier.shape = 16) + theme_minimal() + labs(title = "Outliers in Actday", x = "Actdays", y = "") ``` ## 2.1 Confirm outliers ```{r Confirm-Outliers} # Pharvis mu <- mean(db$pharvis) sigma <- sd(db$pharvis) sum(db$pharvis < mu - 3*sigma | db$pharvis > mu + 3*sigma) QI <- quantile(db$pharvis, 0.25) QS <- quantile(db$pharvis, 0.75) IQR <- QS - QI sum(db$pharvis < QI - 1.5*IQR | db$pharvis > QS + 1.5*IQR) # Illdays mu <- mean(db$illdays) sigma <- sd(db$illdays) sum(db$illdays < mu - 3*sigma | db$illdays > mu + 3*sigma) QI <- quantile(db$illdays, 0.25) QS <- quantile(db$illdays, 0.75) IQR <- QS - QI sum(db$illdays < QI - 1.5*IQR | db$illdays > QS + 1.5*IQR) # Actdays mu <- mean(db$actdays) sigma <- sd(db$actdays) sum(db$actdays < mu - 3*sigma | db$actdays > mu + 3*sigma) QI <- quantile(db$actdays, 0.25) QS <- quantile(db$actdays, 0.75) IQR <- QS - QI sum(db$actdays < QI - 1.5*IQR | db$actdays > QS + 1.5*IQR) # Outliers were detected in pharvis, illdays, and actdays using the 3-sigma rule # and the (IQR) criterion. # However, these outliers correspond to individuals with high medical expenses or # long periods of illness. These outliers could represent realistic cases of people # heavily affected by diseases or illnesses, so we keep them ``` ### 3. Feature engineering # In this point we will create new useful variables form the ones that already exist ```{r Feature-Engineering} # is_ill is a factor variable that works saying TRUE if the person was ill # or FALSE if not db$is_ill <- db$illness > 0 # age_group differentiates by age groups db$age_group <- cut(db$age_real, breaks = c(0, 20, 40, 60, 100), labels = c("Young", "Adult", "Senior", "Elder"), include.lowest = TRUE) # education_level differentiates by education level db$education_level <- cut(db$educ, breaks = c(0, 3, 7, 11), labels = c("Low", "Medium", "High"), include.lowest = TRUE) # People with pharmacy visits weirdly high (detect risk profiles) Q90 <- quantile(db$pharvis, 0.9) db$high_pharvis <- db$pharvis > Q90 # health_index combines the number of sick days and the presence of injuries # into a single measure, where higher values indicate better physical health db$health_index <- pmax(0, (1 - db$illdays/60) + (1 - as.numeric(db$injury))) # wealth_group differentiates by wealth level breaks <- quantile(db$lnhhexp, probs = seq(0, 1, length.out = 6)) db$wealth_group <- cut(db$lnhhexp, breaks = breaks, labels = c("Very Low", "Low", "Medium", "High", "Very High"), include.lowest = TRUE) ``` \# ------------------------------------------------- \# 2. VISUALIZACIÓN \# ------------------------------------------------- # Question 1. Does having insurance reduce pharmacy visits? ```{r Plot-Q1} p = db %>% ggplot(aes(x = as.factor(insurance), y = pharvis, fill = as.factor(insurance))) + geom_boxplot(outlier.color = "red", alpha = 0.7) + scale_y_log10() + scale_fill_manual(values = c("lightblue", "orange")) + theme_minimal() + labs(title = "Pharmacy visits by insurance status (Log Scale)", x = "Insurance (0 = No, 1 = Yes)", y = "Pharmacy visits (log scale)", fill = "Insurance") ggplotly(p) # With this graphic we concluded that having insurance does not change much the # usual number of pharmacy visits. However, people without insurance show more # extreme values, meaning some of them go to the pharmacy much more often. ``` # Question 2. Are insured households economically wealthier? # In this case I use the log to see more clearly the difference between household # total expenses ```{r Plot-Q2} p = db %>% ggplot(aes(x = lnhhexp, fill = as.factor(insurance))) + geom_histogram(position = "identity", bins = 30, alpha = 0.6, color = "white") + scale_fill_manual(values = c("lightblue", "orange")) + theme_minimal() + labs(title = "Distribution of household expenditure by insurance", x = "Household expenditure (lnhhexp)", y = "Count", fill = "Insurance") ggplotly(p) # Both insured and uninsured households show a similar normal distribution of # expenses. Only among the highest spenders there are slightly more insured # households, suggesting that richer families are more likely to have insurance ``` # Question 3. Do pharmacy visits vary across communes (regions)? ```{r Plot-Q3} p = db %>% group_by(commune) %>% summarise(mean_exp = mean(pharvis, na.rm = TRUE)) %>% ggplot(aes(x = reorder(as.factor(commune), mean_exp), y = mean_exp, fill = "orange")) + geom_bar(stat = "identity") + theme_minimal() + labs(title = "Average pharmacy visits by commune", x = "Commune (Region ID)", y = "Average pharmacy visits") + theme(axis.text.x = element_blank(), plot.title = element_text(size = 16, face = "bold"), legend.position = "none") ggplotly(p) # In this plot it is possible to see the differences between the regions, as # pharmacy visits vary a lot between communes. Some regions have very low # averages, while others show much higher values, suggesting big regional # differences in access or health habits ``` # Question 4. Are there gender differences in illness duration? ```{r Plot-Q4} p = db %>% ggplot(aes(x = sex, y = illdays, fill = sex)) + geom_violin(trim = FALSE, alpha = 0.7, color = "gray") + scale_fill_manual(values = c("tomato", "turquoise")) + theme_minimal() + labs(title = "Distribution of Illness Days by Gender", x = "Sex", y = "Illness Days", fill = "Sex") + theme(plot.title = element_text(size = 16, face = "bold")) ggplotly(p) # From this plot, we can see that both women and men have a very similar # distribution of illness days. This suggests that there are no clear gender # differences in illness duration or access to healthcare ``` # Question 5. How does health change with age and gender? ```{r Plot-Q5} p = db %>% group_by(age_group, sex) %>% summarise(mean_health = mean(health_index)) %>% ggplot(aes(x = age_group, y = mean_health, color = sex, group = sex)) + geom_line(size = 1.2) + geom_point(size = 2) + scale_color_manual(values = c("tomato", "turquoise")) + theme_minimal() + labs(title = "Health Index by Age Group and Gender", x = "Age Group", y = "Average Health Index", color = "Sex") + theme(plot.title = element_text(size = 16, face = "bold")) ggplotly(p) # Following the previous plot, both genders show a similar distribution. The # difference is small until old age, when women show better health, which may # explain their longer life expectancy. However, this pattern, where health # decreases with age for both men and women, with women having better health # in older ages, is common worldwide, so nothing unusual appears in the data ``` # Question 6. Does higher income lead to higher pharmacy visits? ```{r Plot-Q6} p = ggplot(db, aes(x = hhexp, y = pharvis, color = as.factor(insurance))) + geom_smooth(method = "loess", se = TRUE, size = 1.2) + scale_color_manual(values = c("red3", "steelblue"), name = "Insurance", labels = c("No", "Yes")) + labs(title = "Trend between household expenditure and pharmacy visits", x = "Household expenditure", y = "Pharmacy visits") + theme_minimal() ggplotly(p) # With this plot you can see that pharmacy visits don’t change much with # household spending. he only thing we can really conclude is that people without # insurance go a bit more, especially when expenses are high ``` # Question 7. How do health, wealth, and illness interact across regions? ```{r Plot-Q7} p = db %>% group_by(wealth_group) %>% ggplot(aes(x = lnhhexp, y = health_index, size = illdays, color = as.factor(wealth_group))) + geom_point(alpha = 0.7) + theme_minimal() + labs(title = "Health and Wealth Relationship by Income Group", x = "Household Expenditure (lnhhexp)", y = "Health Index", size = "Illness Days", color = "Wealth group") ggplotly(p) # From this plot we can conclude that diseases and illness impact very similar # to all regions. So there is not any difference on how people is exposed to # being sick or ill, but as we saw in plot number 3, rich regions handle better the # the situation better after getting ill ``` # Question 8. How do insurance and education influence health and spending? ```{r Plot-Q8} p = db %>% ggplot(aes(x = education_level, y = pharvis)) + geom_boxplot(fill = "turquoise", alpha = 0.7, outlier.color = "red") + facet_wrap(~ insurance, labeller = as_labeller(c(`0` = "No insurance", `1` = "With insurance"))) + theme_minimal() + labs(title = "Pharmacy visits by education level and insurance status", x = "Education level", y = "Pharmacy visits") + theme(plot.title = element_text(size = 16, face = "bold")) ggplotly(p) # With this plot we can see that people with insurance go less to the pharmacy # in all education levels. Education does not change much, but having insurance # helps to go less ``` ## \# ## \# 3. RESAMPLING AND REGRESSION MODELS \# # Firstly, we will select the most important variables, and discard the ones # which contain irrelevant information ```{r Model-Data-Selection} data_model <- db %>% select(lnhhexp, age, educ, married, sex, illness, illdays, pharvis, insurance) # Check the structure to visualize the selected variables str(data_model) ``` ## Train-test ```{r Train-Test} set.seed(123) n <- nrow(data_model) index <- sample(1:n, size = 0.8*n) train <- data_model[index, ] test <- data_model[-index, ] # Simple linear model lm_fit <- lm(lnhhexp ~ age + educ + married + illness + illdays + pharvis + insurance + sex, data = train) summary(lm_fit) # Predictions pred <- predict(lm_fit, newdata = test) # Computing the mean squared error (MSE) to evaluate the performance and # accuracy of the predictions mse <- mean((test$lnhhexp - pred)^2) rmse <- sqrt(mse) rmse ``` ## Cross Validation (K-Fold), in this case we will declare K = 10 ```{r Cross-Validation} glm_fit <- glm(lnhhexp ~ age + educ + married + illness + illdays + pharvis + insurance + sex, data = train, family = gaussian()) cv_error <- cv.glm(data = train, glmfit = glm_fit, K = 10) cv_error$delta # Now we will create two models, one will be simpler than the other. This will # help us to observe the trade-off and balance between bias and variance # (under-fitting/over-fitting) glm_simple <- glm(lnhhexp ~ educ + insurance, data = train, family = gaussian()) glm_full <- glm(lnhhexp ~ ., data = train, family = gaussian()) cv_simple <- cv.glm(train, glm_simple, K = 10)$delta[1] cv_full <- cv.glm(train, glm_full, K = 10)$delta[1] cv_simple cv_full ``` ## Visualization of the bias/variance trade-off ```{r Bias-Variance-Plot} errors <- data.frame( Model = c("Simple", "Complex"), CV_Error = c(cv_simple, cv_full) ) ggplot(errors, aes(x = Model, y = CV_Error, fill = Model)) + geom_bar(stat = "identity", width = 0.5) + labs(title = "Comparison between the error of models with different complexity", y = "Mean Squared Error (MSE)", x = "") # We can observe through the plot that the simpler model produces a greater # error than the more complex one. However, the variance value may be increased # if many irrelevant variables are added. This would break with the balance # between bias and variance. ``` #---------------------------------------------------------------------------- \# 4. PRINCIPAL COMPONENT ANALYSIS (PCA) #---------------------------------------------------------------------------- # First of all, we will select only numeric variables, as those are the only # ones which are useful for PCA, as it is not compatible with categorical # discrete values. ```{r PCA-Data-Selection} data_num <- db %>% select(where(is.numeric)) # Check the columns to visualize the selected variables colnames(data_num) ``` ## Principal Component Analysis # Firstly, we will compute the principal components, which are linear # combinations of the original values ```{r Compute-PCA} pca <- prcomp(data_num, center = TRUE, scale. = TRUE) summary(pca) # Check the variance loading of each principal component and the total variance summary(pca)$importance[2, ] summary(pca)$importance[3, ] round(summary(pca)$importance[2:3, ], 3) # We can observe that the last principal components contain a minuscule # percentage of the variance, which means that they provide very little # new information, and therefore they can be omitted without losing big # amounts of data. Although principal components 4 and 5 contain a good # amount of information, we will stick to the first three components, which # contain almost a 60% of the variance. # Now, we will check which are the most influential variables on each of the # first (and therefore most important), principal components sort(abs(pca$rotation[,1]), decreasing = TRUE) sort(abs(pca$rotation[,2]), decreasing = TRUE) sort(abs(pca$rotation[,3]), decreasing = TRUE) # We can observe that the first two components are more related to health, # while the third one is related to the socioeconomic context of the person. # Also, PC1 seems more related to health and PC2 more related to age. # Now, we will visualize how the variance decreases with each principal # component. plot(pca, type = "l", main = "PCA Scree Plot") summary(pca)$importance pca$rotation ``` ## Visualization # Before visualizing the information with the help of plots, we need to convert # the results of the pca into a new dataset. ```{r PCA-Visualization-Setup} pca_df <- as.data.frame(pca$x) # Classic biplot, simple but not too visually attractive or appealing biplot(pca, scale = 0) ``` # More elaborated versions with ggplot2, easier to visualize the information. # We will add some extra information to the previous plot, which will also # make the graphic more attractive by adding colors. We will divide the data # shown in the previous plot in terms of sex. ```{r PCA-Visualization} pca_df$sex <- db$sex # PC1 and PC2 ggplot(pca_df, aes(x = PC1, y = PC2, color = sex)) + geom_point(alpha = 0.5) + labs(title = "PCA colored by sex") + theme_minimal() # There is a slight ascending inclination, which makes sense as older people # tend to get ill easier than younger people. # PC1 and PC3 ggplot(pca_df, aes(x = PC1, y = PC3, color = sex)) + geom_point(alpha = 0.5) + labs(title = "PCA colored by sex") + theme_minimal() # We can observe that the individuals with lower socioeconomic level do not # present data for health records. This is shown in the left side of the plot, # where the data is almost non-existent. This may seem weird at first, but it # is extremely logical because they do not have access to healthcare. This # means that the less fortunate do get ill but cannot attend hospitals, so the # data does not exist. # Also, we can see that the upper half of the graphic presents more data than # the lower one, it is denser. This also makes sense, as people with a great # social and economic situation regularly attend hospitals and healthcare, # while people from a lower social class find more difficulties. # PC2 and PC3 ggplot(pca_df, aes(x = PC2, y = PC3, color = sex)) + geom_point(alpha = 0.5) + labs(title = "PCA colored by sex") + theme_minimal() # There is not a lot of information to extract out of a dense, uniform, round # cloud of data. As well as in the first plot, the inclination is slightly # ascending. This makes sense as older people tend to be richer, as they have # had more time to earn money than the youth. # After taking a look at the three graphics, we can get some more information. # We can observe that there does not exist a clear distinction between male # and female, what suggests that the sex does not affect the variability of the # dataset. However, it is clear that when looking at the two first plots, # which are PC1/PC2 and PC1/PC3, the right half of the data is more crowded # with females than the left half, which is more crowded with males. This is # perceived because of the majority of pink color on the right side and the # majority of blue color on the left one. This probably means that the # Vietnamese society contains many old rich women and many young poor men. ``` ## \# ## \# 5. FACTOR ANALYSIS (FA) \# # First of all, we will do the exact same thing as in the Principal Component # Analysis. We will select only numeric variables, as those are the only # ones which are useful for FA, as it is not compatible with categorical # discrete values. ```{r FA-Data-Selection} data_num <- db %>% select(where(is.numeric)) # Check the structure to visualize the selected variables str(data_num) # Now, we will compute the correlation matrix, which will allow us to visualize # which variables are related to each other cor_matrix <- cor(data_num) round(cor_matrix, 2) # We need to compute the eigenvalues and use the Kaiser rule to evaluate which # factors are important and which ones can be omitted. Those that are over the # red line (>1) cannot be omitted, while the rest can be considered irrelevant. eigenvalues <- eigen(cor_matrix)$values # Eigenvalues visualization plot(eigenvalues, type = "b", main = "Eigenvalues for FA", xlab = "Factor", ylab = "Eigenvalue") abline(h = 1, col = "red", lty = 2) ``` ## Factor analysis # We will use just three factors, because of the eigenvalues computed before ```{r Factor-Analysis} fa <- fa(data_num, nfactors = 3, rotate = "varimax", fm = "ml", scores = "regression") print(fa, digits = 2, cutoff = 0.3) # Loadings and Visualization loadings <- as.data.frame(unclass(fa$loadings)) names(loadings) ggplot(loadings, aes(x = ML1, y = ML2, label = rownames(loadings))) + geom_point(color = "darkblue") + geom_text(vjust = -0.5) + labs(title = "Loadings plot (Factors 1 and 2)", x = "Factor 1", y = "Factor 2") + theme_minimal() ggplot(loadings, aes(x = ML1, y = ML3, label = rownames(loadings))) + geom_point(color = "darkgreen") + geom_text(vjust = -0.5) + labs(title = "Loadings plot (Factors 1 and 3)", x = "Factor 1", y = "Factor 3") + theme_minimal() ggplot(loadings, aes(x = ML2, y = ML3, label = rownames(loadings))) + geom_point(color = "magenta") + geom_text(vjust = -0.5) + labs(title = "Loadings plot (Factors 2 and 3)", x = "Factor 2", y = "Factor 3") + theme_minimal() # After observing the grouping of the variables, we can extract the following # conclusions out of the plots: # Factor 1: Social and economic context of the individual (variables such as # lnhhexp, insurance and education). A high score in this factor implies a # great social and economic condition (e.g.: great insurance and education) # Factor 2: Health situation of the individual (variables such as # illness, illdays and pharvis). A high score in this factor implies a # great health condition (e.g.: no recent illnesses) # Factor 3: General information about the individual (variables such as age # and sex). A high score in this factor probably implies that the person is old. ``` ## Scores per individual # Now, we will add to the dataset the importance each factor has in each person ```{r FA-Scores-Setup} fa_scores <- as.data.frame(fa$scores) head(fa_scores) data_fa <- cbind(db, fa_scores) ``` ## Visualization of the distribution # As well as in the PCA, we colored the plots by sex so that we can obtain # some extra information while making the graphic more visually appealing. ```{r FA-Scores-Visualization} # Factors 1 and 2 ggplot(data_fa, aes(x = ML1, y = ML2, color = sex)) + geom_point(alpha = 0.5) + labs(title = "Individuals distribution by factors (1 and 2)") + theme_minimal() # We can observe that the individuals with lower socioeconomic level do not # present data for health records. This may seem weird at first, but it is # extremely logical because they do not have access to healthcare. This means # that the less fortunate do get ill but cannot attend hospitals, so the # data does not exist. # Also, we can see that the right half of the graphic presents more data than # the left one, it is denser. This also makes sense, as people with a great # social and economic situation regularly attend hospitals and healthcare, # while people from a lower social class find more difficulties. # We can obtain some more information. At the end of the right side of the # graphic, the health condition presents a decreasing inclination. This is # probably because the richest people are usually older people, as they have # had more time to earn money, an entire lifetime. This is why those with a # better socioeconomic context (probably elders) have worse health. # Factors 1 and 3 ggplot(data_fa, aes(x = ML1, y = ML3, color = sex)) + geom_point(alpha = 0.5) + labs(title = "Individuals distribution by factors (1 and 3)") + theme_minimal() # This graphic does not provide a lot of information. The result, a round # cloud of uniform data, suggests that there are not many socioeconomic # differences in terms of age. We could say that the society in Vietnam is # equally distributed in terms of socioeconomic context. There is young rich # people, old rich people, young poor people and old poor people. # Factors 2 and 3 ggplot(data_fa, aes(x = ML2, y = ML3, color = sex)) + geom_point(alpha = 0.5) + labs(title = "Individuals distribution by factors (2 and 3)") + theme_minimal() # There are two peculiar nuances in this graphic. Firstly, we can notice the # division between the people with the worst health situation and those who # have a slightly better condition, there exists a blank space that breaks # with the tendency of the rest of the plot. This may be because those who # present a serious illness spend way more time at hospital than those who # have a less serious one. # Secondly, we can notice the inclination of the graphic, which looks like the # shape of an arrow pointing to the right side. This shape makes perfect sense, # as both the elders and the kids tend to get ill regularly, while adults are # less likely to get ill. Also, when a kid or an elder gets ill, the hospital # is usually the first destination, while many adults avoid going to the # hospital unless it is a serious illness or injury. # Out of the third graphics relating the different factors, shown above, # we can also observe that, as well as in the Principal Component Analysis, # there does not exist a clear distinction between male and female, what # suggests that the sex does not affect the first two factors. In all of the # plots, the males and females are scrambled throughout the whole plot. # In conclusion, all of the graphics confirm the conclusions made with PCA. # ---------------------------------------------------------------------------- # 6. CLUSTERING TOOLS # ---------------------------------------------------------------------------- ### 1. Global clustering to see general profiles using demographic and health # variables. Prepare the data ## 1.1 Select the variables X = select(db, age_real, illdays, actdays, pharvis, hhexp, education_level, insurance, sex) rownames(X) = rownames(db) X$education_level = as.numeric(factor(db$education_level, levels = c("Low", "Medium", "High"))) X$sex <- as.numeric(factor(db$sex, levels = c("male", "female"))) # 1 = men, 2 = female X$insurance <- as.numeric(factor(db$insurance, levels = c("0", "1"))) # 1 = no insurance, 2 = with insurance ``` ## 1.2 We escalate the data ```{r Scale-Data} X = scale(X) str(X) #------------------------------------------------------------------------------- # IMPORTANT MESSAGE # From now on, we will use a 5,000-observation sample of the dataset to speed up # computations, as the full dataset made some methods take unreasonably long to execute # Already done in the gap statistic method of the PAM's selection number of clusters X_sample <- X[sample(1:nrow(X), 5000), ] #------------------------------------------------------------------------------- ``` ## 1.3 Check structure and NAs ```{r Check-NAs} dim(X) summary(X) sum(is.na(X)) ``` ## 1.4 Look for correlations with GGally ```{r GGally-Correlation} GGally::ggcorr(X, label=T, hjust=0.75) ``` ------------------------------------------------------------------------ ### 2. K-means ## 2.1 Determine the number of cluster to use. Elbow method and silhouette method # 2.1.1 USE fviz_nbclust to choose the number of clusters with the elbow method ```{r Kmeans-WSS} fviz_nbclust(X_sample, kmeans, method = 'wss', k.max = 5, nstart = 8) ``` # 3.1.2 USE fviz_nbclust to choose the number of clusters with the silhouette method ```{r Kmeans-Silhouette} fviz_nbclust(X_sample, kmeans, method = 'silhouette', k.max = 5, nstart = 8) ``` # 2.1.3 USE fviz_nbclust to choose the number of clusters with the gap method ```{r Kmeans-Gap} fviz_nbclust(scale(X_sample), kmeans, method = "gap_stat", k.max = 5, nstart = 8, nboot = 20) ``` # 2.1.4 Conclusion # With the elbow method we were unsure between 3 and 4 clusters, but the silhouette # method clearly indicated that 3 was the optimal choice.The gap statistic suggested # four clusters, but the increase over 3 was small., so we kept 3 for simplicity and interpretability # Now we will run the k-means function with the information we have gathered ```{r Kmeans-Run} set.seed(123) k <- 3 fit <- kmeans(X_sample, centers = k, nstart = 50) groups = fit$cluster barplot(table(groups), col="darkblue") ``` ## 2.2 Interpretation of the centers # Convert cluster information to useful data to ggplot ```{r Kmeans-Centers-Setup} centers = fit$centers tidy = cbind( gather(as.data.frame(t(centers)), "cluster", "coor"), var = rep(colnames(centers), each = nrow(centers)), size = rep(table(fit$cluster), each = ncol(centers)) ) ``` # 1º Graph - How the variables change between cluster ```{r Kmeans-Centers-Plot1} tidy %>% ggplot(aes(x = cluster, y = coor, fill = cluster)) + geom_col() + facet_wrap(~var) + geom_text(aes(label = size), position = position_stack(1.2)) ``` # 2º Graph - How is each clusters in respect to all variables ```{r Kmeans-Centers-Plot2} tidy %>% ggplot(aes(x = var, y = coor, fill = var)) + geom_col() + facet_wrap(~cluster) + theme(axis.text.x = element_text(angle = 45)) # General conclusion # 1 - Represents the average population, with younger individuals and slightly # lower education levels. This group shows overall moderate values across all variables # 2 - Characterized by high pharmaceutical spending (pharvis) and higher # household expenditure. This cluster may include individuals with chronic or # costly illnesses, possibly belonging to wealthier households # 3 - Composed mainly of individuals sharing similar gender characteristics, with # higher insurance coverage and better education levels. This group represents a # more protected and socioeconomically advantaged segment of the population ``` ## 2.3 Clusplot # Use fviz_cluster to visualize the groups ```{r Kmeans-Clusplot} fviz_cluster(fit, data=X_sample, geom=c("point"), ellipse.type = "norm") + theme_minimal() # See what reflects each PCA pca = prcomp(X, scale = TRUE) ``` # PC1 reflects health and wellbeing, contrasting healthy, educated individuals # with those showing illness and higher medical expenses ```{r Kmeans-PC1-Contribution} barplot(pca$rotation[, 1], las = 2, col = "darkblue", main = "Contribution of Variables to PC1") ``` # PC2 represents a socioeconomic dimension, separating individuals by income, # education, and insurance coverage ```{r Kmeans-PC2-Contribution} barplot(pca$rotation[, 2], las = 2, col = "darkblue", main = "Contribution of Variables to PC2") ``` # See if the PCA correlation interpretations correct or not ```{r Kmeans-PCA-Correlation} cor(pca$x[,1], db$health_index) cor(pca$x[,1], db$age_real) cor(pca$x[,1], as.numeric(db$pharvis)) cor(pca$x[,1], as.numeric(db$illdays)) cor(pca$x[,2], db$hhexp) cor(pca$x[,2], as.numeric(db$insurance)) cor(pca$x[,2], as.numeric(db$education_level)) ``` ## 2.4 Now we will use eclust from factorextra to have more clustering tools # Run K-means using eclust (same result, more analysis tools) ```{r Kmeans-Eclust} fit.kmeans = eclust(X_sample, "kmeans", k = 3, stand = TRUE, nstart = 100) ``` ## 2.5 Create the Silhouette Plot to see the quality of the clusters. We will use # directly factorextra ```{r Kmeans-Silhouette-Plot} fviz_silhouette(fit.kmeans) # We see that the first cluster is well defined. The third cluster is moderately # clear. However the second cluster is poorly separated. To solve this we will # use other clustering tools, trying to reach a better approach ``` ------------------------------------------------------------------------ ### 3. PAM ## 3.1. Determine the optimal number of clusters for PAM # 3.1.1 With Within Sum of Squares (WSS) ```{r PAM-WSS} fviz_nbclust(X_sample, pam, method = "wss", k.max = 5) ``` # 3.1.2 With silhouette coefficient ```{r PAM-Silhouette} fviz_nbclust(X_sample, pam, method = "silhouette", k.max = 5) ``` # 3.1.3 With gap statistic ```{r PAM-Gap} fviz_nbclust(X_sample, pam, method = "gap_stat", k.max = 5, nboot = 10) ``` # 3.1.4 Conclusion # The elbow, silhouette and gap methods suggest an optimal range between 3 and 5 # clusters. Since the silhouette and the gap statistics continues to improve slightly at 5 but the # elbow stabilizes around 4, we selected 4 clusters as the most balanced and # interpretable option # We run PAM with eclust using the information we got ```{r PAM-Run} fit.pam = eclust(X_sample, FUNcluster = "pam", stand = TRUE, k = 4, graph = FALSE, nstart = 8) # Visualize clusters fviz_cluster(fit.pam, geom = "point", main = "Clusters with PAM", ellipse.type = "norm", star.plot = TRUE) + theme_minimal() # The PAM clusters overlap considerably, suggesting weak separation between groups # or high similarity among observations across dimensions. We are going to use # other clustering tools trying to achieve a clearer structure ``` ## 3.2 Compare the k-means and the PAM results. # I will repeat fit.kmeans but with X_sample, so both have the same length ```{r PAM-Kmeans-Comparison} fit.kmeans = eclust(X_sample, "kmeans", k = 3, stand = TRUE, nstart = 100) fit.pam = eclust(X_sample, FUNcluster = "pam", stand = TRUE, k = 4, graph = FALSE, nstart = 8) adjustedRandIndex(fit.kmeans$cluster, fit.pam$clustering) # This comparison showed that k-means and PAM clusters are quietly different, following # the conclusion that PAM clusters are more overlapped ``` ------------------------------------------------------------------------ ### 4. Kernel k-means # We will use the same clusters as in k-means ```{r Kernel-Kmeans-Run} fit.ker = kkmeans(as.matrix(X_sample), centers = 3, kernel = "rbfdot") ``` ## 4.1 Inspect the model results # Get the cluster centers (in transformed feature space) ```{r Kernel-Kmeans-Inspect} centers(fit.ker) # Get the within-cluster sum of squares (cluster compactness) withinss(fit.ker) # Get the number of elements in each cluster size(fit.ker) ``` ## 4.2 Visualize the clusters ```{r Kernel-Kmeans-Visualization} object.ker = list(data = X_sample, cluster = fit.ker@.Data) fviz_cluster(object.ker, geom = "point", ellipse = T) + geom_text(aes(label = 1:nrow(X_sample)), check_overlap = T, hjust = 0, vjust = 0) + theme_minimal() + labs(title = "Kernel K-means Clustering (RBF Kernel)") ``` ## 4.3 Conclusion ```{r Kernel-Kmeans-Conclusion} # It looks like the kernel k-means tool also failed in improving the clusters # that k-means already gave. The second cluster is almost totally nested in the # first one and the third cluster overlaps both ``` ------------------------------------------------------------------------ ### 5. Hierarchical clustering ```{r Hierarchical-Clustering-Run} distances = dist(X_sample, method = "euclidean") hc = hclust(distances, method = "ward.D2") ``` ## 5.1 Visualize the cluster's structure # 5.1.1 Dendrogram ```{r Hierarchical-Dendrogram} fviz_dend(x = hc, k = 3, rect = TRUE, rect_fill = TRUE, palette = "Set1", rect_border = "Set1", main = "Hierarchical Clustering - Dendrogram") # The dendrogram visualization allowed us to identify three well-separated groups, # making the use of a phylogenic representation unnecessary in this case. The # clusters are clearly distinguishable, confirming the adequacy of the hierarchical # approach. Making the group profiles much clearer to interpret. ``` ## 5.2 Obtaining the clusters and information about them # Cut the dendrogram into 3 clusters ```{r Hierarchical-Cluster-Info} hc.cluster <- cutree(hc, k = 3) # See how many observations belong to each cluster and what each cluster is made of table(hc.cluster) aggregate(X_sample, by = list(cluster = hc.cluster), mean) # Cluster 1. Older, lower-education, uninsured # Cluster 2. Young and healthy # Cluster 3. Educated, insured, and wealthier ``` ## 5.3 Compare with K-means Results ```{r Hierarchical-Kmeans-Comparison} adjustedRandIndex(fit.kmeans$cluster, hc.cluster) ``` ## 5.4 Conclusion ```{r Hierarchical-Conclusion} # The hierarchical clustering (Ward’s method) produced three clearly separated # groups. # The ARI with K-means suggests partial similarity, but, to us, the hierarchical # model reveals a much clearer and more interpretable cluster structure ``` ------------------------------------------------------------------------ ### 6. Expectation Maximization ```{r EM-Run} res.Mclust = Mclust(X_sample) summary(res.Mclust) # The algorithm selected an EEV model (ellipsoidal, equal volume and shape, # variable orientation). BIC and ICL values are close, indicating moderately # well-defined clusters. Six clusters were identified, most of them balanced # in size, while one very small cluster likely represents outliers. ``` ## 6.2 Inspect the model result # 6.2.1 The probability of each observation to belong to each cluster ```{r EM-Probabilities} head(res.Mclust$z) ``` # 6.2.2 The assigned cluster to each observation based on the highest probability ```{r EM-Classification} head(res.Mclust$classification) ``` ## 6.3 Visualizing the selection process with the BIC model selection criteria ```{r EM-BIC-Plot} fviz_mclust(res.Mclust, what = "BIC") ``` ## 6.4 Visualizing the clusters ```{r EM-Classification-Plot} fviz_mclust(res.Mclust, what = "classification") + theme_minimal() + labs(title = "EM Clustering - Gaussian Mixture Model") # The EM clustering produced six groups, but most of them are highly overlapped, # indicating weak separation between clusters ``` ## 6.5 Comparison with previous models ```{r EM-Comparisons} adjustedRandIndex(res.Mclust$classification, fit.pam$clustering) adjustedRandIndex(res.Mclust$classification, hc.cluster) adjustedRandIndex(res.Mclust$classification, fit.kmeans$cluster) adjustedRandIndex(res.Mclust$classification, fit.ker) # This model have similar structure to the PAM model, a moderate one with the # Hierarchical tool. And low similarities with normal and kernel k-means. # This follows the conclusions that we got from the cluster visualization, # because both, PAM and EM algorithms, show a lot of overlapping between # the clusters ``` ------------------------------------------------------------------------ ### 7. Final conclusion ```{r Final-Conclusion} # Across all methods, K-means and Hierarchical clustering provided the most # interpretable and well-separated group structures, both revealing three distinct # profiles within the data. PAM, Kernel K-means and EM showed excessive overlap. # Overall, K-means and Hierarchical approaches offered the best balance between # simplicity, clarity and interpretability. ```