library (ISLR2)Warning: package 'ISLR2' was built under R version 4.2.3nci.labs <- NCI60$labs
nci.data <- NCI60$dataTALLER DE CLUSTERING-CLASIFICACION
CLUSTERING-CLASIFICACION
Punto 2
Replique en R la sección 12.5.4
Ejemplo de Datos NCI60
Las técnicas no supervisadas se utilizan a menudo en el análisis de datos genómicos. En particular, el PCA (Análisis de Componentes Principales) y el agrupamiento jerárquico son herramientas populares. Ilustramos estas técnicas con los datos de microarreglo de líneas celulares de cáncer NCI60, que consisten en 6,830 mediciones de expresión génica en 64 líneas celulares de cáncer.
Cada línea celular está etiquetada con un tipo de cáncer, dado en nci.labs. No utilizamos los tipos de cáncer al realizar el PCA y el agrupamiento, ya que estas son técnicas no supervisadas. Pero después de realizar el PCA y el agrupamiento, verificaremos hasta qué punto estos tipos de cáncer coinciden con los resultados de estas técnicas no supervisadas. Los datos tienen 64 filas y 6,830 columnas.
dim (nci.data)[1] 64 6830Comenzamos examinando los tipos de cáncer para las líneas celulares.
nci.labs[1:4][1] "CNS" "CNS" "CNS" "RENAL"table (nci.labs)nci.labs
BREAST CNS COLON K562A-repro K562B-repro LEUKEMIA
7 5 7 1 1 6
MCF7A-repro MCF7D-repro MELANOMA NSCLC OVARIAN PROSTATE
1 1 8 9 6 2
RENAL UNKNOWN
9 1 PCA en los Datos NCI60
Primero realizamos el PCA sobre los datos después de escalar las variables (genes) para que tengan una desviación estándar de uno, aunque se podría argumentar razonablemente que es mejor no escalar los genes.
pr.out <- prcomp (nci.data, scale = TRUE)Ahora trazamos los primeros vectores de puntuación de los componentes principales, con el fin de visualizar los datos. Las observaciones (líneas celulares) correspondientes a un tipo de cáncer dado se trazan en el mismo color, para que podamos ver hasta qué punto las observaciones dentro de un tipo de cáncer son similares entre sí. Primero creamos una función simple que asigna un color distinto a cada elemento de un vector numérico. La función se utilizará para asignar un color a cada una de las 64 líneas celulares, según el tipo de cáncer al que corresponde.
Cols <- function(vec) {
cols <- rainbow(length(unique(vec)))
return(cols[as.numeric(as.factor (vec))])
}Tenga en cuenta que la función rainbow() toma como argumento un número entero positivo y devuelve un vector que contiene ese número de colores distintos. Ahora podemos trazar los vectores de puntuación de los componentes principales.
par (mfrow = c(1, 2))
plot (pr.out$x[, 1:2], col = Cols (nci.labs), pch = 19, xlab = "Z1", ylab = "Z2")
plot (pr.out$x[, c(1, 3)], col = Cols (nci.labs), pch = 19, xlab = "Z1", ylab = "Z3")
Los gráficos resultantes se muestran en la Figura 12.17. En general, las líneas celulares correspondientes a un solo tipo de cáncer tienden a tener valores similares en los primeros vectores de puntuación de los componentes principales. Esto indica que las líneas celulares del mismo tipo de cáncer tienden a tener niveles de expresión génica bastante similares.
Podemos obtener un resumen de la proporción de varianza explicada (PVE) de los primeros componentes principales utilizando el método summary() para un objeto prcomp (hemos truncado la salida):
summary(pr.out)$importance[,1:5] PC1 PC2 PC3 PC4 PC5
Standard deviation 27.85347 21.48136 19.82046 17.03256 15.97181
Proportion of Variance 0.11359 0.06756 0.05752 0.04248 0.03735
Cumulative Proportion 0.11359 0.18115 0.23867 0.28115 0.31850Usando la función plot(), también podemos trazar la varianza explicada por los primeros componentes principales.
plot(pr.out)
Tenga en cuenta que la altura de cada barra en el gráfico de barras está dada por el cuadrado del correspondiente elemento de pr.out$sdev. Sin embargo, es más informativo trazar la Proporción de Varianza Explicada (PVE) de cada componente principal (es decir, un gráfico de codo) y la PVE acumulada de cada componente principal. Esto se puede hacer con solo un poco de trabajo.
pve <- 100 * pr.out$sdev^2 / sum (pr.out$sdev^2)
par (mfrow = c(1, 2))
plot (pve, type = "o", ylab = "PVE", xlab = "Componente Principal", col = "blue")
plot (cumsum (pve), type = "o", ylab = "PVE Acumulada", xlab = "Componente Principal", col = "brown3")
(Tenga en cuenta que los elementos de pve también se pueden calcular directamente desde el resumen, summary(pr.out)$importance[2, ], y los elementos de cumsum(pve) se dan por summary(pr.out)$importance[3, ])
Vemos que, en conjunto, los primeros siete componentes principales explican alrededor del 40% de la varianza en los datos. Esto no es una gran cantidad de la varianza. Sin embargo, al observar el gráfico de codo, vemos que, aunque cada uno de los primeros siete componentes principales explica una cantidad sustancial de varianza, hay una marcada disminución en la varianza explicada por los componentes principales posteriores. Es decir, hay un “codo” en el gráfico después de aproximadamente el séptimo componente principal. Esto sugiere que podría haber poco beneficio en examinar más de siete componentes principales (aunque incluso examinar siete componentes principales podría ser difícil).
Clustering de las Observaciones de los Datos NCI60
Ahora procedemos a realizar un agrupamiento jerárquico de las líneas celulares en los datos NCI60, con el objetivo de descubrir si las observaciones se agrupan en tipos distintos de cáncer. Para comenzar, estandarizamos las variables para que tengan una media de cero y una desviación estándar de uno. Como se mencionó anteriormente, este paso es opcional y debe realizarse solo si queremos que cada gen esté en la misma escala.
sd.data <- scale (nci.data)Ahora realizamos el agrupamiento jerárquico de las observaciones utilizando los métodos de enlace completo, enlace simple y enlace promedio. La distancia Euclidiana se utiliza como medida de disimilitud.
data.dist <- dist (sd.data)
plot (hclust (data.dist), xlab = "", sub = "", ylab = "", labels = nci.labs, main = "Enlace Completo",cex = .6) 
plot (hclust (data.dist, method = "average"), labels = nci.labs, main = "Enlace Promedio", xlab = "", sub = "", ylab = "",cex = .6) 
plot (hclust (data.dist, method = "single"), labels = nci.labs, main = "Enlace Simple", xlab = "", sub = "", ylab = "",cex = .6) 
Vemos que la elección del tipo de enlace ciertamente afecta los resultados obtenidos. Típicamente, el enlace simple tiende a generar clústeres prolongados: grandes clústeres a los que las observaciones individuales se van uniendo una por una. Por otro lado, el enlace completo y el enlace promedio tienden a generar clústeres más equilibrados y atractivos. Por esta razón, el enlace completo y el enlace promedio se prefieren generalmente al enlace simple. Claramente, las líneas celulares dentro de un solo tipo de cáncer tienden a agruparse juntas, aunque el agrupamiento no es perfecto. Usaremos el agrupamiento jerárquico con enlace completo para el análisis que sigue.
Podemos cortar el dendrograma en la altura que dará un número particular de clústeres, digamos cuatro:
hc.out <- hclust (dist (sd.data))
hc.clusters <- cutree (hc.out, 4)
table (hc.clusters, nci.labs) nci.labs
hc.clusters BREAST CNS COLON K562A-repro K562B-repro LEUKEMIA MCF7A-repro
1 2 3 2 0 0 0 0
2 3 2 0 0 0 0 0
3 0 0 0 1 1 6 0
4 2 0 5 0 0 0 1
nci.labs
hc.clusters MCF7D-repro MELANOMA NSCLC OVARIAN PROSTATE RENAL UNKNOWN
1 0 8 8 6 2 8 1
2 0 0 1 0 0 1 0
3 0 0 0 0 0 0 0
4 1 0 0 0 0 0 0Existen patrones claros. Todas las líneas celulares de leucemia caen en el clúster 3, mientras que las líneas celulares de cáncer de mama se distribuyen en tres clústeres diferentes. Podemos trazar el corte en el dendrograma que produce estos cuatro clústeres:
par (mfrow = c(1, 1))
plot (hc.out, labels = nci.labs,cex = .6)
abline (h = 139, col = "red")
La función abline() dibuja una línea recta sobre cualquier gráfico existente en R. El argumento h = 139 traza una línea horizontal en la altura 139 del dendrograma; esta es la altura que da como resultado cuatro clústeres distintos. Es fácil verificar que los clústeres resultantes son los mismos que los que obtuvimos usando cutree(hc.out, 4).
Impresión de los resultados de hclust
Imprimir la salida de hclust proporciona un resumen breve pero útil del objeto:
hc.out
Call:
hclust(d = dist(sd.data))
Cluster method : complete
Distance : euclidean
Number of objects: 64 hclust(d = dist(sd.data))
Call:
hclust(d = dist(sd.data))
Cluster method : complete
Distance : euclidean
Number of objects: 64 Como mencionamos anteriormente en la Sección 12.4.2, el agrupamiento K-means y el agrupamiento jerárquico con el corte del dendrograma para obtener el mismo número de clústeres pueden dar resultados muy diferentes. ¿Cómo se comparan estos resultados de agrupamiento jerárquico NCI60 con lo que obtenemos si realizamos un agrupamiento K-means con K = 4?
set.seed (2)
km.out <- kmeans (sd.data, 4, nstart = 20)
km.clusters <- km.out$cluster
table(km.clusters, hc.clusters) hc.clusters
km.clusters 1 2 3 4
1 11 0 0 9
2 20 7 0 0
3 9 0 0 0
4 0 0 8 0Vemos que los cuatro clústeres obtenidos usando el agrupamiento jerárquico y el agrupamiento K-means son algo diferentes. El clúster 4 en el agrupamiento K-means es idéntico al clúster 3 en el agrupamiento jerárquico. Sin embargo, los otros clústeres difieren: por ejemplo, el clúster 2 en el agrupamiento K-means contiene una parte de las observaciones asignadas al clúster 1 por el agrupamiento jerárquico, así como todas las observaciones asignadas al clúster 2 por el agrupamiento jerárquico.
En lugar de realizar el agrupamiento jerárquico sobre toda la matriz de datos, podemos realizarlo simplemente sobre los primeros vectores de puntuación de los componentes principales, como sigue:
hc.out <- hclust (dist (pr.out$x[, 1:5]))
plot (hc.out, labels = nci.labs, main = "Agrupamiento Jerárquico sobre los Primeros Cinco Vectores de Puntuación",cex = .6) 
table (cutree (hc.out, 4), nci.labs) nci.labs
BREAST CNS COLON K562A-repro K562B-repro LEUKEMIA MCF7A-repro MCF7D-repro
1 0 2 7 0 0 2 0 0
2 5 3 0 0 0 0 0 0
3 0 0 0 1 1 4 0 0
4 2 0 0 0 0 0 1 1
nci.labs
MELANOMA NSCLC OVARIAN PROSTATE RENAL UNKNOWN
1 1 8 5 2 7 0
2 7 1 1 0 2 1
3 0 0 0 0 0 0
4 0 0 0 0 0 0No es sorprendente que estos resultados sean diferentes de los que obtuvimos cuando realizamos el agrupamiento jerárquico sobre el conjunto completo de datos. A veces, realizar un agrupamiento sobre los primeros vectores de puntuación de los componentes principales puede dar mejores resultados que realizar el agrupamiento sobre los datos completos. En esta situación, podríamos ver el paso de componentes principales como uno de “eliminación de ruido” de los datos. También podríamos realizar el agrupamiento K-means sobre los primeros vectores de puntuación de los componentes principales en lugar de sobre todo el conjunto de datos.
Punto 3
Solucione el ejercicio 10 del libro (página 551) usando K-means y clustering jerárquico, excluyendo del dataset 10 observaciones.
Punto 3 K-means
(a) Genera un conjunto de datos simulado con 20 observaciones en cada una de tres clases (es decir, 60 observaciones en total), y 50 variables.
# Parámetros
p <- 50
n <- 20
set.seed(123)
# Generación de las matrices y las variables aleatorias
A <- matrix(runif(p*p, -6, 6), nrow = p)
diag(A) <- abs(diag(A)) # Convertir diagonal a valores positivos
Sigma <- t(A) %*% A
A2 <- matrix(runif(p*p, -20, 20), nrow = p)
diag(A2) <- abs(diag(A2)) # Convertir diagonal a valores positivos
Sigma2 <- t(A) %*% A
X_1 <- MASS::mvrnorm(n, mu = rep(0, p), Sigma = Sigma)
X_2 <- MASS::mvrnorm(n, mu = rep(-5, p), Sigma = Sigma2)
X_3 <- MASS::mvrnorm(n, mu = rep(15, p), Sigma = Sigma)
# Crear el dataframe
df_X <- data.frame(rbind(X_1, X_2, X_3)) # Combina X_1, X_2 y X_3
df_X$source <- factor(c(rep(1, n), rep(2, n), rep(3, n))) # Indicador de origen (1, 2, 3)(b) Realiza PCA sobre las 60 observaciones y grafica los dos primeros vectores de puntuación del componente principal. Usa un color diferente para indicar las observaciones en cada una de las tres clases. Si las tres clases aparecen separadas en esta gráfica, continúa con la parte (c). Si no, regresa a la parte (a) y modifica la simulación para que haya mayor separación entre las tres clases. No continúes con la parte (c) hasta que las tres clases muestren al menos algo de separación en los dos primeros vectores de puntuación del componente principal.
pr.out <- prcomp( df_X[,-51], scale = TRUE)
par (mfrow = c(1, 2))
plot (pr.out$x[, 1:2], col = Cols (df_X$source), pch = 19, xlab = "Z1", ylab = "Z2")
plot (pr.out$x[, c(1, 3)], col = Cols (df_X$source), pch = 19, xlab = "Z1", ylab = "Z3")
(c) Realiza el clustering K-means de las observaciones con K = 3. ¿Qué tan bien se comparan los clusters obtenidos en el clustering K-means con las etiquetas de clase reales?
Pista: Puedes usar la función table() en R para comparar las etiquetas de clase reales con las etiquetas de clase obtenidas por clustering. Ten cuidado al interpretar los resultados: el clustering K-means numerará los clusters de manera arbitraria, por lo que no puedes simplemente verificar si las etiquetas de clase reales y las etiquetas de clustering son las mismas.
set.seed(123) # Establece una semilla para reproducibilidad
rows_to_remove <- sample(nrow(df_X), 10) # Selecciona 10 filas aleatorias
df_X_new <- df_X[-rows_to_remove, ] # Elimina las filas seleccionadas
km.out <- kmeans (df_X_new, 3, nstart = 20)
km.clusters <- km.out$cluster
table(km.clusters, df_X_new$source)
km.clusters 1 2 3
1 9 8 0
2 8 10 1
3 0 0 14Se puede observar que el k-means no está logrando clasificar adecuadamente los dos primeros grupos. Esto puede deberse a varias razones, siendo probablemente la principal que media de los dos primeros grupos está más cerca entre sí que del tercer grupo.
(d) Realiza el clustering K-means con K = 2. Describe tus resultados.
set.seed(123)
km.out <- kmeans (df_X_new, 2, nstart = 20)
km.clusters <- km.out$cluster
table(km.clusters, df_X_new$source)
km.clusters 1 2 3
1 1 0 15
2 16 18 0Como se observó en el punto anterior, el k-means tiende a agrupar los dos primeros grupos, ya que en este caso los ha combinado en un solo clúster, dejando el tercer grupo en un clúster aparte.
set.seed(123)
km.out <- kmeans (df_X_new, 4, nstart = 20)
km.clusters <- km.out$cluster
table(km.clusters, df_X_new$source)
km.clusters 1 2 3
1 3 6 0
2 5 4 0
3 0 0 15
4 9 8 0Se puede observar que al hacer 4 clústeres, se mantiene el mismo problema con los grupos 1 y 2. Sin embargo, se nota que el clúster 4 probablemente está agrupando los elementos de ambos grupos que están más cerca entre sí.
(f) Ahora realiza el clustering K-means con K = 3 sobre los dos primeros vectores de puntuación del componente principal, en lugar de sobre los datos crudos. Es decir, realiza el clustering K-means sobre la matriz 60 × 2, cuya primera columna es el primer vector de puntuación del componente principal, y la segunda columna es el segundo vector de puntuación del componente principal. Comenta los resultados.
set.seed (2)
pr.out <- prcomp (df_X_new[,-51], scale = TRUE)
km.out <- kmeans (pr.out$x[,1:5], 3, nstart = 20)
km.clusters <- km.out$cluster
table(km.clusters, df_X_new$source)
km.clusters 1 2 3
1 8 8 0
2 8 10 0
3 1 0 15pve <- 100 * pr.out$sdev^2 / sum (pr.out$sdev^2)
par (mfrow = c(1, 2))
plot (pve, type = "o", ylab = "PVE", xlab = "Componente Principal", col = "blue")
plot (cumsum (pve), type = "o", ylab = "PVE Acumulada", xlab = "Componente Principal", col = "brown3")
Utilizar los primeros 5 componentes parece empeorar el agrupamiento. Esto puede deberse a la baja cantidad de variabilidad que estos explican.
(g) Usando la función scale(), realiza el clustering K-means con K = 3 sobre los datos después de escalar cada variable para que tenga una desviación estándar de uno. ¿Cómo se comparan estos resultados con los obtenidos en (b)? Explica.
set.seed(123)
df_X_news<- scale (df_X_new[,-51])
km.out <- kmeans (df_X_news, 3, nstart = 20)
km.clusters <- km.out$cluster
table(km.clusters, df_X_new$source)
km.clusters 1 2 3
1 9 8 0
2 8 10 0
3 0 0 15Al realizar un escalamiento de los datos, se observa que los resultados son similares a los obtenidos en el primer k-means.
Punto 3 clustering jerárquico
Enlace Completo
data.dist <- dist(df_X_new)
plot (hclust(data.dist), xlab = "", sub = "", ylab = "", labels = df_X_new$source, main = "Enlace Completo",cex = .6)
hc.clusters<-hclust(dist(df_X_new))
hc.clusters <- cutree(hc.clusters, 3)
table(hc.clusters,df_X_new$source)
hc.clusters 1 2 3
1 8 5 0
2 3 1 13
3 6 12 2hc.clusters<-hclust(dist(pr.out$x[,1:5]))
hc.clusters <- cutree(hc.clusters, 3)
table(hc.clusters,df_X_new$source)
hc.clusters 1 2 3
1 4 6 0
2 8 10 2
3 5 2 13hc.clusters<-hclust(dist(df_X_news))
hc.clusters <- cutree(hc.clusters, 3)
table(hc.clusters,df_X_new$source)
hc.clusters 1 2 3
1 15 14 4
2 2 4 1
3 0 0 10hc.clusters<-hclust(dist(df_X_new))
hc.clusters <- cutree(hc.clusters, 2)
table(hc.clusters,df_X_new$source)
hc.clusters 1 2 3
1 11 6 13
2 6 12 2hc.clusters<-hclust(dist(pr.out$x[,1:5]))
hc.clusters <- cutree(hc.clusters, 2)
table(hc.clusters,df_X_new$source)
hc.clusters 1 2 3
1 9 8 13
2 8 10 2hc.clusters<-hclust(dist(df_X_news))
hc.clusters <- cutree(hc.clusters, 2)
table(hc.clusters,df_X_new$source)
hc.clusters 1 2 3
1 15 14 14
2 2 4 1hc.clusters<-hclust(dist(df_X_new))
hc.clusters <- cutree(hc.clusters, 4)
table(hc.clusters,df_X_new$source)
hc.clusters 1 2 3
1 8 5 0
2 3 1 13
3 1 4 0
4 5 8 2hc.clusters<-hclust(dist(pr.out$x[,1:5]))
hc.clusters <- cutree(hc.clusters, 4)
table(hc.clusters,df_X_new$source)
hc.clusters 1 2 3
1 4 6 0
2 7 8 1
3 5 2 13
4 1 2 1hc.clusters<-hclust(dist(df_X_news))
hc.clusters <- cutree(hc.clusters, 4)
table(hc.clusters,df_X_new$source)
hc.clusters 1 2 3
1 14 12 1
2 2 4 1
3 1 2 3
4 0 0 10Enlace Promedio
plot (hclust (data.dist, method = "average"), labels = df_X_new$source, main = "Enlace Promedio", xlab = "", sub = "", ylab = "",cex = .6) 
hc.clusters<-hclust(dist(df_X_new), method = "average")
hc.clusters <- cutree(hc.clusters, 3)
table(hc.clusters,df_X_new$source)
hc.clusters<-hclust(dist(pr.out$x[,1:5]), method = "average")
hc.clusters <- cutree(hc.clusters, 3)
table(hc.clusters,df_X_new$source)
hc.clusters<-hclust(dist(df_X_news), method = "average")
hc.clusters <- cutree(hc.clusters, 3)
table(hc.clusters,df_X_new$source)hc.clusters<-hclust(dist(df_X_new), method = "average")
hc.clusters <- cutree(hc.clusters, 2)
table(hc.clusters,df_X_new$source)
hc.clusters<-hclust(dist(pr.out$x[,1:5]), method = "average")
hc.clusters <- cutree(hc.clusters, 2)
table(hc.clusters,df_X_new$source)
hc.clusters<-hclust(dist(df_X_news), method = "average")
hc.clusters <- cutree(hc.clusters, 2)
table(hc.clusters,df_X_new$source)hc.clusters<-hclust(dist(df_X_new), method = "average")
hc.clusters <- cutree(hc.clusters, 4)
table(hc.clusters,df_X_new$source)
hc.clusters<-hclust(dist(pr.out$x[,1:5]), method = "average")
hc.clusters <- cutree(hc.clusters, 4)
table(hc.clusters,df_X_new$source)
hc.clusters<-hclust(dist(df_X_news), method = "average")
hc.clusters <- cutree(hc.clusters, 4)
table(hc.clusters,df_X_new$source)Enlace Simple
plot (hclust (data.dist, method = "single"), labels = df_X_new$source, main = "Enlace Simple", xlab = "", sub = "", ylab = "",cex = .6) 
hc.clusters<-hclust(dist(df_X_new), method = "single")
hc.clusters <- cutree(hc.clusters, 3)
table(hc.clusters,df_X_new$source)
hc.clusters<-hclust(dist(pr.out$x[,1:5]), method = "single")
hc.clusters <- cutree(hc.clusters, 3)
table(hc.clusters,df_X_new$source)
hc.clusters<-hclust(dist(df_X_news), method = "single")
hc.clusters <- cutree(hc.clusters, 3)
table(hc.clusters,df_X_new$source)hc.clusters<-hclust(dist(df_X_new), method = "single")
hc.clusters <- cutree(hc.clusters, 2)
table(hc.clusters,df_X_new$source)
hc.clusters<-hclust(dist(pr.out$x[,1:5]), method = "single")
hc.clusters <- cutree(hc.clusters, 2)
table(hc.clusters,df_X_new$source)
hc.clusters<-hclust(dist(df_X_news), method = "single")
hc.clusters <- cutree(hc.clusters, 2)
table(hc.clusters,df_X_new$source)hc.clusters<-hclust(dist(df_X_new), method = "single")
hc.clusters <- cutree(hc.clusters, 4)
table(hc.clusters,df_X_new$source)
hc.clusters<-hclust(dist(pr.out$x[,1:5]), method = "single")
hc.clusters <- cutree(hc.clusters, 4)
table(hc.clusters,df_X_new$source)
hc.clusters<-hclust(dist(df_X_news), method = "single")
hc.clusters <- cutree(hc.clusters, 4)
table(hc.clusters,df_X_new$source)Al realizar el agrupamiento jerárquico utilizando diferentes métodos de enlace, distintos números de clústeres y la matriz de distancias de los datos estandarizados, normales y teniendo en cuenta únicamente los 5 componentes principales, se puede observar que el mejor resultado, en términos generales, se obtiene utilizando el enlace completo con la matriz de distancias de los componentes principales.
Punto 3 Clasificación:
Se realiza LDA, QDA y KNN para predecir los puntos que se excluyeron en el análisis anterior, y se elige el mejor modelo según su rendimiento al predecir estos puntos, se ajustan varios modelos con los datos sin estandarizar, los datos estandarizados y los primeros 10 componentes principales
LDA
library (MASS)
Attaching package: 'MASS'The following object is masked from 'package:ISLR2':
Bostonlda.fit <- lda(source ~ ., data = df_X_new)Warning in lda.default(x, grouping, ...): variables are collinearlda.pred <- predict (lda.fit , df_X[rows_to_remove, ])
names(lda.pred)[1] "class" "posterior" "x" lda.class <- lda.pred$class
table (lda.class, df_X[rows_to_remove, ]$source)
lda.class 1 2 3
1 2 1 0
2 1 1 0
3 0 0 5mean (lda.class == df_X[rows_to_remove, ]$source)[1] 0.8rotation <- pr.out$x[, 1:10]
# Extraer la columna 'source' del dataframe original
source_column <- df_X_new$source
# Crear un nuevo dataframe combinando ambos
new_df <- as.data.frame(cbind(rotation, source = source_column))
lda.fit <- lda(source ~ ., data = new_df)df_subset <- df_X[rows_to_remove, -51]
# Proyectar las filas seleccionadas en los ejes principales del PCA realizado sobre df_X_new
projection <- as.matrix(df_subset) %*% pr.out$rotation
lda.pred <- predict (lda.fit , as.data.frame(projection[, 1:20]))
lda.class <- lda.pred$class
table (lda.class, df_X[rows_to_remove, ]$source)
lda.class 1 2 3
1 0 0 0
2 2 1 0
3 1 1 5mean (lda.class == df_X[rows_to_remove, ]$source)[1] 0.6QDA
qda.fit <- qda (source ~ ., data = new_df)
qda.pred <- predict (qda.fit ,as.data.frame(projection[, 1:10]))
names(qda.pred)[1] "class" "posterior"qda.class <- qda.pred$class
table (qda.class, df_X[rows_to_remove, ]$source)
qda.class 1 2 3
1 3 1 1
2 0 0 0
3 0 1 4mean (qda.class == df_X[rows_to_remove, ]$source)[1] 0.7KNN
library (class)Warning: package 'class' was built under R version 4.2.3set.seed (1)
knn.pred <- knn (df_X_new[,-51 ],df_X[rows_to_remove,-51 ] ,df_X_new$source ,k = 1)
table (knn.pred , df_X[rows_to_remove, ]$source)
knn.pred 1 2 3
1 2 1 1
2 1 1 1
3 0 0 3mean(knn.pred== df_X[rows_to_remove, ]$source)[1] 0.6library (class)
set.seed (1)
knn.pred <- knn (df_X_new[,-51 ],df_X[rows_to_remove,-51 ] ,df_X_new$source ,k = 3)
table (knn.pred , df_X[rows_to_remove, ]$source)
knn.pred 1 2 3
1 2 0 1
2 1 2 0
3 0 0 4mean(knn.pred== df_X[rows_to_remove, ]$source)[1] 0.8library (class)
set.seed (1)
knn.pred <- knn (df_X_new[,-51 ],df_X[rows_to_remove,-51 ] ,df_X_new$source ,k = 5)
table (knn.pred , df_X[rows_to_remove, ]$source)
knn.pred 1 2 3
1 1 0 0
2 2 2 0
3 0 0 5mean(knn.pred== df_X[rows_to_remove, ]$source)[1] 0.8Se selecciona KNN debido a sus buenos resultados al predecir los puntos excluidos en el análisis, destacándose frente a otros modelos en términos de precisión y rendimiento en las predicciones.
Punto 4
Solucione el ejercicio 13 del libro (página 551) usando K-means y clustering clustering jerárquico, excluyendo del dataset 10 observaciones.
Punto 4 clustering jerárquico
En el sitio web del libro, www.statlearning.com, hay un conjunto de datos de expresión génica (Ch12Ex13.csv) que consta de 40 muestras de tejidos con mediciones de 1,000 genes. Las primeras 20 muestras son de pacientes sanos, mientras que las siguientes 20 son de un grupo enfermo.
(a) Carga los datos usando read.csv(). Deberás seleccionar header = F.
library(readr)Warning: package 'readr' was built under R version 4.2.3Ch12Ex13 <- read_csv("Ch12Ex13.csv", col_names = FALSE)Rows: 1000 Columns: 40
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (40): X1, X2, X3, X4, X5, X6, X7, X8, X9, X10, X11, X12, X13, X14, X15, ...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.set.seed(123) # Establece una semilla para reproducibilidad
rows_to_remove <- sample(nrow(Ch12Ex13), 10) # Selecciona 10 filas aleatorias
Ch12Ex13_new <- Ch12Ex13[-rows_to_remove, ] # Elimina las filas seleccionadas(b) Aplica el clustering jerárquico a las muestras utilizando una distancia basada en la correlación y dibuja el dendrograma. ¿Se separan las muestras en los dos grupos según los genes? ¿Dependenden tus resultados del tipo de enlace utilizado?
# Estandarizar los datos
Ch12Ex13_scaled <- scale(Ch12Ex13_new)
# Calcular la matriz de covarianza de los datos estandarizados
covarianza <- cov(Ch12Ex13_scaled)
# Calcular la distancia de Mahalanobis para cada fila de la matriz de datos
n <- nrow(Ch12Ex13_scaled)
dist_mahalanobi <- matrix(NA, nrow = n, ncol = n)
# Calcular la distancia de Mahalanobis entre todos los puntos
for (i in 1:n) {
for (j in 1:n) {
dist_mahalanobi[i, j] <- mahalanobis(Ch12Ex13_scaled[i, ], Ch12Ex13_scaled[j, ], covarianza)
}
}d<-as.dist(dist_mahalanobi)
plot (hclust(d), xlab = "", sub = "", ylab = "", main = "Enlace Completo",labels = rep('',1000-10))
plot (hclust(d,method = "average"), xlab = "", sub = "", ylab = "", main = "Enlace Promedio",labels = rep('',1000-10))
plot (hclust((d),method = "single"), xlab = "", sub = "", labels = rep('',1000-10),ylab = "", main = "Simple")
plot (hclust((d),method = "ward.D"), xlab = "", sub = "", labels = rep('',1000-10),ylab = "", main = "Ward")
# Realizar el agrupamiento jerárquico utilizando el método de Ward
hc.clusters <- hclust(d, method = "ward.D")
# Cortar el árbol para obtener el número deseado de clústeres (por ejemplo, 2)
hc.clusters <- cutree(hc.clusters, 2)
# Agregar el resultado del clustering al dataframe original
Ch12Ex13_new <- cbind(Ch12Ex13_new, hc.clusters)Se observa que, al utilizar diferentes enlaces, los resultados varían. Sin embargo, se puede decir que el método de Ward ofrece mejores resultados, como se evidencia en el dendrograma.
k<-kmeans(Ch12Ex13_new[,-41],centers=2)
Ch12Ex13_kmeans <- cbind(Ch12Ex13_new, k$cluster )(c) Tu colaborador quiere saber qué genes presentan las mayores diferencias entre los dos grupos. Sugiere una forma de responder a esta pregunta y aplícalo aquí.
Al analizar los diferentes dendrogramas, se puede observar que los dos grupos de genes que presentan más diferencias son los que se agrupan al hacer un corte del dendrograma generado con el método de Ward. Otra forma de realizar esto es utilizando k-means, como se muestra a continuación.
Punto 4 K-means
library(dplyr)Warning: package 'dplyr' was built under R version 4.2.3
Attaching package: 'dplyr'The following object is masked from 'package:MASS':
selectThe following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, union# Para el cluster k==1
C_k1 <- Ch12Ex13_kmeans %>%
filter(`k$cluster` == 1) %>%
select(-hc.clusters, -`k$cluster`)
C_k1_mean <- colMeans(C_k1)
sum_k1 <- sum((C_k1 - C_k1_mean)^2)
# Para el cluster k==2
C_k2 <- Ch12Ex13_kmeans %>%
filter(`k$cluster` == 2) %>%
select(-hc.clusters, -`k$cluster`)
C_k2_mean <- colMeans(C_k2)
sum_k2 <- sum((C_k2 - C_k2_mean)^2)
# Cálculo de la media global
C_global_k <- colMeans(Ch12Ex13_kmeans %>% select(-hc.clusters, -`k$cluster`))
# Cálculo del error ponderado para los clusters k
error_k <- nrow(C_k1) * sum((C_k1_mean - C_global_k)^2) + nrow(C_k2) * sum((C_k2_mean - C_global_k)^2)
# Para el cluster hc.clusters==1
C_hc1 <- Ch12Ex13_kmeans %>%
filter(hc.clusters == 1) %>%
select(-hc.clusters, -`k$cluster`)
C_hc1_mean <- colMeans(C_hc1)
sum_hc1 <- sum((C_hc1 - C_hc1_mean)^2)
# Para el cluster hc.clusters==2
C_hc2 <- Ch12Ex13_kmeans %>%
filter(hc.clusters == 2) %>%
select(-hc.clusters, -`k$cluster`)
C_hc2_mean <- colMeans(C_hc2)
sum_hc2 <- sum((C_hc2 - C_hc2_mean)^2)
# Cálculo de la media global
C_global_hc <- colMeans(Ch12Ex13_kmeans %>% select(-hc.clusters, -`k$cluster`))
# Cálculo del error ponderado para los clusters hc
error_hc <- nrow(C_hc1) * sum((C_hc1_mean - C_global_hc)^2) + nrow(C_hc2) * sum((C_hc2_mean - C_global_hc)^2)
# Comparación de los resultados
comparison <- data.frame(
Cluster = c("k==1", "k==2", "hc.clusters==1", "hc.clusters==2"),
Dentro = c(sum_k1, sum_k2, sum_hc1, sum_hc2),
Entre = c(error_k, error_k, error_hc, error_hc)
)
comparison Cluster Dentro Entre
1 k==1 35539.890 7942.744
2 k==2 13438.755 7942.744
3 hc.clusters==1 39916.478 4757.354
4 hc.clusters==2 8870.233 4757.354Al realizar K-means y calcular la inercia entre y dentro de los grupos, se puede observar que este método ofrece mejores resultados. La reducción de la inercia sugiere una mejor asignación de los puntos dentro de sus respectivos clústeres, lo que indica un agrupamiento más eficiente.
Punto 4 Clasificación
Se aplican los métodos LDA, QDA y KNN para ajustar modelos de clasificación. El conjunto de datos se divide en un 80% para validación y un 20% para prueba. El modelo que obtiene el mejor resultado en el conjunto de prueba se utiliza para predecir los 10 puntos excluidos del análisis.
# Establecer semilla para reproducibilidad
set.seed(123)
# Seleccionar un 80% de los índices de las filas al azar para el conjunto de entrenamiento
train_size <- round(0.8 * nrow(Ch12Ex13_new)) # Definir el tamaño del conjunto de entrenamiento
train <- sample(1:nrow(Ch12Ex13_new), train_size) # Seleccionar aleatoriamente los índices
# Crear un vector lógico (booleano) para identificar las filas del conjunto de entrenamiento
train_bool <- seq(1, nrow(Ch12Ex13_new)) %in% trainLDA
lda.fit <- lda (hc.clusters ~ ., data = Ch12Ex13_new ,
subset = train_bool)
lda.pred <- predict (lda.fit , Ch12Ex13_new[-train, ])
names(lda.pred)[1] "class" "posterior" "x" lda.class <- lda.pred$class
table (lda.class, Ch12Ex13_new[-train, ]$hc.clusters)
lda.class 1 2
1 175 1
2 7 15mean (lda.class == Ch12Ex13_new[-train, ]$hc.clusters)[1] 0.959596QDA
qda.fit <- qda (hc.clusters ~ ., data = Ch12Ex13_new ,
subset = train)
qda.pred <- predict (qda.fit , Ch12Ex13_new[-train, ])
names(qda.pred)[1] "class" "posterior"qda.class <- qda.pred$class
table (qda.class, Ch12Ex13_new[-train, ]$hc.clusters)
qda.class 1 2
1 181 15
2 1 1mean (qda.class == Ch12Ex13_new[-train, ]$hc.clusters)[1] 0.9191919KNN
library (class)
set.seed (1)
train_df <- Ch12Ex13_new[train, !colnames(Ch12Ex13_new) %in% 'hc.clusters']
test_df<- Ch12Ex13_new[-train,!colnames(Ch12Ex13_new) %in% 'hc.clusters']
set.seed(123)
knn.pred <- knn (train_df, test_df, Ch12Ex13_new[train,]$hc.clusters, k = 4)
table (knn.pred , Ch12Ex13_new[-train, ]$hc.clusters)
knn.pred 1 2
1 178 5
2 4 11mean(knn.pred== Ch12Ex13_new[-train, ]$hc.clusters)[1] 0.9545455test_df<- Ch12Ex13[rows_to_remove,]
knn.pred <- knn (train_df, test_df, Ch12Ex13_new[train,]$hc.clusters, k = 4)
knn.pred [1] 1 1 1 2 1 1 1 1 1 1
Levels: 1 2Se puede observar que el método KNN con k=4 ofrece los mejores resultados.
Punto 2 LDA
library(MASS)
library(ISLR2)
library(dplyr)
data(Smarket)
attach(Smarket)
train <- Smarket$Year < 2005
lda.fit <- lda(Direction ~ Lag1 + Lag2, data = Smarket, subset = train)
lda.fitCall:
lda(Direction ~ Lag1 + Lag2, data = Smarket, subset = train)
Prior probabilities of groups:
Down Up
0.491984 0.508016
Group means:
Lag1 Lag2
Down 0.04279022 0.03389409
Up -0.03954635 -0.03132544
Coefficients of linear discriminants:
LD1
Lag1 -0.6420190
Lag2 -0.5135293Ajustamos un modelo de análisis de discriminante lineal, obteniendo que las probabilidades a priori para cada uno de los grupos (alza y baja) son parecidas, teniendo una probabilidad más alta el alza de retorno. Estas probabilidades fueron tomadas como la proporción de registros dentro de cada grupo.
Las medias de los grupos sugieren que hay una tendencia a la baja del mercado cuando los dos días anteriores tuvieron un alza, y viceversa.
Los coeficientes del discriminante lineal nos permiten construir la regla de decisión para clasificación. Luego, nuestra regla de clasificación será \[-0.642L_1 - 0.513L_2\] La cual en caso de ser mayor a cero, asignará la observación \(x_0\) al primer grupo, es decir, predecirá que ese día irá en alza el mercado, y viceversa.
Observando ahora las predicciones para los datos de 2005 en adelante, tenemos que
Smarket.2005 <- Smarket[!train, ]
Direction.2005 <- Smarket$Direction[!train]
lda.pred <- predict(lda.fit, Smarket.2005)
table(lda.pred$class, Direction.2005) Direction.2005
Down Up
Down 35 35
Up 76 106mean(lda.pred$class == Direction.2005)[1] 0.5595238Donde observamos que clasificó de manera incorrecta 111 observaciones, lo que le da al modelo una presición aproximada del 56 %. Aplicando una regla de clasificación del 50 % sobre las probabilidades posteriores, podemos hallar de manera manual los resultados del objeto “class” resultante del ajuste del modelo.
sum(lda.pred$posterior[, 1] >= 0.5)[1] 70sum(lda.pred$posterior[, 1] < 0.5)[1] 182sum(lda.pred$posterior[, 1] > 0.9)[1] 0lda.pred$posterior[1:20, 1] 999 1000 1001 1002 1003 1004 1005 1006
0.4901792 0.4792185 0.4668185 0.4740011 0.4927877 0.4938562 0.4951016 0.4872861
1007 1008 1009 1010 1011 1012 1013 1014
0.4907013 0.4844026 0.4906963 0.5119988 0.4895152 0.4706761 0.4744593 0.4799583
1015 1016 1017 1018
0.4935775 0.5030894 0.4978806 0.4886331 lda.pred$class[1:20] [1] Up Up Up Up Up Up Up Up Up Up Up Down Up Up Up
[16] Up Up Down Up Up
Levels: Down UpObservamos que las únicas dos observaciones que tienen probabilidades posteriores mayores a 0.5 (1010 y 1016), son las que tienen la predicción de baja dentro del objeto class. \\ Por último, si qusieramos ser muy rigurosos con las predicciones, y clasificar una observación únicamente cuando se tenga una certeza muy grande de que será baja o alza, digamos un 90 %, podemos observar que no hay ninguna observación en 2005 que tenga una probabilidad tan grande, y en particular, la probbailidad más grande que se presenta dentro de este conjunto de prueba es de 54 %.
Para comprobar los resultados obtenidos mediante el comando LDA, realizaremos el cálculo de los coeficientes del discriminante lineal de manera manual. Como habíamos visto en la sección anterior, estos coeficientes corresponden a la combinación lineal dada por el vector \(a^T = (\mu_1 - \mu_2)^T\Sigma^{-1}\), donde estimamos \(\mu_1\) y \(\mu_2\) como los promedios muestrales de los grupos de alza y baja en los datos de entrenamiento, y \(\Sigma\) como la matriz \(S_{pool}\).
ent1 <- Smarket %>%
filter((Year < 2005) & Direction == "Up") %>%
as.data.frame()
ent2 <- Smarket %>%
filter((Year < 2005) & Direction == "Down") %>%
as.data.frame()
# Crear matrices de variables predictoras
x1 <- data.frame(l1 = ent1$Lag1, l2 = ent1$Lag2)
x2 <- data.frame(l1 = ent2$Lag1, l2 = ent2$Lag2)
n1 <- nrow(ent1)
n2 <- nrow(ent2)
spool <- ((n1 - 1) * cov(x1) + (n2 - 1) * cov(x2)) / (n1 + n2 - 2)
m1 <- colMeans(Smarket[train & Direction == "Down", c("Lag1", "Lag2")])
m2 <- colMeans(Smarket[train & Direction == "Up", c("Lag1", "Lag2")])
a <- solve(spool)%*%(m2-m1)
a [,1]
l1 -0.05544078
l2 -0.04434520donde observamos que los coeficientes parecen ser distintos a los obtenidos mediante la función, pero esto se debe a que la función reescala estos coeficientes por un factor, lo que podemos corroborar observando que los coeficientes de la salida de la función, son proporcionales a los obtenidos mediante cálculo manual.
c <- -0.6420190/a[1]
c[1] 11.58027c*a [,1]
l1 -0.6420190
l2 -0.5135292Ajustamos ahora un modelo de análisis de discrimnante cuadrático. Recordemos que la diferencia entre este modelo y el anterior es que para este modelo no se asume la igualdad de estructura de covarianza entre las poblaciones.
qda.fit <- qda(Direction ~ Lag1 + Lag2 , data = Smarket ,
subset = train)
qda.fitCall:
qda(Direction ~ Lag1 + Lag2, data = Smarket, subset = train)
Prior probabilities of groups:
Down Up
0.491984 0.508016
Group means:
Lag1 Lag2
Down 0.04279022 0.03389409
Up -0.03954635 -0.03132544En este caso, la función no proporciona coeficientes de discriminante, pues el vector de coeficientes \((\mu_1 - \mu_2)\Sigma^{-1}\) no es único en este caso, al tener dos matrices de covarianza distintas.
qda.class <- predict(qda.fit , Smarket.2005)$class
table(qda.class , Direction.2005) Direction.2005
qda.class Down Up
Down 30 20
Up 81 121mean(qda.class == Direction.2005)[1] 0.5992063Observando los resultados de este modelo, concluímos que se ajusta de mejor manera que el modelo de discriminante lineal, teniendo casi un 60% de precisión.
Ajustamos ahora un modelo de KNN, el cual se basa en encontrar las k observaciones más cercanas a la observación a clasificar, y asignar esta observación al grupo del que provengan la mayoría de ellas. En caso de existir empates en términos de distancias, R toma aletoriamente una de las observaciones, por lo que también se debe fijar una semilla para reproducibilidad de los resultados.
library(class)
train.X <- cbind(Lag1 , Lag2)[train , ]
test.X <- cbind(Lag1 , Lag2)[!train , ]
train.Direction <- Direction[train]
set.seed (1)
knn.pred <- knn(train.X, test.X, train.Direction , k = 1)
table(knn.pred , Direction.2005) Direction.2005
knn.pred Down Up
Down 43 58
Up 68 83mean(knn.pred == Direction.2005)[1] 0.5Claramente los resultados mejoraron, pero siguen siendo inferiores a los análisis de discrimante, especialmente, el QDA.
Caravan Data
Trabajaremos ahora con los datos de Caravan, el cual cuenta con 85 variables predictoras medidas en 5822 individuos, y en donde nuestra variable respuesta corresponde a la compra de un seguro para caravanas. En este conjunto de datos, únicamente el 6 % de los individuos adquirieron este seguro, por lo que nos encontramos con un conjunto de datos desbalanceado.
data("Caravan")
attach(Caravan)
summary(Purchase) No Yes
5474 348 348/nrow(Caravan)[1] 0.05977327Para el procedimiento de KNN, las unidades de medida o escala son extremadamente importantes, pues variables que tengan escalas mayores incurrirán en medidas de disimilaridad mayores entre observaciones, por lo que se debe en este caso estandarizar el conjunto de variables predictoras para favorecer la comparabilidad y no afectar los resultados.
standardized.X <- scale(Caravan[, -86])
var(Caravan[, 1])[1] 165.0378var(Caravan[, 2])[1] 0.1647078var(standardized.X[, 1])[1] 1var(standardized.X[, 2])[1] 1test <- 1:1000
train.X <- standardized.X[-test, ]
test.X <- standardized.X[test, ]
train.Y <- Purchase[-test]
test.Y <- Purchase[test]
set.seed (1)
knn.pred <- knn(train.X, test.X, train.Y, k = 1)
mean(test.Y != knn.pred)[1] 0.118mean(test.Y != "No")[1] 0.059table(knn.pred , test.Y) test.Y
knn.pred No Yes
No 873 50
Yes 68 9donde se observa que el modelo predijo incorrectamente solo en 12% de las observaciones de prueba, por lo que podría pensarse que es un muy buen modelo. Sin embargo, en este caso debemos tener en cuenta el desbalance de los datos, pues es lógico que la mayoría de las observaciones sean predichas como “No” y sean correctas, pues el 94% de los datos originales tenían esta etiqueta. Además, otro factor importante a tener en cuenta es el interés que se tiene en los resultados del modelo. En este caso no tiene sentido medir la precisión general del modelo, pues el interés se basa únicamente en estudiar a los individuos que sí adquirieron el seguro de caravanas, por lo que nos interesa medir la presición sobre los individuos que el modelo predijo que sí comprarían este seguro. \ Para el caso del modelo para k=1, se tuvo una precisión general del 88 %, pero un precisión para los individuos clasificados como “Sí” del 12%. \\ Ajustamos ahora un modelo con k=3
knn.pred <- knn(train.X, test.X, train.Y, k = 3)
table(knn.pred , test.Y) test.Y
knn.pred No Yes
No 920 54
Yes 21 55 / 26[1] 0.1923077Observando que la precisión de interés aumenta a casi un 20%, y para un modelo con k=5
knn.pred <- knn(train.X, test.X, train.Y, k = 5)
table(knn.pred , test.Y) test.Y
knn.pred No Yes
No 930 55
Yes 11 44/15[1] 0.2666667Obtenemos casi un 27%. \\ Sin embargo, aunque esta estrategia es costo-efectiva, en la práctica, la empresa querría detectar más de 15 clientes potenciales para la adquisición del seguro, por lo que ajustaremos ahora a modo de comparativa, modelos de regresión logística.
glm.fits <- glm(Purchase ~ ., data = Caravan ,
family = binomial , subset = -test)Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredglm.probs <- predict(glm.fits , Caravan[test , ],
type = "response")
glm.pred <- rep("No", 1000)
glm.pred[glm.probs > .5] <- "Yes"
table(glm.pred , test.Y) test.Y
glm.pred No Yes
No 934 59
Yes 7 0Observamos que con un punto de corte de probabilidad del 50%, el modelo predice que únicamente 7 personas comprarán el seguro, y aún peor, los predice de manera incorrecta. La solución a esto, y al problema de querer detectar un mayor número de potenciales compradores, es hacer el modelo más flexible, fijando un punto de corte más pequeño, digamos del 25%.
glm.pred <- rep("No", 1000)
glm.pred[glm.probs > .25] <- "Yes"
table(glm.pred , test.Y) test.Y
glm.pred No Yes
No 919 48
Yes 22 1111/33[1] 0.3333333Para este caso, podemos observar que fijando un punto de corte más bajo, se detectan 33 individuos como potenciales compradores, y el 33% de ellos si están predichos correctamente, por lo que este enfoque resulta más beneficioso que cualquier enfoque de KNN de los que ajustamos.