################################################################################
### TALLER FUNDAMENTOS DE R                                                  ###
### EJERCICIO 9.2: CONTROL DE FLUJO                                          ###
###                                                                          ###
### Center for Conservation and Sustainable Development                      ###
### Missouri Botanical Garden                                                ###
### Website: rbasicsworkshop.weebly.com                                      ###
################################################################################

## OBJETIVO:
## Practicar el uso de de bucles ("loops")

# El análisis de componentes principales (PCA) es una transformación que rota
# un conjunto de variables (posiblemente correlacionadas) para expresarlas en
# términos de un segundo conjunto de variables ortogonales (no correlacionadas)
# llamadas componentes principales. El número de componentes principales es menor
# o igual al número de variables originales. La rotación se define de tal manera
# que el primer componente principal tiene la mayor varianza que puede ser
# representada en un eje lineal. De forma similar, cada componente principal
# subsiguiente tiene la mayor varianza posible, bajo la restricción de ser
# ortogonal a los componentes anteriores.

# Frecuentemente el PCA se utiliza para reducir la dimensionalidad de los datos.
# La idea es que el conjunto de variables originales pueden expresarse resumida
# pero adecuadamente en los primeros componenetes principales, bajo la suposición
# que estos componentes principales capturan gran parte de la variación de los
# datos originales. En esta aplicación del PCA surge entonces la pregunta:
# ¿cuántos componentes principales deben usarse para representar la variación
# original? El código en este ejercicio implementa un método propuesto para
# responder esta pregunta. Este método utiliza un modelo nulo para poner a prueba
# la significancia estadística de cada uno de los componentes principales. 

# A continuación utilizará los datos morfológicos de los Iris de Edgar Anderson
# (vea la página de ayuda para esta base de datos, escriba en su consola 
# "help(Iris)"). Ud implementará un método descrito en el apartado 2.1.2. de
# Peres-Neto et al. 2005 (How many principal components? stopping rules for
# determining the number of non-trivial axes revisited. Computational Statistics
# & Data Analysis 49: 974-997). Este método pone a prueba la significancia
# estadística de los "eigenvalues" de cada componente principal. Estos
# "eigenvalues" miden la varianza de los datos a lo largo de un componente
# principal. Para ser estadísticamente significativo, un componente principal
# debe tener un "eigenvalue" mayor que lo esperado bajo un modelo nulo que
# aleatoriza las relaciones entre las variables originales.

data(iris) # Carga el set de datos Iris del paquete "datasets"

iris[1:5, ] # Imprime las primeras 5 filas del marco de datos


## TAREA 1 ##
## Cree un objeto llamado morfo que contiene sólo las primeras 4 columnas
## del marco de datos "iris", estas columnas continen los datos morfológicos
## de longitud de los sépalos, amplitud de los sépalos, longitud de los pétalos
## y amplitud de los pétalos.

# El siguiente código corre un PCA de los datos morfológicos, usando la función *prcomp*
pca.iris <- prcomp(morpho, scale=TRUE) 

# El siguiente código imprime un resumen de los resultados del PCA
summary(pca.iris) 

# El siguiente código muestra la estructura del onbjeto "pca.iris" que contiene los resultados del PCA
str(pca.iris)

# Este elemento de "pca.iris" es un vector que contiene la raíz cuadrada de los "eigenvalues"
# de los componentes principales 
pca.iris$sdev 

# El siguiente código pone los "eigenvalues" empíricos (de los datos reales) en un objeto. Estos
# valores se compararán con "eigenvalues" obtenidos a partir de datos en los que se han
# aleatorizado las relaciones entre las variables originales. 
empirical.ev <- pca.iris$sdev
names(empirical.ev) <- paste("PC", 1:length(empirical.ev), sep=" ")

# Para cada componente principal existe un "eigenvalue" o "valor propio". Estos
# valores miden la variación en las variables originales (e.g. las variables morfológicas) que está
# representada (o "capturada") en cada componente principal. 


## TAREA 2 ##
## Haga una "barplot" que muestre los "eigenvalues" empíricos

# El siguiente código crea un objeto que define el número de iteraciones del modelo
# nulo (i.e., el número de veces que se aleatorizarán las relaciones entre las variables
# originales
k <- 999

# El siguiente código crea una matriz de NAs que será editada en cada iteración 
# del bucle con "eigenvalues" obtenidos en las k iteraciones del modelo nulo
rand.ev <- matrix(data=NA, nrow=k, ncol=length(empirical.ev))
rownames(rand.ev) <- paste("rand", 1:k, sep="_")


# Tenga en cuenta que la expresión dentro del bucle termina con la TAREA 5
for (i in 1:k)
{
  
  rand.morpho <- morpho
  # Esto copia los datos morfológicos empíricos que luego serán randomizados
  
  rand.morpho[,1] <- sample(rand.morpho[,1])
  # Esto aleatoriza la posición de los valores en la primera variable (columna) 
  # en "rand.morpho". 

  ## TAREA 3 ##
  ## Esta prueba se basa en la comparación de los "eigenvalues" empíricos en 
  ## una PCA con "eigenvalues" del PCA donde las relaciones entre las variables 
  ## han sido aleatorizadas (desruyendo así cualquier correlación entre las variables
  ## originales). Esto requiere que los valores de cada variable sean aleatorizados de
  ## manera independiente. Por lo tanto: aplique el paso anterior de forma independiente
  ## a cada una de las otras columnas del marco de datos morfológico. 
  

  ## TAREA 4 ##
  ## Despues de aleatorizar todas las variables, ejecute un nuevo análisis de 
  ## componentes principales, pero esta vez utilizando los datos morfológicos 
  ## aleatorizados. Guarde los resultados de este PCA en un objeto denominado 
  ## "rand.pca.iris"


  ## TAREA 5 ##
  ## La siguiente línea de código pone los "eigenvalues" del PCA randomizado en la
  ## matriz "rand.ev". Con cada iteración del bucle, se supone que el código debe poner 
  ## estos valores en una fila distinta de esta matriz. Sin embargo hay un error,
  ## y los resultados están siendo colocandos siempre en la última fila. Corrija el error
  ## en la siguiente línea de código.
  rand.ev[k,] <- rand.pca.iris$sdev
  
}
# termina el bucle


# El siguiente código combina por filas ("rbind") el vector de los "eigenvalues" empíricos
# con la matriz de "eigenvalues" aleatorios, poniendo los valores empíricos en la primera
# fila de la matriz "rand.ev":
rand.ev <- rbind(empirical.ev, rand.ev)


## TAREA 6 ##
## Abra una ventana y dividala en 4 páneles con las función *par* o *layout*

# El siguiente código crea un histograma con la distribución del primer "eigenvalue"
# (columna 1), mostrando la variación representada en el primer componente principal
hist(rand.ev[,1]^2, breaks=100, xlab="Eigenvalue", cex.lab=1.5, 
  cex.axis=1.5, main="PC 1")

# El siguiente código resalta con una línea la posición del "eigenvalue" empírico en esa 
# distribución
abline(v=rand.ev[1,1]^2, col="red")

# La siguiente línea de código calcula el valor de p como la proporción de
# "eigenvalues" de la distribución aleatorea que son iguales o mayores que 
# el "eigenvalue" empírico. Tome en cuenta que estamos haciendo la suma de un 
# vector lógico donde los TRUEs se tratan como 1 y los FALSE como 0.
sum(rand.ev[,1] >= rand.ev[1,1]) / (k+1) 


## TAREA 7 ##
## Cree las figuras correspondientes para los "eigenvalues" 2, 3 y 4 en el resto
## de paneles de la figura que está produciendo. Además, calcule también los 
## valores de p para el resto de "eigenvalues".


## TAREA 8 ##
## Elabore una nueva versión de la figura que sea tan profesional como sea posible y
## guardela como un archivo TIFF de alta resolución (600 ppp).




















################################################################################
### SOLUCIONES PARA TAREAS #####################################################
################################################################################

## TAREA 1 ##
morpho <- iris[,1:4]


## TAREA 2 ##
barplot(empirical.ev^2, col="firebrick1", border="firebrick1", 
  xlab="Principal Components", ylab="Eigenvalue")


## TAREA 3 ##
rand.morpho[,2] <- sample(rand.morpho[,2])
rand.morpho[,3] <- sample(rand.morpho[,3])
rand.morpho[,4] <- sample(rand.morpho[,4])
  

## TAREA 4 ##
rand.pca.iris <- prcomp(rand.morpho, sclale=TRUE) 

 
## TAREA 5 ##
rand.ev[i,] <- rand.pca.iris$sdev


## TAREA 6 ##
par(mfrow=c(2,2))
#or
layout(m=matrix(1:4, ncol=2))
layout.show(4)


## TAREA 7 ##
hist(rand.ev[,2]^2, breaks=100, xlab="Eigenvalue", cex.lab=1.5, 
  cex.axis=1.5, main="PC 2")
abline(v=rand.ev[1,2]^2, col="red")
sum(rand.ev[,2] >= rand.ev[1,2]) / (k+1) 

hist(rand.ev[,3]^2, breaks=100, xlab="Eigenvalue", cex.lab=1.5, 
  cex.axis=1.5, main="PC 3")
abline(v=rand.ev[1,3]^2, col="red")
sum(rand.ev[,3] >= rand.ev[1,3]) / (k+1) 

hist(rand.ev[,4]^2, breaks=100, xlab="Eigenvalue", cex.lab=1.5, 
  cex.axis=1.5, main="PC 4")
abline(v=rand.ev[1,4]^2, col="red")
sum(rand.ev[,4] >= rand.ev[1,4]) / (k+1) 


## TAREA 8 ##
getwd() # Revise el directorio de trabajo y cámbielo si es necesario


line.wd <- 2
line.col <- "black"
bar.col <- "firebrick3"
main.size <- 1.75
axis.size <- 1.5
lab.size <- 1.5
breaks.n <- 50


tiff(filename = "PCA_exercise.tiff", width = 25, height = 20, units = "cm", 
  pointsize = 12, res = 600)

  par(mfrow=c(2,2), mar=c(5, 6, 4, 2), mgp=c(3.75, 1, 0))
  
  hist(rand.ev[,1]^2, breaks=breaks.n, xlab="Eigenvalue", las=1, cex.lab=lab.size, 
    cex.axis=axis.size, main="PC 1", cex.main=main.size, col=bar.col, border=bar.col)
  abline(v=rand.ev[1,1]^2, col=line.col, lwd=line.wd)
  
  hist(rand.ev[,2]^2, breaks=breaks.n, xlab="Eigenvalue", las=1, cex.lab=lab.size, 
    cex.axis=axis.size, main="PC 2", cex.main=main.size, col=bar.col, border=bar.col)
  abline(v=rand.ev[1,2]^2, col=line.col, lwd=line.wd)
  
  hist(rand.ev[,3]^2, breaks=breaks.n, xlab="Eigenvalue", las=1, cex.lab=lab.size, 
    cex.axis=axis.size, main="PC 3", cex.main=main.size, col=bar.col, border=bar.col)
  abline(v=rand.ev[1,3]^2, col=line.col, lwd=line.wd)
  
  hist(rand.ev[,4]^2, breaks=breaks.n, xlab="Eigenvalue", las=1, cex.lab=lab.size, 
    cex.axis=axis.size, main="PC 4", cex.main=main.size, col=bar.col, border=bar.col)
  abline(v=rand.ev[1,4]^2, col=line.col, lwd=line.wd)

dev.off()