---
title: "Componentes principales"
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
---  


# Componentes principales

## El problema

- En este tema nos ocupamos de problemas de **reducción de dimensión**.
- ¿Qué significa reducir la dimensión? 
- Trabajando con expresión de genes tenemos tantas filas como
genes y tantas columnas como muestras.
- En resumen miles de filas y decenas o centenares de columnas.
-  En este tema tratamos el tema de cómo reducir el número de columnas o el
número de filas.

## Componentes principales con dimensión 2  

- Expresión de las dos primeras muestras de los datos golub
  para aquellos genes que contienen la expresión **Cyclin**.

```{r}
#| label: pca1
data(golub,package="multtest")
sel = grep("Cyclin",golub.gnames[,2])
golub.red = golub[sel,1:2]
plot(golub.red,xlab="Primera muestra",ylab="Segunda muestra")
```

- Centramos los datos: restamos a cada columna la media de la columna.

```{r}
#| label: pca2
#| echo: false
medias = colMeans(golub.red)
golub.red = sweep(golub.red,2,medias)
plot(golub.red,xlab="Primera muestra",ylab="Segunda muestra")
abline(v=mean(golub.red[,1]),col="red")
abline(h=mean(golub.red[,2]),col="red")
```

- La idea de las componentes principales es considerar una combinación lineal
de los valores originales.

- Se pretende elegir un vector (de dimensión dos) 
$a_{1} = (a_{11},a_{12})$ de modo que en lugar de  
utilizar $x_i$ consideremos (el resumen) 
$u_i = a_{11} x_{i1} + a_{12} x_{i2}$.

- ¿Qué $a_1$ elegimos? 

-  La idea es lograr que los valores $u_i$ tengan
la mayor variabilidad que se pueda con objeto de no perder información.

- En concreto se elige $a_1$ de modo que maximizamos
$$
\frac{1}{n}\sum_{i=1}^n (u_i - \bar{u}) ^2.
$$

## prcomp

- Obtengamos las **componentes principales**. 

```{r}
#| label: pca3
a.pca = prcomp(golub.red)
```
- Los vectores directores de las líneas sobre las que 
  proyectamos aparecen en la siguiente figura.

```{r}
#| label: pca4
#| echo: false
plot(golub.red,xlab="Primera muestra",ylab="Segunda muestra")
abline(h=0,col="red")
abline(v=0,col="red")
arrows(0,0,a.pca$rotation[1,1],a.pca$rotation[2,1],col="blue")
arrows(0,0,a.pca$rotation[1,2],a.pca$rotation[2,2],col="blue")
```

- Estos vectores son 
```{r}
#| label: pca5
a.pca$rotation
```

- Mostrando las líneas sobre las que proyectamos
(en azul) para obtener las dos componentes principales.

```{r}
#| label: pca6
#| echo: false
plot(golub.red,xlab="Primera muestra",ylab="Segunda muestra")
abline(h=0,col="red")
abline(v=0,col="red")
abline(a=0,b = a.pca$rotation[2,1]/a.pca$rotation[1,1],col="blue")
abline(a=0,b = a.pca$rotation[1,2]/a.pca$rotation[2,2],col="blue")
```

-  Mostrando las líneas sobre las que proyectamos
(en azul) para obtener las dos componentes principales y las proyecciones,
esto es, la **primera componente** en verde.

```{r}
#| label: pca7
#| echo: false
plot(golub.red,xlab="Primera muestra",ylab="Segunda muestra")
abline(h=0,col="red")
abline(v=0,col="red")
abline(a=0,b = a.pca$rotation[2,1]/a.pca$rotation[1,1],col="blue")
abline(a=0,b = a.pca$rotation[1,2]/a.pca$rotation[2,2],col="blue")
for(i in 1:nrow(golub.red)){
  x = golub.red[i,]
  y1 = a.pca$rotation[,1]*predict(a.pca)[i,1]
  segments(x[1],x[2],y1[1],y1[2],col="green")
  }
```

-  Y ahora consideremos la segunda componente principal. 

```{r}
#| label: pca8
#| echo: false
plot(golub.red,xlab="Primera muestra",ylab="Segunda muestra")
abline(h=0,col="red")
abline(v=0,col="red")
abline(a=0,b = a.pca$rotation[2,1]/a.pca$rotation[1,1],col="blue")
abline(a=0,b = a.pca$rotation[1,2]/a.pca$rotation[2,2],col="blue")
for(i in 1:nrow(golub.red)){
  x = golub.red[i,]
  y2 = a.pca$rotation[,2]*predict(a.pca)[i,2]
  segments(x[1],x[2],y2[1],y2[2],col="green")
}
```

- Los valores de estas proyecciones (o componentes principales)
  los obtenemos con

```{r}
#| label: pca9
predict(a.pca)
```

## Componentes principales de los datos golub

```{r}
#| label: pca10
golub.pca = prcomp(golub,scale=FALSE,center=TRUE)
```
- El argumento **center=TRUE** centra los datos restando la media de
la columna de modo que las variables tengan medias nulas.
- El argumento **scale=TRUE** hace que las variables originales sean
divididas por su desviación estándar de modo que la varianza (y
la desviación estándar) de las nuevas variables sea la unidad. 

## Criterios de selección del número de componentes

- Uno puede ser la  **proporción total explicada**. Fijar un nivel mínimo
  y quedarnos con el número de componentes necesario para superar este
  valor mínimo.
  
- El segundo puede ser que una componente no puede tener una desviación
  estándar menor que una de las variables originales.
  **Si hemos escalado** cada variable original dividiendo por su desviación
  estándar entonces la   desviación estándar de cada componente ha de ser mayor
  que uno.
  
- Otro criterio puede ser ver en qué momento se produce un descenso 
  de la desviación estándar muy notable. Quedarnos con las componentes 
  previas.


- Un resumen de las componentes nos puede indicar con cuántas nos quedamos.

```{r}
#| label: pca11
summary(golub.pca)
```

-  Si fijamos, por ejemplo un 80% de la variación total a explicar entonces
debemos quedarnos con las cinco primeras componentes.

```{r}
#| label: pca12
a = predict(golub.pca)
a = a[,1:5]
```

- Podemos representar, como es habitual, las dos primeras componentes.

```{r}
#| label: pca13
plot(a[,1],a[,2],xlab="Primera componente",ylab="Segunda componente")
```

## Interpretación

-  Es interesante observar los valores del vector asociado a la primera componente.

```{r}
#| label: pca14
golub.pca$rotation[,1]
```
- Podemos ver que son coeficientes muy parecidos y (casi) todos positivos. 

- Básicamente tenemos la media muestral de todos los niveles de expresión 
en las 38 muestras. 

-  La primera componente es casi media sobre las 38 muestras.

- ¿Y la segunda componente?

```{r}
#| label: pca15
golub.pca$rotation[,2]
```
- Si observamos los coeficientes vemos que las primeros 27 valores son
positivos y los 11 últimos son negativos.

- Además no hay una gran diferencia entre los 27 primeros y tampoco entre 
los 11 últimos.

- Básicamente (aunque no exactamente) estamos comparando, 
para cada gen, la media de los niveles de expresión sobre los datos 
ALL (leucemia linfoblástica aguda) con la media sobre los datos AML 
(leucemia mieloide aguda).

# ¿Componentes principales de las muestras o de los  genes?


##  Genes (características)

- Podemos realizar un análisis de componentes principales de los genes.

```{r}
#| label: pca16
tgolub.pca = prcomp(t(golub),scale=FALSE,center=TRUE)
summary(tgolub.pca)
```

```{r}
#| label: pca17
library(affycoretools)
plotPCA(t(golub),legend = FALSE)
```

## Dos primeras componentes de las muestras

```{r}
#| label: pca18
affycoretools::plotPCA(golub,legend = FALSE)
```


- Cuando realizamos un análisis de los genes sabemos que las observaciones
están clasificadas (las muestras son de dos tipos de leucemia).
- La función **plotPCA** está preparada para mostrarnos si las componentes
nos **reproducen** los grupos que sabemos que previamente tenemos.
- En la siguiente figura mostramos las dos primeras componentes 
pero diferenciando el tipo de muestra (según el tipo de leucemia).
```{r}
#| label: pca19
tipo = factor(golub.cl+1,levels = 1:2,labels = c("ALL","AML")) 
plotPCA(golub,groups = tipo, groupnames = tipo,legend=FALSE)
```

# gse20986

## 
```{r}
#| label: pca20
pacman::p_load("Biobase")
data(gse20986,package="tamidata")
```

```{r}
#| label: pca21
plotPCA(gse20986,groups = pData(gse20986)[,"tissue"],
        groupnames = levels(pData(gse20986)[,"tissue"]))
```


# gse1397 

## Los datos

```{r}
#| label: pca22
data(gse1397,package="tamidata")
pData(gse1397)[,c("tissue","type")]
```

## Tipo

```{r}
#| label: pca23
plotPCA(gse1397, groups = pData(gse1397)[,'type'],
        groupnames = levels(pData(gse1397)[,'type']))
```

## Tejido

```{r}
#| label: pca24
plotPCA(gse1397, groups = pData(gse1397)[,'tissue'],
        groupnames = levels(pData(gse1397)[,'tissue']))
```
