---
title: "Análisis cluster"  
author: "Guillermo Ayala Gallego"  
date: "`r Sys.Date()`"  
address: Departamento de Estadística e I.O. Universidad de Valencia  
format: 
   revealjs:  
      scrollable: true
      echo: true
      embed-resources: true  
   pdf:
       toc: true
       number-sections: false
       colorlinks: true
---

## Un ejemplo artificial: datos Ruspini

```{r}
#| label: cluster1
library(cluster)
data(ruspini)
```

## ...  
```{r}
#| label: cluster2
#| echo: false  
plot(ruspini)
```

## Un ejemplo con muestras

```{r}
#| label: cluster3
library(multtest)
data(golub)
grep("CCND3 Cyclin D3",golub.gnames[,2])
grep("Zyxin",golub.gnames[,2])
cz.data = data.frame(golub[1042,],golub[2124,])
colnames(cz.data) = c("CCND3 Cyclin D3","Zyxin")
```

```{r}
#| label: cluster4
#| echo: false  
plot(cz.data)
```

# Disimilaridades

## Distancias 

- Distancia euclídea

$$
d(x,y) = \sqrt{\sum_{k=1}^d (x_k -y_k)^2}
$$

- Distancia de Manhattan 
$$
d(x,y) = \sum_{k=1}^d |x_k - y_k|
$$


## Correlación de Pearson

- Supongamos que cuando evaluamos si dos genes están **próximos** 
interpretamos que están corregulados.
- Una medida de **co**rrelación es 
$$
r = \frac{\sum_{i=1}^n (x_i - \bar{x}_n) (y_i - \bar{y}_n)}
{\sqrt{\sum_{i=1}^n (x_i - \bar{x}_n)^2} \sqrt{\sum_{i=1}^n
    (y_i - \bar{y}_n)^2}}.
$$

## ...  

```{r}
#| label: cluster5
#| echo: false
x1 = rnorm(135,mean=2.34,sd = 1.34)
y1 = x1 + rnorm(135,mean=0,sd = 1)
plot(x1,y1,xlab="x",ylab="y")
```

## ...  

```{r}
#| label: cluster6
#| echo: false
plot(x1,y1,xlim=c(min(x1),max(x1)),ylim=c(min(y1),max(y1)))
abline(v=mean(x1),col="red",lty=2)
abline(h=mean(y1),col="red",lty=2)
```

## 

```{r}
#| label: cluster7
#| echo: false
crossproduct = (x1-mean(x1)) * (y1 - mean(y1))
plot(x1[crossproduct >=0],y1[crossproduct >=0],
     xlim=c(min(x1),max(x1)),ylim=c(min(y1),max(y1)),col="red",
     xlab="x",ylab="y")
points(x1[crossproduct <0],y1[crossproduct <0],col="blue",pch=2)
abline(v=mean(x1))
abline(h=mean(y1))
```

```{r}
#| label: cluster8
cor(x1,y1)
```

## Una correlación mayor

```{r}
#| label: cluster9
#| echo: false
x1 = rnorm(135,mean=2.34,sd = 1.34)
y1 = 3 * x1 + rnorm(135,mean=0,sd = .56)
plot(x1,y1,xlab="x",ylab="y")
```

## 

```{r}
#| label: cluster10
#| echo: false
crossproduct = (x1-mean(x1)) * (y1 - mean(y1))
plot(x1[crossproduct >=0],y1[crossproduct >=0],
     xlim=c(min(x1),max(x1)),ylim=c(min(y1),max(y1)),col="red",ylab="y",
     xlab="x")
points(x1[crossproduct <0],y1[crossproduct <0],col="blue",pch=2)
abline(v=mean(x1),lty=2)
abline(h=mean(y1),lty=2)
```

Y tenemos una correlación observada de
```{r}
#| label: cluster11
cor(x1,y1)
```

## Correlación muy pequeña

```{r}
#| label: cluster12
#| echo: false
x1 = rnorm(135,mean=2.34,sd = 1.34)
y1 = 0.003*x1 + rnorm(135,mean=0,sd = .56)
plot(x1,y1)
```

## ...  

```{r}
#| label: cluster13
#| echo: false
crossproduct = (x1-mean(x1)) * (y1 - mean(y1))
plot(x1[crossproduct >=0],y1[crossproduct >=0],
     xlim=c(min(x1),max(x1)),ylim=c(min(y1),max(y1)),col="red",
     xlab="x",ylab="y")
points(x1[crossproduct <0],y1[crossproduct <0],col="blue",pch=2)
abline(v=mean(x1),lty=2)
abline(h=mean(y1),lty=2)
```

Y la correlación muestral es

```{r}
#| label: cluster14
cor(x1,y1)
```


## Disimilaridades utilizando coeficientes de correlación

- $d(x,y) = 1 - r_{x,y}$,
- $d(x,y) = (1 - r_{x,y})/2$,
- $d(x,y) = 1 - |r_{x,y}|$,
- $d(x,y) = \sqrt{1 - r_{x,y}^2}$.


## Evaluando las disimilaridades

- Si la disimilaridad que estamos considerando
es simétrica entonces 
$$
d(x_i,x_j) = d(x_j,x_i)
$$
- D es una matriz simétrica
- Además $d(x_i, x_i) = 0$.
- Se almacena una parte de la matriz (superior o inferior).

## Calculando disimilaridades

- Distancia euclídea
```{r}
#| label: cluster15
cz.dist.euclidea = dist(cz.data,method="euclidian")
class(cz.dist.euclidea)
```
- Distancia de Manhattan

```{r}
#| label: cluster16
cz.dist.manhattan = dist(cz.data,method="manhattan")
```
- Basadas en correlación 
```{r}
#| label: cluster17
data(golub,package="multtest")
golub.cor1 = 1-cor(t(golub))
golub.cor2 = (1-cor(t(golub)))/2
golub.cor3 = 1-abs(cor(t(golub)))
golub.cor4 = sqrt(1-cor(t(golub))^2)
```

## Disimilaridades entre grupos de observaciones

- Supongamos que tenemos un banco de datos con $n$ individuos
cuyos índices son $\{1, \ldots, n\}$.
- Sean $A$ y $B$ dos subconjuntos disjuntos del conjunto
de índices de la muestra $\{1, \ldots, n\}$
- Si denotamos la disimilaridad entre $A$ y $B$ como $d(A,B)$ entonces
    + **Enlace simple**: $d(A,B) = \min_{a \in A, b \in B} d(a,b)$.
    + **Enlace completo**: $d(A,B) = \max_{a \in A, b \in B} d(a,b)$.
    + **Promedio**: $d(A,B) = \frac{1}{|A| \times |B|} \sum_{a \in A, b \in B} d(a,b)$ siendo $|A|$ es el cardinal del conjunto $A$.


## ¿Qué es una clasificación?
-  Sea $\{1, \ldots, n\}$ el conjunto de índices que indexan las
distintas observaciones.
- Supongamos que $\{C_1, \ldots, C_r\}$ es una partición de este
conjunto de índices: 
    1. $C_i \subset \{1, \ldots, n\}$;  son disjuntos dos a dos, 
    $C_i \cap C_j= \emptyset$ si $i \neq j$ con $i,j = 1, \ldots, n$ y 
    2. $\cup_{i=}^r C_i = \{1, \ldots, n\}$. 

## Cluster jerárquico aglomerativo
1. Tenemos grupos unitarios formados por cada una
de las observaciones. Tenemos pues una partición inicial $C_i=\{i\}$
con $i=1,\ldots,n$. En un principio, cada dato es un grupo.
2. Calculamos las disimilaridades entre los elementos de la
partición. Para ello utilizamos cualquiera de los procedimientos
antes indicados.
3. Agrupamos los dos conjuntos de la partición más próximos
y dejamos los demás conjuntos igual. Tenemos ahora $C_i$ con
$i=1,\ldots,k$.
4.  Si tenemos un solo conjunto en la partición paramos el
procedimiento.
5. Volvemos al paso **2**.

## agnes: dendograma

Dendograma para un cluster jerárquico aplicado a los datos
ruspini con métrica euclídea y el promedio para la disimilaridad
entre grupos.

```{r}
#| label: cluster18
ruspini.ag = agnes(ruspini,metric = "euclidean",method="average")
```
```{r}
#| label: cluster19
#| eval: false
ruspini.ag = agnes(ruspini,metric = "euclidean",method="average")
```

## ...  

```{r}
#| label: cluster20
#| echo: false
plot(ruspini.ag,which = 2)
```
## Clasificación
- Supongamos que decidimos quedarnos con cuatro grupos. 
- Las clasificaciones de los datos son las siguientes donde la etiqueta
que se asigna a cada grupo es arbitraria. 

```{r}
#| label: cluster21
cutree(ruspini.ag,4)
```


## ...  

Clasificación obtenida para los datos ruspini utilizando
un cluster jerárquico con distancia euclídea,  promedio entre grupos
y cuatro grupos.

```{r}
#| label: cluster22
plot(ruspini,col=cutree(ruspini.ag,4),pch=cutree(ruspini.ag,4))
```

# Métodos de particionamiento

## Características

- Suponemos ahora que tenemos una idea de cuántos grupos hay.
- Posiblemente hemos realizado un análisis jerárquico previo con todos
los datos o, si eran muchos, con una selección aleatoria de los
datos.
-  En principio, vamos a suponer que fijamos el número de grupos a considerar. 
- Suponemos pues que **sabemos** el número de grupos y lo denotamos por $k$.

## Método de las k-medias 

- Supongamos que tenemos $C_1, \ldots, C_k$ una partición de $\{1, \ldots, n\}$.
- Un modo bastante natural de valorar la calidad  del agrupamiento que la 
partición nos indica sería la siguiente función:
$$
\sum_{i=1}^k \sum_{j \in C_i} d_E(x_j,\bar{x}_{C_i})^2,
$$
donde $d_E$ denota aquí la distancia euclídea y
$$
\bar{x}_{C_i} = \frac{1}{|C_i|} \sum_{j \in C_i} x_j,
$$
es el vector de medias del grupo cuyos índices están en $C_i$. 
- Una partición será tanto mejor cuanto menor sea el valor de la función anterior.
- El procedimiento de agrupamiento de las k-medias simplemente se basa en 
elegir como partición de los datos aquella que nos da el **mínimo**
 de la función objetivo.
 
 
## kmeans

```{r}
#| label: cluster23
ruspini.km = kmeans(ruspini,4)
```
La clasificación viene dada por 

```{r}
#| label: cluster24
ruspini.km$cluster
```


## Particionamiento alrededor de los mediodes (PAM)
- Supongamos que tomamos $k$ individuos de la muestra que denotamos
por $m_i$ con $i=1,\ldots,k$. 
- Particionamos la muestra en $k$ grupos de modo que el grupo $C_i$ 
está formado por los individuos más próximos a $m_i$ que a cualquier otro 
$m_j$ con $j \ne i$,
$$
C_i = \{l: d(l,i) = \min_{j \ne i} d(l,j) \}.
$$
- Consideremos la siguiente cantidad:
$$
\sum_{i=1}^k \sum_{j \in C_i} d(j,m_i).
$$
- En el método de particionamiento alrededor de los mediodes nos
planteamos encontrar las observaciones $m_1,\ldots,m_k$ que
minimizan el valor anterior.

## pam

```{r}
#| label: cluster25
ruspini.pam = pam(ruspini,4)
```

## ...  

```{r}
#| label: cluster26
plot(ruspini.pam,which=1)
```

## Silueta  

- Para la observación $i$ y el grupo $C$ consideramos 
$$
\bar{d}(i,C) = \frac{1}{|C|} \sum_{j \in C} d(i,j).
$$ 
- Para cada observación $i$, sea $A$ su cluster
 y $a(i) = \bar{d}(i,A)$. 
- Consideremos 
$$
b(i) = \min_{C \ne A} \bar{d}(i,C).
$$

## 

- Definimos
$$s(i)  =  1 - \frac{a(i)}{b(i)} \text{ si} a(i) < b(i),$$ 
$$s(i)  =  0 \text{ si } a(i) = b(i),$$ 
$$s(i)  =  \frac{b(i)}{a(i)} - 1  \text{ si } a(i) > b(i).$$

- O equivalentemente
$$s(i) = \frac{b(i) - a(i)}{\max\{a(i),b(i)\}}$$



## Anchura media de la silueta  

¿Cómo interpretar $\bar{s} = \sum_{i=1}^n s_i /n$?

------------     ---------------------------
$0.71-1.00$      Fuerte estructura 

$0.51-0.70$      Estructura razonable 

$0.26-0.50$      Estructura débil. Probar otros métodos 

$\leq 0.25$      No se encuentra estructura 

----------------------------------------------


## ...  

```{r}
#| label: cluster27
ruspini.pam = pam(ruspini,4)
summary(silhouette(ruspini.pam))
```

## ...  

```{r}
#| label: cluster28
plot(silhouette(ruspini.pam))
```

## Análisis cluster y datos de expresión de gen

- Hemos visto algunos métodos de clasificación
así como de valoración de hay estructura cluster en un conjunto de
datos dados.
- Estos métodos se aplican (sin modificación) en contextos muy diversos.
- ¿Qué interés particular tiene su aplicación para datos de expresión de gen? 
- Lo primero a considerar es que nos podemos plantear o bien la clasificación
de las filas de la matriz de expresión o bien la clasificación de las
columnas.
- Si clasificamos las filas lo que tenemos en cada fila 
es el perfil de expresión del gen en todas las muestras del estudio.
- Si clasificamos las columnas, las muestras, entonces podemos
estar, por ejemplo, buscando subclases de células tumorales.
- Los resultados de la clasificación de genes y muestras puede utilizarse
como control de calidad. Sobre todo en las muestras tenemos un diseño
experimental previo (habitualmente, el diseño grupo control frente a 
grupo de tratamiento).


## GSE20986

- Cargamos los datos.

```{r}
#| label: cluster29
data(gse20986,package="tamidata")
```

- Seleccionamos genes.

```{r}
#| label: cluster30
library(genefilter)
gse0 = nsFilter(gse20986,var.func = mean, var.cutoff = 0.7)
gse1 = nsFilter(gse0$eset,var.func = IQR, var.cutoff = 0.7)
gse = gse1$eset
```
- ¿Con cuántas filas nos hemos quedado?

```{r}
#| label: cluster31
dim(exprs(gse))
```


## ...  

- Vamos a clasificar las muestras y compararemos con la clasificación 
original de las mismas. 
- Utilizamos distancia Manhattan y promedio con el método PAM.

- Primero calculamos la matriz de distancias.

```{r}
#| label: cluster32
d0 = dist(t(exprs(gse)),method = "manhattan")
```
- Y aplicamos PAM.

```{r}
#| label: cluster33
library(cluster)
gse.pam = pam(d0,diss = TRUE,k=4)
```
- La clasificación original es

```{r}
#| label: cluster34
tipo0 = rep(1:4,rep(3,4))
```
- Comparamos la clasificación original con la obtenida.

```{r}
#| label: cluster35
table(tipo0,gse.pam$clustering)
```

- ¿Qué significa esto?

## Heatmap  

```{r}
#| label: cluster36
heatmap(exprs(gse),hclustfun=agnes)
```
