demogR : A Package for the Construction and Analysis of Age-structured Demographic Models in R

The analysis of matrix population models has become a fundamental tool in ecology, conservation biology, and life history theory. In this paper, I present demogR , a package for analyzing age-structured population models in R . The package includes tools for the construction and analysis of matrix population models. In addition to the standard analyses commonly used in evolutionary demography and conservation biology, demogR contains a variety of tools from classical demography. This includes the construction of period life tables, and the generation of model mortality and fertility schedules for human populations. The tools in demogR are generally applicable to age-structured populations but are particularly useful for analyzing problems in human ecology. I illustrate some of the capabilities of the package by doing an evolutionary demographic analysis of several human populations.


Introduction
demogR is a package for doing demographic analyses of age-structured populations in R (R Development Core Team 2007). By default, it works with the data format typically used by human demographers. However, all of the routines for the analysis of matrix population models will work for any age-structured life cycle. In addition to the functions particularly of interest to ecologists, there are a variety of routines for classical formal demography of human populations.
In this paper, I will focus primarily on describing the functions associated with matrix population models. Matrix population models have proved especially important in two complemen-tary areas of ecology: (1) conservation biology and (2) evolutionary demography, including human evolutionary ecology (Hill and Hurtado 1996;Kaplan, Hill, Lancaster, and Hurtado 2000) 1. 1

. Conservation biology
The judicious use of matrix population models has become an indispensable tool for conservation biology (Lande 1988;Doak, Kareiva, and Klepetka 1994;Beissinger and Westphal 1998;Benton and Grant 1999). The rate of increase of a population is a critical parameter of interest in conservation biology, since a robust growth rate -and low variance in this growth rate -is the best insurance against extinction. The asymptotic growth rate of an age-structured population is the dominant eigenvalue of the projection matrix. This dominant eigenvalue, λ 1 is guaranteed to be unique, strictly greater than all other eigenvalues, and real if standard conditions for irreducibility and aperiodicity are met (Caswell 2001). Sensitivity and elasticity analyses of λ with respect to perturbations of vital events provide important insights into the relative importance of different vital rates for population viability and the potential impact of measurement error on the estimation of a population's rate of increase (Caswell 2001). The assumptions of the perturbation analysis that lie behind demographic elasticity analysis limit the uncritical use of elasticities in conservation biology. For example, real perturbations in a conservation context can be large and affect multiple vital rates simultaneously (Mills, Doak, and Wisdom 1999). In fact, a number of the routines in demogR are specifically designed for just these problems, including those for calculating the second derivatives of λ with respect to perturbations, and the sensitivities of eigenvalue elasticities.

Evolutionary demography
The asymptotic growth rate, λ, is also the fitness measure in a structured population with overlapping generations. Intense interest has focused on understanding how λ responds to perturbations in the elements of the life cycle (Caswell 1978;Benton and Grant 1996;Pfister 1998;de Kroon, van Groenendael, and Ehrlen 2000). Lande (1982) showed that the expected change in the mean life history phenotype (or any other phenotype) ∆z is given by where ∇λ is the fitness gradient (which is simply a vector of the sensitivities of λ) and G. is the additive genetic covariance matrix. The G matrix describes how a multivariate phenotype can respond to selection. Fisher (1958) showed that the rate of evolution is proportional to the variance in fitness. The G matrix generalizes this idea to multiple fitness components. These components all may have additive genetic variance, but they are also likely to covary with each other. In the simplest case where there are not (non-zero) off-diagonal elements and the diagonal elements of G are greater than zero, the life cycle traits will move in the direction of the gradient ∇λ which is the vector of sensitivities. When the covariances (i.e., off-diagonal elements) are non-zero, the response to selection will be more complex. The methods encompassed by matrix population models are therefore fundamental for understanding evolutionary change in life histories. There is currently some excellent work being done linking evolutionary quantitative genetics studies with ecology (e.g., Coulson, Benton, Lundberg, Dall, and Kendall 2006;Pelletier, Clutton-Brock, Pemberton, Tuljapurkar, and Coulson 2007).
In addition, Hamilton's perturbation analysis of the discrete-time characteristic equation for the human life history provided the canonical theory of senescence for the past 40 years (Hamilton 1966;Baudisch 2005;Steinsaltz, Evans, and Wachter 2005). The force of selection against deleterious mutations with age-specific effects declines with age and this decline is measured by the sensitivities of λ. The demographic methods used by demogR are therefore of central importance to the evolution of senescence and life histories in general and such analyses are increasingly being carried out in an ecological context on wild populations (Charmantier and Garant 2005;Nussey, Clutton-Brock, Albon, Pemberton, and Kruuk 2005;Charmantier, Perrins, McCleery, and Sheldon 2006;Nussey, Kruuk, Donald, Fowlie, and Clutton-Brock 2006)

Human evolutionary ecology
The observation that the asymptotic rate of increase of the population projection matrix is the proper fitness measure for an age-structured population places demography at the center of studies of human evolutionary ecology. The human life cycle is highly age-structured and there is little value to studies of human adaptation that fail to account for the variable force of selection on traits with respect to their age expression. Work in human evolutionary ecology has tended to focus on the demography of small-scale populations such as hunter-gatherers or horticulturalists (e.g., Howell 1979;Early and Peters 1990;Hill and Hurtado 1996). Understanding the evolution of human life histories, characterized as they are by long lifespan, modest fertility, extensive biparental care, overlapping periods of juvenile dependency, and long post-reproductive lifespan, remains a major challenge in human evolutionary ecology. An important area where human evolutionary ecology and conservation biology intersect is the analysis of subsistence and resource management in small populations (e.g. Alvard 1998;Smith and Wishnie 2000;Bird, Bird, and Parker 2005). Humans often play an instrumental role in regulating (and disturbing) communities and ecosystem processes and understanding the role that changes in the size and composition of human populations play is essential.

The need for a package for demographic analysis
The methods of analysis for matrix population models are widely used and a number of excellent texts are available to describe their use. Caswell (2001) provides Matlab code fragments and pseudo-code for many of the analyses described in his comprehensive book, while Morris and Doak (2002) provide extensive Matlab code to support their text. However, these functions are written for fairly expensive proprietary software, limiting their availability to all users. This is particularly problematic for scientists in lesser developed countries. Having a single package that collects many of the most commonly used analyses for age-structured population models, written for freely available software and distributed under the GPL solves this problem. demogR also allows population biologists and demographers to enjoy the many other advantages of the R programming environment including a more natural environment for handling complex data, advanced tools for statistical inference and superior graphics.

Life table methods
I will first discuss the construction of life tables from the census counts and death tabulations. Throughout the rest of the paper, I will use the vital event data from three countries: Venezuela (1965), Madagascar (1966, USA (1967). These are the three populations used by Goodman, Keyfitz, and Pullum (1974) in their classic paper on the demographic theory of kinship and they illustrate differing combinations of extremely high and low fertility and mortality. The dataset goodman is included in the package and consists of raw numbers used to construct the life tables and age-specific fertility schedules analyzed in the Goodman et al. (1974) paper. This dataset is a 19 × 10 data frame with columns representing age followed by the mid-year population size ( n K x ), number of deaths ( n D x ) and births (B x ) for Venezuela, Madagascar, and the USA for ages x ∈ 0, 1, 5, 10, . . . 80, 85. Following standard demographic notation, x is the age index and n indicates the width of the age interval. Thus, n K x is the midyear population for individuals age x to x + n years.
Construct a period life table with the routine life.table, which uses the Greville equation to convert observed death rates ( n M x = n D x / n K x ) into interval probabilities n q x . A period life table is a model for a hypothetical population. Through it, we apply period measures (i.e., cross-sectional measures taken at a specific time) and create a hypothetical cohort whose mortality experience we summarize. In effect, the period life table asks questions like: if current mortality conditions pertained to a particular population, what would the average age of death of that population be? If we had an actual cohort, then we would observe (cohort) age-specific death rates n m x = n d x / n L x , the observed number of deaths divided by the number of person-years in an age class. For the synthetic cohort of the period life table, we substitute the period death rates n M x , taking the mid-interval population size in age class x as our estimate for the number of person-years lived in the interval ( n L x ). As Greville showed, the only information needed to move between n M x and n q x is the average person-years lived by individuals dying in the interval n a x (Greville 1943). Given an estimate of the central death rate n M x , the conversion to n q x is given by: There are a variety of methods available for specifying the n a x schedule (Keyfitz 1977;Preston, Heuveline, and Guillot 2001). Probably the most common approach uses regression on the observed n M x values. In terms of mortality, the most rapidly changing portion of the human life cycle are the ages under five, and especially from birth to the first birthday. This general pattern of high and declining early mortality probably pertains to most mammals (Charnov 1993). For ages above five, we can typically assume that deaths are randomly distributed across the 5-year age interval, so that the average number of person-years lived by individuals dying at ages x > 5 is n a x = 2.5. Thus, the regression approach typically focuses on the first two n a x values (i.e., for age classes 0-1 and 1-4). The most common regression-based approach to graduating the n a x schedule is attributable to Keyfitz and Flieger (1990) and this is the default value for life.table.
In using regression-based graduation for the n a x schedules, human demographers capitalize on the remarkably robust nature of the human life cycle. Despite the great diversity of eco-nomic, social, and demographic environments that people experience, there are many general features of the human mortality experience. For many ecological applications, the detailed demographic information that allows human demographers to specify an n a x schedule by regression or borrowing from a similar population may be absent. It is common in ecological applications to simply construct a cohort life table (Deevey 1947;Caughley 1966). The assumptions underlying a cohort life table are not unreasonable if we assume that a population has lived in a relatively constant (average) environment. In the presence of non-stationary demographic change (e.g., in the case of demographic transitions in recent human populations), the cohort and period life tables will give very different pictures of the mortality experience of the population. Thus, for most applications to human ecology, the options type="cd" or type="kf" should be used, while for most non-human ecological applications, type="cohort" is probably the best option.
First, load the data and then calculate a period life tables using the default Keyfitz-Flieger graduation for the three populations. For the purposes of this paper, I will also limit the number of significant digits that are printed.
R> library("demogR") R> options(digits = 3, scipen = 1) R> data("goodman") R> names(goodman) [1] "age" "ven.nKx" "ven.nDx" "ven.bx" "mad.nKx" "mad.nDx" "mad.bx" [8] "usa.nKx" "usa.nDx" "usa.bx" R> ven <-with(goodman, life.table(x = age, nKx = ven.nKx, nDx = ven.nDx)) R> mad <-with(goodman, life.table(x = age, nKx = mad.nKx, nDx = mad.nDx)) R> usa <-with(goodman, life. The life table contains nine columns: age (x), n a x , the period central death rate ( n M x ), the interval mortality probability ( n q x ), survivorship (l x ), the proportion of deaths in the interval ( n d x ), person-years lived in the interval ( n L x ), person-years of remaining life in the cohort (T x ) and life expectancy (e x ). In most human demography -particularly for national-level data -the radix of the life table is typically set to l 0 = 100, 000, leading the the interpretation of l x as the number of survivors in the artificial cohort at age x and n d x as the number of deaths in the interval. Following ecological and evolutionary conventions, the default of life.table is set for a radix of l 0 = 1, making l x the probability of surviving to exact age x and n d x the proportion of deaths between x and x + n (= n m x l x ).
life.table can also use the Coale-Demeny approach to graduating the n a x schedule, which can be more useful in anthropological applications (i.e., those of small populations, frequently experiencing high mortality) that are more likely to be used in ecological or evolutionary analyses. The Coale-Demeny approach uses different values in the regression depending on the severity of early mortality (Coale, Demeny, and Vaughn 1983). In practice, the two different methods produce very similar life tables for the three national-level populations.
To demonstrate the usage of the cohort life table, I use Caughley's survival data on Himalayan thar (Hemitargus jemlahicus). These data are included in the dataframe thar.

Projection matrix construction
Having constructed the life tables, we are in a position to construct and analyze demographic projection matrices for the age-structured population. The routine leslie.matrix takes as arguments a measure of age-specific survivorship (either l x or n L x ) and the age-specific fertility rates, m x . life.table assumes that you want to construct a life table from data tabulated in a manner similar to goodman. That is, there are age-specific tabulations of mid-year populations size (or its equivalent), deaths, and births for ages x = 0, 1, 5, 10, . . .
The inputs for the Leslie matrix are a measure of survivorship (either l x or n L x ) and agespecific fertility m x . The default value is to use n L x , with corresponding argument L=TRUE. For national-level demographic data frequently employed in human demography, it is conventional to use person-years ( n L x ) to construct a projection matrix. This is because projections for national (or other large) populations are typically done with fairly coarse age intervals (e.g., 5 or 10 years). Since n L x = x+n x l(x)dx, n L x will give better estimates of survival for such large age-classes. It is assumed that both input vectors will be of the same length and that they will both contain a value for one year-olds which will need to be removed from the m x and l x vectors when L=FALSE and will need to be combined with 4 L 1 when L=TRUE.
Births are assumed to be continuously distributed throughout the year. That is, leslie.matrix assumes a birth-flow population. Following standard demographic conventions 2 (Keyfitz 1977;Preston et al. 2001), the Leslie matrix entry for the ith element of the first row is given by When the option L=FALSE is used, l(0.5) · n is substituted for n L x in equation 2, where l(0.5) = l(2) (Caswell 2001).
If data are not in the format expected by leslie.matrix, the utility function odiag can facilitate constructing a Leslie matrix. odiag takes two arguments: x and d. For vector argument x, odiag produces a square matrix with x along the off-diagonal given by the integer d. For matrix argument x, odiag extracts as a vector the dth off-diagonal. To place values along the sub-diagonal of a matrix, use d=-1.

Projection
Caswell (2001) notes that projection is the simplest form of analysis. Having constructed the Leslie matrix, and given a starting population vector, we can now project it forward some number of years and observe its behavior. Here we construct the Leslie matrix for the USA in 1967 and project it forward 100 years using the routine project.leslie, which takes as arguments the Leslie matrix, a starting population vector and a time horizon.
We now project the Leslie matrix with project.leslie using the census population of 1967 as the starting population size and plot the trajectory of the 10 age classes over the 100 year projection.
R> no <-goodman$usa.nKx The call to project.leslie take tmax=20 as an argument. However, since the width of the age classes in the Leslie matrix is five years, the projection actually spans 100 years, as shown

Eigenvalue analysis
The most common forms of analysis of the projection matrix are contained in the function eigen.analysis, which takes the projection matrix as its sole argument. The return value of eigen.analysis is a list containing: (1) lambda1: the dominant eigenvalue λ of the projection matrix, (2) rho: the damping ratio ρ = λ 1 /|λ 2 |, which is a measure of the rate of convergence to the stable age distribution, (3) sensitivities: a matrix of eigenvalue sensitivities, (4) elasticities: a matrix of eigenvalue elasticities, (5) stable.age: the stable age distribution u normalized so that u = 1, and (6) repro.value: reproductive value v normalized such that v 1 = 1. We see that the population of Venezuela in 1965 was increasing very rapidly: r = log(λ)/n = 0.0382, or nearly 4% annual increase. This is what happens when a population with a life expectancy at birth of nearly 70 also has a total fertility rate of of 6.4. Using the projection matrix A, we can calculate both the net reproduction number R 0 and the generation time T . Following Caswell (2001), we can decompose A into two components, A = T + F, where T is a matrix of transitions (simply the subdiagonal of an age-classified matrix) and F contains the fertilities. Using this decomposition, Caswell (2001) shows how to convert the projection matrix into a matrix of Markov transition probabilities for the states of the life cycle. Define the fundamental matrix of the Markov chain as N = (I − T) −1 , where I is an identity matrix of the same rank as A. The fundamental matrix N contains the expected number of visits to a transient state i prior to absorption (i.e., death) conditional on starting state j and is an object of interest in its own right in addition to its utility for calculating other quantities. Caswell shows that R 0 is the dominant eigenvalue of the matrix R = FN. For the ageclassified model, this is simply the r 11 th entry of R. The generation time is the amount of time it takes for a population growing at period rate λ to increase by a factor R 0 : λ T = R 0 , and can be found easily by solving this equation for T .

R> gen.time(A)
[1] 27.6 The average woman in Venezuela in 1965 replaced herself with almost three daughters and took approximately 27.6 years to do so.
The sensitivities of the growth rate, ∂λ/∂a ij , provide important demographic and evolutionary information. One possibility for calculating sensitivities is implicit differentiation of the characteristic equation (Hamilton 1966;Demetrius 1969). Caswell (1978) derived the sensitivities for a demographic projection matrix using a perturbation analysis of the matrix equation Au = λu. The calculation of sensitivities in this way only requires suitably scaled dominant left and right eigenvectors of the projection matrix. Caswell shows that ∂λ/∂a ij = v i u j , where v, u = 1, v is the dominant left eigenvector and u is the dominant right eigenvector. Elasticities are simply proportional sensitivities e ij = (a ij /λ)s ij and have the convenient property that i,j e ij = 1. We can compare the both the sensitivities and elasticities of the three populations in a convenient graphical manner using the plotting function plot.leslie.matrix. plot.leslie.matrix takes a projection matrix as its primary argument. It then plots the subdiagonal and first row of a square matrix on common axes. The figure uses the generic function plot to plot the fitness sensitivities (a-c) and elasticities (d-f) from the leslie.matrix objects.
R> par(mfrow = c(2, 3)) R> plot(Aea$sensitivities) R> title("(a)") R> plot(Bea$sensitivities) R> title("(b)") R> plot(Cea$sensitivities) R> title("(c)") R> plot(Aea$elasticities) R> title("(d)") R> plot(Bea$elasticities) R> title("(e)") R> plot(Cea$elasticities) R> title("(f)") The results are plotted in figure 2. Overall, sensitivities to pre-reproductive survival are highest, but there is a great deal of diversity in these three populations in the magnitude of the sensitivities to fertility. When the sensitivities are scaled to be proportional, we find that the general pattern is much more uniform across the three populations despite the tremendous diversity in vital rates.
As Stearns (1992) noted, sensitivities are "situational." As partial derivatives of the growth rate, they are conditional on all other values of the projection matrix remaining constant. However, there are at least three important contexts where the vital rates are likely to change. First in an evolutionary context, when the structure of G, the additive genetic covariance matrix for the life cycle, permits it, the life cycle transitions will respond to selection and change. Second in a conservation context, we may undertake specific policies to improve the survival and/or fertility of specific age classes (Crouse, Crowder, and Caswell 1987;Doak et al. 1994;Morris, Tuljapurkar, Haridas, Menges, Horvitz, and Pfister 2006). A successful management intervention will change vital rates. Third in the case of culture change, vital rates may change in response to exogenous drivers such as economic development or endogenous drivers such as ideational change (Borgerhoff Mulder 1998).
Given this interest, we can ask the question how will a small change to some matrix element a kl affect the other sensitivities? Sensitivities are partial derivatives. To investigate the change in sensitivities, we need to calculate eigenvalue second derivatives. (Caswell 1996) derived a formula for the second derivative of λ: where the superscript (i) indicates which eigenvalue/eigenvector (assuming they have been sorted in descending order of the real parts of the eigenvalues).
Calculate the second derivatives of λ to a perturbation of the klth element of A using secder, which takes as arguments the projection matrix A followed by the row index k and column index l of the perturbation.
We can see that selection of the survival of 0-4 year-olds is highly concave, which is a multivariate generalization of stabilizing selection (Caswell 1996). Since the second derivative of λ with respect to P 1 is negative (a standard definition of a concave function), selection acting on P 1 will reduce variance in it (see also Arnold and Wade (1984) for more discussion of concavity vs. convex selection). In contrast, selection on fertility prior to the mean age of childbearing T is convex, which is the multivariate generalization of directional selection. Selection on early fertility will increase variance in fertility. The convexity of selection on survival of ages after 0-4 is more complex. There is a sign change in age 20-24, which happens to be the age-class with peak fertility.
We may also want to ask the question for elasticities. We use the eigenvalue second derivatives to calculate the sensitivities of the elasticities according to the relationship: To do this, we use the function elassens, which takes the same arguments as secder.
[10,] 0.000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0. The function fullsecder computes all the second derivative of the dominant eigenvalue and returns them in a single matrix showing only the second derivatives with respect to non-zero matrix elements. The rows of this matrix, which will have a maximum size of 2k − 1 × 2k − 1 (where k is the rank of the Leslie matrix) are sorted such that all the fertilities are at the top. The matrix dimensions are labeled to facilitate easy access of the desired value of ∂ 2 λ/∂a ij ∂a kl . This application reveals a clear advantage to R over Matlab, as such dimension labeling is not possible in Matlab. This may seem, at first glance, like a small point but it greatly simplifies the interpretation of a very complex analysis. We can use this sorting to conveniently extract second derivatives for plotting. For instance, suppose that we want to know ∂ 2 λ/∂P i ∂P j for Madagascar and plot the 10 lines. R> cols <-rgb(0, (1:10)/10, (10:1)/10) R> plot (age,SD[11,11:20], type = "l", col = cols [1], xlab = "Age", + ylab = "secder") R> for (i in 12:20) lines (age,SD[i,11:20]

, col = cols[i -11])
The plot (figure 3) reveals that the value of ∂ 2 λ/∂P 2 i is negative for all ages, indicating concave selection, though for ages greater than 25, it is very nearly zero. ∂ 2 λ/∂P i ∂P j > 0 for j = i − 5, indicating convex selection. This indicates that selection on survival at any age class will reduce variance in it. A detailed discussion of these results for human life cycles is beyond the scope of the present paper, but is taken up in a forthcoming paper (Jones, in prep).
It is a general feature of the human life cycle that the largest fraction of total selection occurs on pre-reproductive survival (Jones 2005). However, elasticities of any of the life cycle transitions are conditional upon the other values of the life cycle. Survival is important for fitness only to the extent that it leads eventually to reproduction. van Groenendael, De Kroon, Kalisz, and Tuljapurkar (1994) presented a novel decomposition of the life-cycle using elasticities. They decompose the life cycle into its constituent loops (illustrated in figure 4). Using the properties of elasticities, they show that the total elasticity of the life cycle can be usefully decomposed into its loop elasticities. The algorithm presented by van Groenendael et al. (1994) is to identify all the loops of the life cycle, find the unique arc of each loop and multiply its elasticity, known as the loop's "characteristic elasticity," by the number of arcs in the loop. Wardle (1998) noted that identifying all the loops of a complex life cycle is not always a trivial task and presents a graph theoretic algorithm for enumerating all life cycle graphs. Fortunately, in the age classified life cycle, the unique arc of any loop will be the fertility arc of the loop. Thus the characteristic elasticity of any loop in an age-classified life cycle is simply the elasticity of the fertility arc in each loop.
With the function loop.elas, we can calculate and plot the loop elasticities of an age-classified projection matrix. loop.elas takes as arguments the projection matrix and the two optional arguments draw.plot (default=TRUE) and peryear (default=5).   presents the results. All three populations have a peak loop elasticity coinciding with the mean age of childbearing. The extent of this peak seems to be related to how spread out reproduction is around this mean. A full analysis of this pattern is again beyond the scope of this paper but is taken up in another forthcoming paper (Jones, in prep b).

Model life tables
Many parts of the world still lack vital event registration. This is particularly true for populations in the developing world that are of special interest to anthropologists and human behavioral ecologists. In addition to the lack of vital event registration, many populations of anthropological interest are quite small and modeling age-patterns of vital events (gathered, for example through genealogical surveys or archival reconstructions) can be challenging because the small number of events means that sampling variance can swamp out the demographic signal. Model life tables have been used to address both of these problems. The functions cdmltw, cdmlts, cdmltn, and cdmlte provide the Coale-Demeny families of model life tables (West, South, North, and East respectively). These four functions take one optional argument, sex, which defaults to sex=F. The output of the functions is a list containing the ages, age-class widths, and a series of 25 × 21 arrays of the various life table functions (25 levels of e 10 by 21 age classes): survivorship l x , interval mortality probability n q x , person-years lived by those dying in the interval n a x , deaths in the interval n d x , person-years lived in the interval n L x , central death rate n m x , person-years of remaining life T x , and life expectancy e x .
The code underlying these functions uses the coefficients provided in Coale et al. (1983)  Coale and Demeny's method for constructing the regional model life tables involved regression the observed age-specific probabilities of mortality in an interval x to x + n ( n q x ) on the observed life expectancy at age 10 (e 10 ) for the 192 input tables. These regressions were done by ordinary least-squares on both log 10 transformed and untransformed values of n q x . Plots of the predicted values of both these regressions against age will intersect twice, and in constructing the model life tables, Coale and Demeny used predicted n q x values from the simple regression for to the left of the first intersection and from the logarithmic regression to the right of the second intersection. Between the two, they used the arithmetic mean of the two predicted values. The n q x values are extended to ages 105 by assuming Gompertz mortality above age 80: l x = l 80 exp{−(µ(80)/k)(e k(x−80) − 1)}, x = 85, 90, . . . , 105 where k = log(µ(105) − µ(77.5))/27.5, µ(105) = 0.551 + 1.75( 5 q 75 ) for males, µ(105) = 0.613 + 1.75( 5 q 75 ) for females, µ(77.5) = l 75 ( 5 q 75 / 5 L 75 ), and µ(80) = µ(77.5) exp(2.5k). In all these equations µ(·) is the mortality hazard for the specified age The values of e 10 were chosen by an iterative procedure in such as way as to yield even 2.5year intervals of e 0 from 20 to 80 years. Male values of e 10 for a given model lifetable are a linear function of the female e 10 . Unfortunately, Coale and Demeny do not provide the coefficients for this relationship, but they can be recovered by copying the published e 0 values and regressing the values of e 10 for males onto the values of e 10 for females.
We can demonstrate the use of the Coale-Demeny model life tables with an example from the demography of hunter-gatherers. Howell (1979) pioneered the application of formal demographic methods to anthropological populations with her analysis of the Dobe !Kung of Botswana. Faced with a very small number of events and a population that did not keep track of their ages, Howell dealt with some serious challenges in constructing age-specific schedules of mortality and reproduction. She determined relative ages of the extant population and then assigned ages to as many people as possible using commonly known unusual events in the recent history of the people to triangulate individuals births and deaths with calendar years. She then used the Coale-Demeny West family model life table to smooth her estimated age distribution, which allowed her to estimate the age-specific mortality rate of the total population.
In their demographic study of the Ache of eastern Paraguay, Hill and Hurtado (1996) developed field techniques to allow them to avoid relying on model life tables for smoothing the age distribution. We can compare the !Kung and Ache l x curves to model model l x curves generated from the West family of life tables.
While this section is more specifically applicable to human demography, model life tables have been applied in nonhuman ecological contexts. For example, Dyke and colleagues have constructed a series of model life tables for nonhuman primates (Gage and Dyke 1988;Dyke, Gage, Ballou, Petto, Tardiff, and Williams 1993;Dyke, Gage, Alford, Swenson, and Williams-Blangero 1995). Caswell, Brault, Read, and Smith (1998)

Future work
Future development of demogR will focus on two primary areas: (1) incorporation of stochastic matrix models, and (2) curve-fitting to model age-schedules. For stochastic population models, I plan to write functions to estimate stochastic growth rate under independent and identically distributed (iid) and Markovian environments and calculate the sensitivities of these stochastic growth rates to perturbations of both the vital rates and, in the case of Markovian environments, the environmental transition probabilities . The curve-fitting for model schedules is more specifically related to human demography. Smoothing of empirical age distributions (particularly for small populations) using model life tables and observed growth rates is on particularly important application. A second will be to fit observed marriage transitions to a Coale-McNeil model nuptiality schedule (Coale and McNeil 1972). Another extension planned for the near future is a series of functions that provide other model life tables (e.g., UN, INDPETH) (United Nations 1981;Network 2002).