################################################################################
### R BASICS WORKSHOP ###
### EXERCISE 7.1: Generation of regular sequences (and a bit on indexing and ###
### graphics) ###
### ###
### Center for Conservation and Sustainable Development ###
### Missouri Botanical Garden ###
### Website: rbasicsworkshop.weebly.com ###
################################################################################
### INTRODUCTION ###############################################################
# This exercise is based on a macroecological model: the second version of the
# "interim general model", hereafter IGM2. This model attempts to predict
# geographic variation in broad-scale species richness of woody plants,
# as measured in areas of about 100 km x 100 km (Field et al. 2005, GLOBAL MODELS
# FOR PREDICTING WOODY PLANT RICHNESS FROM CLIMATE: DEVELOPMENT AND EVALUATION,
# Ecology 86: 2263–2277, available at the workshop's website under "PAPERS").
# The IGM2 proposes that species richness of woody plants (S) in any
# given area can be predicted from three properties of that area:
# annual rainfall (Ran), minimum monthly potential evapotranspiration
# (PETmin), and elevation range (TOPOG). The IGM2 is expressed in this equation:
# S = b0 + b1*Ran + b2*PETmin + b3*PETmin^2 + b4*(log(TOPOG))
# where b0, b1, b2, b3, and b4 are coefficients that define the effect of each
# predictor variable (Ran, PETmin y TOPOG) on woody plant richness. Field et al.
# estimated the values of these coefficients to predict woody plant richness
# anywhere in the world: b0 = -371, b1 = 0.2987, b2 = 5.1186, b3 = -0.0257, and
# b4 = 42.7155 (Table 1 in Field et al. 2005, available in the workshop page). To
# interpret the IGM2, it is critical to keep in mind that Ran y PETmin must be
# expressed in milimeters, and TOPOG in meters (see Table 1 en Field et al. 2005).
# In this exercise you will explore predictions derived from the IGM2 for the
# Neotropics.
## TASK 1 ##
## Use operator "<-" to assign to b0, b1, b2, b3 and b4 the values
## estimated by Field et al. (2005). Freebie: for the first coefficient type:
b0 <- -371
# this verifies that the value assigned to b0 is correct:
b0
## TASK 2 ##
## Explore how S changes with PETmin, when Ran is 0 mm and TOPOG is 1 m (i.e.,
## conditional on Ran = 0 mm and TOPOG = 1 m). Start by using operator "<-" to assign
## values to Ran and TOPOG. Then, use function "seq" to create a series of values
## for PETmin, from 0 mm to 160 mm (this is approximately the range of values in
## figure 5b of Field et al. (2005)), in intervals of 10 mm. Go and see the help
## page for function "seq" by typing "?seq" or "help(seq)". Last, use operator
## "<-" to assign values to S using the equation described in the introduction to
## Exercise 5: S = b0 + b1*Ran + b2*PETmin + b3*PETmin^2 + b4*(log(TOPOG)).
## TASK 3 ##
## Type "PETmin" and "S" in the console to see these two objects you just created.
## Then use the function "length" to see the length of S and PETmin. Use "?length"
## or "help(length)" to see the help page for function "length". Objects S and
## PETmin should have the same length.
## TASK 4 ##
## Use function "plot" to see the relationship between S and PETmin. Use "?plot"
## or "help(plot)" to see the help page for function "plot". Ensure that S is in
## the ordinate (vertical or y axis) and PETmin in the abscissa (horizontal axis).
## Compare your results with Figure 5b in Field et al. (2005). Do you have an
## opinion about the sign of the S values? Field et al. (2005) discuss the
## biological meaning of negative S values in the first complete paragraph in the
## right column of page 2266, and in the legend of Figure 5b. According to Field
## et al. (2005), negative S values inform about the increase in Ran or PETmin
## required for there to be at least one woody plant species in a given area.
## Recall that the S values you calculated and plotted are conditional on Ran =
## 0 mm (see TASK 3). In other words, you calculated S values for
## places similar to the Atacama Desert.
## TASK 5 ##
## Next examine predictions from the IGM2 for Ran values higher than zero. First
## use function "seq" to create a series of Ran values, from 0 mm to 5000 mm (this
## is the range of values used by Field et al. (2005)) in intervals of 1000 mm.
## Use function "length" to examine the length of Ran, and type "Ran" in the
## console to examine object Ran. Now create all possible combinations between
## the values of PETmin and Ran. There are several ways to achieve this. Try
## first this approach: use function "rep" to repeat PETmin a number of times
## equal to the length of Ran, assigning the result to an object called PETmin_r
## (using operator "<-"). If needed, go and see the help for function "rep" by
## typing "?rep" or "help(rep)". The next step is to use function "rep" to repeat
## each of the elements of Ran a number of times equal to the length of PETmin,
## and assign the result to an object named Ran_r. Note that, in Ran_r, each
## element of Ran must be repeated a number of times equal to the length of
## PETmin before the next element of Ran shows up for the first time. Make sure
## to use argument "each" of the function "rep". See the help page for "rep". Last,
## use operator "<-" to assign values to S using the equation described in the
## introduction to this exercise. Make sure to use Ran_r rather than Ran; and use
## PETmin_r instead of PETmin when you calculate S.
## TASK 6 ##
## Use functions "plot" and "points" to visualize the relationship between S,
## PETmin_r and Ran_r. First use "plot" to see the relationship between S and
## PETmin_r. Ensure that S is in the ordinate and PETmin_r in the abscissa. Next
## use function "points" to distinguish points that correspond to different values
## of Ran_r. Go and see the help page for function "points" by typing "?points"
## or "help(points)" in the console. Use function "points" to distinguish points
## that correspond to the relationship between S and PETmin_r conditional on
## Ran_r = 1000 mm. To select S values that correspond to Ran_r = 1000 mm, use
## square brackets (for indexing) this way:
S[Ran_r==1000]
## and to select PETmin_r values corresponding to Ran_r = 1000 mm use square
## brackets (for indexing) this way:
PETmin_r[Ran_r==1000]
## This use of square brackets for indexing is extremely useful, as we will see in
## some detail when we discuss "indexing" (see also indexing in "R for beginners"
## (section 3.5.4, page 26). Note you should use square brackets. Using round
## brackets for indexing is a mistake.
## TASK 7 ##
## In TASK 5, you created all possible combinations of PETmin and Ran values
## using function "rep". Now use function "expand.grid" to do the same task,
## and assign the result (which will be a "data frame") to an object named "PRcomb":
PRcomb <- expand.grid(PETmin, Ran)
## to see the result type in the console:
PRcomb
## use square brackets (for indexing) to see the first five rows of PRcomb:
PRcomb[1:5,]
## for objects that have two dimensions, like PRcomb, the numbers within the
## square brackets before the comma correspond to the rows to be extracted by
## indexing, and the numbers after the square brackets correspond to the columns
## to be extracted. Absence of numbers after the comma means that all columns
## will be extracted. The code above extracts rows 1 to 5 for all columns of
## PRcomb. The columns of PRcomb have names "Var1" and "Var2". Assign more
## informative titles using this code:
names(PRcomb)<-c("PETmin", "Ran")
## examine the result:
PRcomb[1:5,]
## to extract the first column of PRcomb type:
PRcomb[,1]
## or you could also use column names with the operator "$" to extract a particular
## column in a data fame (we will see this in more detail later today):
PRcomb$Ran
## to extract the second column of PRcomb write:
PRcomb[,2]
## or:
PRcomb$PETmin
## now, to calculate S use this code:
S <- b0 + b1*PRcomb[,2] + b2*PRcomb[,1] + b3*(PRcomb[,1])^2 + b4*(log(TOPOG))
## or:
S <- b0 + b1*PRcomb$Ran + b2*PRcomb$PETmin + b3*(PRcomb$PETmin)^2 + b4*(log(TOPOG))
## visually examine the relationship between S and PETmin using function "plot":
plot(PRcomb$PETmin, S, xlab="Minimum monthly potential evapotranspiration (mm)",
ylab="Woody plant species richness in 100 x 100 km")
## note the use of arguments "xlab" y "ylab" and their effect. Open the help page
## for function "plot" (type "?plot" in the console) and read the section on
## arguments "xlab" and "ylab". Now use square brackets and the function points to
## distinguish points in the relationship between PETmin and S that correspond to
## different values of 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)
## note the use and effect of arguments "col" and "pch". Open the help page for
## function "plot" (type "?plot" in the console) and read the section on arguments
## "col" and "pch". If you want to change colors write "colours()" in the console
## to see a list of more than 600 available colors. You can also change symbols
## using argument "pch".
## TASK 8 ##
# Examine predictions of the IGM2 for the relationship between S and PETmin for
# various TOPOG values, conditional on Ran = 0. Use function "seq" to create
# TOPOG values from 100 m to 3100 m (this is the approximate range of values
# used by Field et al. (2005)) in intervals of 500 m. Then use function
# "expand.grid" to create all combinations of PETmin and TOPOG values. Last,
# calculate S and visually examine the relationship between S and PETmin for
# different values of TOPOG using functions "plot" and "points". You may use
# TASK 7 as a guide. You need to create each object keeping in mind that
# you previously created a variety of objects. To see a list of the existing
# objects in the R session you can type in the console "objects()". You can
# remove objects from your session using function "remove" or "rm".
################################################################################
################################################################################
################################################################################
################################################################################
### TASK ANSWERS ###############################################################
################################################################################
################################################################################
################################################################################
################################################################################
## TASK 1 ##
b1 <- 0.2987
b1
b2 <- 5.1186
b2
b3 <- -0.0257
b3
b4 <- 42.7155
b4
## TASK 2 ##
Ran <- 0
TOPOG <- 1
PETmin <- seq(0, 160, 10)
S <- b0 + b1*Ran + b2*PETmin + b3*PETmin^2 + b4*(log(TOPOG))
S
## TASK 3 ##
PETmin
length(PETmin)
S
length(S)
## TASK 4 ##
plot(PETmin, S)
## TASK 5 ##
Ran <- seq(0, 5000, 1000)
length(Ran)
PETmin_r <- rep(PETmin, times=length(Ran))
Ran_r <- rep(Ran, each=length(PETmin))
#examine results
PETmin_r
Ran_r
cbind(PETmin_r, Ran_r) #function "cbind" may be used here to examine the two vectors
S <- b0 + b1*Ran_r + b2*PETmin_r + b3*PETmin_r^2 + b4*(log(TOPOG))
#examine result
S
## TASK 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)
## TASK 7 ##
# you already have all answers
## TASK 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]))
#or
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)