Introduction

Here are a few additional exercises to accompany the Rcompadre package vignettes found here. We will use the Rcompadre package in R for the following items.

1. Simple outputs from many matrices
   • Calculate population growth rate and damping ratio for COMPADRE/COMADRE data
2. Geographic distribution
   • Plot a subset of COMPADRE/COMADRE data on a map

Video tutorial

A video tutorial has been created to accompany this page and demonstrate these excercises.

Check it out here on the COM(P)ADRE Learn YouTube channel!

Preliminaries

For these exercises, we’ll use the Rcompadre, popbio, and maps packages. The packages popbio and maps are available throught CRAN. Instructions for Rcompadre installation can be found here.

library(Rcompadre)
library(popbio)
library(maps)

Walkthrough and R code

The R code for these exercises can be downloaded as an R file or PDF in addition to being displayed on this page.

Download as R file
Download as PDF

Simple outputs from many matrices

This example produces some basic outputs such as the population growth rate (λ) and damping ratio (Caswell 2001) for a subset of species and populations given some selection criteria.

Begin by loading the COMADRE database. The latest version can be found here or pulled using the Rcompadre package with the code below.

comadre <- cdb_fetch("comadre")

Important note: If using the cdb_fetch function, the database is automatically converted to an S4 object.

Once the database is loaded, we’ll use as_cdb (from the Rcompadre package) to make the loaded database into an S4 COMPADRE object. Then, we’ll subset the database to our data of interest: only mean matrices for bony fish from studies of 3 years in duration or longer, and with a matrix dimension of 3 or greater.

comadre <- as_cdb(comadre) # Converts the S3 database format from the website into S4

x <- subset(comadre, MatrixComposite == "Mean" & # Subsets the database based on variables in the metadata
              Class == "Actinopterygii" & 
              StudyDuration >= 3 &
              MatrixDimension > 3)

The object x is now a version of the COMADRE database object that contains only the matrices that match the search criteria.

These matrices can now be analyzed by applying functions in a loop, or by using lapply. For example, to calculate population growth rate and damping ratio for the subset of matrices, first create an empty data.frame to accommodate the output using popbio.

output <- data.frame(lambdas = rep(NA, length(x$mat)), damps = rep(NA, length(x$mat)))
for (i in 1:length(x$mat)){
  output$lambdas[i] <- Re(eigen(x$mat[[i]]@matA)$value)[1]
  output$damps[i] <- damping.ratio(x$mat[[i]]@matA)
}

Let’s examine the output.

output
#>     lambdas    damps
#> 1 1.5061504 1.303947
#> 2 0.7459078 1.432266
#> 3 0.9176391 1.397521
#> 4 1.3501445 1.697755
#> 5 2.3156105 1.303005
#> 6 0.9091956 1.031997
#> 7 0.9027406 1.027850
#> 8 0.9140835 1.035914
#> 9 0.9999592 1.005268

These parameters are presented in the same order as the metadata so we could also add the species name to this data.frame.

data.frame(Species = x@data$SpeciesAccepted, output)
#>                      Species   lambdas    damps
#> 1       Ammocrypta pellucida 1.5061504 1.303947
#> 2   Catostomus platyrhynchus 0.7459078 1.432266
#> 3                 Cottus sp. 0.9176391 1.397521
#> 4           Erimyzon sucetta 1.3501445 1.697755
#> 5      Hybognathus argyritis 2.3156105 1.303005
#> 6   Oncorhynchus tshawytscha 0.9091956 1.031997
#> 7   Oncorhynchus tshawytscha 0.9027406 1.027850
#> 8   Oncorhynchus tshawytscha 0.9140835 1.035914
#> 9 Sprattus sprattus balticus 0.9999592 1.005268

We can also plot the population growth rates and damping ratios derived from these matrices. In the plot below, the vertical, dashed red line indicates population growth rate = 1 (or log(λ) = 0).

par(mfrow = c(1,2))
  hist(log(output$lambdas),xlab = "Log population growth rate", col = "gold", main = "")
  abline(v=0,col = "red", lwd = 4, lty = 3)
  hist(output$damps, xlab = "Damping ratio", col = "brown", main = "")

Geographic distribution

This example produces on a world map with the viability (population growth rate λ > 1, λ = 1, λ < 1) of a subset of studied populations given some selection criteria and color-codes the location of each population according to the value of λ.

First, we’ll subset mean matrices for all Carnivora in the wild in the northern hemisphere, with no issues for survival > 1, for which matrices have been split into A = U + F + C, and for which reproduction was explicitly modeled.

x <- subset(comadre, MatrixComposite == "Mean" &
              Order == "Carnivora" &
              MatrixCaptivity == "W" &
              Lat > 0 &
              SurvivalIssue < 1 &
              MatrixSplit == "Divided" &
              MatrixFec == "Yes")

The object x is now a version of the COMADRE database object that contains only the matrices that match the search criteria. To calculate population growth rate for the subset matrices, we’ll first create an empty variable to accommodate the output from lambda calculations.

x@data$lambdas <- NA

Then, we’ll create a for loop to examine each matrix in turn. Here, it may be advisable to use the function tryCatch as a wrapper to cope with the situation if/when the function in the loop fails.

for (i in 1:length(x$mat)){
  tryCatch({
    x@data$lambdas[i] <- Re(eigen(x$mat[[i]]@matA)$value)[1]
    }, error = function(e){})
}

Now, we’ll create a vector of color hex codes that can be applied according to the estimate of λ. This is done using the colorRampPalette function to go from green for high values of λ to red for low values of λ. Here, paste is used to append a value of 90 to the hex codes to allow transparency for aesthetic reasons.

rampfunc <- colorRampPalette(c("green", "red"))
colVect <- rampfunc(100)
colVect <- paste(colVect, "90", sep = "")
s1 <- seq(min(x@data$lambdas, na.rm=TRUE),max(x@data$lambdas, na.rm=TRUE),
          length.out = 100)

Using the maps package, we’ll plot the world map and overlay the points from our data, color coded by value of λ. In this case, the points are jittered slightly to improve visibility of nearby populations.

par(mfrow = c(1,1))
{map("world", col = "gray", fill = TRUE, bg = "light blue",
    xlim = c(-175, 176), ylim = c(-60, 85), border = "white")
points(jitter(x@data$Lon, 0.6), jitter(x@data$Lat, 0.6),
       col = colVect[findInterval(x@data$lambdas, s1)], cex = 2, pch = 16)
}