################################################################################
### TALLER FUNDAMENTOS DE R                                                  ###
### EJERCICIO 10.1: ESCRIBIR FUNCIONES                                       ### 
###                                                                          ###
### Center for Conservation and Sustainable Development                      ###
### Missouri Botanical Garden                                                ###
### Website: rbasicsworkshop.weebly.com                                      ### 
################################################################################


## TAREA 1 ##
# Use el siguiente código para crear la lista llamada "L1"
L1 <- list(c(0.01, 3.1), c(0.02, 4.0, 2.6, 7.01, 0.2), c(0.05, 3.5), c(0.01, 2.9, 0.44), c(0.03, 3.1), c(0.04, 3.4)) 
# examine las propiedades de la lista 
L1
class(L1)
mode(L1)
length(L1)
str(L1)

# visite y estudie la página de ayuda de la función *lapply*
?lapply

# ahora utilice la función *lapply* para determinar si cada elemento de L1 tiene
# una longitud mayor a 2. Para hacerlo, escriba una función y usela como valor
# del argumento "FUN" de la función *lapply*. El resultado debe ser una lista 
# de longitud igual a L1.Esta es una forma de escribir la función:  

my.task1.fun <- function(x)
{
  result <- length(x)>2
  result
}

# otra opción:

my.task1.fun <- function(x)
{
  length(x)>2
}

# y una tercera opción:

my.task1.fun <- function(x) length(x)>2

# ahora use la función que acaba de escribir como valor del argumento "FUN"
# en la función "lapply" 

lapply(L1, FUN = my.task1.fun)

# Listo! Asegúrese que la respuesta es correcta, examine la longitud de cada
# elemento de L1.


## TAREA 2 ##
# Defina su directorio de trabajo con la función *setwd*, y verifíque lo con
# la función *getwd*.

setwd("C:/_transfer/R_Basics_Workshop/St_Louis_May_2017/Datasets") #this is the working directory at Ivan's laptop
getwd()

# Use la función *read.csv* para crear un marco de datos llamado "Nicaragua.data" 
# a partir del archivo "nicaraguadata.csv" (disponible en la página del taller).
 
Nicaragua.data <- read.csv("nicaraguadata.csv", sep = ",", header=T)
#examine las propiedades del marco de datos 
head(Nicaragua.data)
str(Nicaragua.data)

# En el marco de datos "Nicaragua.data" cada fila representa un especímen de herbario
# collectado en Nicaragua, y diferentes columnas proveen información sobre cada
# especímen, incluyendo la familia, el género y la especie, así como información
# sobre la localidad de colecta. Toda esta información fue obtenida de la base de
# datos Tropicos, del Missouri Botanical Garden.

# Use la función *tapply* para determinar el número de especies por género en el marco
# de datos "Nicaragua.data""Nicaragua.data". Para ello, escriba una función y use esa
# función como valor del argumento "FUN" de la función *tapply*. Recuerde que el
# ejercicio es determinar el número de especies en cada género (y NO el número de 
# especímenes en cada género). Para empezar, visite la página de ayuda para la
# función *tapply*.

## TAREA 3 ##
# Vea el Ejemplo 3 en la presentación sobre "escribir funciones". Modique la función
# "my.twosam" para obtener un valor de significancia estadística (es decir, un valor
# de p), además del valor del estadístico t. Use la nueva versión de la función para
# calcular el estadístico t y su significancia estadística (valor de p) para la comparación
# de la longitud de los pétalos entre Iris versicolor and Iris setosa, basandose en los
# datos del marco de datos "Iris".

## TAREA 4 ##
# Estudie cuidadosamente y utilice el siguiente código. Este código define una función
# que incorpora la estocasticidad ambiental en el modelo de tamaño poblacional de 
# Ricker. La estocasticidad ambiental es la vairabilidad en la tasa crecimiento del
# tamaño poblacional que se debe a variación temporal en las condiciones ambientales
# (y no debe confundirse con estocasticidad demográfica, que se refiere a la variación
# estocástica en el orden de nacimientos y muertes). La siguiente función asume que
# la estocasticidad ambiental puede describirse aproximadamente con una distribución
# normal de la tasa de crecimiento poblacional. El modelo incluye los siguientes
# parámetros:
# "nzero" es el tamaño inicial de la población
# "r" es la media de la tasa de crecimiento poblacional
# "es" es la desciación estandar de la tasa de crecimiento poblacional y, por lo tanto,
#      representa la maginitud de la estocasticidad ambiental
# "K" es la capacidad de carga
# "time" es el número de unidades de tiempo en las cuales se modelará la dinámica de la población.
# Cree la función llamada *my.second.ricker.fun*:

my.second.ricker.fun <- function(nzero, r, es, K, time)
{
	N <- numeric(time+1)
	N[1] <- nzero
	r.es <- rnorm(time, r, es)
	for(i in 1:time)
	{
		N[i+1] <- N[i]*exp(r.es[i]*(1 - (N[i]/K)))
	}
	Time <- 0:time
	plot(Time, N, type="l", xlim=c(0, time), ylim=c(0,K+10), cex.axis=1.5, 
	    cex.lab=1.5, bty="n", lwd=2)
	abline(h=K, lty=3, col="red")
}

# use la función *my.second.ricker.fun"*

for(j in 1:1000)
{
	my.second.ricker.fun(nzero=1, r=0.1, es=0.5, K=30, time=100)
	Sys.sleep(0.1)
}


## TAREA 5 ##
# Modifique la función anterior (*my.second.ricker.fun*) para que el resultado de la
# función sea un vector con el tamaño poblacional en cada tiempo. Visite la página
# de ayuda para la functión *function* y lea detenidamente la sección sobre "Details".

## TAREA 6 ##
# Escriba un bucle que i) realice 1000 iteraciones de la función que usted
# creó en la Tarea 5, ii) capture los 1000 vectores producidos (cada vector contiene
# el tamaño poblacional en cada punto en el tiempo en cada iteración). Construya
# una gráfica que muestre la dinámica de la población en las 1000 iteraciones.

## TAREA 7 ##
# El bucle que escribió en la Tarea 6 simula 1000 iteraciones de la dinámica  
# de la población. Escriba código que determine el número de iteraciones en que
# la población se extinguió (i.e., alcanzó un tamaño poblacional cercano a cero).
# Note que la función produce números reales, que podrían hacer sentido si el tamaño
# poblacional se mide en ciertas unidades (e.g., biomasa), pero que pueden ser arbitrariamente
# pequeños. Por lo tanto, deberá definir un umbral de tamaño poblacional bajo el cual
# la extinción ocurre. 

## TASK 8 ##
# Modifique la función *my.second.ricker.fun* (vea la Tarea 4) para incorporar
# variación temporal estocástica en la capacidad de carga "K", suponiendo
# que tal variación puede representarse aproximadamente con una distribución
# normal. Escriba un bucle para realizar 1000 iteraciones de la dinámica de la
# población y grafique cada iteración.





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

## TAREA 1 ##

# ya tiene la respuesta.

## TAREA 2 ##

# esta es una forma de escribir la función:

my.task2.fun <- function(x)
{
  result <- length(unique(x))
  result
}

# otra forma:

my.task2.fun <- function(x)
{
  length(unique(x))
}

# una tercera forma:

my.task2.fun <- function(x) length(unique(x))

# ahora utilice la función *my.task2.fun* en la función *tapply*:
 
tapply(Nicaragua.data$Species, Nicaragua.data$Genus, FUN = my.task2.fun)

# voila!


## TAREA 3 ##

my.twosam.version2 <- function(y1, y2) # two arguments without pre-determined values
{
  # Sample Size
  n1 <- length(y1)
  n2 <- length(y2)

  # Averages
  b1 <- mean(y1)
  b2 <- mean(y2)

  # Variances
  s1 <- var(y1)
  s2 <- var(y2)

  # Joint Variance
  s <- ((n1-1) * s1 + (n2-1) * s2) / (n1 + n2 - 2)

  # t statistic
  tst <- (b1 - b2) / sqrt(s * (1 / n1 + 1 / n2))

  # p value
  p.1 <- pt(tst, df=n1 + n2 - 2, lower.tail = TRUE) #starting at the lower tail
  p.2 <- pt(tst, df=n1 + n2 - 2, lower.tail = F) #starting at the upper tail
  p.final <- min(p.1, p.2) #the p-value is the smallest of p.1 and p.2	

  # Result: t-value and p-value
  list(t.statistic = tst, p.value = p.final) 
}

# Ahora utilice la función:

data(iris)
PL.versi <- iris$Petal.Length[which(iris$Species == "versicolor")]
PL.seto <- iris$Petal.Length[which(iris $ Species == "setosa")]
boxplot(PL.versi, PL.seto, names=expression(italic("Iris versicolor"), italic("Iris setosa")), 
    ylab="Petal length (cm)", cex.axis=1.5, cex.lab=1.5)
tstat <- my.twosam.version2(PL.versi, PL.seto)
tstat


## TAREA 4 ##

# ya tiene la respuesta

## TAREA 5 ##

my.second.ricker.fun <- function(nzero, r, es, K, time)
{
	N <- numeric(time+1)
	N[1] <- nzero
	r.es <- rnorm(time, r, es)
	for(i in 1:time)
	{
		N[i+1] <- N[i]*exp(r.es[i]*(1 - (N[i]/K)))
	}
	Time <- 0:time
	plot(Time, N, type="l", xlim=c(0, time), ylim=c(0,K+10), cex.axis=1.5, 
	    cex.lab=1.5, bty="n", lwd=2)
	abline(h=K, lty=3, col="red")
	return(N)
}


## TAREA 6 ##

my.results.N <- matrix(NA, 1000,101)

for(j in 1:1000)
{
	my.results.N[j,] <- my.second.ricker.fun(nzero=1, r=0.1, es=0.5, K=30, time=100)
	#Sys.sleep(0.1)
}

#examine los resultados
my.results.N[1:5, 1:5]

plot(0:100, my.results.N[1,], type="n", ylim=range(my.results.N), xlab="Time", ylab="Population size", cex.axis=1.5, cex.lab=1.5, bty="n")
for(i in 1:1000)
{
	points(0:100, my.results.N[i,], type="l", ylim=range(my.results.N))
}
abline(h=30, lty=3, col="red")


## TAREA 7 ##

#defina el umbral poblacional bajo el cual la extinción ocurre:
extinction.threshold <- 0.001

#cree un vector que tomará un valor de 1 cuando la población se extingue y 0 cuando sobrevive
extinct <- rep(NA, times=nrow(my.results.N))

#use un bucle para determinar cuales poblaciones se extinguieron
for(i in 1:nrow(my.results.N))
{
	extinct[i] <- (sum(my.results.N[i,]<=extinction.threshold))>0
}

#examine los resultados
summary(extinct)
sum(extinct)

#determine qué filas de "my.results.N" corresponden a las poblaciones que se extinguieron
which(extinct)
#plot the population size through time for the cases in which the population went extinct
plot(0:100, my.results.N[which(extinct)[1],], type="l", xlab="Time", ylab="Population size", cex.axis=1.5, cex.lab=1.5, bty="n")
abline(h=extinction.threshold, lty=3, col="red")
plot(0:100, my.results.N[which(extinct)[2],], type="l", xlab="Time", ylab="Population size", cex.axis=1.5, cex.lab=1.5, bty="n")
abline(h=extinction.threshold, lty=3, col="red")


## TAREA 8 ##

# El modelo incluye los siguientes parámetros:
# "nzero" es el tamaño inicial de la población
# "r" es la media de la tasa de crecimiento poblacional
# "es" es la desciación estandar de la tasa de crecimiento poblacional y, por lo tanto,
#      representa la maginitud de la estocasticidad ambiental
# "K" es la capacidad de carga
# "Ks" es variación temporal estocástica en "K"
# "time" es el número de unidades de tiempo en las cuales se modelará la dinámica de la población.
# Cree la función *my.third.ricker.fun*:

my.third.ricker.fun <- function(nzero, r, es, K, Ks, time)
{
	N <- numeric(time+1)
	N[1] <- nzero
	r.es <- rnorm(time, r, es)
	K.es <- rnorm(time, K, Ks)
	for(i in 1:time)
	{
		N[i+1] <- N[i]*exp(r.es[i]*(1 - (N[i]/K.es)))
	}
	Time <- 0:time
	plot(Time, N, type="l", xlim=c(0, time), ylim=c(0,K+10), cex.axis=1.5, 
	    cex.lab=1.5, bty="n", lwd=2)
	abline(h=K, lty=3, col="red")
	return(N)
}

#use la función
my.third.ricker.fun(nzero=1, r=0.1, es=0.5, K=30, Ks=2, time=100)

#realice 1000 iterciones 
R.N <- matrix(NA, 1000,101)
for(j in 1:1000)
{
	R.N[j,] <- my.third.ricker.fun(nzero=1, r=0.1, es=0.5, K=30, Ks=2, time=100)
	#Sys.sleep(0.1)
}

#examine los resultados
R.N[1:5,1:5]

#grafique la dinámica de la población en cada iteración
plot(0:100, R.N[1,], type="n", ylim=range(R.N), xlab="Time", ylab="Population size", cex.axis=1.5, cex.lab=1.5, bty="n")
for(i in 1:1000)
{
	points(0:100, R.N[i,], type="l", ylim=range(my.results.N))
}
abline(h=30, lty=3, col="red")