index.qmd

---
title: "Casos y Ejemplos de Análisis Multivariante con R"
author: "Alex Sánchez y Francesc Carmona"
date: "`r Sys.Date()`"
format:
    html: 
      toc: true
      toc-depth: 3
      code-fold: show
      code-summary: "Muestra el código"
      code-tools: true
      fig-width: 8
      fig-height: 6
    pdf:
      toc: true
      number-sections: true
      colorlinks: true
      geometry:
      - top=20mm
      - left=15mm
      papersize: A4
quarto:
    chunk_options:
      echo: true
      cache: false
      prompt: false
      tidy: true
      comment: NA
      message: false
      warning: false
    knit_options:
      width: 75
reference-location: margin
execute:
    echo: true
    message: false
    warning: false
    cache: true
# bibliography: bibliography: AMVRef.bib
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, eval=FALSE, echo=FALSE}
quarto::quarto_render("index.Rmd", output_format="all")
```


# Presentación 

Este documento contiene ejemplos de análisis multivariante con R.

Los ejemplos están pensados para ilustrar los capítulos del documento "Notas de Análisis Multivariante" pero se mantienen en un documento aparte para poder ser usados autónomamente.

Algunos de los ejemplos se han inspirado en un documento no publicado del Dr. Carles Cuadras nuestro mentor en Análisis Multivariante por lo que deseamos hacer agradecimiento explícito.

## Instalación de paquetes utilizados

```{r}
if (!require(FactoMineR)) install.packages("FactoMineR", dep=TRUE)
if (!require(factoextra)) install.packages("factoextra", dep=TRUE)
if (!require(ade4))
  install.packages("ade4", dep=TRUE)
if (!require(BiocManager))
  install.packages("BiocManager", dep=TRUE)
if (!require(made4))
  BiocManager::install("made4")
```


# Análisis de componentes principales

## Ejemplo PCA-1: Búsqueda de factores latentes en datos ecológicos[^1]

[^1]:El ejemplo de los gorriones se ha tomado del documento *Exemples d'Anàlisi Multivariant* del  Dr. Carles Cuadras disponible en su página web.

Un estudio ecológico recogió datos de 49 gorriones hembras que fueron recogidos casi moribundos después de un temporal. 
Los investigadores midieron 5 magnitudes de los animales ("length","wing","head","humerus","sternum") y anotaron si los animales sobrevivían o no a los pocos días de ser recogidos ("survival").
Entre otros objetivos, el estudio perseguía encontrar *factores*, es decir algunas variables latentes, no medibles pero presentes en los datos, que explicaran la supervivencia de los mismos ante tempestades como la registrada.
Como veremos, las dos primeras componentes principales resultan ser estos factores.

```{r}
load("datos/gorriones.RData")
colnames(gorriones) <-c("length","wing","head","humerus","sternum","survival")
gorrionesNum<-as.matrix(gorriones[,1:5])
```

Empezamos con una exploración de los datos. Aunque existen muchos paquetes que permiten hacerla nos restringiremos a las funciones de R Base (o de tidyverse) para facilitar la reproducibilidad.

```{r}
str(gorriones)
summary(gorriones)
```

```{r out.width="90%"}
f<- function(x){
  ifelse (is.numeric(x), 
          hist(x, breaks=5),
          barplot(table(x))
  )
}
par(mfrow=c(3,2))
apply(gorriones,2,f)
par(mfrow=c(1,1))
```

La escala de las variables numéricas es heterogénea, pero los datos son del mismo tipo (magnitudes biométricas) por lo que nos basaremos en la matriz de covarianzas de los datos *centrados*.

```{r}
gorrionesNum <- scale(gorrionesNum, center = TRUE, scale=FALSE)
apply(gorrionesNum,2, mean)
```

Calculamos la matriz de varianzas ajustando a dividir por n en vez de por (n-1) para que los resultados sean comparables en los dos casos que haremos.

```{r}
n<- dim(gorriones)[1]
S<-cov(gorrionesNum)*(n-1)/n
show(S)
```
Calculamos la matriz de correlaciones.

```{r}
R<-cor(gorrionesNum)
show(R)
```

A modo de ilustración empezamos calculando las componentes principales con funciones básicas:

### Mediante diagonalización de la matriz de covarianzas

```{r}
EIG <- eigen(S)
show(EIG)
```
Los vectores propios, es decir las columnas de la matriz "eigen$vectors" son las coordenadas de las cinco componentes principales.

La transformación de los datos asociada a las componentes principales se obtiene multiplicando la matriz original por la matriz de vectores propios

```{r}
eigenVecs1 <- EIG$vectors
PCAS1 <- gorrionesNum %*% eigenVecs1
head(PCAS1)
```

Con esto podemos hacer una primera representación.

```{r}
plot(PCAS1[,1], PCAS1[,2], main = "Gorriones. 2 primeras PCs")
```

Los valores propios representan la varianza de las componentes por lo que cada valor dividido por la suma de éstos es el porcentaje de variabilidad explicado.

```{r}
vars1<- EIG$values/sum(EIG$values)
round(vars1,3)
```

Podemos usar estos porcentajes para etiquetar los ejes indicando la importancia de cada componente.

```{r}
xlabel <- paste("PCA1 ", round(vars1[1]*100, 2),"%" )
ylabel <- paste("PCA2 ", round(vars1[2]*100,2),"%" )
plot(PCAS1[,1], PCAS1[,2], main = "Gorriones. 2 primeras PCs",
     xlab=xlabel, ylab=ylabel)

```

Finalmente podemos visualizar si el gorrión ha sobrevivido o no representándolo en el gráfico con colores distintos.

```{r}
bgSurv<- colSurv <- ifelse(gorriones$survival=="N", "red", "blue")
pchSurv <- ifelse(gorriones$survival=="N",1,2)
plot(PCAS1[,1], PCAS1[,2], main = "Gorriones. 2 primeras PCs",
     xlab=xlabel, ylab=ylabel, 
     pch=pchSurv, col=colSurv, bg=bgSurv)
legend( "topright"  , inset = c(0.01,0.01), cex =1, 
        bty = "y", legend = c("N", "S"), 
        text.col = c("red", "blue"),
        col = c("red", "blue"), 
        pt.bg = c("red","blue"), 
        pch = c(1,2)
        )
```

### Calculo de las componentes principales mediante las funcionesn `princomp` y `prcomp`

La función `princomp` calcula las componentes principales recurriendo básicamente al mismo método, a diferencia de la función `prcomp` que calcula las componentes principales mediante la descomposición en valores singulares de la matriz de datos.

```{r}
PCAS2 <- princomp(gorrionesNum)
names(PCAS2)
PCAS3 <- prcomp(gorrionesNum)
names(PCAS3)
```

Observemos como algunos resultados coinciden mientras que otros tendremos que ajustarlos.

EL argumento `sdev` contiene las desviaciones estándard es decir las raíces cuadradas de los valores propios.

```{r}
PCAS2$sdev
PCAS3$sdev
sqrt(EIG$values)
```

Hay una ligera diferencia en los resultados calculados por  `prcomp` que podemos atribuir a la forma de cálculo.

El argumento `loadings` contiene los vectores propios en el objeto calculado por `princomp. 

En el caso de `prcomp` estan en el objeto `rotation`.

```{r}
PCAS2$loadings
PCAS3$rotation[,1]
EIG$vectors
```

Observamos como en este caso coinciden todos los valores, salvo por el signo del primer valor propio del cálculo basado en la diagonalización. Esto no es un error sino que se debe a que el vector propio no es único y está indeterminado por un producto por -1.

Finalmente el argumento `scores` contiene las componentes principales calculadas por `princomp`. en el caso de `prcomp` se encuentran en el campo `x`. Para hacerlas comparables centraremos las componentes que hemos calculado nosotros.

```{r}
head(PCAS2$scores)
head(PCAS3$x)
head(PCAS1)
```

De nuevo los valores coinciden aunque los dos métodos `prcomp` y `princomp` han centrado las componentes principales antes de almacenarlas.

### Interpretación de las componentes

A partir de los resultados obtenidos podemos escribir así las dos primeras componentes:

$Y_1 = 0.5365\times length + 0.8290\times wing +  0.0965\times head + 0.0743\times humerus + 0.1003\times sternum$

$Y_2 = -0.8281\times length + 0.5505\times wing -0.0336\times head + 0.0146\times humerus -0.0992\times sternum$

La primera componente tiene todos los coeficientes positivos del mismo signo y se puede interpretar como el **tamaño** del gorrión. 

El peso principal de las variables originales se da en las dos primeras. 
El mayor peso se da a la extensión de las alas, con un coeficiente de 0.83. Le sigue la longitud del pájaro con un coeficiente de 0.54. Las demás variables presentan una contribución notablemente menor a la PC1. Esta primera componente explica un 82% de la variabilidad.

La segunda componente tiene coeficientes positivos y negativos y se puede interpretar como un factor **forma**. En este caso por contraposición de las dos primeras variables. En el eje PC2 las variables con mayor peso en valor absoluto son las mismas pero aquí una aparece con signo positivo y la otra con signo negativo: 0.55 para las alas, −0.83 para la longitud. Esta segunda componente explica un 11.2% de la variabilidad.

Valores positivos de PC2 corresponderán a pájaros cortos con gran envergadura de alas, y valores negativos de PC2 corresponderán a pájaros más proporcionados o incluso de mayor longitud que envergadura.

Tamaño y forma son pues dos componentes interrelacionadas que, entre ambas explican el 97.51% de la variabilidad. Si atendemos a la distribución de los "S" y "N" en el gráfico podríamos sugerir que los gorriones "extremos" ya sea en tamaño o en forma tienen menos posibilidades de sobrevivir.


## Ejemplo PCA-2: Detección del efecto batch en datos de microarrays.

Muchos datos de alto rendimiento como los microarrays o los de proteómica presentan un tipo particular de complicación técnica que se conoce como efecto \emph{batch}. Dicho efecto consiste en que muestras producidas en el mismo ``lote'' (mismo día, misma tanda de hibridación, mismo operario, ...) se parecen más entre ellas que muestras producidas en lotes distintos pudiendo llegar a causar confusión en estudios comparativos en los que el efecto batch enmascare el efecto de los tratamientos en estudio.

Si la presencia del efecto batch se conoce o espera \emph{a priori} puede, en principio, ser controlado o cancelado mediante un adecuado diseño experimental, en el que el lote se tratará usualmente como un bloque. En este caso un adecuado balanceo entre bloques y tratamientos puede cancelar el efecto batch. Alternativamente, dicho efecto puede incluirse como un  factor en el modelo de análisis de la varianza, lo que eliminara del análisis la variación atribuible al lote.

En muchas (demasiadas) ocasiones el efecto batch no ha sido considerado de antemano, o incluso, cuando sí lo ha sido, puede que no se hayan tenido en cuenta todos los posible efectos (a lo mejor se ha tenido en cuenta el día pero no que los reactivos utilizados provenían de dos lotes o que las piezas se han procesado en grupos, o ...). 

En prevención de los problemas que esto puede generar es conveniente realizar algún tipo de análisis exploratorio que permita detectar la posible presencia de estos efectos. En caso de detectarse efectos indeseados puede plantearse su eliminación mediante alguna de las técnicas disponibles para ello.

### Los datos para el análisis

Los datos para este ejemplo consisten en datos de microarrays de expresión génica, tipo hgu95-a utilizados para analizar un estudio de cáncer de mama en el que se analiza el efecto de distintos tratamientos y distintos tiempos de exposición a éstos, sobre la expresión de los genes. 

Los datos y la información del estudio se hallan disponibles en la base de datos Gene Expression Omnibus (GEO, 
[http://www.ncbi.nlm.nih.gov/geo](http://www.ncbi.nlm.nih.gov/geo) con número de acceso  ``GSE848''. 

Para este ejemplo hemos utilizado un subconjunto de estos datos basado en 18 muestras, que hemos guardado en un objeto binario `datosMicroarrays.Rda` con el fin de simplificar el proceso y centrar la atención en el análisis.

```{r}
load(file.path("datos", "datosMicroarrays.Rda"))
head(x)
rownames(targets) <- targets$ShortName
```

La tabla `targets`contiene información sobre el tratamiento y el lote de cada muestra analizada. los individuos analizados. 

```{r}
kableExtra::kable(targets)
```
Obsérvese que los nombres de las columnas de los datos coinciden con los de las filas de `targets`

```{r}
rownames(targets)==colnames(x)
```

### Análisis del efecto batch

La detección de efectos batch puede hacerse de diversas formas pero esencialmente, de lo que se trata es de ver si las muestras se agrupan por causas distintas a las que se podría esperar, por ejemplo si en vez de parecerse más entre si muestras que han recibido un mismo tratamiento, lo hacen las que han recibido algún tratamiento el mismo día.

La detección de estos efectos puede hacerse mediante distintas técnicas de las que, la más popular es el análisis de componentes principales. Una vez detectado si existe efecto batch es posible mirar de eliminarlo. Si cada tratamiento en estudio está presente en cada lote suele ser posible separar ambos efectos. Si, no es así no suele ser posible evitar cierto grado de \emph{confusión} tratamiento--lote.

En este ejercicio nos centraremos únicamente en la detección del efecto, no en su eliminación.

Cabe destacar una diferencia entre esta aplicación del PCA y la anterior. 

- En el caso de los gorriones hemos utilizado el PCA para buscar y explicar variables latentes que nos ayuden a interpretar los datos.
- En este caso no pretendemos detectar variables asociadas al efecto batch, es decir no buscaremos interpretarlos sino tan sólo ponerlo en evidencia. 

Esto no es porque no sea importante saber que causa, sino porque muchas veces puede tener orígenes desconocidos y sobretodo es plausible que no hayamos recogido información de las variables que lo explican. Esto ha llevado al desarrollo de métodos multivariantes, como el análisis de variables subrogadas del paquete SVA [https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3307112/](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3307112/) orientadas a modelizar y quitar el efecto batch globalmente sin buscar su interpretación.


### Detectando el efecto batch

La detección de efectos batch puede hacerse de diversas formas pero esencialmente, de lo que se trata es de ver *si las muestras se agrupan por causas distintas a las que se podría esperar*, por ejemplo si en vez de parecerse más entre si muestras que han recibido un mismo tratamiento, lo hacen las que han recibido cualquier tratamiento el mismo día.

La detección de estos efectos puede hacerse mediante distintas técnicas de las que, la más popular es el análisis de componentes principales. 

Una vez detectado si existe efecto batch es posible mirar de eliminarlo. Si cada tratamiento en estudio está presente en cada lote suele ser posible separar ambos efectos usando el análisis de la varianza. Si, no es así no suele ser posible evitar cierto grado de \emph{confusión} tratamiento-lote.

```{r}
colores <- c(rep("black", 4), rep("blue", 4), rep("red", 4), rep("green", 4))
```

Un diagrama de cajas de todos los datos no muestra ninguna tendencia clara que separe los lotes A y B

```{r}
sampleNames <- targets$ShortName
boxplot(as.data.frame(x), cex.axis=0.6, col=colores, las=2, names=sampleNames, 
        main="Signal distribution for selected chips")
```

Para simplificar la visualización, la función siguiente combina la realización de un PCA, la extracción de las dos primeras componentes y su visualización, debidamente etiquetada, por tratamientos y por lotes.

```{r}
plotPCA <- function ( X, labels=NULL, colors=NULL, dataDesc="", scale=FALSE, cex4text=0.8)
{
  pcX<-prcomp(X, scale=scale) # o prcomp(t(X))
  loads<- round(pcX$sdev^2/sum(pcX$sdev^2)*100,1)
  xlab<-c(paste("PC1",loads[1],"%"))
  ylab<-c(paste("PC2",loads[2],"%"))
  if (is.null(colors)) colors=1
  plot(pcX$x[,1:2],xlab=xlab,ylab=ylab, col=colors,  
       xlim=c(min(pcX$x[,1])-10, max(pcX$x[,1])+10),
       ylim=c(min(pcX$x[,2])-10, max(pcX$x[,2])+10))
  text(pcX$x[,1],pcX$x[,2], labels, pos=3, cex=cex4text)
  title(paste("Plot of first 2 PCs for expressions in", dataDesc, sep=" "), cex=0.8)
}
```

Con esta función podemos realizar y visualizar de forma sencilla el PCA de los datos. Obsérvese que realizamos el PCA de los datos traspuestos puesto que los datos de microarrays se presentan traspuestos a la forma habitual, es decir las filas son las variables y las columnas los individuos.


```{r}
sampleNames <- targets$ShortName
plotPCA(t(x), labels=sampleNames, colors=colores, dataDesc="selected samples", cex4text=0.6)
```

Las dos primeras componentes explican muy poca variabilidad, pero es imediato ver que

- La primera componente se asocia bien al factor "time" (valor 8 en general a la derechay valor 48 a la izquierda)
- La segunda componente parece asociada al factor "Batch", las muestras del Batch A, arriba y las del batch "B" abajo.

Esto sugiere que el efcto "Batch" puede enmascarar los efectos de los tratamientos por lo qyue será recomendable eliminarlo para poder estudiar bien el efecto Tratamiento.

En este caso esto será posible, al menos en parte, porque el diseño se encuentra balanceado entre Tratamiento, Tiempo y Batch, por lo que un modelo de Análisis de la Varianza nos permitirá descomponer la variabilidad asociada a cada causa y eliminar la debida al efecto batch.

### Análisis basado en distancias 

Un enfoque alternativo aunque relacionado es realizar un análisis basado en distancias. 
Podemos hacerlo calculando la matriz de distancias y visualizándola mediante un mapa de colores o usando escalamiento multidimensional que se discute más adelante.

```{r}
manDist <- dist(t(x))
heatmap (as.matrix(manDist), col=heat.colors(16))
```

```{r}
require(MASS)
sam1<-sammon (manDist, trace=FALSE)
plot(sam1$points)
text(sam1$points, targets$Batch, pos=4)
```
 
Todas las visualizaciones coinciden en mostrar una separación asociada al factor batch A o B.

# Análisis de proximidades

## Breve repaso sobre medidas de similaridad

Supongamos que tenemos $n$ variables  *dicotómicas*, basadas en ausencia (-) o presencia (+) de caracteres cualitativos. Un individuo queda entonces caracterizado por la presencia o la ausencias de los caracteres, y dos individuos seran tano más similares cuanto más parecido sea su perfil de presencias-ausencias.

|   |   |   |   |   |   |
|---|---    |---    |---    |---       |---    |
|       | $x_1$ | $x_2$ | $x_3$ | $\cdots$ | $x_n$ |
| $i$   | $+$   | $-$   | $+$   | $\cdots$ | $+$ |

Podemos determinar a asociación entre dos individuos $i$, $j$ a partir de la tabla de frecuencias

|   |   |   |   |   |
|---|---|---|---|---|
|   |   |   |  $i$  |   |
|   |   | $+$ |   | $-$  |
|   | $+$  | $a$  |  | $b$  |
|$j$  |   |   |   |
|   | $-$  | $c$  |   |  $d$  | 

donde: $n = a + b + c + d$, que contiene:

- el número de caracteres *comunes* a $i$ y a $j$,  $a$, 
- el número de caracteres *presentes en $i$ pero no en $j$*, $b$, 
- el número de caracteres *presentes en $j$ pero no en $i$*, $c$,
- el número de caracteres *no comunes* a $i$ y a $j$,  $d$,

La asociación se mide mediante un *coeficiente de similaridad* $s_{ij}$ que
verifica, en general:

- $0 \leq s_{ij} \leq 1$, 
- $s_{ij} = 0$ si $c+b = n$ **(discrepancia total)**
- $s_{ij} = 1$ si $a+d = n$ **(concordancia total)**

La similaridad $s_{ij}$ da el grado de semejanza entre $i$, $j$ en relación a los $n$ caracteres. Ejemplos de coeficientes de similaridad son

\begin{itemize}
  \item [*Sokal y Michener*:] {$\frac{a+d}{n}$}
  \item [*Sokal y Sneath*:] {$\frac{a}{a+2(b+c)}$}
  \item [*Jaccard*:] {$\frac{a}{a+b+c}$}
  \item [*Russel y Rao*:]{$\frac{a}{n}$}
\end{itemize}

## Ejemplo 1: Distancias geográficas[^2]

[^2]: El ejemplo de la visualización de distancias entre ciudades se ha tomado del libro de @Everitt2011 *An introduction to aplied multivariate analysis using R* y del repositorio de github en donde se encuentra el código R de dicho libro: [https://rdrr.io/cran/MVA/](https://rdrr.io/cran/MVA/)

En muchas ocasiones se dispone de datos en forma de una matriz de _distancias_, de _similaridades_ o de _disimilariddes_. En algún sentido, que puede variar de un caso a otro, estas matrices cuantifican el parecido, o la diferencia, o incluso la distancia geográfica por parejas de un grupo de variables. Nuestro objetivo, en estos casos, suele ser intentar recuperar, a partrir de las maytrices anteriores, las coordenadas de los puntos, de forma que podamos visualizarlos mostrando las relaciones entre ellos.

Este ejemplo utiliza las distancias *aéreas*, en millas, entre 10 ciudades de Estados Unidos. A partir de éstas deseamos recuperar las coordenadas de cada dciudad de forma que podamos representarlas en un plano que refleje lo mejor posible las distancias entre ellas. Hemos de tener en cuenta que puesto que se trata de duistancias aéreas, no son distancias euclídeas (los aviones no se desplazan en línea recta) por lo que la recuperación tan ssólo será aproximada.

```{r}
library(readxl)
DistSim <- as.data.frame(read_excel("datos/Dist_Sim_Disim.xlsx", sheet = "GeoDistances"))
rownames(DistSim) <- DistSim[,1]
DistSim <-as.matrix(DistSim [,-1])
```

```{r}
library(magrittr)
DistSim %>%
  kableExtra::kbl() %>%
  kableExtra::kable_paper("hover", full_width = F, font_size = 9)
```

Nuestro objetivo es lograr una visualización de los datos que represente lo mejor posible las distancias entre éstas, aunque no conozcamos, como en este caso, las coordenadas de cada ciudad. Además, puesto que se trata de distancias aéreas, no será posible recuperar directamente las distancias euclídeas entre las ciudades, por lo que procederemos a realizar un escalamiento multidimensional clásico

Siguiendo a @Everitt2011, capítulo 4, procederemos a:

- Diagonalizar la matriz de distancias.
- Obtener las coordenadas de cada individuo (ciudad).

$$
X = V_1\cdot \Lambda_1^{\frac 12} ,
$$
donde $V_1$ contiene los vectores propios no nulos de $B$ y 
$$
\Lambda_1^{\frac 12}= \mbox{diag}(\lambda_1^{\frac 12}, \lambda_2^{\frac 12},..., \lambda_q^{\frac 12}),
$$ 
contiene la raíz cuadrada de los primeros $q$ valores propios no nulos de $B$.

Para realizar un multidimensonal scaling clásico nos basaremos en la instrucción `cmdscale` del paquete `stats`:

```{r}
airdist <- as.dist(as.matrix(DistSim))
airline_mds <- cmdscale(airdist, k = 9, eig = TRUE)
airline_mds$points
```

Si inspeccionamos los valores propios veremos que algunos son negativos:

```{r}
lam <- airline_mds$eig
show(lam)
```

Esto se explica porque, tal como se ha comentado, las distancias no son euclídeas. Para determinar si podemos prescindir de algunas componentes calculamos, como hacíamos en el análisis de componentes principales, el porcentaje de la suma de los valores propios (en valor absoluto o al cuadrado) que representa cada valor propio individual.

```{r}
cumsum(abs(lam)) / sum(abs(lam))
cumsum(lam^2) / sum(lam^2)
```

El elevado valor del primer término en ambos casos, muestra como, bastará con las dos primeras coordenadas para obtener una buena representación.

```{r}
lim <- range(airline_mds$points[,1] * (-1)) * 1.2
plot(airline_mds$points[,1] * (-1), airline_mds$points[,2] * (-1),
     type = "n", xlab = "Coordinate 1", ylab = "Coordinate 2",
     xlim = lim, ylim = lim)
text(airline_mds$points[,1] *(-1), airline_mds$points[,2] * (-1), 
     labels(airdist), cex = 0.7)
```

Como puede verse en el gráfico la primera coordenada representa la dirección Este Oeste, y la seguna la Norte-Sur. Las posiciones de los puntos reproducen relativamente bien las posiciones relativas de las ciudades en el mapa.


## Ejemplo 2: Similaridades entre partidos políticos[^3]

[^3]:El ejemplo de la relación entre partidos políticos se ha tomado del documento @Cuadras2007, *Exemples d'Anàlisi Multivariant* del  Dr. Carles Cuadras disponible en su página web.

La tabla siguiente contiene una matriz de distancias entre partidos políticos españoles en la primera década del siglo XXI obtenida a partir de 11 características sociológicas y económicas (sobre constitución del estado, poder judicial y ejecutivo, financiación autonómica, etc.) aplicando la similaridad de Jaccard. Las distancias estan elevadas al cuadrado.


```{r}
require(readxl)
simPartidos <- as.data.frame(read_excel("datos/Dist_Sim_Disim.xlsx", sheet = "SimilarJaccard-Partidos"))
rownames(simPartidos) <- simPartidos[,1]
simPartidos <-as.matrix(simPartidos [,-1])
```

```{r}
library(magrittr)
simPartidos %>%
  kableExtra::kbl() %>%
  kableExtra::kable_paper("hover", full_width = F, font_size = 9)
```

Para trabajar con los datos los convertimos en una matriz simétrica:

```{r}
makeSymm <- function(m) {
   m[upper.tri(m)] <- t(m)[upper.tri(m)]
   return(m)
}
distPartidos <- makeSymm(as.matrix(simPartidos))
show(distPartidos)
```

En este caso tenemos una matriz de distancias $\Delta=\delta_{ij}$. Siguiendo a @Cuadras2019, cap 8, podemos transformarla en una matriz de similaridad mediante la transformación: $S_{ij} = -\frac 12 \delta_{ij}^2$.

En este caso se nos indica que las distancias ya estan al cuadrado por lo que tan sólo multiplicaremos por $-1/2$:
```{r}
simPartidos <-  S<- -0.5 *distPartidos
```

Una vez tenemos la matriz de similaridades $S$ la transformaremos en 
$$
S^*=s_{ij}^* = s_{ij} - \overline{s}_{i}- 
 \overline{s}_{j} +  \overline{s},
$$

Esto puede hacerse multiplicando por la izquierda y poir la derecha por la matriz $H$, es decir: 
$$
S*= HSH, \mbox{ donde $H$ es la matriz de centrado de datos: } H = I_n- 1_n 1_n'.
$$

Aquí:

```{r}
n = nrow(S)
ones = rep(1, n)
H = diag(ones) - (1/n) * (ones %*% t(ones))
Sprime <- H %*% S %*% H
```

A continuación diagonalizamos esta matriz:

```{r}
diagS <- eigen(Sprime)
diagS
```

Las coordenadas principales seran:

```{r}
U <- diagS$vectors
Dhalf <- sqrt(diagS$values)
coords<- U %*% diag(Dhalf)[,1:2]
```

Y la representación gráfica de los partidos según las dos primeras CP es:

```{r}
xlim <- range(coords[,1] * (-1)) * 1.2
ylim <- range(coords[,2] * (-1)) * 1.2
plot(coords[,1] * (-1), coords[,2] * (-1),
     type = "n", xlab = "Coordinate 1", ylab = "Coordinate 2",
     xlim = xlim, ylim = ylim)
text(coords[,1] * (-1), coords[,2] * (-1), 
     rownames(simPartidos), cex = 0.7)
```

Como puede verse surge una interpretación bastante natural, en tanto que el primer eje representa si el partido es (pretendidamente) de izquierdas o de derecha, ientrtas que el segundo eje alinea partidos nacionalistas en la parte superior y partidos de ámbito más estatal en la inferior. Como era de esperar, en el centro no hay nadie.

# Análisis de correspondencias

- El análisis de correspondencias es una técnica que resulta especialmente adecuada para analizar y visualizar datos de tablas de contingencia con datos de frecuencias numéricas.

- Como en el caso del PCA y el MDS este análisis proporciona una representación de datos en dimensión reducida que facilita la interpretación y la comprensión de los datos.

- En este caso, además, es posible una representación simultánea muy intuitiva de los atributos de las filas y las columnas, lo que, en las circunstancias adecuadas, permite visualizar el grado de asociación entre unas y otras.


## Ejemplo 1: Representación de mutaciones cromosómicos entre poblaciones

Este primer ejemplo es un análisis clásico basado en estudios de genética de poblaciones de los años 70. En primer lugar se presenta la realización del análisis "empaquetado" para después ilustrarl su realización paso a paso.

*@Alonso1975* hace un amplio estudio de la distribución geográfica del polimorfismo cromosómico de _Drosophila subobscura_ utilizando análisis factorial de correspondencias. 

Sobre una tabla de frecuencias de 66 poblaciones y $8 + 3 + 7 + 23 + 11$ ordenaciones de los $5$ cromosomas $A$, $J$, $E$, $U$ y $O$ respectivamente, hace un análisis de correspondencias global, y varios análisis parciales.
Utilizaremos, como ilustración uno de los análisis parciales, concretamente la que coge $13$ poblaciones y $3$ inversiones del cromosoma $A$. Los datos, se dan en forma de porcentaje en la **Tabla 1**.

```{r ejemplo1CA, echo=FALSE}
polimorph <- c("A-ST", "A-1", "A-2")
populat <- c("HELSINKI","DROBACK","HARIOT","DALKEITH","GRONINGEN","FONTAINEBLEAU", "VIENA","ZURICH","FRUSKA-GORA","LAGRASSE","MONTPELLIER","CARASCO")
datos <-  matrix(byrow=TRUE, ncol=3, data=c(96.0 , 4.0  , 0.0 ,78.4 , 16.2 , 5.4 ,100.0, 0.0  , 0.0 , 100.0, 0.0  , 0.0 ,80.0 , 16.0 , 4.0 ,88.5 , 7.7  , 3.8 ,56.9 , 35.8 , 7.4 ,67.8 , 24.4 , 7.8 ,36.1 , 55.6 , 8.3 , 72.5 , 17.5 , 10.0, 60.2 , 24.3 , 15.5, 50.0 , 31.8 , 18.2 ))
rownames(datos)=populat
colnames(datos)=polimorph
```


**TABLA**.- _Porcentajes de frecuencias de tres inversiones del cromosoma A para 13 poblaciones de Drosophila subobscura_:

```{r}
kableExtra::kable(datos)
```

Obsérvese como un gráfico sencillo puede dar buena idea de qué relaciones entre filas y columnas esperamos encontrar:

```{r}
library("gplots")
# 1. convert the data as a table
dt <- as.table(datos)
# 2. Graph
balloonplot(t(dt), main ="Polimorfismos en Europa", xlab ="Cromosoma A", ylab="Ciudades",
            label = FALSE, show.margins = FALSE)
```

Si realizamos el análisis usando la función `CA` del paquete `FactoMineR` obtenemos los siguientes resultados:

```{r}
library(FactoMineR)
datos.ca <- CA(datos, graph = FALSE)
plot.CA(datos.ca)
```

<!-- **_Figura _** - Análisis de correspondencias sobre 13 poblaciones europeas de _**Drosophila subobscura**_, en relación a tres inversiones del cromosoma A. Las inversiones están también representadas; reflejan las incidencia que tienen en cada población. -->

Observamos que _Heriot_ y _Dalkeith_ quedan representadas en un mismo punto, dada su distribución idéntica. También queda reflejada la similaridad entre _Drobak_ y _Fontainebleau_. La población _Fruska-Gora_  queda apartada del resto, debido a la influencia de la ordenación $A-1$, menos frecuente en las demás poblaciones. Las proporciones de las 3 ordenaciones quedan muy bien reflejadas en la gráfica; se ve que las poblaciones quedan más cercanas de las ordenaciones que se presentan más.

### Realización del análisis paso por paso

Aunque, en la práctica, solemos utilizar librerías de R que "empaquetan" en funciones propias las transformaciones y análisis de los datos, resulta instructivo realizar algún análisis paso a paso, ya sea para tener una mejor comprensión del análisis realizado o, incluso, para saber cuando puede ser preciso buscar alternativas porque se dan condiciones particulares.

En esta sección realizaremos paso a paso el análisis de los datos anteriores comparando los resultados que obtenemos con los proporcionados por el paquete utilizado.

Empezamos transformando los datos en frecuencias relativas y calculando las frecuencias relativas marginales.

```{r}
tabla.N <- datos
N <- sum(tabla.N)
```

El análisis de esta tabla de frecuencias mediante análisis de correspondencias puede plantearse de dos formas:

1- Como un escalado multidimensional de filas y de columnas con la distancia ji-cuadrado.
2- Como un análisis de componentes principales sobre la matriz Z estandarizada.

#### Como un escalado multidimensional

Calculamos las tablas de frecuencias relativas y marginales:

```{r}
tabla.F <-   tabla.N/N
margin.f <- margin.table(tabla.F,1)
margin.c <- margin.table(tabla.F,2)
show(round(addmargins(tabla.F),4))
```

La solución mediante escalado multidimensional pasa por calcular la matriz de distancias jui-cuadrado por filas y/o por columnas y,  sobre esta matriz, realizamos un escalado multidimensional.
Es decir realizamos dos cálculos separados que _podremos superponer después_

**Análisis de la tabla basado en los perfiles "fila"**

A partir de la tabla de frecuencias relativas podemos calcular la tabla de frecuencias relativas condicionada por filas (haciendo que éstas sumen 1).

```{r}
Y<- tabla.P <- diag(1/margin.f) %*% tabla.F
show(round(addmargins(Y),4))
```

Calculamos la matriz de distancias ji cuadrado entre las filas.
Obsérvese que para el cálculo de la distancia entre dos perfiles filas aplicamos la fórmula:

\begin{equation}
D^{2}\left(\mathbf{p}_{u}, \mathbf{p}_{v}\right)=\left(\mathbf{p}_{u}-\mathbf{p}_{v}\right)^{\prime} \mathbf{D}_{c}^{-1}\left(\mathbf{p}_{u}-\mathbf{p}_{v}\right)
\end{equation}

```{r}
nf <- nrow(tabla.F)
D2.chisq.f <- matrix(0,nf,nf)
for(i in 1:(nf-1))
  for(j in i:nf)
    D2.chisq.f[i,j] <-  t(Y[i,]-Y[j,]) %*% diag(1/margin.c) %*% (Y[i,]-Y[j,])
   D2.chisq.f <- D2.chisq.f + t(D2.chisq.f)
rownames(D2.chisq.f) <- colnames(D2.chisq.f) <- rownames(tabla.F)
show(round(D2.chisq.f,6))
```

Finalmente realizamos escalado multidimensional sobre la matriz de distancias _entre filas_. 

```{r}
mds.f <- cmdscale(sqrt(D2.chisq.f),eig=TRUE)
mds.f
```

El resultado del MDS se puede representar para visualizar la relación entre las filas.

```{r}
plot(mds.f$points,ty="n",xlab="PC1",ylab="PC2")
abline(v=0,h=0, col="gray",lty=4)
text(mds.f$points[,1],mds.f$points[,2],labels=substr(rownames(tabla.N),1,6),cex=0.6)
```

**Análisis de la tabla basado en los perfiles "columna"**
```{r}
Y<- tabla.C <-   tabla.F %*% diag(1/margin.c)
show(round(addmargins(Y),4))
```

A continuación calculamos la matriz de distancias ji cuadrado entre las columnas.

```{r}
nc <- ncol(tabla.F)
D2.chisq.c <- matrix(0,nc,nc)
for(i in 1:(nc-1))
  for(j in i:nc)
    D2.chisq.c[i,j] <-  t(Y[,i]-Y[,j]) %*% diag(1/margin.f) %*% (Y[,i]-Y[,j])
   D2.chisq.c <- D2.chisq.c + t(D2.chisq.c)
rownames(D2.chisq.c) <- colnames(D2.chisq.c) <- colnames(tabla.F)
show(round(D2.chisq.c,6))
```

Finalmente realizamos escalado multidimensional sobre la matriz de distancias _entre columnas_. 

```{r}
mds.c <- cmdscale(sqrt(D2.chisq.c),eig=TRUE)
mds.c
```

El resultado del MDS se puede representar para visualizar la relación entre las filas.

```{r}
plot(mds.c$points,ty="n",xlab="PC1",ylab="PC2")
abline(v=0,h=0, col="gray",lty=4)
text(mds.c$points[,1],mds.c$points[,2],labels=substr(colnames(tabla.F),1,6),cex=0.6)
```

Finalmente observemos que podemos superponer los dos gráficos:

```{r}
xinf<- min(min(mds.c$points[,1]), min(mds.f$points[,1]))
xsup <- max(max(mds.c$points[,1]), max(mds.f$points[,1]))
yinf<- min(min(mds.c$points[,2]), min(mds.f$points[,2]))
ysup <- max(max(mds.c$points[,2]), max(mds.f$points[,2]))
xdelta<- (xsup-xinf)/10
ydelta<- (ysup-yinf)/10
xLim <- c(xinf-xdelta, xsup+xdelta)
yLim <- c(yinf-ydelta, ysup+ydelta)

plot(mds.c$points,ty="n",xlab="PC1",ylab="PC2", xlim=xLim, ylim=yLim)
abline(v=0,h=0, col="gray",lty=4)
text(mds.f$points[,1],mds.f$points[,2],labels=substr(rownames(tabla.F),1,6),cex=0.6, col="blue")
text(mds.c$points[,1],mds.c$points[,2],labels=substr(colnames(tabla.F),1,6),cex=0.6, col="red")
```

#### AC como un ACP sobre la matriz estandarizada

El análisis se realiza sobre la matriz $\mathbf{Z}$:

\begin{equation}
\mathbf{Z}=\mathbf{D}_{f}^{1 / 2} \mathbf{Y}=\mathbf{D}_{f}^{-1 / 2} \mathbf{F} \mathbf{D}_{c}^{-1 / 2}=\left(\frac{f_{i j}}{\sqrt{f_{i} f_{\cdot j}}}\right)
\end{equation}

Sin embargo y para evitar que tengamos un valor propio 1, inútil ya que es una solución trivial,
es mejor estandarizar la matriz de correspondencias y utilizar la matriz

\begin{equation}
\mathbf{Z}=\mathbf{D}_{f}^{-1 / 2}\left(\mathbf{F}-\mathbf{f} \mathbf{c}^{\prime}\right) \mathbf{D}_{c}^{-1 / 2}
\end{equation}


```{r}
Dfmh <- diag(1/sqrt(margin.f))
Dcmh <- diag(1/sqrt(margin.c))
Z <-  Dfmh %*% (tabla.F -margin.f%*% t(margin.c)) %*% Dcmh
rownames(Z) <- rownames(tabla.F)
colnames(Z) <- colnames(tabla.F)
round(Z,4)
```

Realizaremos el ACP mediante la descomposición en valores singulares de Z:

```{r}
Z.svd <- svd(Z)
Z.svd
```

De aquí se obtienen las coordenadas principales y estándar:

```{r}
f.sc <- diag(1/sqrt(margin.f)) %*% Z.svd$u
c.sc <- diag(1/sqrt(margin.c)) %*% Z.svd$v
f.pc <- f.sc %*% diag(Z.svd$d)
c.pc <- c.sc %*% diag(Z.svd$d)
```

Podemos finalmente representar las distancias con filas y columnas de **Z**

```{r}
library(MASS)
# solución simétrica
eqscplot(f.pc[,1:2],type="n",xlab="PC1",ylab="PC2")
abline(v=0,h=0, col="gray",lty=4)
text(f.pc[,1],f.pc[,2],labels=rownames(tabla.F),cex=0.8,font=2,col="blue")
text(c.pc[,1],c.pc[,2],labels=colnames(tabla.F),cex=0.8,font=2,col= "red")
title(main="Solución simétrica",line=1)
#
# solución asimétrica
eqscplot(c.sc[,1:2],type="n",xlab="PC1",ylab="PC2")
abline(v=0,h=0, col="gray",lty=4)
text(f.pc[,1],f.pc[,2],labels=rownames(tabla.F),cex=0.8,font=2,col="blue")
text(c.sc[,1],c.sc[,2],labels=colnames(tabla.F),cex=0.8,font=2,col= "red")
title(main="Solución asimétrica",line=1)
```

## Ejemplo 2: Análisis de correspondencias de datos de microarrays

El AC se presenta tradicionalmente como una técnica para el análisis de tablas de frecuencias, pero estríctamente hablando, si podemos expresar unos datos cuantitativos como perfiles fila y perfiles columna, tambien se puede aplicar a otros tipos de datos.

Un ejemplo de esto fue la aplicación a principios de este siglo a datos de microarrays por @Fellenberg2001. El interes de esta aplicación reside sobretodo en la posibilidad de visualizar simultaneamente expresión de genes e individuos/condiciones experimentales. 

El ejemplo se lleva a cabo usando el paquete de R/Bioconductor `made4` que es una extensión, para ser usada en datos ómicos, del paquete `ade4`, un programa muy popular en francia para el análisis exploratorio de datos.

El paquete `made4` (@Culhane2005) contiene algunos _datasets_ de ejemplo. Usaremos el subconjunto de datos `data.train` del dataset `khan`. La función `overview` permite una rápida visualización de un dataset dado.

```{r}
library(made4)
data("khan")
names(khan)
k.data<-khan$train
k.class<-khan$train.classes
overview(k.data)
```

EL paquete contiene la función `ord` que permite realizar distintos métodos de reducción de la dimensión cambiando un argumento. Para análisis de correspondencias es `coa`. Para saber más puede hacerse `help (ord)`.

```{r}
k.coa<- ord(k.data, type="coa")
show(k.coa)
```

La forma más sencilla de ver los resultados producidos por `ord` es utilizar `plot`. `plot (k.ord)` dibujará un gráfico de los valores propios, junto con los gráficos de las variables (genes) y un gráfico de los casos (muestras de microarrays). En este ejemplo, las muestras de microarrays están codificadas por colores utilizando la variable `k.classes` que se guarda como `k.class`.

```{r}
plot(k.coa, classvec=k.class, genecol="grey3")
```

- La imagen muestra arriba a la izquierda (A) el  gráfico de los valores propios, lo que permite hacerse una idea de cuantas dimensiones retener.
- Arriba a la derecha (B) se muestra la proyección de las muestras de microarrays de pacientes coloreadas por tipos de tumores EWS (rojo), BL (azul), NB (verde) o RMS (marrón) 
- Abajo a la izquierda (C) se muestra la proyección de los genes (círculos grises rellenos). 
- Finalmente, abajo a la derecha (D) se muestra el biplot con los genes superpuestos con la muestras. 

Las muestras y genes con una fuerte asociado se proyectan en la misma dirección desde el origen. Cuanto mayor sea la distancia desde el origen, más fuerte será la asociación.


# Análisis de conglomerados ("Cluster Analysis")

El análisis de conglomerados pertenece al conjunto de técnicas conocidas com "análisis no supervisado" porque su objetivo es descubrir grupos en los datos, es decir grupos de individuos (o eventualmente de variables) que pueden considerarse más similares entre ellos que a los individuos de los otros grupos, según el criterio de similaridad que se decida utilizar.

Presentaremos dos ejemplos distintos:

- Un estudio sobre como se agrupan los países europeos por su consumo de carne
- Un análisis sobre la clasificación de los estados del ciclo celular y los genes que influyen en cada uno de ellos basado en un estudio de expresión génica.

## Agrupación de genes y de muestras en el análisis del ciclo celular

## Clustering y visualización

En este caso se presentan distintos métodos de agrupación (\"clustering\") disponibles en R y Bioconductor. Se utilizará el conjunto de datos <tt>yeast</tt> generado por el proyecto  \"Yeast Cell Cycle Project\" (Proyecto de estudio del ciclo celular de la levadura).

Estos datos ya no se encuentran disponibles en la web, pero se han descargado para su uso, al subidrectorio datos de este proyecto en github [https://github.com/ASPteaching/AMVCasos](https://github.com/ASPteaching/AMVCasos) 


### Pre-procesamiento de los datos

La agrupación de datos se suele aplicar después del pre-procesamiento y filtrado de los mismos. El conjunto de datos de ejemplo ya ha sido normalizado y filtrado por lo que no nos ocuparemos de este aspecto.

El archivo de datos contiene información de varios experimentos (<tt>factor alfa, cdc15 y cinéticas de decantación</tt>).
Para simplificar en este ejercicio nos ocuparemos únicamente de los datos denominados \"cdc15\" (\"cell division control protein 15\").

EL primer paso consiste en extraer los valores del grupo <tt>cdc15</tt> del conjunto de datos pre-normalizados y eliminar los genes con valores faltantes:

```{r readData, eval=TRUE}
workingDir <- getwd()
dataDir <- file.path(workingDir, "datos")
d <- read.table(file.path(dataDir, "combined.txt"), sep="\t", header=T)
names(d)
cdc15 <- which(substr(names(d),1,6)=="cdc15_")
dat <- d[,cdc15]
names(dat)
rownames(dat) <- d$X
dat <- na.omit(dat)
colnames(dat)
```

Los datos descargados de la web habían sido previamente normalizados. Para comprobarlo podemos realizar un boxplot y comprobar que los datos presentan una distribución similar y centrada en el cero como sería de esperar con valores de expresión relativa normalizados.

```{r drawboxplot}
boxplot(dat, las=2, cex.axis=0.7)
```

A continuación seleccionaremos sólo aquellos genes que se encuentran entre los que poseen el 1\% de las desviaciones estándar más altas. Esto es una práctica habitual en estudios de clasificación puesto que se considera que solo aquellos individuos con algo variabilidad resultan relevantes para la construcción o la distinción de grupos.

```{r filterData, eval=TRUE}
percentage <- c(0.975)
sds <- apply(dat, MARGIN=1, FUN="sd")
sel <- (sds>quantile(sds,percentage))
dat.sel <- dat[sel, ]
dim(dat.sel)
```

### Agrupación jerárquica de las muestras

La primera aproximación que suele hacerse para establecer la relación entre entre muestras es la aplicación de algun método de agrupamiento jerárquico. En este caso se realizará un cluster jerárquico basado en distancias euclídeas y enlaces promedio ("average linkage'').

```{r hclust, eval=TRUE}
distmeth <- c("euclidian")
Distan <- dist(t(dat.sel), method=distmeth)
treemeth <- c("average")
hc <- hclust(Distan, method=treemeth)
```

Para representar la estructura jerárquica de un agrupamiento se utiliza el dendrograma, un gráfico que tiene forma de árbol invertido, donde los nombres de los genes equivaldrían a las hojas.

```{r plotclust1,fig=T, echo=F, eval=T}
library(factoextra)
fviz_dend(hc, cex = 0.5)
```

Esta no es la única forma de realizar este tipo de agrupación y , en general se recomienda probar con distintos métodos \"razonables\" y ver hasta que punto cambian los resultados que se obtiene.

### Ejercicio

Dibujar un dendrograma utilizando el método de distancias mínimas (\"single linkage\") y el de distancias máximas (\"complete linkage\").


\paragraph{Visualización un dendrograma a partir de las correlación entre genes}

A pesar que en los ejemplos anteriores se han utilizado distancias euclídeas, los perfiles de expresión génica acostumbran a agruparse en función de los coeficientes de correlación.

<<<<<<< HEAD:index.Rmd
Para ello es necesario calcular la correlación entre los genes y cambiar el formato de la matriz de correlación para pueda contener las distancias. 

```{r}
cor.pe <- cor(as.matrix(dat.sel), method=c("pearson"))
cor.pe.rows <- cor(t(as.matrix(dat.sel)), method=c("pearson"))
cor.sp <- cor(as.matrix(dat.sel), method=c("spearman"))
dist.pe <- as.dist(1-cor.pe)
dist.pe.rows <- as.dist(1-cor.pe.rows)
dist.sp <- as.dist(1-cor.sp)
hc.cor <- hclust(dist.pe, method=treemeth)
hc.cor.rows <- hclust(dist.pe.rows, method=treemeth)
```

Para ello es necesario calcular la correlación entre los genes y cambiar el formato de la matriz de correlación para pueda contener las distancias. A continuación, el árbol puede dibujarse normalmente.

```{r plotclust2,fig=T, echo=F, eval=T}
#par(mfrow=c(2,1))
fviz_dend(hc.cor, cex = 0.5, main="Agrupacion de genes\n basada en la correlacion de Pearson")
fviz_dend(hc.cor.rows, cex = 0.5,main="Agrupacion de muestras\n basada en la correlacion de Pearson")
#par(oldpar)
```

\paragraph{Ejercicio}
Repetir el ejercicio anterior cambiando el coeficiente de Pearson por el de Spearman. ?Los resultados son diferentes?


### Agrupación de genes por el método de las <tt>k-means</tt>
Si en lugar de muestras deseamos agrupar por genes es mejor usar métodos no jerárquicos partitivos como <tt>k-means</tt>.

Un inconveniente de l'agrupación por <tt>k-means</tt> es que hace falta escoger un número el número de clusters ($k$) antes de realizar la agrupación.

Por ejemplo para producir un agrupamiento k-means con 5 grupos haremos:

```{r Kmeans1, eval=T}
k <- c(4)
km <- kmeans(dat.sel, k, iter.max=1000)
```

El método de las <tt>k-means</tt> permite estudiar la consistencia de la agrupación a partir del cálculo de las sumas de cuadrados
dentro de cada grupo (_Within SS_). El promedio de estas sumas de cuadrados para todos los grupos creados (variable  \"withinss\") es una medida
de la variabilidad  dentro de los grupos, es decir  qué tan parecidos son los genes dentro de los grupos.

```{r withinness, eval=FALSE}
mean(km$withinss)
```

Una manera de evaluar las agrupaciones generadas mediante \"k-means\" es ejecutar el mismo análisis varias
veces -- guardando cada vez el resultado como un nuevo objeto -- y seleccionar aquella agrupación que ofrece
un menor valor de  \"withinss\".


#### Dibujo de una agrupación <tt>k-means</tt>

Los resultados de los métodos divisivos como k-means<tt>k-means</tt> no pueden visualizarse mediante un dendrograma por lo que es preciso recurrir a otras representaciones gráficas.

Podemos, por ejemplo, pintar los miembros de cada grupo con un color y un símbolo distintos.

```{r KmeansAssess, fig=T, echo=F, eval=T}
cl <- km$cluster
plot(dat.sel[,1], dat.sel[,2], col=cl)
points(km$centers, col = 1:4, pch = 8)
```


Una alternativa es utilizar de nuevo el paquete factoextra 
```{r}
fviz_cluster(km, data = dat.sel,
palette = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
ellipse.type = "euclid", # Concentration ellipse
star.plot = TRUE, # Add segments from centroids to items
repel = TRUE, # Avoid label overplotting (slow)
ggtheme = theme_minimal()
)

```

Para visualizar un sólo grupo, y sus valores de expresión, seleccionaremos los genes que lo constituyen y a continuación,
dibujaremos los perfiles de expresión de estos genes utilizando diferentes colores.

```{r plotKMeans,fig=T, echo=F, eval=T}
set3 <- data.frame(dat.sel, cl)
cl1 <- dat.sel[which(cl==1),] # Here we have 3 genes
plotcol <- colorRampPalette(c("Grey", "Black"))(5)
plot(t(cl1[1,]), type="l", col=plotcol[1])
lines(t(cl1[2,]), type="l", col=plotcol[3])
lines(t(cl1[3,]), type="l", col=plotcol[5])
```

Si hay más de un grupo podemos visualizarlos todos en el mismo gráfico utilizando un peque\~no bucle <tt>for</tt>

```{r fourgroups, eval=TRUE}
km <- kmeans(dat.sel, 4, iter.max=1000)
```

```{r fourplots,fig=T, echo=F, eval=T}
par(mfrow=c(2,2))
for(i in 1:4) {
  matplot(t(dat.sel[km$cluster==i,]), type="l",
  main=paste("cluster:", i), ylab="log expression", xlab="time")
}
```

### Anotación de resultados

Una vez que se han identificado algunos grupos de interés, se puede comprobar que genes constituyen de cada uno de los \"clusters\".
Esto puede hacerse con el paquete de anotaciones para levadura de donde podremos extraer información sobre nombres, descripciones y otras anotaciones en múltiples bases de datos:

```{r yeastAnots1}
if(!require(org.Sc.sgd.db)){
 BiocManager::install("org.Sc.sgd.db")
}
require("org.Sc.sgd.db")
anotData <- capture.output(org.Sc.sgd())
print(anotData[1:15])
cat ("... output continues until ", length(anotData), " lines.\n")
```

Podemos utilizar las funciones de anotacion para convertir los identificadores propios a identificadores más generals ("Entrezs") enriqueciéndolos con información sobre dichos genes.

```{r}
genes<-rownames(dat.sel) 
columns(org.Sc.sgd.db)
geneAnots <-select(org.Sc.sgd.db, genes,c("SGD", "ENTREZID", "GENENAME"))
head(geneAnots)

###################################################
### code chunk number 15: genesInClusters
###################################################
genesInClusters <- cbind(geneAnots, cl)
kableExtra::kable(genesInClusters[order(genesInClusters$cl),])
```

## Algunas fuentes con ejemplos adicionales sobre análisis de conglomerados:

[Ejemplos de cluster analysis con R del profesor Gabriel Martos](https://rpubs.com/gabrielmartos/ClusterAnalysis)
[El paquete cluster.datasets, con los datasets de un clásico sobre análisisd e conglomerados (Hartigan)](https://cran.r-project.org/web/packages/cluster.datasets/cluster.datasets.pdf)
[Materiales de Thomas Girke, bioinformatico de la UCR, sobre Clustering](http://girke.bioinformatics.ucr.edu/GEN242/pages/mydoc/Rclustering.html)






# Referencias