Adding variables to a PresenceAbsence object

Adding variables in raster format

To add variables in raster format to a PresenceAbsence object we can use the function lets.addvar from the letsR package. This function takes a raster object with any resolution and extent, and transform it to match the information in your PresenceAbsence object. Once this is done, the variables are included as additional columns containing the aggregate/summarize value of the variable(s) in each cell of the presence-absence matrix. Let’s see an example using the bioclimatic data from WorldClim that can be downloaded with the package raster.

library(letsR)

First we get the data.

r <- getData("worldclim", var = "bio", 
             res = 10)
plot(r)

Here I will use the PresenceAbsence object for Phyllomedusa species generated previously.

data(PAM)
plot(PAM, main = "Phyllomedusa\nRichness")

Let’s check the projection differences.

projection(r)

## [1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

projection(PAM$Richness_Raster)

## [1] "+proj=longlat +datum=WGS84"

Since the CRS descriptions differ, we have to fix them first. In this case the projections are actually the same, but the character describing them differ, so we can just change the name directly. But if the projections are actually different, you may want to use the function raster::projectRaster.

projection(PAM$Richness_Raster) <- projection(r)

We can now run the lets.addvar function.

PAM_env <- lets.addvar(PAM, r, fun = mean)

The climatic data have a higher resolution than our PAM. In this case, we could choose a function to aggregate the values with the argument fun. In most of the situations, people will be interested in averaging values to aggregate multiple cells, but in some specific cases you may want to sum them, or get the standard deviation, or use any another function.

The result is a presence absence matrix, where the last columns now include the raster values. Check the table:

head(PAM_env)

If you do not want the coordinates and species included you can set the argument onlyvar = TRUE.

climate <- lets.addvar(PAM, r, fun = mean, onlyvar = TRUE)

head(climate)

Now that we have the variables, we can use it to relate to our species data in many ways. For example, you could graph the relationship between precipitation/temperature and species richness.

library(ggplot2)

rich <- rowSums(PAM$P[, -(1:2)])

mpg1 <- data.frame("Temperature" = climate[, 1]/10,
                   "Precipitation" = climate[, 12],
                   "Richness" = rich)
f1 <- ggplot(mpg1, aes(Temperature, Richness))
f1 <- f1 + geom_smooth(model = lm) + 
  geom_point(col = rgb(0, 0, 0, .6)) + 
  theme_bw()

f2 <- ggplot(mpg1, aes(Precipitation, Richness))
f2 <- f2 + geom_smooth(model = lm) + 
  geom_point(col = rgb(0, 0, 0, .6)) + 
  theme_bw()

f1

f2

Adding variables in polygon format

Data in shapefile format like ecorregions, conservation units or countries, can be added to a PAM using the function lets.addpoly. This function adds polygons’ attributes as columns at the right-end of the matrix. The values represent the percentage of the cell covered by the polygon attribute used. As an example, we can use the South America countries map available in the package maptools.

library(maptools)

data("wrld_simpl")
SA <- c("Brazil", "Colombia",  "Argentina",
        "Peru", "Venezuela", "Chile",
        "Ecuador", "Bolivia", "Paraguay",
        "Uruguay", "Guyana", "Suriname",
        "French Guiana")
south_ame <- wrld_simpl[wrld_simpl@data$NAME %in% SA, ]
plot(south_ame)

Now we can add this variables to our PAM.

PAM_pol <- lets.addpoly(PAM, south_ame, "NAME")

head(PAM_pol)

This information can be used to calculate the number of species per country for example.

vars_col <- (ncol(PAM$P) + 1):ncol(PAM_pol)
n <- length(vars_col)
rich_count <- numeric(n)
for (i in 1:n) {
  rich_count[i] <- sum(colSums(PAM$P[PAM_pol[, vars_col[i]] > 0,
                                     -(1:2)]) > 0)
}
labs <- as.factor(colnames(PAM_pol)[vars_col])
names(rich_count) <- labs

mpg <- data.frame("Richness" = rich_count, "Country" = as.factor(labs))
g <- ggplot(mpg, aes(labs, Richness))
g + geom_bar(stat = "identity") + labs(x = "") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

To cite letsR in publications use: Bruno Vilela and Fabricio Villalobos (2015). letsR: a new R package for data handling and analysis in macroecology. Methods in Ecology and Evolution. DOI: 10.1111/2041-210X.12401