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

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


################################################################################
### PARTE 1: Usar bucles para ejemplos simples de "bootstrap" no paramétrico ###
################################################################################

# En estadística, el término "bootstrap" no paramétrico se refiere a métodos que
# utilizan muestreos aleatorios de datos con remplazo. Estos métodos sirven para
# medir la exactitud (definida en términos de sesgo, varianza, intervalos de
# confianza, error de prediction, entre otras medidas) de los estimados de un
# parámetro.

##################################
### Ejemplo 1 de "bootstrap"   ###
##################################

## TAREA 1 ##
## en esta sección utilizaremos los datos de Blonder et al (2014). Lea el archivo
## "Blonder_etal_2014_PlosBiology_VeinDensityData.csv" y asignelo a un marco de
## datos llamado "VD.data". El archivo está disponible en la página del taller,
## y tiene datos separados por comas.


# Asegúrese de que las primeras filas del marco de datos parezcan estar bien:
head(VD.data)

# El marco de datos debe tener 468 filas y 8 columnas. Verifique que así es:
dim(VD.data)

# Examine los nombres de la variables:
colnames(VD.data)


# Las variables de mayor interés son "VD" y "depth". "VD" es la densidad de las
# venas medida en hojas fosilizadas, "depth" es profundidad estratigráfica del
# fósil y es una medida de el tiempo en el que se fosilizó la hoja. Un valor
# para "depth" de cero corresponde a la frontera K-P (i.e., la frontera entre
# el Cretácico y el Paleógeno), hace al rededor de 65.5 millones de años,
# cuando un asteroide se impactó la tierra causando una extinción en masa.


# Primero crearemos un objeto "time.group" para separar fósiles de plantas que
# se encuentran a profundidades estratigráficas anteriores y posteriores a la
# fromtera K-P:

# Inicialmente, todos los valores son "before": 
time.group <- rep("before", times=nrow(VD.data))

# Después, edite el objeto "time.group", utilizando indexación, para asignar "after"
# a todos los fósiles con valores de "depth" mayores a cero:
time.group[VD.data$depth>0] <- "after" 


# El siguiente código produce una gráfica similar a el panel B en la Figura 2 de
# Blonder et al. (2014):
plot(VD.data$VD ~ VD.data$depth, type="n", xlab="Stratigraphic Depth", 
  ylab="Vein Density")
  
points(VD.data$VD[time.group=="before"] ~ VD.data$depth[time.group=="before"], 
  pch=21, cex=1.5, col="white", bg="firebrick2")
  
points(VD.data$VD[time.group=="after"] ~ VD.data$depth[time.group=="after"], 
  pch=21, cex=1.5, col="white", bg="dodgerblue2")
  
abline(v=0, lwd=2)

# La media de la densidad de venas para fosiles anteriores a la frontera K-P es:
mean.before <- mean(VD.data$VD[time.group=="before"])
mean.before

# Ahora podemos utilizar "bootstrap" para construir intervalos de confianza alrededor
# de esta media. Hay varios métodos para hacer esto, aquí utilizaremos el método de
# percentiles. empezamos por crear una muestra "boostrap" muestreando los datos
# de densidad de venas con remplazo (note el argumento "replace=TRUE"):
boot.sample.before <- sample(VD.data$VD[time.group=="before"], replace=TRUE)

# Ahora calcule el un estimado "bootstrap" de la media:
boot.mean.before <- mean(boot.sample.before)

# Este estimado "bootstrap" es, nuy probablemente, diferente del estimado original de
# la media:
boot.mean.before
mean.before

# A continuación vamos a utilizar un bucle para calcular 'n' estimados bootstrap
# de la media:
n <- 10000

# cree un vector con NAs que será editado para asignar estimados bootstrap de la
# media a medida que se obtienen
boot.mean.before <- rep(NA, times=n)

## TAREA 2 ##
## El bucle a continuación está casi listo. Cada iteración del bucle debe calcular 
## un estimado bootstrap de la media y asignarlo al un sitio correspondiente del
## vector 'boot.mean.before'. Sin embargo, hay un error en el código. La media
## calculada en cada iteración se asigna al primer valor del vector 'boot.mean.before'.
## Corrija el error.
for(i in 1:n)
{
  boot.sample.before <- sample(VD.data$VD[time.group=="before"], replace=TRUE)
  boot.mean.before[1] <- mean(boot.sample.before)
}


# El siguiente código grafica el histograma de estimados bootstrap de la media
# de la densidad de venas de las hojas fosilizadas antes de la frontera K-P y
# lo compara con el estimado original de la media:
hist(boot.mean.before, breaks=100, col="firebrick1", border="firebrick1", 
  main="Mean before KP boundary")
abline(v=mean.before, lwd=2)

# El método de percentiles utiliza los percentiles 2.5 y 97.5 de la distriución
# bootstrap para calcular el intervalo del 95% de confianza para la media: 
CIs.before <- quantile(boot.mean.before, prob=c(0.025, 0.975))


## TAREA 3 ##
## Modifique el código anterior para calcular la distribución bootstrap y el
## intervalo de confianza del 95% para la media de la densidad de venas en
## hojas fosilizadas después de K-P boundary. Use nombres similares a los
## utilizados en el códugo anterior, por ejemplo, 'CIs.after' en vez de
## 'CIs.before'.

# Si siguió las convenciones sugeridas para nombrear objetos, puede utilizar
# el siguiente código para crear una gráfica de barras que muestre las medias
# y sus intervalos de confianza correspondientes:
bp <- barplot(c(mean.before, mean.after), names=c("Before", "After"), 
  ylim=range(c(0, CIs.before, CIs.after)), ylab="Vein Density")
lines(x=bp[c(1,1),1], y=CIs.before, lwd=3)
lines(x=bp[c(2,2),1], y=CIs.after, lwd=3)


##################################
### Ejemplo 2 de bootstrapping ###
##################################
# Hay una relación lineal entre "VD" y "depth":
VD.depth.lm <- lm(VD.data$VD ~ VD.data$depth)
summary(lm(VD.depth.lm)) # Note la significancia estadística de esta regresión

# Podemos utilizar un gráfico para examinar esta relación:
plot(VD.data$VD ~ VD.data$depth, pch=21, cex=1.5, col="white", bg="grey30", 
  lwd=1.5, xlab="Stratigraphic Depth", ylab="Vein Density")
abline(VD.depth.lm, lwd=3, col="goldenrod2") #adiciona la linea predecida por la regresión

# El método bootstrap puede utilizarse para construir intervalos de confianza de los
# parámetros de la regresión. Estos son los estimados de los coeficientes obtenidos con
# con los datos originales:
reg.coeff <- coefficients(VD.depth.lm)
reg.coeff

# Para este tipo de bootstrap, es necesario separar los valores predecidos por el
# modelo lineal de los residuos del modelo:
pred.VD <- fitted(VD.depth.lm)
resid.VD <- residuals(VD.depth.lm) 

# El siguiente código estandariza los residuos:
resid.VD <- ( (length(resid.VD)/(length(resid.VD)-1))^(1/2) ) * resid.VD

# Ahora vamos a crear una muestra bootstrap de los residuos:
boot.resid.VD <- sample(resid.VD, replace=TRUE)

# Y utilizamos la muestra bootstrap para crear una nueva variable, sumando la muestra
# bootstrap a los valores predecidos:
boot.VD <- pred.VD + boot.resid.VD

# Estos valores de VD creados a partir de una muestra bootstrap se usan ahora
# en un modelo linear de regresión para calcular los coeficientes de regresión:
boot.VD.depth.lm <- lm(boot.VD ~ VD.data$depth)
boot.reg.coeff <- coefficients(boot.VD.depth.lm)
boot.reg.coeff

# Este procedimiento se repite 'n' veces para obtener una distribución bootstrap
# de los coeficientes de regresión.


### TAREA 4 ##
## Siguiendo los pasos descritos arriba, complete el siguiente código para
## producir estimados bootstrap para los coeficientes de la regresión de
## "VD" en función de "depth". Note que los residuos y valores predecidos
## ya se calcularon, no es necesario re-calcularlos dentro del bucle. Sólo
## es necesario 1) muestrear con remplazo los residuos, 2) sumar la muestra
## bootstrap de los residuos a los valores predecidos para construir una
## variable de respuesta, 3) calcular la regresión con esta variable de
## respuesta, y 4) guardar los coeficientes.

n <- 5000

# El siguiente código crea una matriz con valores NA y que será editada 
# para alojar los valores de los coeficientes de las regresiones en
# cada iteración del bucle:
boot.reg.coeff <- matrix(NA, nrow=n, ncol=length(reg.coeff))
colnames(boot.reg.coeff) <- names(reg.coeff)

head(boot.reg.coeff)

# El código a continuación es el bucle para construir una distribución
# bootstrap de los coeficientes de regresión:
for(i in 1:n)
{

  ## ESCRIBA AQUÍ SU CÓDIGO ##
  
  # El siguiente código asigna valores de los coeficientes de regresión
  # en la fila "i" de la matriz "boot.reg.coeff":
  boot.reg.coeff[i,] <- ## SU CÓDIGO VA AQUÍ ##

}


# Usando el mismo método de percentiles, calcule el intervalo de confianza
# del 95% para cada coeficiente de regresión: 
CIs.intercept <- quantile(boot.reg.coeff[,1], prob=c(0.025, 0.975))
CIs.slope <- quantile(boot.reg.coeff[,2], prob=c(0.025, 0.975))

CIs.intercept
CIs.slope

# El intervalo de confianza para la pendiente no incluye cero y, por lo tanto
# es estadísticamente significatica.

# Este es el histograma de la distribución bootstrap de la pendiente:
hist(boot.reg.coeff[,2], breaks=100, col="firebrick1", border="firebrick1", 
  main="", xlab="Slope")
  
# Este polígono delimita intervalo del 95% de confianza:  
polygon(x=CIs.slope[c(1,1,2,2)], y=c(0, n, n, 0), col=gray(0.5, 0.5), 
  border=gray(0.5, 0.5))

# Esta linea muestra la posición del estimado original de la pendiente:  
abline(v=reg.coeff[2], lwd=2)



################################################################################
### PARTE 2: Uso de bucles en ejemplos simples de pruebas de aleatorización  ###
################################################################################

# Las pruebas realizadas en aleatorización tienen como objetico calcular la
# significancia estadística, en términos de valores de p, y el tamaño de los
# efectos. Esto se logra comparando los estadísticos empíricos a la
# distribución esperada bajo una hipótesis nula particular. Esta hipótesis
# nula se define aleatorizando los datos. 


###################################
### Ejemplo 1 de aleatorización ###
###################################

# La hipótesis clave que estudian Blonder et al. es que la extinción masiva de
# la frontera K-P constituyó un evento de selección marcado. Esto
# significaría que las plantas que sobrevivieron la extinción en masa
# no serían una muestra aleatoria de los linajes que existían antes de tal
# evento. Para poner a prueba esta hipótesis, Blonder et al. derivaron de la
# hipótesis la predicción según la cual la densidad media de las venas de las
# hojas antes y después de la frontera K-P son estadísticamente diferentes.
# Blonder et al. utilizaron análisis paramétricos para hacer esta comparación
# de medias. Aquí utilizaremos una prueba basada en aleatorización para comparar
# las medias.
# 

# Este es un gráfico de cajas (boxplot) similar al del panel B de la Figure 2
# de Blonder et al. 
boxplot(VD.data$VD[time.group=="before"], VD.data$VD[time.group=="after"], 
  ylab="Vein Density", names=c("Before", "After"), 
  col=c("firebrick2", "dodgerblue2"))

# Este código realiza la prueba paramétrica, una prueba t. Como puede verse
# esta prueba indica que la diferencia es estadísticametne significativa:
t.test.res <- t.test(VD.data$VD[time.group=="after"], 
  VD.data$VD[time.group=="before"])
t.test.res

# Esto aisla el valor del estadístico de la prueba t:
emp.t <- t.test.res$statistic
emp.t


# If the prediction is that the two groups ("before" and "after") come from two
# different populations with different means, then the null hypothesis in that
# both groups come from the same population. Thus, an algorithm whereby values
# of vein density are randomly assigned among the two groups is a reasonable way
# to produce random (i.e. null) values of the statistic.

# La forma más fácil de hacer una prueba basada en un modelo nulo que aleatoriza
# los datos es asignar aleatoriamente las hojas fósiles a los grupos "before" y
# "after":
rand.time.group <- sample(time.group, replace=FALSE)

# Note que al argumento 'replace' se le da el valor FALSE, de forma que el
# vector es re-ordenado aleatoriamente (no muestreado con remplazo como en
# bootstrap).

# el siguiente gráfico de cajas (boxplot) muestra la aleatorización de los datos
# de densidad de venas de las hojas:
boxplot(VD.data$VD[rand.time.group=="before"],  
  VD.data$VD[rand.time.group=="after"], ylab="Vein Density", 
  names=c("Before", "After"))

# Ahora el nuevo vector de asignación de grupo se utiliza para calcular el estadístico t:
rand.t.test.res <- t.test(VD.data$VD[rand.time.group=="before"], 
  VD.data$VD[rand.time.group=="after"])
rand.t <- rand.t.test.res$statistic


# Debemos llevar a cabo este procedimiento un número grande de veces:
n <- 9999

rand.t <- rep(NA, times=n)

for(i in 1:n)
{
  rand.time.group <- sample(time.group, replace=FALSE)
  
  rand.t.test.res <- t.test(VD.data$VD[rand.time.group=="after"], 
    VD.data$VD[rand.time.group=="before"])

  rand.t[i] <- rand.t.test.res$statistic
}

# Finalmente, el siguiente código combina el valor empírico del estadístico t 
# con los valores aleatorizados:
rand.t <- c(emp.t, rand.t)

# El siguiente código ilustra la distribución de t obtenida de los datos aleatorizados
# y la compara al valor del estadístico t obtenido a partir de los datos originales:
hist(rand.t, breaks=100, col="firebrick1", border="firebrick1", 
  main="", xlab="t-statistic")
abline(v=emp.t, lwd=2)
 

# Ahora podemos calcular el valor de p de una cola:
p.value <- length(which(rand.t>=emp.t)) / length(rand.t)
p.value

# La diferencia es estadísticamente significativa. 



###################################
### Ejemplo 2 de aleatorización ###
###################################

# Los resultados obtenidos arriba son iguales a los obtenidos por Blonder et al.;
# sin embargo, Blonder et al. no consideraron el posible efecto de la relación
# positiva entre la densidad de las venas y la profundidad estratigráfica.
plot(VD.data$VD ~ VD.data$depth, pch=21, cex=1.5, col="white", bg="grey30", 
  lwd=1.5, xlab="Stratigraphic Depth", ylab="Vein Density")
abline(VD.depth.lm, lwd=3, col="goldenrod2") #esta es la línea de regresión
abline(v=0, lwd=3, col="black") #esta línea representa la frontera K-P

summary(lm(VD.depth.lm)) # este es el resumen de la regresión

# Dado esta relación positiva, la división de los datos en dos partes
# estratigráficas necesariamente resulta en medias distintas. Por lo tanto,
# la hipótesis nula apropiada debe considerar esta relación positiva. Como
# antes, empezamos por calcular un estadístico observado:
t.test.res <- t.test(VD.data$VD[time.group=="after"], 
  VD.data$VD[time.group=="before"])
  
# Esta es el estadístico t que mide la diferencia entre las medias antes y
# después de la frontera K-P:
emp.t <- t.test.res$statistic
emp.t


# Ahora escoja al azar una profundidad estratigráfica:
rand.depth <- runif(n=1, min=min(VD.data$depth), max(VD.data$depth))

# Adicionémos este valor de profundidad estratigráfica al gráfico:
abline(v=rand.depth, lwd=2, col="red")

# Usando el valor estratigráfico escojido al azar, divida los datos de densidad
# de venas en dos grupos:
rand.time.group <- rep("before", times=nrow(VD.data))
rand.time.group[VD.data$depth>rand.depth] <- "after"

# Finalmente, calculamos el estadístico t para la división aleatoria de los
# datos:
rand.t.test.res <- t.test(VD.data$VD[rand.time.group=="after"], 
  VD.data$VD[rand.time.group=="before"])
  
rand.t <- rand.t.test.res$statistic
rand.t  


## TAREA 5 ##
## Siguiendo los pasos descritos arriba, complete el código a continuación
## para producir una distribución aleatorizada del estadístico t. Esta
## distribución representa la hipótesis nula. 
n <- 99999

rand.t <- rep(NA, times=n)

# This is the bootstrapping loop where 'i' will take values from 1 to 'n':
for(## ESCRIBA AQUÍ SU CÓDIGO ##)
{
  ## ESCRIBA AQUÍ SU CÓDIGO ##
}

  
rand.t <- c(emp.t, rand.t) # combina los valores observados y aquellos obtenidos por medio de aleatorización

## TAREA 6 ##
## Haga un histograma que muestre los valores de la distribución nula del estadístico
## t y el valor observado de este estadístico.


## TAREA 7 ##
## Calcule la significancia estadística, o valor de p, que es la probabilidad de que emp.t >= rand.t. Note que 
## el resultado es solo marginalmente significativo.




















################################################################################
### RESPUESTAS #################################################################
################################################################################

## TAREA 1 ##
VD.data <- read.table(file=file.choose(), header=TRUE, sep=",")


## TAREA 2 ##
boot.mean.before[i] <- mean(boot.sample.before)


## TAREA 3 ##
mean.after <- mean(VD.data$VD[time.group=="after"])

n <- 10000

boot.mean.after <- rep(NA, times=n)

for(i in 1:n)
{
  boot.sample.after <- sample(VD.data$VD[time.group=="after"], replace=TRUE)
  boot.mean.after[i] <- mean(boot.sample.after)
}
 
CIs.after <- quantile(boot.mean.after, prob=c(0.025, 0.975))


## TAREA 4 ##
n <- 10000

boot.reg.coeff <- matrix(NA, nrow=n, ncol=length(reg.coeff))
colnames(boot.reg.coeff) <- names(reg.coeff)

head(boot.reg.coeff)

for(i in 1:n)
{
  boot.resid.VD <- sample(resid.VD, replace=TRUE)

  boot.VD <- pred.VD + boot.resid.VD

  boot.VD.depth.lm <- lm(boot.VD ~ VD.data$depth)

  boot.reg.coeff[i,] <- coefficients(boot.VD.depth.lm)
}


  
## TAREA 5 ##
for(i in 1:n)
{
  rand.depth <- runif(n=1, min=min(VD.data$depth), max(VD.data$depth))

  rand.time.group <- rep("before", times=nrow(VD.data))
  rand.time.group[VD.data$depth>rand.depth] <- "after"

  rand.t.test.res <- t.test(VD.data$VD[rand.time.group=="after"], 
    VD.data$VD[rand.time.group=="before"])
  
  rand.t[i] <- rand.t.test.res$statistic
}

  

## TAREA 6 ##
hist(rand.t, breaks=30, col="firebrick1", border="firebrick1", 
  main="", xlab="t-statistic")
abline(v=emp.t, lwd=2)


## TAREA 7 ##
p.value <- length(which(rand.t>=emp.t)) / length(rand.t)
p.value