Introduction

COM(P)ADRE databases contain thousands of matrix population models. If we want to derive traits from a large set of these matrices, we’ll need to use either loops or vectorization.

Vectorizing means applying a function to each element of a vector. Here vector is defined broadly — it could be a sequence of character strings, a column of a data frame, a list of matrices, etc. Vectorized code generally runs faster than loops, and many R users find that vectorized code is easier to write and understand.

Preliminaries

We’ll start by loading a few packages and a dataset that we’ll be using throughout this vignette. The dataset Compadre is a subset of a recent COMPADRE release that’s built into Rcompadre.

library(Rcompadre)
library(popbio)
data(Compadre)

Introduction to vectorization

To understand vectorization, we first need a vector. For this purpose, we’ll extract a list of A matrices from the mat column of Compadre, and add this list of matrices to Compadre as a new column.

This new column, Compadre$matA, is both a vector and a list.

Let’s say we want to calculate the dimension of every matrix in Compadre$matA. In fact, Compadre already has this data in the column ‘MatrixDimension’, but let’s say we want to double-check it. We’ll use the function nrow(), and assume that the number of rows and columns are equal. But we can’t use nrow() directly on Compadre$matA, because the function nrow() isn’t vectorized. It can only take one object at a time.

Manual approach

We could do something like this…

Compadre$dim <- numeric(nrow(Compadre))  # create empty vector to store output
Compadre$dim[1] <- nrow(Compadre$matA[[1]])  # nrow matrix 1
Compadre$dim[2] <- nrow(Compadre$matA[[2]])  # nrow matrix 2
Compadre$dim[3] <- nrow(Compadre$matA[[3]])  # nrow matrix 3
# ... all the way to 150

But that’s not very efficient for 150 matrices.

Loop approach

A loop would be much more efficient here.

Vectorized approach

An even nicer approach is to vectorize the function nrow() over the vector Compadre$matA using sapply().

Compadre$dim <- sapply(Compadre$matA, nrow)

sapply() applies the function specified in the 2nd argument (nrow()) to every element of the vector in the first argument (Compadre$matA), and returns a vector of the results. The advantage of sapply() over the loop is that we don’t need to pre-define an object to store the results.

Vectorizing custom functions

We can also vectorize with a custom function. Let’s say we want to know, for every matrix, whether there are stages with no transitions (i.e. any column sums equal to zero). Here’s a vectorized approach.

The key to vectorizing is to make sure the function works on individual elements of the vector.

Accessor functions and vectorization

Note that, in the example above, it wasn’t necessary to create the column Compadre$matA before vectorizing over the A matrices. We could have simply used the matA() accessor within sapply().

Compadre$null_stages <- sapply(matA(Compadre), NullStages)

Using cdb_unnest() to avoid accessors

That said, using accessor funtions can get tedious. Rather than constantly using accessor functions to extract components of the mat column, we could use the function cdb_unnest() to extract separate columns for all matrix components at the start of our analysis.

Then we can refer to any component using $, e.g.

Other apply functions

vapply() is similar to sapply(), except that the output type is specified as an argument.

sapply(CompUnnest$matA[1:6], nrow)
#> [1] 4 6 3 5 5 3
vapply(CompUnnest$matA[1:6], nrow, numeric(1)) # must specify output type
#> [1] 4 6 3 5 5 3

lapply() always returns a list, so it’s useful if our output is more complex than a single value for each input. For example, we could use lapply to calculate vectors of stage-specific survival (column sums of matU).

mapply() is for vectorizing over multiple arguments. For example, the lifeExpectancy() function below (taken from the package Rage) calculates life expectancy given two arguments: a U matrix, and an integer indicator for the stage class reflecting the ‘start of life’. The start of life is often defined as the first ‘active’ stage class (i.e. not propagule or dormant), the index of which will vary from row to row.

When functions fail

Just like loops, vectorization fails if the function being vectorized throws an error on any element of the vector. Here’s an example. The lambda() function from popbio calculates the expected population growth rate given a projection matrix. It’ll work on most of the A matrices in Compadre, but fails on matrices that contain missing values (NA).

There are two basic approaches to overcoming this:

1. Remove or skip problem elements

If lambda() doesn’t work on matrices with missing values, one approach is to simply remove matrices with missing values. The cdb_flag() function is an easy way to check for missing values, and other common issues that might hinder our analyses.

Alternatively, if we want to avoid subsetting, we could pre-define a placeholder column for the result, and then selectively apply the lambda() function to only those matrices that don’t contain missing values.

Rows where there were missing values in matA retain the original placeholder value of NA.

2. Modify the function

A second approach is to modify the function we want to vectorize with so that it can natively handle special cases. For example, we might modify lambda() so that it returns NA if a matrix contains missing values.

If the special cases are harder to test for, we could use R’s condition-handling functions like try() or tryCatch(). Here’s an example.

This latter approach requires caution, as we’ll get an NA for any error. Some errors might reflect problems with our data or code that are fixable, in which case an NA may be misleading.