learn.R

6. Constructing a Leslie Matrix

# If we want to project an age-structured population forward, given
# age-specific mortality probabilities and fertility rates, we need to
# construct a Leslie matrix. This is quite straightforward in R.

# for inputs, we take the lifetable lx or nLx values and the
# age-specific birth rates (typically in dughters -- two sex model
# later!)

# If your projection interval (=age class width) is 5 years or more,
# you probably want to use nLx. If it's 1 year, then lx values will
# probably suffice (this is typically how biologists do it)

# say the width of the projection interval is 5 (note I will continue
# to call the vector I use to calculate age-specific survival
# probabilities lx, even if it actually contains nLx values)

n <- 5

# To ensure that our projection matrix is irreducible, we include only
# those age classes younger than the age of last reproduction
# (=omega). We'll use as input lx and mx vectors of the same length.
# If you want to project the post-reproductive ages, this is a simple
# matter of taking the output from the projection and applying the
# survival probabilities to all ages > omega. The dynamics of these
# age classes will be completely determined by the dynamics of the
# connected parts of the life cycle together with the decrements of
# the unconnected parts.

# we have lx, but want the interval-specific survival probabilities =
# l(x+n)/l(x). We can use a trick of vectorized calculations to do
# this without a loop:

px <- exp(diff(log(lx)))

# make the first-row entries for the Leslie matrix. These
# "fertilities" are typically calculated as the survival-weighted
# average of fertility in age class i and i+1 reduced by the
# probability that an infant born at the outset of the interval will
# survive long enough to be counted as a live birth for that
# interval. The obvious thing to do is to use a for-loop to calculate
# these Fx values. We could do this as follows:

Fx <- NULL
for(i in 1:k-1) {
Fx[i] <- sqrt(lx[2]) * n*(mx[i] + px[i]*mx[i+1])/2
}
Fx <- c(Fx,mx[k])

# A more elegant way of doing this is to use a vectorized calculation:

Fx1 <- Fx # store the old Fx as something else, so you can compare it
mmax <- length(mx)
Fx <- sqrt(lx[2]) * n*(mx[mmax] + px*mx[-1])/2
Fx <- c(Fx,mx[k])

# Fx1 == Fx. Cool, eh?

# Note that the formula I use to calculate the Fx values is different
# from that used in Preston et al. The one I use comes from Caswell
# (2000). In practice, the two are very similar.

# You can usually safely assume that the sex ratio at birth (for
# humans) is approximately 1.05, meaning 105 male births per female
# birth. If you have information to the contrary for your population
# (e.g., China [Tuljapurkar et al. 1995]), use it!

SRB <- 1.05

# one-sex population

Fx <- Fx/(1+SRB)

# Now to assemble our leslie matrix. We first initialize a kxk matrix
# of zeros. Here we capitalize on R's recyling rule to make a matrix
# of zeros

A <- matrix(0, nrow=k, ncol=k)

# We want the survival probabilities on the sub-diagonal. Here is a
# nice way of using logical indexing to do that, again without a loop.
# The function row(A) returns a matrix of the same dimension as A with
# the ijth element equal to the the row number i. You can probably
# gues what col(A) does...

# The subdiagonal of a square matrix is where the column index is one
# higher than the row index

A[row(A) == col(A)+1] <- px

# stick in the first row

A[1,] <- Fx

# You now have a Leslie matrix!

# Now, to make thes commands useful beyond this one session, you'd
# want to assemble them into an R function which you could call
# anytime you wanted to create a Leslie matrix from given lx and mx
# schedules. If you're going to the trouble to create a Leslie
# matrix, chances are you'll need to do it again. Having a function
# will save you a lot of work in the final analysis.

Posted by jhj1 at August 19, 2005 10:01 PM

Comments

Post a comment