##############################
# ANÁLISIS CLUSTER #
# ANÁLISIS JERÁRQUICO #
# EJEMPLO Pancreatitis #
# Mg. Jesús Salinas Flores #
# [email protected] #
##############################
#--------------------------------------------------------------
# Para limpiar el workspace, por si hubiera algun dataset
# o información cargada
rm(list = ls())
dev.off()
#--------------------------------------------------------------
# Cambiar el directorio de trabajo
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
getwd()
#--------------------------------------------------------------
# Otras opciones
options(scipen=999) # Eliminar la notación científica
options(digits = 3) # Número de decimales
#--------------------------------------------------------------
# Paquetes
library(factoextra)
library(cluster)
library(StatMatch)
library(corrplot)
library(ade4)
library(fpc)
library(dplyr)
library(DataExplorer)
library(funModeling)
library(aplpack)
library(TeachingDemos)
library(DescTools)
###############################
# I. PREPARACIÓN DE LOS DATOS #
###############################
# Titulo: "UTILIDAD DE BISAP Y PCR COMO PREDICTORES DE
# PANCREATITIS GRAVE"
# Presentada por: JOSÉ MARIANO PINZÓN FERNÁNDEZ
# Tesis presentada para obtener el grado de Maestro en
# Ciencias Médicas con Especialidad en Medicina Interna
# Universidad de San Carlos de Guatemala
#
# Se incluyeron 88 pacientes durante los meses de Enero de
# 2013 a marzo de 2014 en el departamento de medicina
# interna del HR de Guatemala con pancreatitis aguda a
# quienes se les calcularon índices pronósticos de
# severidad de ingreso incluyendo BISAP APACHE II y PCR a
# las 48 hrs así como estudios bioquímicos y de imágenes
# en búsqueda de la causa etiológica y/o complicaciones.
#--------------------------------------------------------------
# Cargando los datos
datos <- read.csv("database_bisap.csv", header = TRUE)
attach(datos)
str(datos)
# Vamos a utilizar únicamente variables independientes
# numéricas.
# pcr, ldh, bun, creatinina, ph, haco3, pco2, amilasa,
# lac, sodio, calcio, potasio, globulos_blancos, sp02,
# temperatura, presión arterial sistólica, presión arterial
# diastólica y frecuencia cardiaca.
my_data <- data.frame(pcr,ldh,bun,creatitina,ph,hco3,pco2,
amilasa,lac,sodio,calcio,potasio,
globulos_blancos,spo2,temperatura,
presion_arterial_sistolica,
presion_arterial_diastolica,
plaquetas)
# Visualizando la estructura de los datos
str(my_data)
# Visualizando los primeros 5 registros
head(my_data, n = 5)
# Visualizando los últimos 5 registros
tail(my_data, n = 5)
#############################
# II. ANÁLISIS EXPLORATORIO #
#############################
summary(my_data)
#--------------------------------------------------------------
library(DataExplorer)
# Estructura de los datos
plot_str(my_data)
# Detectando y graficando los % de datos perdidos
plot_missing(my_data)
# Histograma para las variables numéricas
plot_histogram(my_data)
create_report(my_data)
#--------------------------------------------------------------
library(funModeling)
# Descripción de los datos
df_status(my_data)
# Gráfico de variables numéricas
plot_num(my_data)
# Descripción de las variables numéricas
profiling_num(my_data)
####################################
# III. USANDO MEDIDAS DE DISTANCIA #
####################################
# Principales funciones: dist(), get_dist() y daisy()
# Paquetes stats, factoextra y cluster
# Para estandarizar los datos (center y scale)
datos.e <- as.data.frame(scale(my_data))
str(datos.e)
#--------------------------------------------------------------
# 1.a Calculando la matriz de distancia euclidiana
# con la funcion dist()
dist.eucl <- dist(datos.e, method = "euclidean")
# Visualizando un subconjunto de la matriz de distancia
round(as.matrix(dist.eucl)[1:6, 1:6], 1)
# Otras distancias: "maximum", "manhattan", "minkowski"
#--------------------------------------------------------------
# 1.b Calculando la matriz de distancia euclidiana
# con la funcion daisy()
library(cluster)
dist.eucl2 <- daisy(datos.e, metric= "euclidean")
# Visualizando un subconjunto de la matriz de distancia
round(as.matrix(dist.eucl2)[1:6, 1:6], 1)
# Otras distancias: "gower", "manhattan"
# 1.c Calculando la matriz de distancia euclidiana con la
# funcion get_dist()
library(factoextra)
res.dist <- get_dist(my_data, stand = TRUE,
method = "euclidean")
# Visualizando un subconjunto de la matriz de distancia
round(as.matrix(res.dist)[1:6, 1:6], 1)
# Otras distancias: "pearson", "kendall", "spearman"
#--------------------------------------------------------------
# 2.a Visualizando la matriz de distancia con fviz_dist()
fviz_dist(res.dist)
fviz_dist(res.dist,
gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"))
# 2.b Visualizando la matriz de distancia con corrplot()
library(corrplot)
corrplot(as.matrix(dist.eucl),
is.corr = FALSE,
method = "color")
# 2.c Visualizando solo el triangulo superior
corrplot(as.matrix(dist.eucl), is.corr = FALSE,
method = "color",
order="hclust", type = "upper")
#--------------------------------------------------------------
# * Calculando distancias con datos mixtos con la distancia
# de Gower
library(cluster)
data(flower)
View(flower)
str(flower)
head(flower,3)
str(flower)
dd <- daisy(flower)
round(as.matrix(dd)[1:3, 1:3], 2)
# Ejemplo de distancia Gower
library(StatMatch)
datos.go <- read.csv("Ejemplo-Gower.csv",sep=";")
str(datos.go)
variables <- c("X3","X4","X5","X6")
datos.go[,variables] <- lapply(datos.go[,variables] , factor)
str(datos.go)
gower.dist(datos.go)
gower.dist(datos.go)[2,7]
##############################################
# IV. CLUSTER JERÁRQUICO AGLOMERATIVO: AGNES #
##############################################
#--------------------------------------------------------------
# Calculando la matriz de disimilaridad
# usando la distancia euclidiana
d <- dist(datos.e, method = "euclidean")
# method = euclidean, maximum, manhattan, canberra,
# binary, minkowski
as.matrix(d)[1:6, 1:6]
#--------------------------------------------------------------
# Cluster Jerarquico usando el método de enlace Ward
res.hc <- hclust(d, method = "ward.D2" )
# method = ward.D, ward.D2, single, complete, average, median
print(res.hc)
#--------------------------------------------------------------
# Visualizar el dendrograma
plot(res.hc, cex = 0.6, hang=-1)
rect.hclust(res.hc, k = 3, border = 2:5)
# Dendograma con el paquete factoextra
library(factoextra)
fviz_dend(res.hc, cex = 0.5)
# Dendrograma sin etiquetas
hcd <- as.dendrogram(res.hc)
nodePar <- list(lab.cex = 0.6, pch = c(NA, 19),
cex = 0.7, col = "blue")
plot(hcd, ylab = "Height",
nodePar = nodePar,
leaflab = "none")
# Dendrograma Horizontal
plot(hcd, xlab = "Height",
nodePar = nodePar, horiz = TRUE)
# Cambiando los colores
plot(hcd, xlab = "Height", nodePar = nodePar,
edgePar = list(col = 2:3, lwd = 2:1))
# Dendrograma Triangular
plot(hcd, type = "triangle", ylab = "Height")
#--------------------------------------------------------------
# Cophenetic Correlation Coefficient
names <- c("A","B","C","D","E","F")
x1 <- c(1,1.5,5,3,4,3)
x2 <- c(1,1.5,5,4,4,3.5)
data.cc <- data.frame(names,x1,x2)
rownames(data.cc) <- data.cc$names
data.cc$names <- NULL
dist.orig <- dist(data.cc, method = "euclidean")
dist.orig
as.matrix(dist.orig)[1:6, 1:6]
hc.single <- hclust(dist.orig, method = "single" )
plot(hc.single, cex = 0.6, hang=-1)
dist.coph <- cophenetic(hc.single)
dist.coph
# Correlacion entre la distancia cophenetic
# y la distancia original
cor(dist.orig, dist.coph)
#--------------------------------------------------------------
# Verificando el cluster tree
# Calcular la distancia cophentic
res.coph <- cophenetic(res.hc)
# Correlacion entre la distancia cophenetic distance
# y la distancia original
cor(d, res.coph)
res.hc2 <- hclust(d, method = "average")
cor(d, cophenetic(res.hc2))
#--------------------------------------------------------------
# Dividir en 3 clusters
grp <- cutree(res.hc, k = 3)
library(factoextra)
fviz_dend(res.hc, k=3, cex = 0.5,
k_colors = rainbow(3),
color_labels_by_k = TRUE,
rect=TRUE)
fviz_cluster(list(data = datos.e, cluster = grp),
palette = c("#2E9FDF", "#E7B800", "#FC4E07"),
ellipse.type = "convex", # Concentration ellipse
repel = F, # Avoid label overplotting (slow)
show.clust.cent = FALSE, ggtheme = theme_minimal())
# Junta el archivo de datos con la columna de cluster
datos.j <- cbind(datos.e,grp)
str(datos.j)
datos.j$grp <- factor(datos.j$grp)
str(datos.j)
# write.csv(datos.j,"Pancreatitis con Jerarquico Aglomerativo.csv")
#--------------------------------------------------------------
# Cluster jerarquico aglomerativo con el paquete cluster
library(cluster)
res.agnes <- agnes(x = my_data, # matriz de datos
stand = TRUE, # estandariza los datos
metric = "euclidean", # métrica de distancia
method = "ward" # método de enlace
)
fviz_dend(res.agnes, cex = 0.6, k = 3)
grp <- cutree(res.agnes, k = 3)
fviz_cluster(list(data = datos.e, cluster = grp),
palette = c("#2E9FDF", "#E7B800", "#FC4E07"),
ellipse.type = "convex", # Concentration ellipse
repel = F, # Avoid label overplotting (slow)
show.clust.cent = FALSE, ggtheme = theme_minimal())
##########################################
# V. CLUSTER JERARQUICO DIVISIVO : DIANA #
##########################################
#--------------------------------------------------------------
# Cluster jerárquico divisivo con el paquete cluster
library(cluster)
res.diana <- diana(x = my_data, # matriz de datos
stand = TRUE, # estandariza los datos
metric = "euclidean") # métrica de distancia
#--------------------------------------------------------------
# Visualizar el Dendograma
pltree(res.diana, cex = 0.6,
hang = -1,
main = "Dendrograma de Diana")
# Coeficiente divisivos
res.diana$dc
# Dendograma con la funcion plot.hclust()
plot(as.hclust(res.diana), cex = 0.6, hang = -1)
# Dendograma Horizontal con plot.dendrogram()
plot(as.dendrogram(res.diana), cex = 0.6, horiz = TRUE)
# Dendograma con factoextra
library(factoextra)
fviz_dend(res.diana, cex = 0.6, k = 3)
grp <- cutree(res.agnes, k = 3)
#####################################
# VI. CARACTERIZANDO A LOS CLUSTERS #
#####################################
# Caracterizando al cluster k-means
str(datos.j)
library(compareGroups)
group <- compareGroups(grp ~.,data=datos.j)
clustab <- createTable(group,
digits=3,
show.p.overall=FALSE)
clustab
table(datos.j$grp)
# Parallel coordinates plot using GGally
library(GGally)
#--------------------------------------------------------------
# Diagrama de lineas de promedios por cluster
library(dplyr)
datos.j %>%
group_by(grp) %>%
summarise_all(list(mean)) -> medias
medias
datos.j %>% summarise_if(is.numeric,mean) %>% round(4) -> general
general
general <- cbind(grp="general",general)
general
medias <- as.data.frame(rbind(medias,general))
medias
# Convirtiendo la data formato tidy
library(tidyr)
gathered_datos.j <- gather(data = medias,
-grp,
key = variable,
value = valor)
gathered_datos.j
# Otra forma con pivot_longer
gathered_datos.j <- pivot_longer(data = medias,
-grp,
names_to = "variable",
values_to = "valor")
#install.packages("gganimate")
library(gganimate)
ggplot(gathered_datos.j,aes(x=variable,y=valor,color=grp)) +
geom_point() + geom_line(aes(group = grp)) +
theme_bw() +
theme(legend.position = "bottom",legend.title=element_blank()) +
labs(title="Diagrama de líneas de Cluster por Variable",
x="Variable",y="") +
coord_flip() + facet_grid(grp ~.)
#+ transition_reveal(as.numeric(grp))+ enter_fade() + exit_fade()
###################################
# VII. COMPARANDO DOS DENDOGRAMAS #
###################################
library(dendextend)
#--------------------------------------------------------------
# Calculando la matriz de distancia
res.dist <- dist(datos.e, method = "euclidean")
#--------------------------------------------------------------
# Calculando 2 clusters jerárquicos aglomerativos
hc1 <- hclust(res.dist, method = "average")
hc2 <- hclust(res.dist, method = "ward.D2")
#--------------------------------------------------------------
# Creando dos Dendrogramas
dend1 <- as.dendrogram (hc1)
plot(dend1)
dend2 <- as.dendrogram (hc2)
plot(dend2)
#--------------------------------------------------------------
# Creando una lista de Dendrogramas
dend_list <- dendlist(dend1, dend2)
tanglegram(dend1, dend2)
tanglegram(dend1, dend2,
highlight_distinct_edges = FALSE, # Turn-off dashed lines
common_subtrees_color_lines = FALSE, # Turn-off line colors
common_subtrees_color_branches = TRUE, # Color common branches
main = paste("entanglement =", round(entanglement(dend_list), 2))
)
# Entanglement es una medida entre 1 (full entanglement)
# y 0 (no entanglement).
# Un bajo coeficiente de entanglement corresponde a
# un buen alineamiento
#--------------------------------------------------------------
# Verificando el cluster
# Calcular la distancia cophentic
res.coph <- cophenetic(hc1)
# Correlación entre la distancia cophenetic
# y la distancia original
cor(res.dist, res.coph) # Enlace average
cor(res.dist, cophenetic(hc2)) # Enlace ward
#--------------------------------------------------------------
# Matriz de Correlación entre una lista de dendogramas
# La función cor.denlist() es usada para calcular la matriz de
# correlación entre una lista de gráficos "Baker" o "Cophenetic".
# El valor puede caer entre -1 o 1
# Un valor cercano a 0 indica que los dos gráficos no son
# estadísticamente similares
#--------------------------------------------------------------
# Matriz de Correlación Cophenetic
cor.dendlist(dend_list, method = "cophenetic")
# Coeficiente de Correlación Cophenetic
cor_cophenetic(dend1, dend2)
#--------------------------------------------------------------
# Matriz de Correlación Baker
cor.dendlist(dend_list, method = "baker")
# Coeficiente de Correlación Baker
cor_bakers_gamma(dend1, dend2)
#--------------------------------------------------------------
# Creando múltiples dendogramas por encadenamiento
dend1 <- datos.e %>% dist("euclidean") %>% hclust("complete") %>% as.dendrogram
dend2 <- datos.e %>% dist("euclidean") %>% hclust("single") %>% as.dendrogram
dend3 <- datos.e %>% dist("euclidean") %>% hclust("average") %>% as.dendrogram
dend4 <- datos.e %>% dist("euclidean") %>% hclust("centroid") %>% as.dendrogram
#--------------------------------------------------------------
# Calcular la Matriz de Correlación
dend_list <- dendlist("Complete" = dend1, "Single" = dend2,
"Average" = dend3, "Centroid" = dend4)
cors <- cor.dendlist(dend_list)
# Imprimir la Matriz de Correlacion
round(cors, 2)
# Visualizar la Matriz de Correlacion usando el
# paquete corrplot
library(corrplot)
corrplot(cors, method="circle", type= "lower")
# method = c("circle", "square", "ellipse", "number", "shade",
# "color", "pie")
# type = c("full", "lower", "upper")
#############################
# VIII. PHYLOGENETICS TREES #
#############################
# La libreria ape (Analyses of Phylogenetics and Evolution)
# puede ser usada para producir un dendrograma más
# sofisticado.
# La función plot.phylo() puede ser usada para graficar un
# dendrograma
#--------------------------------------------------------------
# Cladogram
library(ape)
d <- dist(datos.e, method = "euclidean")
hc <- hclust(d, method = "ward.D2")
plot(as.phylo(hc), type = "cladogram", cex = 0.6,
label.offset = 0.5)
#--------------------------------------------------------------
# Unrooted
plot(as.phylo(hc), type = "unrooted", cex = 0.6,
no.margin = TRUE)
library(igraph)
fviz_dend(hc,k=3,k_colors="jco",
type="phylogenic", repel=TRUE)
#--------------------------------------------------------------
# Fan
plot(as.phylo(hc), type = "fan")
fviz_dend(hc, cex= 0.5, k=3, k_colors="jco",
type="circular")
#--------------------------------------------------------------
# Radial
plot(as.phylo(hc), type = "radial")
# Dividir el dendrograma en 3 clusters
colors = c("red", "blue", "green", "black")
clus4 = cutree(hc, 3)
plot(as.phylo(hc), type = "fan", tip.color = colors[clus4],
label.offset = 1, cex = 0.7)
# Cambiando la apariencia
plot(as.phylo(hc), type = "cladogram", cex = 0.6,
edge.color = "steelblue", edge.width = 2, edge.lty = 2,
tip.color = "steelblue")
#####################################
# IX. CLUSTER CON MÉTODOS GRÁFICOS #
# Caras de Chernoff #
#####################################
# The Use of Faces to Represent Points in K-Dimensional
# Space Graphically Herman Chernoff. Journal of the
# American Statistical Association, Vol. 68, No. 342.
# (Jun., 1973), pp. 361-368.
library(aplpack)
aplpack::faces(my_data)
aplpack::faces(my_data,face.type=0)
aplpack::faces(my_data,face.type=1)
aplpack::faces(my_data,face.type=2)
library(TeachingDemos)
TeachingDemos::faces(my_data)
library(DescTools)
DescTools::PlotFaces(my_data,scale=T)
