################################################################################
### 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, b1, b2, b3, b4 y son coeficientes que definen el efecto de cada
# variable predictora (Ran, PETmin y TOPOG) sobre la riqueza de plantas leñosas. 
# Field et al. calculó 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 el ejercicio 5 usted 
# explorará las predicciones derivadas de la IGM2 para el Neotrópico.


## TAREA 1 ##
## Use el operador "<-" para asignar a objetos con los nombres b0, b1, b2, b3 y 
## b4 los valores estimada por Field et al. (2005). Freebie: 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 de 0 mm y TOPOG es de 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. A continuación, utilice la 
## función "seq" para crear una serie de valores ara 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)), . Vaya a ver la ayuda para la función "seq". Por 
## último, el use el operador "<-" para asignar valores a S utilizando la 
## ecuación descrita en la introducción de 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. Use "?length" o "help(length)" para ver la ayuda de 
## 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 Campo et al. (2005), los 
## valores negativos 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 (véase la tarea 3). En otras palabras, se calcularon los valores de S
## en lugares similares a el 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. Hay varias maneras de lograr esto. Primero,
## trate este enfoque: use la función "rep" para repetir PETmin un número de 
## veces igual a la longitud de Ran, asignando el resultado a un objeto llamado 
## "PETmin_r" (usando el operador "<-"). Si es necesario, vea la ayuda para 
## "rep". El siguiente paso es utilizar la función de "rep" para repetir CADA UNO 
## de los elementos de Ran un número de veces igual a la longitud de PETmin, y 
## asignar el resultado a un objeto denominado "Ran_r". Tenga en cuenta que en 
## "Ran_r" cada elemento de Ran debe repetirse un número de veces igual a la 
## longitud de PETmin antes de que el siguiente elemento de Ran aparezca por 
## primera vez. Asegurarse de usar el argumento "each" de la función "rep". 
## Ver la página de ayuda para "rep". Finalmente, use "<-" para asignar valores 
## a S utilizando la ecuación descrita en el Introducción. Asegúrese de utilizar 
## "Ran_r" en lugar de Ran; y use "PETmin_r" en lugar de PETmin cuando.


## 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]

## Este uso de corchetes para la indexación es extremadamente útil, como veremos 
## en más detalle cuando hablemos de "indexación" (véase también indexación en 
## "R para principiantes", Sección 3.5.4, página 26). Tenga en cuenta que debe 
## utilizar corchetes "[]". El uso de parentesis redondos "()" para indexación 
## es un error.


## 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

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

## Para objetos que tienen dos dimensiones, como PRcomb, los números dentro de 
## los corchetes antes de la coma corresponden a las filas que quieren ser extraídas 
## por indexación, y los números después de la coma corresponden a las columnas
## a ser extraidas. La ausencia de números después de la coma significa que 
## todas las columnas seran extraídas. El código anterior extrae las filas 1 a 5 
## para todas las columnas de PRcomb. 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,]

## Para extraer la primera columna de PRcomb:
PRcomb[,1]

## O también puede utilizar los nombres de columna con el operador "$" para 
## extraer una determinada columna en un margo de datos (veremos esto con más 
## detalle más adelante):
PRcomb$PETmin

## Para extraer la segunda columna de PRcomb de escritura:
PRcomb[,2]

## O:
PRcomb$Ran

## 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 examinan 
# 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 ##
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 ##
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)