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 artículo evalúa variables relacionadas a niñas y niños menores de 1 a 3 años de la Región de Ayacucho - Perú, que fueron atendidos en los establecimientos de salud de Dirección de Regional de Salud (Diresa) respectiva, 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: AGNES, K-Means y BIRCH; obteniendo que el mejor índice de silueta lo brinda la metodología …… con (n) grupos debido a que el perfilado resultante permite observar las diferencias entre infantes, principalmente por las variables …… y ……

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 y 2 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 ““6 047”” registros de menores de 1 y 2 años y fue tomada de la Plataforma Nacional de Datos Abiertos .

Unidad de análisis Niñas y niños menores de 1 y 2 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 método BIRCH (Balanced iterative reducing and clustering using hierarchies) para determinar grupos de niños de acuerdo con sus características nutricionales. El objetivo principal de esta metodología es trabajar la agrupación con grandes cantidades de datos [1] 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.

Algoritmo Considerando los siguientes valores:

N Puntos 1 22 2 9 3 12 4 15

1. Se elige el Treshold (5).
2. Se elige el Branching factor (3).
3. Se toma el primer valor y se eleva al cuadrado.
4. Se eleva al cuadrado ambos valores y se suma.

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

5. El centroide es igual a la suma de los valores sobre N, para el primer valor el centroide sería: Centroide = 22/1 = 22
6. Se aplica la fórmula

………………………..

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

7. Trabajando el siguiente valor:

Valores N Suma Cuadrado 1 9 81 12 2 21 225

8. Se vuelve a aplicar la fórmula, y se obtiene D = 3, por lo tanto, el tercer valor se designa al clúster con el centroide más cercano, el cual es el segundo clúster ya que 12 está más cerca de 9 que de 22.

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.

5 Configuración inicial

#Configuración del entorno
rm(list = ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
graphics.off()
options(scipen = 999)  
options(digits = 3) 
path_python="C:/Users/adwin/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")

6 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)

7 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

7.1 Datos Perdidos

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

7.2 Analisis de Outliers

gg_box_density(aya)

7.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)

7.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)

7.5 Estandarización de los datos

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

8 Determinando factibilidad del cluster

8.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

8.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)
 # }

8.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

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

9.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 

9.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

9.3 Obteniendo el cluster con el metodo: Ward.D

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

9.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)

9.5 Determinado el número de clusters optimo

re_process=FALSE
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 .

9.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")

9.7 Caracterízando los cluster

9.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)

9.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))

9.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

9.8 Obteniendo el cluster con el metodo: average

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

9.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)

9.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 .

9.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")

9.12 Caracterízando los cluster

9.12.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)

9.12.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))

9.12.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

9.12.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.

10 Análisis cluster: K – Means

10.1 Determinando número óptimo de clusters

10.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

10.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")

10.3 Caracterízando los cluster

10.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())

10.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))

10.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

11 Análisis cluster Birch (reticulate)

11.1 Configurando entorno de desarrollo

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

11.2 Cargando librerias de python en R

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

11.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)

11.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]

11.5 Proporciones por cada cluster

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

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

11.6 Guardamos el objeto con Cluster Birch

Para luego subirlo al R

data_total.to_excel('birch.xlsx')

11.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

11.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%")

11.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))

11.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

12 Validando los clusters

12.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

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 

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

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

13 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