if(!require("pacman")) install.packages("pacman")
## Cargando paquete requerido: pacman
pacman::p_load(tidyverse,
               tidymodels,
               janitor,
               purrr,
               stringr,
               highlight,
               corrplot,
               RColorBrewer,
               tidytext,
               FactoMineR, caret, factoextra,
               fpc)

1 Diccionario

1.0.0.0.3 Base de datos
hospital <- read.csv("HCAHPS-Hospital.csv") %>% 
  janitor::clean_names()
1.0.0.0.4 Grafico

Cómo estan evaluadas las encuestas, cuáles reciben mayor calificación.

hospital %>%
  mutate(patient_survey_star_rating = as.numeric(gsub("[^0-9.]", "", patient_survey_star_rating))) %>%
  drop_na() %>% 
  #agrupar y calcular la media
  group_by(hcahps_answer_description) %>%
  summarise(avg_rating = mean(patient_survey_star_rating, na.rm = TRUE)) %>%
  #gráfico de frecuencia
  ggplot(aes(x = reorder(hcahps_answer_description, avg_rating),
             y = avg_rating)) +
  geom_col(fill = "#BBAE5F", color = "black") +
  geom_text(aes(label = round(avg_rating, 2)), hjust = -0.1, size = 3) +
  labs(title = "Calificaciones promedio por tema",
       x = NULL,
       y = "Calificación promedio",
       caption = "HCAHPS - 2024",
       tag = "Fig. 1") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1)) +
  coord_flip() +
  #tema
  theme_minimal() +
  #expandir el límite y para que quepan las etiquetas
  scale_y_continuous(expand = expansion(mult = c(0, 0.1)))

1.0.0.0.4.1 Objetivo
  • Obtener un data frame que indetifique por columna las variables del cuestionario.
  • Frame con columnas de interes
  • Separar variables en distintos frames
  • Limpiar cada frame convirtiendo a númerico
  • Unir cada frame mediante un join
  • Verificar si se puede unir filar y solo pivotear
  • Obtener la variable summary rating como variable dependiente o explicada
  • Realizar un análisis exploratorio de datos
  • Reducir la dimensionalidad de las variables
  • Realizar una clasificación

Nota: post-edición

En un principio se espera separar las variables de calificación final contra las evaluaciones por encuesta, pero es un asco la base de datos, o al menos no pense en un método para limpiarla mejor, así que solo evaluare las calificaciones finales por encuesta vs la calificación final del usuario. Triste, así que si alguien que esta por acá leyendo, se agradece los comentarios y las sugerencias, dejo los links para que puedan consultar la información y replicar el código, gg.

1.0.0.1 2. Variables de interes

Se han identificado las principales variables para trabajar con la base:

  • Identificadores:
  • ID, city, name
  • Variables de interes
  • name, hcahps_question, star rating
  • Variables “predeterminadas”
  • hcahps_question, star rating.

Se debe contemplar que dejar star rating podría generar multicolinealidad ya que es la medida “agrupada” de la evaluación para cada tipo de cuestionario. Por lo que se puede agrupar en una base diferente.

1.0.0.1.1 2.1 Calificaciones por pregunta, sin evaluar rating
reduce.data <- hospital %>%
  select(facility_id, city_town, facility_name, hcahps_measure_id, hcahps_answer_percent)

1.1 Función

#library(purrr)
#library(stringr)

#Creamos función que:

 # Convierte el input a carácter.
 # Elimina todos los caracteres excepto números, puntos y guiones.
 # Reemplaza cadenas vacías con "0".
 # Convierte el resultado a numérico.

clean_numeric <- function(x) {
  x <- as.character(x) #cinvertir a caracter
  x <- str_replace_all(x, "[^0-9.]", "")  # Elimina todo excepto números, puntos y guiones
  x <- str_replace_all(x, "^\\s*$", "0")   # Reemplaza cadenas vacías con "0"
  as.numeric(x)
}
1.1.0.0.1 2.1.2: quitamos ratings y mean
reduce.data <- reduce.data %>%
  #Función a columnas numericas
  mutate(across(
    #Seleccionar columnas
    .cols = !c(facility_id, city_town, facility_name,hcahps_measure_id),
    #Aplicar función a columnas seleccionadas
    .fns = ~clean_numeric(.)
  )) %>%
  filter(hcahps_answer_percent != 0)

Filas con califiación en ratings

Debemos quitar el linear mean score y tomar de footnote las calificaciones que no estan en rating.

El mal diseño de la base hace que haya filas intercaladas con datos en diferentes columnas como se puede ver arriba, es por eso que no se puede hacer el pivot directamente, o al menos tome la primera estrategia para hacerlo.

1.1.0.0.2 2.1.2: Solo ratings y mean
#Base con solo calificaciones finales = "star rating"
data_ratings <- hospital %>%
  select(facility_id, city_town, facility_name, hcahps_measure_id, patient_survey_star_rating) %>%
  filter(str_detect(tolower(hcahps_measure_id), "star_rating")) %>% 
  clean_names()

data_ratings <- data_ratings %>%
  #Función a columnas numericas
  mutate(across(
    #Seleccionar columnas
    .cols = !c(facility_id, city_town, facility_name,hcahps_measure_id),
    #Aplicar función a columnas seleccionadas
    .fns = ~clean_numeric(.)
  ))%>%
  filter(patient_survey_star_rating != 0)
1.1.0.0.3 Paso Tres: Pivoteamos
#Base reducida con evaluación por pregunta
reduce.data <- reduce.data %>%
  pivot_wider(names_from = hcahps_measure_id,
              values_from = hcahps_answer_percent) %>% 
  janitor::clean_names()

#Base reducida solo con evaluación final
data_ratings <- data_ratings %>%
  pivot_wider(names_from = hcahps_measure_id,
              values_from = patient_survey_star_rating) %>% 
  janitor::clean_names()
1.1.0.0.4 Paso cuatro: Unimos
#Solo necesitamos columna summay de la base ratings.
dta.ratings2 <- data_ratings %>%
  select("h_star_rating", "facility_id", "city_town", "facility_name") 

# Unión de reduce.data y data_ratings
data_merge <- merge(reduce.data, 
                    dta.ratings2,
                    by = c("facility_id", "city_town", "facility_name"))

data_merge <- rename(data_merge,
       "summary_star_rating" = h_star_rating)

data_ratings <- rename(data_ratings,
       "summary_star_rating" = h_star_rating)

1.1.1 Base

Tenemos dos bases:

  1. star_ratings, donde encontramos la calificación fincal para cada tipo de cuestionario, con la variable summary_star_rating que toma el total de las evalauaciones aplicadas para cada usuario.

  2. data_merge, en esta se encuentran la evaluación por pregunta de cada tipo de encuesta, junto con la variable summary_star_rating.

En ambas bases también se ubican las variables de indentificación, como id, town, y name que en este caso, name toma el nombre del centro de atención.

1.1.2 Visualización

Caracteristicas Generales

Centros de atención evaluados: 3255 Usuarios evaluados: 3186 Ciudades registradas: 1980

1.1.2.0.1 Principales calificaciones por número de centros de atención
ggplot(data_ratings, aes(x = summary_star_rating)) +
  geom_histogram(binwidth = 0.5, fill = "skyblue", color = "black") +
  labs(title = "Distribución de Calificaciones",
       x = "Calificación",
       y = "Frecuencia") +
  theme_minimal()

Se observa que la calificación promedio en el total de centros de atención se encuentra en 3 sobre la escala de 1 a 5.

# Supongamos que tus datos están en un dataframe llamado 'data_ratings'

# Seleccionar cinco centros de atención al azar
set.seed(123)  # Para reproducibilidad
selected_centers <- sample(unique(data_ratings$facility_name), 5)

# Filtrar los datos para estos centros seleccionados
filtered_data <- data_ratings[data_ratings$facility_name %in% selected_centers, ]

# Crear el histograma para los centros seleccionados

ggplot(filtered_data, aes(y = summary_star_rating, x=facility_name, color = facility_name)) +
  geom_point() +
  #titulos
  labs(title = "Distribución de Calificaciones en Centros Seleccionados",
       x = "Centro de Atención",
       y = "Calificación") +
  #tema del gráfico
  theme_minimal() +
  #paleta de colores
  scale_color_brewer(palette = "Set3") +
  theme(
    axis.text.x = element_blank(),          # Eliminar etiquetas del eje x
    axis.ticks.x = element_blank(),         # Eliminar marcas del eje x
    legend.position = "bottom",             # Mover leyenda abajo
    legend.title = element_blank(),         # Quitar título de la leyenda
    plot.title = element_text(hjust = 0.5)  # Centrar el título del gráfico
  ) +
  guides(color = guide_legend(ncol = 3, nrow = 2))  # Hacer la leyenda de dos filas

1.2 Modelos

  1. Se realiza un ACP con los datos a nivel de pregunta en cada encuesta, la variable dependiente será el rating final. Se aconseja retomar los resultados para ratings de 5, para evaluar cada componente sobre la mejor calificación, sin embargo, al tener una media sobre 3, se decide dejar todos los niveles de ratings.

  2. Posteriormente se raliza el ACP con la evaluación final de cada encuesta, siguiendo el mismo método que el primer análisis.

##A) Tidymodels*

Usamos el componente basico de tidymodels para realizare el análisis PCA, posteriormente usamos otras librerias. Hasta el momento desconozco como aplicar PCA de otras librerias dentro de tidymodels.

###Análisis de Componentes principales

1.2.0.0.0.1 Descriptivo

Limpieza

# Renombrando las columnas
data_ratings <- rename(data_ratings,
       "Nurse communication" = h_comp_1_star_rating,
       "Doctor communication"= h_comp_2_star_rating,
       "Staff responsiveness"= h_comp_3_star_rating,
       "Communication about medicines"= h_comp_5_star_rating,
       "Discharge information" = h_comp_6_star_rating,
       "Care transition" = h_comp_7_star_rating,
       "Cleanliness" = h_clean_star_rating,
       "Quietness" = h_quiet_star_rating,
       "Overall hospital rating" = h_hsp_rating_star_rating,
       "Recommend hospital" = h_recmnd_star_rating)

Distribución

#library(corrplot)
#library(RColorBrewer)
#Variables Númericas
col.numeric <- data_ratings %>% select(where(is.numeric))
#Matriz de Correlación
cor_matrix <- cor(col.numeric, use = "complete.obs")
#Gráfico de correlación
corrplot(
  cor_matrix, 
  tl.col = "black", # Color de las etiquetas de las variables
  col = colorRampPalette(brewer.pal(11, "PuOr"))(200)  # Paleta de colores
)

#Gráficos de distribución de las variables
col.numeric %>% 
  # Convertir la base de datos en formato largo
  pivot_longer(cols = everything(),#Mover todas las variables
               names_to = "variable",#Pasarlas a una sola
               values_to = "valor") %>% #Pasar las observaciones a una sola
  ggplot(aes(x = valor)) +
  geom_histogram() +
  facet_wrap(~ variable, scales = "free") +
  theme_bw()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

1.2.0.0.1 1. Receta

Primero, debemos indicar en recipe() como debemos realizar con nuestro conjunto de datos.

#library(tidymodels)
pca_rec <- recipe(summary_star_rating ~.,data = data_merge) %>%
  #Variable que no usare
  update_role(facility_id,city_town, facility_name, new_role = "Id") %>%
  #Seleccionas solo numericas
  step_select(where(is.numeric)) %>%
  #imputa con la media valores faltantes
  step_impute_mean(all_predictors()) %>%
  #Eliminamos variables sin varianza
  step_zv(all_predictors()) %>%
  step_nzv(all_predictors()) %>%
  #Normalizamos datos
  step_normalize(all_predictors()) %>%
  #Convertinos en Componentes Principales
  step_pca(all_predictors())
1.2.0.0.2 2. Aplicar la receta al conjunto de datos

Con prep() realizamos las instrucciones antes indicadas

pca_prep <- prep(pca_rec)
pca_prep
## 
## ── Recipe ──────────────────────────────────────────────────────────────────────
## 
## ── Inputs
## Number of variables by role
## outcome:    1
## predictor: 72
## Id:         3
## 
## ── Training information
## Training data contained 3255 data points and 128 incomplete rows.
## 
## ── Operations
## • Variables selected: h_comp_1_a_p and h_comp_1_sn_p, ... | Trained
## • Mean imputation for: h_comp_1_a_p and h_comp_1_sn_p, ... | Trained
## • Zero variance filter removed: <none> | Trained
## • Sparse, unbalanced variable filter removed: <none> | Trained
## • Centering and scaling for: h_comp_1_a_p and h_comp_1_sn_p, ... | Trained
## • PCA extraction with: h_comp_1_a_p and h_comp_1_sn_p, ... | Trained
1.2.0.0.3 3. Obtener resultados

Bake nos proporciona los resultados de la transformación de datos a través del PCA.

resultados_pca <- bake(pca_prep, new_data = NULL)
head(resultados_pca)

Interpretación:

  • El primer centro de salud, con summary_star_rating de 4 tiene valores positivos altos en PC2 y PC3, sugiriendo que se desempeña bien en los aspectos principales que estos componentes representan.
  • El segundo centro de salud con summary_star_rating 3 tiene valores negativos en todos los PC, indicando un desempeño más bajo en estos aspectos.
  • Principalmente para ratings mayores a 4, los elementos del primer PC son positivos, mientras que para todos los menores a 4 son negativos.
  • Para la mayoria de los ratigns los elementos influyen de manera significativa en el PC2.
tidy(pca_prep)

Extraer las cargas (loadings) de cada variable original en cada componente principal para los valores en la columna value que indican cuánto contribuye cada variable original a cada componente principal.

tidy(pca_prep, number = 2, type = "coef")
#Gráfico de variables en los CP
tidy(pca_prep, 6) %>% 
  #Elegir rango de PC
  filter(component %in% paste0("PC", 1:4)) %>%
  #Ordenar de PC1 a PC4
  mutate(component = fct_inorder(component)) %>%
  #Gráfico componentes por terminos (variables)
  ggplot(aes(value, terms, fill = terms)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~component, nrow = 1) +
  labs(y = NULL)

Gráfico de los primeros 10 terminos en cada componente

#library(tidytext)
tidy(pca_prep, 6) %>%
  filter(component %in% paste0("PC", 1:4)) %>%
  group_by(component) %>%
  top_n(10, abs(value)) %>%
  ungroup() %>%
  mutate(terms = reorder_within(terms, abs(value), component)) %>%
  ggplot(aes(abs(value), terms, fill = value > 0)) +
  geom_col() +
  facet_wrap(~component, scales = "free_y") +
  scale_y_reordered() +
  labs(
    x = "Absolute value of contribution",
    y = NULL, fill = "Positive?"
  )

1.3 B) FactoMiner

1.3.0.0.0.1 1. Preparar datos
#library(FactoMineR)
#library(caret)
#library(factoextra)

#Restamos media y dividimos entre varianza = Valores Z para normalizar
#train_data.standard <- preProcess(datos, method = c("center", "scale"))
#train_data.norm <- predict(train_data.standard, train_data)

#Seleccionar variables númericas
datos <- data_merge %>% 
  #slecciono columnas numericas
  select(which(sapply(data_merge,is.numeric))) %>% 
  #Quito la variable dependiente
  select(-summary_star_rating) %>% 
  #Quito Faltantes
  drop_na() %>%
  select(where(~ var(.) > 0)) %>% 
# Normalizar los datos
  scale()
1.3.0.0.0.2 2. Aplicar PCA
# Realizar el análisis PCA con FactoMineR
pca1 <- PCA(datos, graph = FALSE)

#Corremos PCA nativo
acp <- prcomp(datos)
1.3.0.0.1 2.1 Análisis Resultados

Evaluamos:

  1. CP con mayor varianza explicada
  2. Contribución al CP por variable
  3. Variables de mayor varianza explicada al CP1
1.3.0.0.1.1 i. Varianza explicada
head(pca1$eig)
##        eigenvalue percentage of variance cumulative percentage of variance
## comp 1  44.542739              61.864915                          61.86492
## comp 2   6.708206               9.316953                          71.18187
## comp 3   3.061410               4.251959                          75.43383
## comp 4   2.581420               3.585306                          79.01913
## comp 5   1.854886               2.576230                          81.59536
## comp 6   1.432225               1.989202                          83.58456

Como se puede esperar, el primer componente alberga la mayor varianza explicada (eigenvalor), siendo el 61% del total.

1.3.0.0.1.2 ii. ScreePlot: Numero de componentes a retener:
# Scree plot para determinar el número de componentes a retener
fviz_eig(pca1, addlabels = TRUE)

Respecto al Screeplot, se obtiene un codo en el segundo componente y se ilustra que el primer componente explica mayormente la varianza.

Representación de observaciones sobre componentes principales.

• fviz_pca_ind(): Representación de observaciones sobre componentes principales.

# Graficar los individuos (primeros dos componentes principales)
fviz_pca_ind(pca1,
             geom.ind = "point", # Mostrar los individuos como puntos
             col.ind = "cos2",   # Color por calidad de representación
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE)       # Evitar la sobreposición de etiquetas

#Scores después de la rotación
#Nos dice como contribuyen de forma positiva o negativa las variables
#En donde estan las variables respecto al PCA
fviz_pca_var(pca1, 
             axes = c(1, 2),
             col.var="contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07", "#008f39"),
             repel = TRUE # Avoid text overlapping
             )

La calidad de esta representación se mide por el valor al cuadrado del coseno (cos2) del ángulo del triángulo formado por el punto del origen, la observación y su proyección sobre el componente. Para una variable dada, la suma del cos2 sobre todos los componentes principales será igual a 1, y si además la variable es perfectamente representable por solo los dos primeros componentes principales, la suma de cos2 sobre estos dos será igual a 1. Variables posicionadas cerca del origen puede ser un indicativo de que serían necesarios más de dos componentes principales para su representación.

Top 10 variables que más contribuyen a PC1

#Contribucaciones del Componente Principal 1,  por Variable
fviz_contrib(pca1, choice = "var", axes = 1, top = 10)

Principales variables en el PCA1

  • Las enfermeras “siempre” se comunicaron bien
  • Las enfermeras “siempre” escucharon atentamente
  • Las enfermeras “siempre” explicaban las cosas para que pudieran entender
  • Las enfermeras “siempre” los trataron con cortesía y respeto.
  • Pacientes que “totalmente de acuerdo” en que entendieron su atención cuando salieron del hospital.
  • Pacientes que dieron una calificación de “9” o “10” (alta)
  • El personal explicó “a veces” o “nunca”.
  • Los médicos “siempre” se comunicaron bien
  • Pacientes que dieron una calificación de “6” o menos (baja)
  • Pacientes que “totalmente de acuerdo” en que el personal tuvo en cuenta sus preferencias
hospital %>% 
  select(hcahps_measure_id,hcahps_answer_description) %>% 
  head()
1.3.0.0.1.3 2. Elección de CP

Criterio valores > 1

fviz_eig(pca1, choice = "eigenvalue")

Criterio valores > 1 y % de varianza

fviz_screeplot(pca1, addlabels = TRUE)

proporción de varianza explicada (PVE) y acumulada

#Contribuciones por dimensión a los CP
evi.val <- get_eigenvalue(acp)
head(evi.val)
par(mfrow = c(1,2))

plot(evi.val$variance.percent, type = "o", 
     ylab = "PVE", 
     xlab = "Componente principal", 
     col = "blue")
plot(cumsum(evi.val$variance.percent), type = "o", 
     ylab = "PVE acumulada", 
     xlab = "Componente principal", 
     col = "brown3")

De manera conjunta, los primeros 5 componentes principales explican en torno al 80% de la varianza de los datos, lo cual es una buena cantidad.

1.4 C) Cluster

#Análisis de Clusters
#library(fpc)
set.seed(123)
#Intentamos con dos grupos, dada la conclusión obtenida en PCA.
k1 <- kmeansruns(datos, #datos centrados
                 k = 5,#2 CP
                 runs = 100)#Al menos 100 iteraciones

#Visualizamos
fviz_cluster(k1, data = datos)

También podemos trabajar con otra función

# Calcule k-medias
set.seed(123)
km.2 <- kmeans(datos,#Base normalizada
               5,#número de cp
               nstart = 25)

#Visualizamos
fviz_cluster(
   km.2,
   data = datos,
   ellipse.type = "euclid",
   # Concentration ellipse
   star.plot = TRUE,
   # Avoid label overplotting (slow)
   ggtheme = theme_minimal()
)

# #Agregar los clusters a la base de datos inicial: datos
# aggregate(datos, #base normalizada
#           by=list(cluster=km.2$cluster),#cluster
#           mean) #media de las variables
1.4.0.0.0.1 Clusters óptimos
#Evaluamos el número de cluster óptimo con la base datos, es es la que se trabaja y se encuentra estandarizada.
set.seed(123)

#Método de la suma de cuadrados dentro de los clusters
fviz_nbclust(datos,
             kmeans, #Con KMedias
             method = "wss")+#método de separación de medias al cuadrado
  labs(subtitle = "Elbow method")#Gráfico de codo

#Método de interpretación y validación de la coherencia dentro del análisis de grupos
fviz_nbclust(datos,
             kmeans, #Con KMedias
             method = "silhouette")+#meétodo de separación de medias al cuadrado
  labs(subtitle = "Silhouette method")

Ambos métodos arrojan un K óptimo de 2.

fviz_nbclust(datos, kmeans, nstart = 25, method = "gap_stat", nboot = 50)+
   labs(subtitle = "Gap Statistic Method")

1.4.0.0.0.2 Estabilidad
#Dado que son aleatorios, veamos que son estables
cf1 <- clusterboot(datos, 
                   #numero de remuestreos
                   B=100,
                   #metodo de remuestreo
                   bootmethod = c("jitter", "boot"),
                   #Metodo de clustering
                   clustermethod = kmeansCBI,
                   #Numero de clusters a evaluar
                   krange=2, seed = 123)

cf1
  • Los clusters son extremadamente estables. Tanto el jittering (pequeñas perturbaciones) como el bootstrap (remuestreo con omisión) muestran índices de Jaccard muy cercanos a 1.
  • Ningún cluster se disolvió en ninguna de las 100 ejecuciones, y ambos se recuperaron en todas las ejecuciones.
  • Esto sugiere que la estructura de 2 clusters es muy robusta y probablemente representa una división natural y significativa en tus datos.

Realizamos de nuevo el agrupamiento en la base de datos.

1.4.0.0.0.3 Agrupamiento
# Calcule k-medias
set.seed(123)
km.3 <- kmeans(datos,#Base normalizada
               2,#número de cp
               nstart = 25)

#Visualizamos
fviz_cluster(
   km.3,
   data = datos,
   ellipse.type = "euclid",
   # Concentration ellipse
   star.plot = TRUE,
   # Avoid label overplotting (slow)
   ggtheme = theme_minimal()
)

#Agregar los clusters a la base de datos inicial: datos
aggregate(datos, #base normalizada
          by=list(cluster=km.3$cluster),#cluster
          mean) #media de las variables

1.5 Conclusiones

La relación de los pacientes con las enfermeras es el principal factor en determinar la experiencia en el servicio de salud. Asimismo, la experiencia del cliente esta determinada en buena parte por el entendimiento de la atención y el diagnóstico que recibio derivado de la consulta.

Finalmente, queda pendiente el ejercicio para evaluar de forma más oportuna cada tema por tipo de evaluación a los centros de atención, sin embargo, estos resultados nos pueden dar una primera señal.