Overview

Session 2, “Population Dynamics” suggests ‘practice’ generating life tables and proffers some “sample” data. We use a couple of references to establish the metrics:

Additionally, the presentation Population parameters by Amber Sherman is helpful.

We expand the female antelope data given the number of female antelope surviving up to age thirteen, the number of female individuals and their age-specific fecundity schedule into a “full” life table and generate a number of plots to visualise the data.

Life Table Variables

From the references we have collated the following list of variables:

  1. \(x\): A categorical age factor (sometimes specified as a range).
  2. \(n_x\): Number of female antelope alive at start of age \(x\)
  3. \(m_x\): Fecundity rate - reproductive output of cohort at each age \(x\)
  4. \(l_x\): Standardised proportion of survivors at the start of age \(x\)
  5. \(d_x\): Number dying during the age interval \((x,x+1]\)
  6. \(q_x\): Mortality rate - the intensity of deaths at each age \(x\)
  7. \(L_x\): Mean number of survivors at age \(x\)
  8. \(e_x\): Average remaining life expectancy at age \(x\)
  9. \(\bar{e}_x\): Average life expectancy of individuals of age \(x\)
  10. \(k_x\): The “Killing Power” - the standardised (that is, log10 converted) intensity of mortality
  11. \(l_x \times m_x\): Mean number of female offspring produced by females at age \(x\)
The code for producing these values is listed at the end of this document. With these metrics, we then work with the table:
x nx mx lx dx qx Lx ex avgex kx lxmx
0 498 0.00 1.000 164 0.329 416 2.761 2.761 0.173 0.000
1 334 0.10 0.671 186 0.557 241 2.871 3.871 0.353 0.067
2 148 0.15 0.297 53 0.358 121 4.851 6.851 0.193 0.045
3 95 0.30 0.191 7 0.074 91 6.284 9.284 0.033 0.057
4 88 0.44 0.177 2 0.023 87 5.750 9.750 0.010 0.078
5 86 0.46 0.173 2 0.023 85 4.872 9.872 0.010 0.080
6 84 0.42 0.169 5 0.060 81 3.976 9.976 0.027 0.071
7 79 0.43 0.159 2 0.025 78 3.203 10.203 0.011 0.068
8 77 0.39 0.155 4 0.052 75 2.273 10.273 0.023 0.060
9 73 0.35 0.147 20 0.274 63 1.370 10.370 0.139 0.051
10 53 0.25 0.106 42 0.792 32 0.698 10.698 0.683 0.026
11 11 0.14 0.022 11 1.000 5 0.455 11.455 Inf 0.003
12 0 0.00 0.000 0 NaN 0 NaN NaN NaN 0.000

Additional variables require post-processing against the table - these, with their values for our sample data, are:

  1. \(R_0=0.606\): Net reproductive rate per generation - mean number of female offspring produced per female over their lifetime.
  2. \(G=5.26\): Average age a female gives birth.
  3. \(\lambda=0.909\): Rate of population change - logistic curve rate parameter.
  4. \(r=-0.095\): Change in population size per individual per unit of time.
  5. \(c_x\): Stable age distribution - the proportion of the population in each age class over time. This is a vector - we plot it in a chart shortly.

Interpretation

Clearly, from these values (\(R_0<1\) and \(r<0\)), the cohort population is in decline! Although we can evaluate \(\lambda\) and thus \(c_x\), we cannot plot a logistic population growth curve since we are not given the base data regarding the cohort’s population. Moreover, as we have no data regarding the gender split, we cannot produce a population pyramid. As for the balance of the variables, we can construct a number of plots:

Proportional Mortality, \(q_x\), Survival, \(l_x\) and killing power, \(k_x\) for the Fecundity of \(m_x\)

For reference, the chart shows \(G\) - the mean age at which females first give birth. Technically, it’s not good practice to mix log10 with other data values on the \(y\)-axis, but we show \(k_x\) nonetheless. Given the difference in scale it may not be appropriate to strongly infer that killing power, \(k_x\) is correlated to mortality, \(q_x\) despite the proximity of the chart but, bear in mind, they are both covariates dependent up \(n_x\). Other than that, the curves, as one would expect from the derivation of the metric as fairly self explanatory - once visualised!

Life Expectancy

The trace for remaining life expectancy appears to be more obvious in its inference; that for average life expectancy at a particular age is, to my mind, less ‘useful’ - especially in stating, essentially, that your average life expectancy is is observed at your death!

Age distribution and Survivor count

I understand that aspects of the stable age distribution are prerequisites of producing a population (logistic) growth plot - however, as we don’t have the data to construct such a plot, I do not otherwise see any benefit in plotting this metric. The survivor count in this plot is “non-standardised” - that is, the trace cannot be compared direly against that for another cohort or species. Consequently, this plot also appears of little use - rather, we produce the standardised (log10 converted) survivorship in the following chart:

Standardised Survival Curve

In this plot the survivorship, \(L_x\), is plotted using a log10 transformation which enables its cross-cohort or cross-species comparison. For completeness I’ve added curves representing survivorship for a Type I, II and III species; this is rather contrived as, to evaluate the points for a Type I species, I have simply taken the logarithm of survivorship binned equally (I used 1000 data points for a smooth curve) across the range of our logged \(L_x\) value. The Type II curve is just a straight line joining the endpoints of the Type I curve. The Type II curve is obtained by the rotation of the Type I curve through \(\pi\) radian, then a linear translation to the endpoints of the Type I curve.

Regardless of the contrivance, it does show lower and upper bounds for our actual data. Moreover, we may infer that female antelopes are more of a Type I species than Type II.

Code Listing

For reference, the base data evaluations and ggplot2 functions for producing this R markdown document are as follows:

library(tidyverse)
library(data.table)
library(dplyr)
library(latex2exp)
library(knitr)
library(kableExtra)
library(grid)
library(gridExtra)

survivalRates = function(data){
  with(data,{
    # Net reproductive rate per generation - mean number of female offspring produced per female lifetime
    R0 = sum(lxmx)/head(lx[1],1)
    # Average age a female gives birth
    G = sum(lx*mx*x)/R0
    # Rate of population change - logistic curve rate parameter
    lambda = R0^(1/G)
    # Change in population size per individual per unit of time
    r = log(lambda)
    # Stable age distribution - the proportion of the population in each age class over time
    cx = data.table(x=data$x, cx=lambda^(-x)*lx/(sum(lambda^(-x)*lx)))
  list(
    R0 = R0, G = round(G,2), lambda = round(lambda,3),
    r = round(r,3), cx = round(cx,3)) }) }

# Female Antelopes
Antelope <- data.table(
  # Categorised Age
  x = seq(0,12, length.out=13),
  # Number alive at start of age x
  nx = c(498,334,148,95,88,86,84,79,77,73,53,11,0),
  # Fecundity rate - reproductive output of cohort at each age x
  mx = c(0.00,0.10,0.15,0.30,0.44,0.46,0.42,0.43,
         0.39,0.35,0.25,0.14,0.00) )  %>%
  dplyr::mutate(
    # Standardised proportion of survivors at the start of age x
    lx = round(nx/nx[1],3),
    # Number dying during the age interval (x,x+1]
    dx = nx - append(tail(nx,-1),0),
    # Mortality rate - the intensity of deaths at each age x
    qx = round(dx/nx,3),
    # Mean number of survivors at age x
    Lx = as.integer(ceiling(nx - append(tail(nx,-1),0))/2 +
                      append(tail(nx,-1),0)),
    # Average remaining life expectancy at age x
    ex = round(rev(cumsum(rev(Lx)))/nx,3),
    # Average life expectancy of individuals of age x
    avgex = round(rev(cumsum(rev(Lx)))/nx + x,3),
    # Killing power = standardised intensity of mortality (log10)
    kx = round(log10(nx)-log10(append(tail(nx,-1),0)),3),
    # Mean number of female offspring produced by females at age x
    lxmx = round(lx*mx,3)
  )

sr <- survivalRates(Antelope)

# Charts
# 1. Proportional Mortality qx, Survival lx and killing power kx for the Fecundity of mx
basePlot <- function(){
  g <- ggplot(data=Antelope) +
    ylab(TeX(paste("Proportion of Cohort $(n_{0}=", max(Antelope$nx), ")$"))) +
    xlab("Age in Years") +
    labs(colour="Metric",
         title = "Female Antelope Life Table Metrics",
         subtitle = TeX(paste("\\textit{Proprtional evolution of mortality, $q_x$, killing power, ",
                              "(a $log_{10}$ value) ",
                              "$k_x$ and (non-standardised) survival, $l_x$ for a fecundity of $m_x$}"))) +
    geom_point(aes(y=mx,x=x,colour="Fecundity"))  +
    geom_point(aes(y=lx,x=x,colour="Survivorship"))  +
    geom_point(aes(y=qx,x=x,colour="Mortality"))  +
    geom_point(aes(y=kx,x=x,colour="Killing Power"))  +
    geom_smooth(aes(y=mx,x=x,colour="Fecundity"), stat="smooth",method="loess",se=F) +
    geom_smooth(aes(y=lx,x=x,colour="Survivorship"), stat="smooth",method="loess",se=F) +
    geom_smooth(aes(y=qx,x=x,colour="Mortality"), stat="smooth",method="loess",se=F) +
    geom_smooth(aes(y=kx,x=x,colour="Killing Power"), stat="smooth",method="loess",se=F) +
    geom_vline(aes(xintercept=sr$G),col="pink",lwd=1.1) +
    geom_text(aes(sr$G,0.6,label="G"),fontface="italic") +
    theme_bw() +
    theme(plot.title =
            element_text(hjust = 0.5,size=12,colour="black",face="bold"),
          plot.subtitle =
            element_text(hjust = 0,size=10,colour="black",face="italic"),
          axis.title =
            element_text(hjust = 0.5,size=10,colour="black"),
          axis.text =
            element_text(hjust = 0.5,size=10,colour="black"),
          legend.position = c(0.6,.83),
          legend.background=element_blank())  +
    scale_colour_hue( l=40) +
    scale_x_continuous(labels = function(x){sprintf("%.0f",x)}) +
    scale_y_continuous(labels = function(y){sprintf("%.2f",y)}) }


# 2. Life Expectancy
lePlot <- function(){
  g <- ggplot(data=Antelope) +
    ylab("Life Expectancy in Years") + xlab("Age in Years") +
    labs(colour="Life Expectancy Metric",
         title = "Female Antelope Life Expectancy Metrics",
         subtitle = TeX(paste("\\textit{Evolution of remaing life expectancy, $e_x$ and ",
                              " average life expectancy, $\\bar{e}_x$ for individuals at age $x$}"))) +
    geom_point(aes(y=ex,x=x,colour="Remaining"))  +
    geom_point(aes(y=avgex,x=x,colour="Average"))  +
    geom_smooth(aes(y=ex,x=x,colour="Remaining"), stat="smooth",method="loess",se=F) +
    geom_smooth(aes(y=avgex,x=x,colour="Average"), stat="smooth",method="loess",se=F) +
    geom_vline(aes(xintercept=sr$G),col="pink",lwd=1.1) +
    geom_text(aes(sr$G,8,label="G"),fontface="italic") +
    theme_bw() +
    theme(plot.title =
            element_text(hjust = 0.5,size=12,colour="black",face="bold"),
          plot.subtitle =
            element_text(hjust = 0,size=10,colour="black",face="italic"),
          axis.title =
            element_text(hjust = 0.5,size=10,colour="black"),
          axis.text =
            element_text(hjust = 0.5,size=10,colour="black"),
          legend.position = c(0.8,.5),
          legend.background=element_blank())  +
    scale_colour_hue( l=40) +
    scale_x_continuous(labels = function(x){sprintf("%.0f",x)}) +
    scale_y_continuous(labels = function(y){sprintf("%.2f",y)}) }

ltPlot <- function(data,vars,labs, tit,st,ypos){
  g <- ggplot(data=data,vars) +
    ylab(labs[1]) + xlab(labs[2]) + labs(title=tit,subtitle=st) +
    geom_point()  +
    geom_smooth(stat="smooth",method="loess",se=F) +
    geom_vline(aes(xintercept=sr$G),col="pink",lwd=1.1) +
    geom_text(aes(sr$G,ypos,label="G"),fontface="italic") +
    theme_bw() +
    theme(plot.title =
            element_text(hjust = 0.5,size=12,colour="black",face="bold"),
          plot.subtitle =
            element_text(hjust = 0,size=10,colour="black",face="italic"),
          axis.title =
            element_text(hjust = 0.5,size=10,colour="black"),
          axis.text =
            element_text(hjust = 0.5,size=10,colour="black"))  +
    scale_colour_hue( l=40) +
    scale_x_continuous(labels = function(x){sprintf("%.0f",x)}) +
    scale_y_continuous(labels = function(y){sprintf("%.2f",y)}) }

# 3. Standardised survival curve
sCurves <- function(data){
  xa <- data$x
  actual = log10(data$Lx+1)
  Actuals <- data.table(x=xa,actual=actual)
  xt <- seq(0,max(xa),length.out=1000)
  tI <- log10(rev(seq(0,data$Lx[1],length.out=1000))+1)
  m <- -actual[1]/max(xa); a <- actual[1]
  tII <- m*xt+a
  alp <- pi
  tIII <- tI * cos(alp) + xt * sin(alp) + max(tI)
  xtIII <- xt * cos(alp) - tI * sin(alp) + max(xt)
  Bounds <- data.table(x=xt,tI=tI,tII=tII,xtIII = xtIII,tIII=tIII)
  ggplot() +
    ylab(TeX(paste("$log_{10}$ of Mean Number of Survivors"))) + xlab("Age in Years") +
    labs(colour="Survival Curve Type",
         title = "Female Antelope Survivorship",
         subtitle = TeX(paste("\\textit{Evolution of standardised, $log_{10}$ survivorship ",
                              "for individuals at age $x$}"))) +
    geom_line(data=Bounds,aes(y=tI,x=x,colour="Type I")) +
    geom_line(data=Bounds,aes(y=tII,x=x,colour="Type II")) +
    geom_line(data=Bounds,aes(y=tIII,x=xtIII,colour="Type III")) +
    geom_point(data=Actuals,aes(y=actual,x=x)) +
    geom_smooth(data=Actuals,aes(y=actual,x=x),col="purple",method="loess",se=F) +
    geom_vline(aes(xintercept=sr$G),col="pink",lwd=1.1) +
    geom_text(aes(sr$G,0.6,label="G"),fontface="italic") +
    theme_bw() +
    theme(plot.title =
            element_text(hjust = 0.5,size=12,colour="black",face="bold"),
          plot.subtitle =
            element_text(hjust = 0,size=10,colour="black",face="italic"),
          axis.title =
            element_text(hjust = 0.5,size=10,colour="black"),
          axis.text =
            element_text(hjust = 0.5,size=10,colour="black"),
          legend.position = c(0.15,0.15),
          legend.background=element_blank())  +
    scale_colour_hue( l=40) +
    scale_x_continuous(labels = function(x){sprintf("%.0f",x)}) +
    scale_y_continuous(labels = function(y){sprintf("%.2f",y)}) }