Algorimo Birtch: Estado nutricional de niñas y niños menores de 1 a 3 años de la Región de Ayacucho - Perú

1 Resumen

El presente proyecto evalúa variables relacionadas a niñas y niños menores de 1 a 3 años de la Región de Ayacucho, que fueron atendidos en los establecimientos de salud de Dirección de Regional de Salud (Diresa) respectiva, en el año 2021, con el objeto de obtener perfiles nutricionales que permitan focalizar el servicio de los centros de salud respectivos. Se aplicaron tres (3) metodologías de agrupamiento: Cluster Jerarquico, cluster basado en particiones y algoritmo BIRCH.

2 Introducción

La desnutrición crónica en la infancia es uno de los principales problemas de salud pública y un indicador de desarrollo social del país; tiene efectos irreversibles en el desarrollo de habilidades y capacidades en la niña y el niño.

El INEI en su informe “Perú: Indicadores de resultados de los programas presupuestales, 2015-2020” señala que el 12.1% de la población menor de cinco años de edad del país sufrió desnutrición crónica, con base al patrón de la OMS, siendo 0.1% menor que el 2019 y -2.3% respecto al 2015. A su vez, el mismo informe señala que el 40% de niños y niñas entre 6 a 35 meses en el 2020 presentaron prevalencia de anemia.

El diseño de políticas públicas, y en específico, el desarrollo de programas y proyectos sociales, para mejorar el estado nutricional de los menores de cinco años es una prioridad para el desarrollo nacional; siendo importante la generación de investigaciones que permitan conocer sus perfiles nutricionales

En ese sentido el presente trabajo tiene como objetivo la obtención de clústeres para conocer los perfiles nutricionales de niños y niñas menores de 1 a 3 años que fueron atendidos en establecimientos de salud de primer nivel de la Diresa Ayacucho el año 2021, con fines de identificar patrones que permitan la planificación de intervenciones focalizadas que contribuyan a mejorar su calidad de vida.

3 Descripción de los datos

Los datos corresponden a información de historias clínicas de los establecimientos de salud del Ministerio de Salud del Perú (Minsa), y que es registrada en el Sistema de Información del Estado Nutricional (SIEN), administrado por el Centro Nacional de Alimentación y Nutrición (CENAN), órgano de línea técnico normativo del Instituto Nacional de Salud (INS). La base de datos que usamos cuenta con un total de ““7,599”” registros de menores de 1 a 3 años y fue tomada de la Plataforma Nacional de Datos Abiertos.

Fuente

Unidad de análisis Niñas y niños menores de 1 a 3 años que fueron atendidos en los establecimientos de salud del primer nivel de atención de la Diresa Ayacucho el año 2021.

Variables Las variables utilizadas en el análisis son las siguientes:

Variables predictoras: Edad/meses: Edad en meses del menor atendido. Peso (kg) Talla (cm) PTZ: Índice que compara el peso del menor con el peso esperado para su talla y permite establecer si ha ocurrido una pérdida/ganancia de peso corporal. ZTE: Índice que compara la talla del menor con la talla esperada para su edad y permite establecer si está ocurriendo un retraso en el crecimiento. ZPE: Índice que compara el peso del menor con el peso esperado para su edad y permite establecer si está ocurriendo desnutrición. IMC: Índice de masa corporal Hemoglobina Hemoglobina ajustada: Hemoglobina ajustada según la altura de la localidad de residencia del menor. AlturaREN: Metros sobre el nivel del mar de la localidad de residencia del menor.

4 Metodología

Utilizamos el algoritmo BIRCH (Balanced iterative reducing and clustering using hierarchies) y lo comparamos con los algoritmos K-meansy AGNES para determinar grupos de niños de acuerdo con sus características nutricionales.

El objetivo principal del algoritmo BIRCH es trabajar la clusterizacion con un gran volumen de datos a través del siguiente proceso:

• A cada uno de los clústers se le asigna un vector de valores, el CF (Clustering Feature).

• Clustering Feature: Si se tiene un set de N datos, el clustering feature se define como:

\(CF=(N,\bar{LS},SS)\)

Donde:

\(\bar{LS}=\sum_{i=0}^n\bar{X_i}\)

\(SS=\sum_{i=0}^n\bar{X_i}^2\)

• Si se tienen los datos \((x_1);(x_2);(x_3);(x_4)\)

  • El valor N será igual a 4

  • El valor de \(\bar{LS}\) es igual a \((x_1+x_2+x_3+x_4)\)

  • El valor de \(SS\) será igual a \((x_1^2+x_2^2+x_3^2+x_4^2)\)

  • De esta forma al final los puntos de muestra están representados como \(CF=(N,\bar{LS},SS)\)

• Luego de construye un árbol de clústers en orden jerárquico de tal modo que las hojas previas representan clústers de mayor tamaño, y que la división de estos produzcan a su vez nuevos nodos (clústers) de menor tamaño.

• El procedimiento culmina una vez que los clústers alcancen un tamaño inferior a un parámetro dado o se logre el número de clústers deseado [2].

• El algoritmo Birch cuenta con dos principales parámetros

  • Branching factor (B): Cantidad máxima de subclústers en cada nodo, en caso de superar este factor el nodo se divide en dos nodos con los subclústers redistribuidos.

  • Threshold (T): Es el umbral establecido para la creación de un nuevo subclúster. Establecer este valor bajo genera la creación de más clústers.

Ejemplo del Algoritmo

Considerando los siguientes valores:

N	Puntos
1	22
2	9
3	12
4	15

Paso 1: Se elige el Treshold (5).

Paso 2: Se elige el Branching factor (3).

Paso 3: Se toma el primer valor y se eleva al cuadrado.

Paso 4: Se eleva al cuadrado ambos valores y se suma.

Valores	N	Suma	Cuadrado
22	1	22	484
9	2	31	565

Paso 5: El centroide es igual a la suma de los valores sobre N, para el primer valor el centroide sería:

\(Centroide_0 = 22/1 = 22\)

Paso 6: Se aplica la fórmula

\(D=\sqrt{\frac{\sum_{i=1}^n\sum_{j=1}^n(x_i-x_j)^2}{n(n-1)}}\)

\(D=\sqrt{\frac{2n(SS)-2(LS)^2}{n(n-1)}}\)

\(D=\sqrt{\frac{2*2(565)-2(31)^2}{2(2-1)}}\)

Entonces \(D=13\), que es superior al treshold, por lo tanto, se crea un nuevo clúster.

Paso 7: Trabajando el siguiente valor: 12

\(Centroide_0 = 22/1 = 22\)

\(Centroide_1 = 9/1 = 9\)

El punto 12 esta más cerca a 9.

Valores	N	Suma	Cuadrado
9	1	9	81
12	2	21	225

\(D=\sqrt{\frac{2*2(225)-2(21)^2}{2(2-1)}}\)

\(D=3\)

Como \(D=3<5\), se juntan los puntos 9 y 12.

Paso 8: Trabajando el siguiente valor: 15

\(Centroide_2 = (9+12)/2 = 10.5\)

\(Centroide_0 = (22)/1 = 22\)

El punto 12 esta más cerca al \(Centroide_2=10.5\).

Valores	N	Suma	Cuadrado
9, 12	2	21	225
9, 12, 15	3	36	450

\(D=\sqrt{\frac{2*3(450)-2(36)^2}{2(2-1)}}\)

\(D=4.24\)

Como \(D=4.24<5\), se juntan los puntos 9, 12 y 15.

\(Centroide_3 = (9+12+15)/3 = 12\)

\(Centroide_0 = (22)/1 = 22\)

Ejemplo en Python

#Configuración del entorno
rm(list = ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
graphics.off()
options(scipen = 999)  
options(digits = 3) 
path_python="C:/Users/msarsozaa/Anaconda3"
#Cargando paquetes necesarios
library(pacman)
p_load(reticulate, PerformanceAnalytics, purrr, skimr, corrplot, cluster, psych, ggplot2, 
       stream, ellipse, tictoc,factoextra, NbClust, BiocManager, naniar, DataExplorer, 
       tidyverse, purrr, dplyr, readxl, readr,stats, DescTools, class,devtools,imager,knitr,kableExtra)#source,compareGroups 
#Cargando funciones de usuario
source("funciones.R")

library(reticulate)
use_python(path_python)

import pandas as pd
from sklearn.cluster import Birch
df_ejemplo = pd.DataFrame([22,9,12,15],columns=['Puntos'])

brc0 = Birch(threshold=5,
            branching_factor=3
            )
brc0.fit(df_ejemplo)

Birch(branching_factor=3, threshold=5)

C:\Users\MSARSO~1\ANACON~1\lib\site-packages\sklearn\cluster\_birch.py:717: ConvergenceWarning: Number of subclusters found (2) by BIRCH is less than (3). Decrease the threshold.
  warnings.warn(

pd.concat([df_ejemplo,pd.DataFrame(brc0.predict(df_ejemplo),columns=['Cluster'])],axis=1)

   Puntos  Cluster
0      22        0
1       9        1
2      12        1
3      15        1

5 Lectura de datos

# Cargamos la data
ayacucho <- read.csv("Data.csv")
head_5=head(ayacucho)
kable(head_5,caption = "Ayacucho") %>% kable_styling("striped") %>% scroll_box(width = "100%")

Ayacucho

Diresa	Microred	EESS	Dpto_EESS	Prov_EESS	Dist_EESS	Renipress	FechaAtencion	Sexo	FechaNacimiento	Juntos	SIS	Qaliwarma	EdadMeses	Peso	Talla	IMC	PTZ	ZTE	ZPE	Hemoglobina	AlturaREN	Hemoglobinaajustada	Cred	Suplementacion	Consejeria	Sesion
AYACUCHO	SANTA ROSA	I-1 - 00003765 - PUESTO DE SALUD COMUNPIARI	AYACUCHO	LA MAR	SANTA ROSA	3765	12/31/2021	F	12/31/2019	NA	NA	NA	24	12.0	76.0	15.8	2.53	-3.02	0.35	10.7	330	10.7	1	0	0	0
AYACUCHO	SAN MARTIN	I-3 - 00003762 - CENTRO DE SALUD SAN MARTIN	AYACUCHO	LA MAR	ANCO	3762	10/06/2021	F	3/24/2020	NA	NA	NA	19	11.8	79.7	14.9	1.80	-0.49	1.09	11.5	3215	9.4	1	0	0	0
AYACUCHO	PAMPA CANGALLO	I-4 - 00003507 - CENTRO DE SALUD PAMPA CANGALLO	AYACUCHO	CANGALLO	LOS MOROCHUCOS	3507	6/14/2021	M	05/05/2020	NA	NA	NA	13	8.9	73.3	12.1	-0.34	-1.63	-1.01	14.4	3330	12.1	1	1	0	0
AYACUCHO	NO PERTENECE A NINGUNA MICRORED	II-E - 00003575 - HOSPITAL JESUS NAZARENO	AYACUCHO	HUAMANGA	JESUS NAZARENO	3575	8/27/2021	M	5/27/2020	NA	NA	NA	15	8.4	74.7	11.2	-1.44	-1.76	-1.86	12.4	2780	10.8	0	1	0	0
AYACUCHO	SAN MARTIN	I-3 - 00003762 - CENTRO DE SALUD SAN MARTIN	AYACUCHO	LA MAR	ANCO	3762	10/30/2021	M	7/25/2020	NA	NA	NA	15	10.5	76.4	13.7	0.82	-1.16	0.11	11.7	3215	9.6	0	1	0	0
AYACUCHO	SANTA ROSA	I-4 - 00003764 - SANTA ROSA	AYACUCHO	LA MAR	SANTA ROSA	3764	9/20/2021	F	1/20/2019	NA	NA	NA	32	13.5	88.8	15.2	0.99	-0.94	0.24	11.6	330	11.6	1	1	0	0

str(ayacucho)

'data.frame':   7599 obs. of  27 variables:
 $ Diresa             : chr  "AYACUCHO" "AYACUCHO" "AYACUCHO" "AYACUCHO" ...
 $ Microred           : chr  "SANTA ROSA" "SAN MARTIN" "PAMPA CANGALLO" "NO PERTENECE A NINGUNA MICRORED" ...
 $ EESS               : chr  "I-1   - 00003765 - PUESTO DE SALUD COMUNPIARI" "I-3   - 00003762 - CENTRO DE SALUD SAN MARTIN" "I-4   - 00003507 - CENTRO DE SALUD PAMPA CANGALLO" "II-E  - 00003575 - HOSPITAL  JESUS NAZARENO" ...
 $ Dpto_EESS          : chr  "AYACUCHO" "AYACUCHO" "AYACUCHO" "AYACUCHO" ...
 $ Prov_EESS          : chr  "LA MAR" "LA MAR" "CANGALLO" "HUAMANGA" ...
 $ Dist_EESS          : chr  "SANTA ROSA" "ANCO" "LOS MOROCHUCOS" "JESUS NAZARENO" ...
 $ Renipress          : int  3765 3762 3507 3575 3762 3764 3600 3782 3783 3603 ...
 $ FechaAtencion      : chr  "12/31/2021" "10/06/2021" "6/14/2021" "8/27/2021" ...
 $ Sexo               : chr  "F" "F" "M" "M" ...
 $ FechaNacimiento    : chr  "12/31/2019" "3/24/2020" "05/05/2020" "5/27/2020" ...
 $ Juntos             : int  NA NA NA NA NA NA NA NA NA NA ...
 $ SIS                : int  NA NA NA NA NA NA NA NA NA NA ...
 $ Qaliwarma          : int  NA NA NA NA NA NA NA NA NA NA ...
 $ EdadMeses          : int  24 19 13 15 15 32 19 27 19 30 ...
 $ Peso               : num  12 11.8 8.9 8.4 10.5 ...
 $ Talla              : num  76 79.7 73.3 74.7 76.4 88.8 77 80 77 90.5 ...
 $ IMC                : num  15.8 14.9 12.1 11.2 13.7 ...
 $ PTZ                : num  2.53 1.8 -0.34 -1.44 0.82 0.99 -0.5 0.07 -0.63 0.82 ...
 $ ZTE                : num  -3.02 -0.49 -1.63 -1.76 -1.16 -0.94 -2.12 -2.45 -1.54 -0.31 ...
 $ ZPE                : num  0.35 1.09 -1.01 -1.86 0.11 0.24 -1.35 -1.33 -1.2 0.44 ...
 $ Hemoglobina        : num  10.7 11.5 14.4 12.4 11.7 ...
 $ AlturaREN          : int  330 3215 3330 2780 3215 330 2800 3499 3499 2734 ...
 $ Hemoglobinaajustada: num  10.7 9.4 12.1 10.8 9.6 ...
 $ Cred               : int  1 1 1 0 0 1 1 1 0 1 ...
 $ Suplementacion     : int  0 0 1 1 1 1 0 0 0 1 ...
 $ Consejeria         : int  0 0 0 0 0 0 0 0 0 1 ...
 $ Sesion             : int  0 0 0 0 0 0 0 0 0 0 ...

#skim(ayacucho)
#describe(ayacucho)

6 Preprocesamiento de datos

# Factorización de variables no númericas
variables_fac <- c("Diresa", "Microred", "EESS", "Dpto_EESS", "Prov_EESS", "Dist_EESS", "Renipress", "Sexo", "Juntos", "SIS", "Qaliwarma", "Cred", "Suplementacion", "Consejeria", "Sesion")
ayacucho[,variables_fac] <- lapply(ayacucho[,variables_fac], factor)
#str(ayacucho) cambio

6.1 Datos Perdidos

aya <-ayacucho[,14:23]
plot_missing(aya)

El análisis muestra que las variables seleccionadas no presentan valores perdidos o faltantes, por lo tanto, no se requiere realizar un procedimiento de imputación.

6.2 Analisis de Outliers

gg_box_density(aya)

A partir de los gráficos anteriores se observa que todas las variables a excepción de EdadMeses poseen valores outliers y extremos en algunos casos, como las variables hemoglobina y AlturaRen y Hemoglobinaajustada.

6.3 Tratamiento de datos Outliers

cols=c('Peso','Talla','IMC','PTZ','ZTE','ZPE','Hemoglobina','AlturaREN','Hemoglobinaajustada')
aya=fn_outliers(aya,cols,2,1.5)

gg_box_density_2(aya)

El procedimiento para tratar los valores outliers y extremos consistió en acotar los valores superiores y inferiores de las variables mediante el principio de 1.5 veces el rango intercuartílico, sumando este valor al percentil 75 para definir el valor superior de la variable y restando el valor indicado al percentil 25 para definir el valor inferior máximo aceptado.

Los valores que fueron superiores o inferiores a los valores en el rango definidos fueron reemplazados por el máximo y mínimo calculado de dicho rango.

6.4 Correlación

#cor <- ayacucho[,14:23]%>% chart.Correlation(histogram=TRUE, pch=15) 
correlacion<-round(cor(ayacucho[,14:23]), 2)
corrplot(correlacion, method="number", type="upper",number.cex = 0.72,tl.cex = 0.8,addCoef.col = 0.5)

De acuerdo con el análisis de correlación de Pearson se interpreta que no existen relación lineal negativa entre las variables.

Existe una correlación lineal positiva fuerte entre las variables:

EdadMeses y Talla: 0.86 Peso y Talla: 0.84 Peso e IMC: 0.92 IMC y PTZ: 0.82 IMC y ZPE: 0.78 PTZ y ZPE: 0.84 Hemoglobina y Hemoglobinaajustada: 0.84

Existe una correlación lineal positiva moderada entre las variables: EdadMeses y Peso: 0.69 EdadMeses y IMC: 0.44 Peso y PTZ: 0.55 Peso y ZTE: 0.42 Peso y ZPE: 0.61 Talla y ZPE: 0.42 Hemoglobina y AlturaRen: 0.44

6.5 Estandarización de los datos

aya <- as.data.frame(scale(aya))

7 Determinando factibilidad del cluster

7.1 Matriz de distancia euclidiana

dis.Data <- dist(aya, metric  = c("euclidean")) # Matriz de distancia #cambio
#dis.Data    <- daisy(aya, metric= "euclidean",stand = TRUE) # Matriz de distancia

7.2 Visualizando la matriz de distancia con fviz_dist()

 library(imager)
 if  (file.exists("graf_max.dist.png")) {      #Cambio  
     img <- load.image('graf_max.dist.png')
     plot(img)
  }

#else{
  # fviz_dist(dis.Data)
 # }

Se realizó el cálculo y visualización la matriz de distancia euclidiana utilizando las funciones fviz_dist () en el paquete factoextra r. Se observan áreas sombreadas de color rojo que indican la cercanía de las observaciones en base al cálculo de distancia indicado. Esto permite inferir que el set de datos es agrupable.

7.3 Estadístico Hopkins

re_process=F
if  (file.exists("res.Hopkins.RData") & !re_process) {      #Cambio  
     load('res.Hopkins.RData') 
}else{
  res.Hopkins <- get_clust_tendency(aya, 
                          n = nrow(aya) - 1, 
                          graph = FALSE, seed = 2022)
  save(res.Hopkins,file='res.Hopkins.RData')
}
res.Hopkins$hopkins_stat

[1] 0.853

Se realizó la prueba estadística de Hopkins, utilizando 0,5 como el umbral para determinar que es poco probable que el set de datos tenga conglomerados estadísticamente significativos. Observamos que el valor resultante de 0.853 se acerca a 1, entonces podemos concluir que el conjunto de datos es significativamente agrupable.

8 Análisis cluster jerarquico (hclust {stats})

8.1 Analizando el metodo de enlace optimo (coeficiente de aglomeración)

metodos= c("single", "complete", "average", "ward.D", "ward.D2")
COEF.AGL(dis.Data,metodos)

  single complete  average   ward.D  ward.D2 
   0.788    0.954    0.888    1.000    0.997 

Hacemos uso de la comparación de los coeficientes de aglomeración para determinar el número de conglomerados y proponer el mejor esquema de agrupación a partir de los diferentes resultados obtenidos al variar todas las combinaciones de métodos de enlace: “single”, “complete”, “average”, “ward.D”, “ward.D2”. Con base en el resultado en la gráfica se puede indicar que el enlace optimo es el de: Ward.D

8.2 Analizando el metodo de enlace optimo (matriz cofenetica)

#Debido a que graficamente no estan muy alejados con coeficientes de aglomeracion, probaremos el metodo de
#correlacion cofenetica para definir el mejor metodo de enlace
re_process=FALSE
if(file.exists("df_cof.RData") & !re_process) {      #Cambio  
     load('df_cof.RData') 
}else{
  metodos=c("ward.D", "ward.D2","complete", "average", "single")
  df_cof=fn_confenetico(dis.Data,metodos)
  save(df_cof,file='df_cof.RData')
}
kable(df_cof,caption = "Correlacion_cofenetica")%>% kable_styling("striped")

Correlacion_cofenetica

Metodo	Correl
average	0.551
single	0.464
ward.D2	0.393
ward.D	0.360
complete	0.350

Realizamos el análisis de la matriz cofenetica para determinar la similaridad entre la matriz de distancia original (euclidiana) y la matriz de las uniones de las observaciones en los cluster jerárquicos efectuados con los métodos de enlaces iniciales. Como resultado se observa que no hay una correlación fuerte pero el método average presenta la mayor correlación. Se decide evaluar el análisis de cluster jerárquico con ambos métodos: Ward.D y average.

8.3 Obteniendo el cluster con el metodo: Ward.D

Clus_AG_W <- hclust(dis.Data,method="ward.D")

8.4 Determinando de manera grafica el numero de clusters

library(gghighlight)
longitud=length(Clus_AG_W$height)
alturas <- data.frame(etapa = 1:longitud, distancia = Clus_AG_W$height)
ggplot(alturas) + aes(x = etapa, y = distancia)  +
  geom_point() + geom_line()  + 
  scale_x_continuous(breaks = seq(1, longitud, 1000)) + 
  geom_vline(xintercept = 7597, col = "red", lty = 2) + 
  geom_text(aes(label  = round(distancia,1)),
            size = 3, hjust= +1, vjust= -1) +
  theme_classic() #+ gghighlight(distancia > 11)

El método grafico nos sugiere la existencia de 2 a 5 agrupaciones.

8.5 Determinado el número de clusters optimo

re_process=F
if  (file.exists("res.nbclustW.RData")& !re_process) {      #Cambio  
     load('res.nbclustW.RData') 
}else{
  seed = 2022
  res.nbclustW <- NbClust(aya, distance ="euclidean", 
                         min.nc = 2, max.nc = 8, method = "ward.D", index ="all") #Cambio
  save(res.nbclustW,file='res.nbclustW.RData')
}
par(mfrow=c(1,1))  
fviz_nbclust_x(res.nbclustW)

Among all indices: 
===================
* 2 proposed  0 as the best number of clusters
* 1 proposed  1 as the best number of clusters
* 6 proposed  2 as the best number of clusters
* 8 proposed  3 as the best number of clusters
* 1 proposed  4 as the best number of clusters
* 1 proposed  5 as the best number of clusters
* 3 proposed  6 as the best number of clusters
* 2 proposed  8 as the best number of clusters
* 2 proposed  NA's as the best number of clusters

Conclusion
=========================
* According to the majority rule, the best number of clusters is  3 .

Mediante el uso de paquete NbClust podemos determinar el mejor número de agrupaciones tomando en cuenta los diferentes resultados obtenidos al variar todas las combinaciones la cantidad de clúster deseados para el método de enlace elegido: Ward.D.

8.6 Realizando clustering con K=3

# Cortando en 2 cluester
grp_WR=cutree(Clus_AG_W, k = 3)
# Number de casos en cada cluster
table(grp_WR)

grp_WR
   1    2    3 
3361 2919 1319 

# Descripción de cada cluster
med<-aggregate(aya_bkp, by=list(cluster=grp_WR), mean) #medias
kable(med) %>% kable_styling("striped") %>% scroll_box(width = "100%")

cluster	EdadMeses	Peso	Talla	IMC	PTZ	ZTE	ZPE	Hemoglobina	AlturaREN	Hemoglobinaajustada
1	26.5	12.66	86.2	14.7	0.785	-0.673	0.240	13.3	2797	11.8
2	19.9	9.99	78.7	12.7	-0.149	-1.399	-0.796	12.5	2827	11.0
3	30.9	11.32	85.6	13.2	-0.400	-1.803	-1.273	13.7	2986	11.9

#knitr::kable(med, format = "markdown")

Se procedió a cortar la agrupación general (dendograma) en 03 clústers tomando como base el resultado óptimo de Nbclust y obtenemos la cantidad de observaciones por cada uno de ellos.

8.7 Caracterízando los cluster

8.7.1 Graficando los clusters.

#fviz_dend(Clus_AG, k = 2, cex = 0.7, horiz = FALSE, k_colors = "jco",
          #rect = TRUE, rect_border = "jco", rect_fill = TRUE)

8.7.2 Diagrama de caracterización - lineas

data_plot=scale(aya_bkp) #Scale
#data_plot=apply(aya_bkp, 2, normalize) #Max-Min
M<-as.data.frame(t(rbind(aggregate(data_plot, by=list(cluster=grp_WR), mean)[,-1])))
a=as.vector(colMeans(data_plot))
fin=data.frame(M,a,names(aya_bkp));names(fin)<-c("Clus1","clus2","clus3","Media","var")
#fin=data.frame(M,a,names(aya_bkp));names(fin)<-c("Clus1","clus2","Media","var")
ali=melt(fin,id.vars = "var")
ggplot(ali, aes(x=var,y=round(value,1),group=variable,colour=variable)) +
  geom_point()+ geom_line(aes(lty=variable))+ expand_limits(y = c(-1.9, 1.9))+
   theme(axis.text.x = element_text(angle = 60, vjust = 0.5, hjust=1))

Se calcula los promedios de las observaciones con la base de datos original agrupado por cada clúster. Se observa que el clúster 01 tiene mayor cantidad de individuos y estos se caracterizan por tener mayor peso, talla, el mínimo retraso en su crecimiento, un índice nutricional positivo, mayor índice de masa corporal y viven a una altura menor. Por otro lado, se puede observar que el clúster 03 posee los individuos de mayor edad, pero, el peso, el índice de masa corporal, talla es menor al clúster 01 y poseen mayor desnutrición.

8.7.3 Diagrama de caracterización - boxplot

dd <- cbind(aya_bkp, cluster =grp_WR )
dd$cluster<-as.factor(dd$cluster)
df.m <- melt(dd, id.var = "cluster")
p <- ggplot(data = df.m, aes(x=variable, y=value)) + 
  geom_boxplot(aes(fill=cluster))+ facet_wrap( ~ variable, scales="free") 
p

Las gráficas permiten analizar los clústeres de manera individual a través de sus variables. Por ejemplo, el clúster con mayor altura de residencia es el número 03. El clúster con menor talla 03, el clúster con mayor índice de masa corporal es el clúster 01.

8.8 Obteniendo el cluster con el metodo: average

Clus_AG_AV <- hclust(dis.Data,method="average")

8.9 Determinando de manera grafica el numero de clusters

library(gghighlight)
longitud=length(Clus_AG_AV$height)
alturas <- data.frame(etapa = 1:longitud, distancia = Clus_AG_AV$height)
ggplot(alturas) + aes(x = etapa, y = distancia)  +
  geom_point() + geom_line()  + 
  scale_x_continuous(breaks = seq(1, longitud, 1000)) + 
  geom_vline(xintercept = 7597, col = "red", lty = 2) + 
  geom_text(aes(label  = round(distancia,1)),
            size = 3, hjust= +1, vjust= -1) +
  theme_classic() #+ gghighlight(distancia > 11)

El método grafico nos sugiere la existencia de 2 a 5 agrupaciones.

8.10 Determinado el número de clusters optimo

re_process=FALSE
if  (file.exists("res.nbclustAV.RData")& !re_process) {      #Cambio  
     load('res.nbclustAV.RData') 
}else{
  seed = 2022
  res.nbclustAV <- NbClust(aya, distance ="euclidean", 
                         min.nc = 2, max.nc = 8, method = "average", index ="all") #Cambio
  save(res.nbclustAV,file='res.nbclustAV.RData')
}
par(mfrow=c(1,1))  
fviz_nbclust_x(res.nbclustAV)

Among all indices: 
===================
* 2 proposed  0 as the best number of clusters
* 1 proposed  1 as the best number of clusters
* 9 proposed  2 as the best number of clusters
* 8 proposed  3 as the best number of clusters
* 2 proposed  7 as the best number of clusters
* 4 proposed  8 as the best number of clusters

Conclusion
=========================
* According to the majority rule, the best number of clusters is  2 .

los diferentes resultados obtenidos al variar todas las combinaciones la cantidad de clúster deseados para el método de enlace elegido: average.

8.11 Realizando clustering con K=2

# Cortando en 2 cluester
grp_AV=cutree(Clus_AG_AV, k = 2)
# Number de casos en cada cluster
table(grp_AV)

grp_AV
   1    2 
7218  381 

# Descripción de cada cluster
med<-aggregate(aya_bkp, by=list(cluster=grp_AV), mean) #medias
kable(med) %>% kable_styling("striped") %>% scroll_box(width = "100%")

cluster	EdadMeses	Peso	Talla	IMC	PTZ	ZTE	ZPE	Hemoglobina	AlturaREN	Hemoglobinaajustada
1	24.6	11.2	83.0	13.5	0.133	-1.192	-0.507	13	2841	11.5
2	27.2	14.4	87.9	16.3	1.875	-0.321	1.211	13	2844	11.5

#knitr::kable(med, format = "markdown")

Se realiza el corte de la agrupación general (dendograma) en 02 clústeres, tomando como base el resultado óptimo de Nbclust y obtenemos la cantidad de observaciones por cada uno de ellos. Se puede observar que el resultado presenta dos grupos, con la mayor cantidad de individuos agrupados en clúster 01. ## Caracterízando los cluster

8.11.1 Graficando los clusters.

#fviz_dend(Clus_AG, k = 2, cex = 0.7, horiz = FALSE, k_colors = "jco",
          #rect = TRUE, rect_border = "jco", rect_fill = TRUE)

8.11.2 Diagrama de caracterización - lineas

data_plot=scale(aya_bkp) #Scale
#data_plot=apply(aya_bkp, 2, normalize) #Max-Min
M<-as.data.frame(t(rbind(aggregate(data_plot, by=list(cluster=grp_AV), mean)[,-1])))
a=as.vector(colMeans(data_plot))
#fin=data.frame(M,a,names(aya_bkp));names(fin)<-c("Clus1","clus2","clus3","Media","var")
fin=data.frame(M,a,names(aya_bkp));names(fin)<-c("Clus1","clus2","Media","var")
ali=melt(fin,id.vars = "var")
ggplot(ali, aes(x=var,y=round(value,1),group=variable,colour=variable)) +
  geom_point()+ geom_line(aes(lty=variable))+ expand_limits(y = c(-1.9, 1.9))+
   theme(axis.text.x = element_text(angle = 60, vjust = 0.5, hjust=1))

El cálculo del promedio por cada grupo no permite diferenciarlos correctamente con una marcada diferencia en muchas variables.

8.11.3 Diagrama de caracterización - boxplot

dd <- cbind(aya_bkp, cluster =grp_AV )
dd$cluster<-as.factor(dd$cluster)
df.m <- melt(dd, id.var = "cluster")
p <- ggplot(data = df.m, aes(x=variable, y=value)) + 
  geom_boxplot(aes(fill=cluster))+ facet_wrap( ~ variable, scales="free") 
p

8.11.4 Usando el indice de rand para definir el cluster optimo

library(fossil)
rand.index(grp_WR,grp_AV)

[1] 0.39

Segun el indice de Rand, no se podria comprobar que las soluciones (ward.D y average) clusters son parecidos, por ese motivos elegiremos el que tiene mayor serparacón a nivel grafico del lineas.

9 Análisis cluster: K – Means

9.1 Determinando número óptimo de clusters

9.1.1 Silhouette method

set.seed(123)
re_process=FALSE
if  (file.exists("gg_sil.RData") & !re_process) {      #Cambio  
     load('gg_sil.RData') 
}else{
   gg_sil=fviz_nbclust(aya, kmeans, method = "silhouette")+
          labs(subtitle = "Silhouette method")
   save(gg_sil,file='gg_sil.RData')
}
gg_sil

El resultado del índice de silueta establece que para K-means es óptimo utilizar 3 clústeres, de este modo se generaron 3 clústeres con 2962, 2210 y 2427 respectivamente.

9.1.2 Elbow method (WSS)

set.seed(123)
re_process=FALSE
if  (file.exists("gg_WSS.RData") & !re_process) {      #Cambio  
     load('gg_WSS.RData') 
}else{
   gg_WSS=fviz_nbclust(aya, kmeans, method = "wss") + geom_vline(xintercept = 3, linetype = 2)+
          labs(subtitle = "Elbow method")
   save(gg_WSS,file='gg_WSS.RData')
}
gg_WSS

9.2 Realizando clustering con K=3

km.res <- kmeans(aya, centers=3,nstart = 25)
grp_km=km.res$cluster
table(grp_km)

grp_km
   1    2    3 
2962 2210 2427 

# Descripción de cada cluster
med<-aggregate(aya_bkp, by=list(cluster=grp_km), mean)
kable(med) %>% kable_styling("striped") %>% scroll_box(width = "100%")

cluster	EdadMeses	Peso	Talla	IMC	PTZ	ZTE	ZPE	Hemoglobina	AlturaREN	Hemoglobinaajustada
1	18.6	9.77	77.6	12.6	-0.089	-1.436	-0.765	12.6	2841	11.1
2	26.7	13.19	86.9	15.1	1.089	-0.496	0.547	13.0	2769	11.6
3	30.5	11.77	86.7	13.6	-0.191	-1.390	-0.882	13.6	2909	12.0

#knitr::kable(med, format = "markdown")

9.3 Caracterízando los cluster

9.3.1 Graficando los clusters.

#fviz_cluster(km.res, data = aya,
            # palette = "jco",
             #ellipse.type = "euclid", # Concentration ellipse
             #star.plot = TRUE, # Add segments from centroids to items
             #repel = TRUE, # Avoid label overplotting (slow)
             #ggtheme = theme_minimal())

9.3.2 Diagrama de caracterización - lineas

data_plot=scale(aya_bkp) #Scale
#data_plot=apply(aya_bkp, 2, normalize) #Max-Min
M<-as.data.frame(t(rbind(aggregate(data_plot, by=list(cluster=grp_km), mean)[,-1])))
a=as.vector(colMeans(data_plot))
fin=data.frame(M,a,names(aya_bkp));names(fin)<-c("Clus1","clus2","clus3","Media","var")
#fin=data.frame(M,a,names(aya_bkp));names(fin)<-c("Clus1","clus2","Media","var")
ali=melt(fin,id.vars = "var")
ggplot(ali, aes(x=var,y=round(value,1),group=variable,colour=variable)) +
  geom_point()+ geom_line(aes(lty=variable))+ expand_limits(y = c(-1.9, 1.9))+
   theme(axis.text.x = element_text(angle = 60, vjust = 0.5, hjust=1))

Con respecto a la caracterización de los clústeres, se desprende del gráfico que el valor más cercano entre clústeres es la altura, mientras que el valor más alejado es el peso.

9.3.3 Diagrama de caracterización - boxplot

dd <- cbind(aya_bkp, cluster =grp_km )
dd$cluster<-as.factor(dd$cluster)
df.m <- melt(dd, id.var = "cluster")
p <- ggplot(data = df.m, aes(x=variable, y=value)) + 
  geom_boxplot(aes(fill=cluster))+ facet_wrap( ~ variable, scales="free") 
p

El gráfico de cajas indica que la variable más cercana entre clústers es la alturaREN, y la más lejana es el IMC.

Para la clusterización utilizando el algoritmo Birch, se hizo uso de un Threshold igual a 0.42 y un Branching Factor igual a 50.

10 Análisis cluster Birch (reticulate)

10.1 Configurando entorno de desarrollo

library(reticulate)
use_python(path_python)
aya_p = r_to_py(aya) #r.aya

10.2 Cargando librerias de python en R

import pandas as pd 
from sklearn.preprocessing import StandardScaler 
from sklearn.cluster import Birch 

10.3 Definiendo Birch con: threshold=1, branching_factor=2

df = r.aya_p# Obteniendo objeto R to python
brc = Birch(threshold=0.42,branching_factor=50,n_clusters=3)
brc.fit(df)

Birch(threshold=0.42)

10.4 Realizando predicción

data_total = pd.concat([df,pd.DataFrame(brc.predict(df))],axis=1)
data_total.head()

   EdadMeses      Peso     Talla  ...  AlturaREN  Hemoglobinaajustada  0
0  -0.103463  0.338491 -1.171154  ...  -1.921536            -0.782780  2
1  -0.794094  0.247707 -0.569719  ...   1.022897            -2.052721  2
2  -1.622850 -1.420439 -1.610039  ...   1.337883             0.584849  1
3  -1.346598 -1.704137 -1.382469  ...  -0.168571            -0.685092  1
4  -1.346598 -0.523952 -1.106134  ...   1.022897            -1.857346  2

[5 rows x 11 columns]

10.5 Proporciones por cada cluster

data_total.groupby(0)['Peso'].count()

0
0    2788
1    2404
2    2407
Name: Peso, dtype: int64

La proporción para los 3 clústeres es de 2788, 2404 y 2407 respectivamente.

10.6 Guardamos el objeto con Cluster Birch

Para luego subirlo al R

data_total.to_excel('birch.xlsx')

10.7 Leemos el objeto con Cluster Birch en R

bi <- read_excel("birch.xlsx")
birch <- as.data.frame(bi[,2:12])
names (birch)[11] = "cluster"
birch$cluster <- factor(birch$cluster, levels = c(0,1,2), labels = c(1,2,3))

head(birch)

  EdadMeses   Peso  Talla    IMC    PTZ     ZTE    ZPE Hemoglobina AlturaREN
1    -0.103  0.338 -1.171  1.593  2.522 -2.0201  0.848      -1.968    -1.922
2    -0.794  0.248 -0.570  0.899  1.725  0.7099  1.663      -1.296     1.023
3    -1.623 -1.420 -1.610 -1.116 -0.612 -0.5202 -0.649       1.138     1.338
4    -1.347 -1.704 -1.382 -1.782 -1.813 -0.6605 -1.585      -0.541    -0.169
5    -1.347 -0.524 -1.106  0.054  0.654 -0.0131  0.584      -1.128     1.023
6     1.002  1.190  0.909  1.157  0.840  0.2243  0.727      -1.212    -1.922
  Hemoglobinaajustada cluster
1             -0.7828       3
2             -2.0527       3
3              0.5848       2
4             -0.6851       2
5             -1.8573       3
6              0.0964       3

10.8 Caracterízando los cluster

data_plot=birch[1:10]
cluster=birch$cluster
table(cluster)

cluster
   1    2    3 
2788 2404 2407 

# Descripción de cada cluster
med<-aggregate(aya, by=list(cluster=cluster), mean)

med

  cluster EdadMeses    Peso   Talla    IMC    PTZ    ZTE    ZPE Hemoglobina
1       1     0.658  0.7457  0.7675  0.590  0.213  0.346  0.321      0.6561
2       2    -0.618 -0.9573 -0.8110 -0.913 -0.605 -0.546 -0.711      0.0541
3       3    -0.144  0.0924 -0.0789  0.229  0.358  0.144  0.338     -0.8140
  AlturaREN Hemoglobinaajustada
1     0.122               0.643
2     0.270              -0.146
3    -0.411              -0.599

# knitr::kable(med) %>% kable_styling("striped") %>% scroll_box(width = "100%")

10.9 Diagrama de caracterización - lineas

M<-as.data.frame(t(rbind(aggregate(data_plot, by=list(cluster=cluster), mean)[,-1])))

a=as.vector(colMeans(data_plot))
fin=data.frame(M,a,names(aya_bkp));names(fin)<-c("Clus1","Clus2","Clus3","Media","var")

ali=melt(fin,id.vars = "var")
ggplot(ali, aes(x=var,y=round(value,1),group=variable,colour=variable)) +
  geom_point()+ geom_line(aes(lty=variable))+ expand_limits(y = c(-1.9, 1.9))+
   theme(axis.text.x = element_text(angle = 60, vjust = 0.5, hjust=1))

No existen diferencias tan marcadas en AlturaREN entre el cluster 1 y 2. En PTZ, ZPE y ZTE tampoco hay diferencias marcadas entre el cluster 1 y 2. Mientras que diferencias más marcadas entre Hemoglobina, HemoglobinaAjustada, IMC, Peso y Talla

10.9.1 Diagrama de caracterización - boxplot

dd <- cbind(aya_bkp, cluster =cluster )
dd$cluster<-as.factor(dd$cluster)
df.m <- melt(dd, id.var = "cluster")
p <- ggplot(data = df.m, aes(x=variable, y=value)) + 
  geom_boxplot(aes(fill=cluster))+ facet_wrap( ~ variable, scales="free") 
p

El diagrama de cajas indica que la variable con los valores más cercanos entre clústers es la AlturaREN, mientras que el más alejado es la Edad en Meses.

11 Validando los clusters

11.1 Interna (cohesión y separación)

re_process=TRUE
# if  (file.exists("df_indexs.RData") & !re_process) {      #Cambio  
     # load('df_indexs.RData') 
# }else{
    library(clusterSim)
    library(clValid)
    #Indice de Davies-Bouldin: buscamos el valor mas alto posible
    DBkm_ <- index.DB(aya, km.res$cluster, centrotypes = "centroids")$DB #kmeans
    DBHc_ <-index.DB(aya, grp_WR, d=dis.Data,centrotypes="centroids")$DB #hclust
    DBirch_ <-index.DB(aya, as.numeric(cluster), centrotypes="centroids")$DB #Birch
    
    #Indice de dunn: buscamos el valor mas bajo posible
    Dnkm_ <- dunn(Data = aya, clusters = km.res$cluster, distance = NULL)#kmeans
    DnHc_ <- dunn(dis.Data, grp_WR) #hclust
    DnBirch_ <- dunn(Data = aya, clusters = as.numeric(cluster), distance = NULL)#Birch
    
    tipo_=c('kmeans','hclust','Birch')
    davies.bouldin_=c(DBkm_,DBHc_,DBirch_)
    dunn_=c(Dnkm_,DnHc_,DnBirch_)
    df_indexs=data.frame(tipo_,davies.bouldin_,dunn_)
    colnames(df_indexs)=c('Cluster','Davies.Bouldin','Dunn')
    save(df_indexs,file='df_indexs.RData')
# }
    # df_indexs
# kable(df_indexs) %>% kable_styling("striped") %>% scroll_box(width = "100%")

silhouette<-rbind(
mean(silhouette(as.numeric(km.res$cluster) ,dis.Data)[,3]), #kmeans
mean(silhouette(as.numeric(grp_WR) ,dis.Data)[,3]), #hclust
mean(silhouette(as.numeric(birch$cluster) ,dis.Data)[,3]) #Birch
)

cbind(df_indexs, silhouette)

  Cluster Davies.Bouldin   Dunn silhouette
1  kmeans           1.74 0.0100      0.184
2  hclust           1.81 0.0287      0.125
3   Birch           2.17 0.0174      0.134

De acuerdo Davies.Bouldin el mayor valor es método es kmeans, respecto a Dunn el menor valor es hclust y con respecto a silhouette el mejor es k-means.

Creamos una nueva variable llamada Anemia en base a Hemoglobinaajustada para analizarlos por grupos de resultados de K-means.

ayacucho_final<-ayacucho
ayacucho_final$Anemia<-if_else(ayacucho$Hemoglobinaajustada*10>=110, 'Sin anemia', 
        if_else(ayacucho$Hemoglobinaajustada*10>=100,'Leve',
                if_else(ayacucho$Hemoglobinaajustada*10>=70,'Moderada',
                        'Grave')))
        
ayacucho_final=cbind(ayacucho_final, cluster=km.res$cluster)

table(ayacucho_final$Anemia)

     Grave       Leve   Moderada Sin anemia 
        11       1186        645       5757 

Tomando los resultados de K-means

g1_F <-
  ggplot(mutate(ayacucho_final, cluster = factor(cluster))) +
  aes(cluster, fill =Anemia ) +
  geom_bar(cluster = position_fill()) +
  labs(title="cluster según Anemia",
            x = NULL, y = "Proporción") +
  theme_bw()
g1_F

Vemos que el mayor grupo de niños con desnutrición moderada y leve esta en el cluster 1, seguido del cluster 2 y por ultimo el cluster 3. El mayor número de niños sin anemia esta en el cluster 3.

Del mismo modo transponiendo los ejes del grafico anterior, el cluster 1 tiene la mayor proporción de niños con anemia leve, lo mismo ocurre en anemia moderada.

g2_F <-
  ggplot(mutate(ayacucho_final, cluster = factor(cluster))) +
  aes(Anemia, fill =cluster ) +
  geom_bar(Anemia = position_fill()) +
  labs(title="Anemia según cluster",
            x = NULL, y = "Proporción") +
  theme_linedraw()
g2_F

12 Conclusiones

• La evaluación tanto por K-prototypes como con la función NbClust nos indica que debemos trabajar con 3 grupos.

• La metodología AGNES y Fuzzy C-Means en las variables numéricas generan clústeres con un comportamiento muy similar, sin embargo, en el caso de las variables categóricas se tienen diferencias significativas.

• Se determina mediante el gráfico de perfiles que la edad gestacional no genera un mayor aporte al comportamiento de los grupos, ya que en los tres los valores

• Finalmente considerando AGNES, se tiene: Grupo 1: Gestantes en su mayoría de embarazo simple, de Regiones Apurímac, Huancavelica y Pasco, por lo que tienen un nivel de altitud mayor y un nivel de hemoglobina superior a los demás grupos. Grupo 2: Gestantes en su mayoría de embarazo múltiple, de la región de Ayacucho, por lo que tienen un nivel de altitud intermedia y con hemoglobina IMC, Peso, PPG y Talla en un nivel medio. Grupo 3: Gestantes con un embarazo simple de la región de Junín con el

  GRUPO	ATRIBUTO	AGNES	K-Means	BIRCH
  Grupo 1	Cantidad	3,361	2,962	2,788
  	%	44.2%	39.0%	36.7%
  Grupo 2	Cantidad	2,919	2,210	2,404
  	%	38.4%	29.1%	31.6%
  Grupo 3	Cantidad	1,319	2,427	2,407
  	%	17.4%	31.9%	31.7%

13 Referencias

• Fanny Ramadhani et al 2020 IOP Conf. Ser.: Mater. Sci. Eng. 725 012090

• Alonso del Saso, Javier “Métodos de detección de anomalías y clustering en series temporales” 2020, Santander, España.