################################################################################
### TALLER FUNDAMENTOS DE R                                                  ###
### EJERCICIO 7.1: Sequencias regulares (y un poco sobre graficos)           ###
###                                                                          ###
### Center for Conservation and Sustainable Development                      ###
### Missouri Botanical Garden                                                ###
### Website: rbasicsworkshop.weebly.com                                      ###
################################################################################

### INTRODUCCION ###############################################################
# La segunda versión del "interim general model" (de aquí en adelante "IGM2"),  
# es un modelo de variación geográfica en la riqueza de especies de plantas 
# leñosas a escalas grandes (medido en áreas de alrededor de 100 km x 100 km; 
# Field et al. 2005, GLOBAL MODELS FOR PREDICTING WOODY PLANT RICHNESS FROM 
# CLIMATE: DEVELOPMENT AND EVALUATION, Ecology 86). El IGM2 propone que la 
# riqueza de especies de plantas leñosas (S) en cualquier zona determinada se 
# puede predecir a partir de tres propiedades de esa zona: precipitación anual 
# (Ran), mínimo mensual de evapotranspiración potencial (PETmin), y el rango de 
# elevación (TOPOG). El IGM2 se expresa en la siguiente ecuación:

# S = b0 + b1*Ran + b2*PETmin + b3*PETmin^2 + b4*(log(TOPOG))

# donde b0 es el intercepto, y b1, b2, b3, y b4 son coeficientes que definen el
# efecto de cada variable predictora (Ran, PETmin y TOPOG) sobre la riqueza de
# plantas leñosas. Field et al. calcularon los valores de estos coeficientes para
# predecir la riqueza de especies en cualquier parte del mundo: b0 = -371, 
# b1 = 0,2987, b2 = 5,1186, b3 = -0.0257, y b4 = 42,7155 (Tabla 1 en Field et al.,
#  2005, disponible en la página del taller). Para interpretar el IGM2, es
# fundamental tener en cuenta que Ran y PETmin debe ser expresadas en milímetros,
# y TOPOG en metros (véase la Tabla 1 en Field et al., 2005). En este ejercicio
# usted explorará las predicciones derivadas del IGM2 para el Neotrópico.


## TAREA 1 ##
## Use el operador "<-" para asignar a los objetos de nombre b0, b1, b2, b3 y 
## b4 los valores estimados por Field et al. (2005). Por ejemplo, para el primer 
## coeficiente:
b0 <- -371

# Esto verifica que el valor asignado a b0 es correcto:
b0


## TAREA 2 ##
## Explore como S cambia con PETmin, cuando Ran es 0 mm y TOPOG es 1 m (es 
## decir, condicionado a que Ran = 0 mm y TOPOG = 1 m). Empieze por usar el 
## operador "<-" para asignar valores a Ran y TOPOG. Después utilice la 
## función "seq" para crear una serie de valores para PETmin, de 0 mm a 160 mm en 
## intervalos de 10 mm (esto es aproximadamente el rango de valores en la figura 
## 5b de Field et al. (2005)). Vea la ayuda para la función "seq". Por 
## último, use el operador "<-" para asignar valores a S utilizando la 
## ecuación descrita en la introducción de este ejercicio: 

## S = b0 + b1*Ran + b2*PETmin + b3*PETmin^2 + b4*(log(TOPOG))


## TAREA 3 ##
## Escriba "PETmin" y "S" en la consola para ver el contenido de estos objetos 
## que acaba de crear. A continuación, utilice la función "length" para ver la 
## longitud de S y PETmin. Los objetos S y PETmin deben tener la misma longitud.


## TAREA 4 ##
## Utilice la función "plot" para ver la relación entre S y PETmin. Asegúrese de 
## que S esté el eje de ordenadas (eje vertical o y) y PETmin en el eje de las 
## abscisas (eje horizontal o x). Compare sus resultados con la figura 5b en 
## Field et al. (2005). Tiene una opinión sobre el signo de los valores de S? 
## Field et al. (2005) discuten el significado biológico de los valores 
## negativos de S en el primer párrafo completo en la columna izquierda de la 
## página 2266, y en la leyenda de la figura 5b. Según Field et al. (2005), los 
## valores negativos de S informan sobre el aumento en Ran o PETmin necesario 
## para que haya al menos una especie de plantas leñosas en un área dada.
## Recuerde que los valores de S que usted calculó están condicionados en Ran =
## 0 mm. En otras palabras, se calcularon los valores de S en lugares similares
## al desierto de Atacama.


## TAREA 5 ##
## A continuación, examine las predicciones de la IGM2 para valores de Ran 
## superiores a cero. Primero, use la función "seq" para crear una serie de 
## valores para Ran de 0 mm a 5000 mm en intervalos de 1000 mm (este es el rango
## de valores utilizados por Field et al. (2005)) . Utilize la función "length" 
## para examinar la longitud de Ran, y escriba "Ran" en la consola para examinar 
## el contenido de este objeto. Ahora cree todas las combinaciones posibles 
## entre los valores de PETmin y Ran.

Ran <- seq(0, 5000, 1000)
length(Ran)

PETmin_r <- rep(PETmin, times=length(Ran))  
Ran_r <- rep(Ran, each=length(PETmin))
PETmin_r
Ran_r
cbind(PETmin_r, Ran_r) # La función "cbind" puede usarse para ver los vectores juntos

S <- b0 + b1*Ran_r + b2*PETmin_r + b3*PETmin_r^2 + b4*(log(TOPOG))
S


## TAREA 6 ##
## Use las funciones "plot" y "points" para visualizar la relación entre S,
## PETmin_r y Ran_r. Use primer "plot" para ver la relación entre S y
## PETmin_r. Asegúrese de que S se encuentra en el eje de ordenadas y PETmin_r 
## en la abscisa. Luego, use la función "points" para distinguir los puntos que 
## corresponden a diferentes valores de Ran_r. Vea la ayuda para la función de 
## "points". Utilize la función "points" para distinguir puntos que corresponde 
## a la relación entre S y PETmin_r condicional en Ran_r = 1,000 mm. Para 
## seleccionar los valores de S que corresponden a Ran_r = 1000 mm, use
## corchetes (para la indexación) de esta manera:

S[Ran_r==1000] # Selecciona los valores en "S" que corresponden a "Ran_r" igual a 1000

## Y para seleccionar los valores correspondientes a PETmin_r donde Ran_r = 1000 mm 
## use los corchetes cuadrados (para la indexación) de esta manera:

PETmin_r[Ran_r==1000]

## TAREA 7 ##
## En la tarea 5, usted creó todas las combinaciones posibles entre los valores 
## PETmin y Ran utilizando la función de "rep". Ahora, use la función 
## "expand.grid" para hacer la misma tarea, y asigne el resultado a un objeto 
## denominado "PRcomb" (que será un marco de datos o data.frame) 


## Para ver el tipo de resultado, escriba en la consola:
PRcomb <- expand.grid(PETmin, Ran)

## Utilice corchetes (para indexación) para ver las primeras cinco filas de 
## "PRcomb":
PRcomb[1:5,]

## Las columnas de PRcomb tienen nombres "Var1" y "var2". Asigne nombres de
## columnas más informativos utilizando este código:
colnames(PRcomb)<-c("PETmin", "Ran")

## Examine el resultado:
PRcomb[1:5,]

## Ahora, para calcular S utilizar este código:
S <- b0 + b1*PRcomb[,2] + b2*PRcomb[,1] + b3*(PRcomb[,1])^2 + b4*(log(TOPOG))

## O:
S <- b0 + b1*PRcomb$Ran + b2*PRcomb$PETmin + b3*(PRcomb$PETmin)^2 + b4*(log(TOPOG))

## Examine visualmente la relación entre S y PETmin usando la función "plot":
plot(PRcomb$PETmin, S, xlab="Minimum monthly potential evapotranspiration (mm)",
    ylab="Woody plant species richness in 100 x 100 km")

## Tenga en cuenta el uso de los argumentos "xlab" "ylab" y, y su efecto. Abra 
## la página de ayuda para la función de "plot" (escriba "?plot" en la consola) 
## y lea la sección sobre argumentos "xlab" y "ylab". Ahora utilice corchetes y 
## la función "points" para distinguir la relación entre PETmin y S que 
## corresponden a diferentes valores de PRcomb$Ran:

points(PRcomb$PETmin[PRcomb$Ran==0], S[PRcomb$Ran==0], col="red", pch=19)
points(PRcomb$PETmin[PRcomb$Ran==1000], S[PRcomb$Ran==1000], col="mistyrose4", pch=19)
points(PRcomb$PETmin[PRcomb$Ran==2000], S[PRcomb$Ran==2000], col="lightseagreen", pch=19)
points(PRcomb$PETmin[PRcomb$Ran==3000], S[PRcomb$Ran==3000], col="mediumpurple", pch=19)
points(PRcomb$PETmin[PRcomb$Ran==4000], S[PRcomb$Ran==4000], col="springgreen3", pch=19)
points(PRcomb$PETmin[PRcomb$Ran==5000], S[PRcomb$Ran==5000], col="springgreen", pch=19)

## Tenga en cuenta el uso y efecto de argumentos "col" y "pch". Abra la página 
## de ayuda para la función "plot" y lea la sección sobre estos argumentos. Si 
## desea cambiar los colores escriben "colors()" en la consola para ver una 
## lista de más de 600 colores disponibles. En la página web del taller bajo
## Resources también hay una lista de colorres. También puede cambiar los 
## símbolos en la figura utilizando el argumento "pch".


## TAREA 8 ##
# Examine las predicciones de IGM2 para la relación entre S y PETmin para 
# diversos valores TOPOG condicionadas a Ran = 0. Utilice la función "seq" para 
# crear valores TOPOG desde 100 m hasta 3100 mm en intervalos de 500 m (este es 
# el rango aproximado de los valores utilizado por Field et al. (2005)). Luego, 
# utilice la función "expand.grid" para crear todas las combinaciones posibles 
# entre los valores PETmin y TOPOG. Finalmente, calcule S y visualmente examine 
# la relación entre S y PETmin para diferentes valores de TOPOG usando las 
# funciones "plot" y "points". Puede utilizar la TAREA 7 como una guía. Es 
# necesario crear cada objeto teniendo en cuenta que creó anteriormente una 
# variedad de objetos. Para ver una lista de los objetos existentes en la sesión 
# de R se puede escribir en la consola "objects()". Usted puede eliminar objetos 
# de la sesión utilizando las funciones "remove" o "rm".







































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

## TAREA 1 ##
b1 <- 0.2987
b1
b2 <- 5.1186
b2
b3 <- -0.0257 
b3
b4 <- 42.7155
b4

## TAREA 2 ##
Ran <- 0
TOPOG <- 1
PETmin <- seq(0, 160, 10)

S <- b0 + b1*Ran + b2*PETmin + b3*PETmin^2 + b4*(log(TOPOG))
S


## TAREA 3 ##
PETmin
length(PETmin)

S
length(S)


## TAREA 4 ##
plot(PETmin, S)


## TAREA 5 ##
#ya tiene la respuesta


## TAREA 6 ##
plot(PETmin_r, S)
points(PETmin_r[Ran_r==1000], S[Ran_r==1000], col="orange", pch=19)
points(PETmin_r[Ran_r==0], S[Ran_r==0], col="red", pch=19)
points(PETmin_r[Ran_r==2000], S[Ran_r==2000], col="orange2", pch=19)
points(PETmin_r[Ran_r==3000], S[Ran_r==3000], col="green", pch=19)
points(PETmin_r[Ran_r==4000], S[Ran_r==4000], col="green2", pch=19)
points(PETmin_r[Ran_r==5000], S[Ran_r==5000], col="blue", pch=19)


## TAREA 7 ##
# Ya tiene las respustas


## TAREA 8 ##
objects()

Ran <- 0
TOPOG <- seq(100, 3100, 500)

PTcomb <- expand.grid(PETmin, TOPOG)
PTcomb
PTcomb[1:5,]

names(PTcomb) <- c("PETmin", "TOPOG")
PTcomb[1:5,]

S <- b0 + b1*Ran + b2*PTcomb[,1] + b3*(PTcomb[,1])^2 + b4*(log(PTcomb[,2]))
# o
S <- b0 + b1*Ran + b2*PTcomb$PETmin + b3*(PTcomb$PETmin)^2 + b4*(log(PTcomb$TOPOG))

plot(PTcomb$PETmin, S, xlab="Minimum monthly potential evapotranspiration (mm)", 
    ylab="Woody plant species richness in 100 x 100 km")
points(PTcomb$PETmin[PTcomb$TOPOG==100], S[PTcomb$TOPOG==100], col="lightseagreen", pch=19)
points(PTcomb$PETmin[PTcomb$TOPOG==600], S[PTcomb$TOPOG==600], col="mistyrose4", pch=19)
points(PTcomb$PETmin[PTcomb$TOPOG==1100], S[PTcomb$TOPOG==1100], col="red", pch=19)
points(PTcomb$PETmin[PTcomb$TOPOG==1600], S[PTcomb$TOPOG==1600], col="purple", pch=19)
points(PTcomb$PETmin[PTcomb$TOPOG==2100], S[PTcomb$TOPOG==2100], col="springgreen", pch=19)
points(PTcomb$PETmin[PTcomb$TOPOG==2600], S[PTcomb$TOPOG==2600], col="yellow", pch=19)
points(PTcomb$PETmin[PTcomb$TOPOG==3100], S[PTcomb$TOPOG==3100], col="cyan", pch=19)