This is a working document describing how we might build on the age structured model that was proposed by Schwacke et al. (2017):
Schwacke et al. 2017 Quantifying injury to common bottlenose dolphins from the Deepwater Horizon oil spill using an age-, sex- and class-structured population model. Endangered Species Research, 33: 265-279.
for Bay Sound and Estuary bottlenose dolphins (BSE BND) stocks, and extend it for species (and stocks) that were originally modeled by a stage structured model. Those correspond to all the other species and stocks in the GOM besides BSE BND. The idea was to try to use the Siler model results for BSE BND and “expand” them to other species.
While the document was originally intended as a where-we-are-today document on the 30th March 2020, aiming to set the stage to decide where to go next, it since then evolved into its own working document. It now fully describes the thought process and different considerations about scaling the BSE BND survival information used in the age structured model of Schwacke et al. (2017) to the different stocks.
The bulk of this document starts with prior references and thoughts about how one might scale survival across different species that might have different longevities.
This is followed by a quick exploration of the Siler model related material sent to TAM by LS via email on the 18th March 2020, including the model considered, the R functions used and the posterior distribution of the model paramenters considered in Schwacke et al. (2017).
This is followed by some comments about how to extend a function in general, building up to what we will have to do next.
We then build on it to investigate the Siler model and extending dolphins to other species and how it might be used to estimate survival probabilities per age by extending the BND model to the stocks originally addressed by the stage structured model. As a reality check we try to do this for the population of Elephants reported in Lahdenpera et al. (2018)
Lahdenpera et al 2018. Differences in age-specific mortality between wild-caught and captive-born Asian elephants. Nature Communications 9: 3023 DOI: http://doi.org/10.1038/s41467-018-05515-8.
This is followed by a side step on the mean age of an animal sampled at random, a quantity that would be involved in the recovery of the survival and fecundity after exposure to oil. This is addressed in a separate document (How to land it).
We conclude with some pending questions that need solving before all this is put to rest. This is essentially a section with hiperlinks to points in the text that require further attention.
The idea of using information about survival from one species to simulate/estimate/model/draw insights about survival in other species is not new. In the following we list some references that have looked into the problem before.
The idea is that we have age indexed survivals from a Siler model, but we want to extend or shrink these to match a species with a diferent longevity. Thinking about a “simple” scaling, we could expand or contract a function such that the resulting survivals match a given target statistic. Different suggestions have been put forward for that target statistic, but each of these have advantages and disadvantages. These target statistics include:
Screenshot 1 of Krementz et al 1989.
Screenshot 1 of Eakin et al. 1994.
Screenshot 1 of Caswell et al. 1998.
A few other references on scaling survival functions and patterns in survival across different species, not so relevant and borderline useless, but here for future reference if needed:
Trends in scale and shape of survival curves https://www.nature.com/articles/srep00504
Survivorship Curves and Their Impact on the Estimation of Maximum Population Growth Rates https://www.jstor.org/stable/25592598?seq=1
Ageing and typical survivorship curves result from optimal resource allocation http://www.cyfronet.krakow.pl/~uxcichon/pub/eer.pdf
Unfortunately, at least for TAM, only late in the game he realized that Barlow & Boveng (1991)
Barlow & Boveng 1991 Modeling age-specific mortality for marine mammal populations. Marine Mammal Science 7: 50-65
was a very relevant source to this problem. Wording might be at times slightly inconsistent and implying we were not aware of this work originally. That paper is really good. Being aware of it would have helped solving many questions listed below, and certainly faster, inclusively as an example the reparametrization issue!
Here we provide some relevant screenshots from that paper, to keep in the back of our minds and inform the discussion.
Screenshot 1 of BBarlow and Boveng 1991.
So essentially, this is what we did to get the Siler model to begin with (Ask Len - right? Were you just following Jay’s footsteps?)
Here’s the basic notation (see how this would have solved the reparametrization that I had to go at by trial and error?)
Screenshot between 1 and 2 of BBarlow and Boveng 1991.
and the notion that functions can be re-scaled:
Screenshot 2 of BBarlow and Boveng 1991.
A detail about how females are the focus here
Screenshot 3 of BBarlow and Boveng 1991.
The species chosen as a model: the northern fur seal.
Screenshot 4 of BBarlow and Boveng 1991.
The lack of information for cetaceans in general. At the date, but still overall true.
Screenshot 5 of BBarlow and Boveng 1991.
The necessary caveats. Have I said I find Jay an amazing scientist before? Not screenshoted here, but the way he criticizes his own work should be an example to us all.
Screenshot 6 of BBarlow and Boveng 1991.
Insights about issues of how hard it might be about finding a representative general model for all cetaceans and the idea of scaling such a model.
Screenshot 7 of BBarlow and Boveng 1991.
A key comment about the why they think mean age is not a good idea. Something we need to get to grasps with
I actually have a feeling that scaling with respect to a mean value makes more sense than do it with respect to an extreme quantile: extremes are harder to estimate. We can estimate them, but they have necessarily much more variance associated with that say mean, or a non-extreme quartile, as they are much more dependent on any given sample not providing an outlier, say. I wonder what Jay would say about this today?
Screenshot 8 of BBarlow and Boveng 1991.
Making some suggestions for when the life table data is not available, i.e., for people like us…
Screenshot 9 of BBarlow and Boveng 1991.
option 1 - combine info from a life table and a model
Screenshot 10 of BBarlow and Boveng 1991.
option 2 - fix some of the parameters of the Siler model and estimate others
Screenshot 11 of BBarlow and Boveng 1991.
option 3 - consider empirical Bayes with priors based on data
Screenshot 12 of BBarlow and Boveng 1991.
Yes, here we are… still :)
Screenshot 13 of BBarlow and Boveng 1991.
This material corresponds to two emails:
Hence, we have the values (in fact, we have samples from the posterior distribution for each of the 5 model parameters, and the growth rate - for the estimated parameters for the Siler model) and the code that allows to evaluate all the relevant functions given parameter values.
Here we explore the different datasets that Lori sent TAM
(note: these will not be used here, so just checking)
Reading the female life table
FLT <- read.table("Female Life Table.txt", quote="\"", comment.char="")
MLT <- read.table("Male Life Table.txt", quote="\"", comment.char="")
This is weird. The data corresponds to a single column with 976061 lines, and there are 16001 (females) and 16001 (males) values of 0 in these. A detective look tells me the 0’s appear every 61 (female) and 61 (males) values in both files, and that can’t be a coincidence… these must correspond to the age classes, and the age 60 probability of survival of 0. After discussion with LS we can ignore everything after the first 61 values. But then again, these will not be used here.
Check this looks about right by plotting the first 61 values for each file
par(mfrow=c(1,1))
plot(0:60,MLT$V1[1:61],type="l",lty=1,ylab="Survival probability",xlab="Age",ylim=c(0,1))
lines(0:60,FLT$V1[1:61],lty=2)
legend("bottomleft",legend=c("Males","Females"),lty=1:2, inset=0.05)
Now look at the parameter values. Read the data
pf <- read.csv("var_female.csv")
pm <- read.csv("var_male.csv")
#remove useless ID first column
pf<-pf[,-1]
pm<-pm[,-1]
names(pf) <- names(pm) <- c("a1","a2","a3","b1","b3","rmean")
and check there are consistent.
First for females
Then for males
Just to get a feeling for the relevant functions, use the code provided by LS to plot these. We consider just the first value from the posterior distribution as an illustrative example.
The required functions used by LS were sent to TAM by LS on the 18th March 2020
#funtions sent via email by LS to TAM on the 18th March 2020
#source("usefulfunctions.R")
# Proportion surviving to age x
lx <- function(x,a1,a2,a3,b1,b3) {
exp(-a2*x) * exp((-a1)*(1-exp(-b1*x))) * exp(a3*(1-exp(b3*x)))
}
# Age specific mortality rate
px <- function(x,a1,a2,a3,b1,b3) {
lx(x+1,a1,a2,a3,b1,b3)/lx(x,a1,a2,a3,b1,b3)
}
qx <- function(x,a1,a2,a3,b1,b3) {
1-px(x,a1,a2,a3,b1,b3);
}
# Number dying in age class x
dx <- function(x,a1,a2,a3,b1,b3,r) {
exp(-r*x)*lx(x,a1,a2,a3,b1,b3)*qx(x,a1,a2,a3,b1,b3)
}
nx <- function(x,a1,a2,a3,b1,b3,r) {
exp(-r*x)*lx(x,a1,a2,a3,b1,b3)
}
ages <- 0:60
plot(ages,lx(ages,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5]),type="l",lty=1,main="Proportion surviving to age x",ylab="Proportion",xlab="Age")
lines(ages,lx(ages,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5]),lty=2)
legend("topright",legend=c("Males","Females"),lty=2:1, inset=0.05)
ages <- 0:60
plot(ages,px(ages,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5]),type="l",lty=1,main="Age specific survival rate",ylab="Survival",xlab="Age",ylim=c(0,1))
lines(ages,px(ages,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5]),lty=2)
legend("topright",legend=c("Males","Females"),lty=2:1, inset=0.05)
(see this is just another realization of the first plot above!)
ages <- 0:60
plot(ages,qx(ages,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5]),type="l",lty=1,main="Age specific mortality rate",ylab="Mortality",xlab="Age",ylim=c(0,1))
lines(ages,qx(ages,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5]),lty=2)
legend("topright",legend=c("Males","Females"),lty=2:1, inset=0.05)
ages <- 0:60
plot(ages,dx(ages,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5],pf[1,6]),type="l",lty=1,main="Number dying in age class",ylab="Number dying",xlab="Age",ylim=c(0,1))
lines(ages,dx(ages,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5],pm[1,6]),lty=2)
legend("topright",legend=c("Males","Females"),lty=2:1, inset=0.05)
ages <- 0:60
plot(ages,nx(ages,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5],pf[1,6]),type="l",lty=1,main="Not sure yet",ylab="Not sure yet",xlab="Age",ylim=c(0,1))
lines(ages,nx(ages,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5],pm[1,6]),lty=2)
legend("topright",legend=c("Males","Females"),lty=1:2, inset=0.05)
We refer loosely to the notion of extending a function to the otherwise loose concept of “going from one function to the other” by changing some feature(s) on the original function to reach a target statistic (or set of statistics) via the “extended” function. This stretching might correspond to stretching, shrinking, increasing, decreasing, time warping, etc, and sometimes these different expressions might be used depending on context.
Specifically, we are here referring to “extending” the Siler model estimated on BND BSE to be a reasonable approximation to a description of survival probabilities for other species.
As LT correctly pointed out, extending the function lx
is straightforward, but the function px
is a bit harder, because that is a quantity derived from the lx
applied in two subsequent years.
In other words, px
is the survival in year \(x\), and that is obtained by the quotient of the proportion of animals that reach age \(x+1\) and those that reach age \(x\).
\[p(x)=\frac{l(x+1)}{l(x)}\]
So when one scales the function over time, one needs also to correspondingly stretch what the comparison between a year and the next is.
The best way to cope for this is perhaps to re-code the functions px
, such that it has explicitly a scaling, that under no stretching is 1 (and is 1 by default).
# Age specific mortality rate
px <- function(x,a1,a2,a3,b1,b3,scaling=1) {
lx(x+scaling,a1,a2,a3,b1,b3)/lx(x,a1,a2,a3,b1,b3)
}
We want to extend this function, that was estimated for a given species, to another species. Assume (these are quick and dirty approximations, but ballpark correct) that the mean age of a dolphin at death is 14 years. For a sperm whale that value increases to about 25 years (24.85 years: see separate R code document “FromStage2AgeModel.R” on how we obtained this number. This is hosted in the folder with the code for the sperm whale analysis included in the GOMOSES 2020 poster.). Then we could conceivable visualize the effect of extending the functions by a factor of 25/14. Lets see what that corresponds to for both the probability of being alive at a given age (lx
) and the survival at year x (px
).
First, we do it by changing the time explicitly:
par(mfrow=c(1,2))
#the scaling factor, as explained above (this is an aproximation based on a gross integers rounding, just for the purposes of this working document)
#original age value/target age value
sca1 <- 14/25
#proportion alive at age
plot(ages/sca1,lx(ages,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5]),type="l",lty=1,main="Proportion surviving to age x",ylab="Proportion",xlab="Age",xlim=c(0,120))
lines(ages/sca1,lx(ages,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5]),lty=2)
legend("topright",legend=c("Males","Females"),lty=1:2, inset=0.02)
#survival px - note the use of the scaling argument in "px"
plot(ages/sca1,px(ages,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5],scaling=sca1),type="l",lty=1,main="Survival probability at age x",ylab="Survival",xlab="Age",xlim=c(0,120),ylim=c(0,1))
lines(ages/sca1,px(ages,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5],scaling=sca1),lty=2)
legend("bottomleft",legend=c("Males","Females"),lty=1:2, inset=0.02)
Or equivalently, not changing the data on the plot but the x argument in the functions themselves, which means that for expanding by a factor c
one needs to re-scale the function f(x)
to f(x*K)
. Note that If c
is greater than one the function will undergo horizontal shrinking, and if c
is less than one the function will undergo horizontal stretching. (see also https://courses.lumenlearning.com/boundless-algebra/chapter/transformations/ for extra details and thoughts on function stretching etc)
#define the new maximum age for the extended function
nma <- max(ages)/sca1
#the range of ages
ager <- 0:nma
par(mfrow=c(1,2))
# Proportion surviving to age x
plot(ager,lx(ager*sca1,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5]),type="l",lty=1,main="Proportion surviving to age x",ylab="Proportion",xlab="Age",xlim=c(0,120))
lines(ager,lx(ager*sca1,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5]),lty=2)
legend("topright",legend=c("Males","Females"),lty=1:2, inset=0.02)
# Survival at age x
plot(ager,px(ager*sca1,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5],scaling=sca1),type="l",lty=1,main="Survival at age x",ylab="Survival",xlab="Age",ylim=c(0,1),xlim=c(0,120))
lines(ager,px(ager*sca1,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5],scaling=sca1),lty=2)
legend("bottomleft",legend=c("Males","Females"),lty=1:2, inset=0.02)
What we want is to find what is the correct value of c
that leads from a mean age at death of 14 years to 24 years.
Note that one can easily set up a large number of age classes and then use the same approach to evaluate what the corresponding functions would look like for different degrees of stretching. In this way one could easily “get ready” for dealing with different species with different mean ages at death. The grey lines in the plot below represent intermediate degrees of stretching, from no stretching (\(c=1\)) to a \(c\) of 0.5 (solid black lines)
stretching <- seq(0.51,0.99,by=0.01)
#survival
#females
par(mfrow=c(1,2))
plot(0:150,lx((0:150)*0.5,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5]),type="l",lty=1,main="Females surviving to age x",ylab="Proportion",xlab="Age")
lines(0:150,lx((0:150),pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5]))
for(i in 1:length(stretching)){
lines(0:150,lx((0:150)*stretching[i],pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5]),col="grey")
}
#males
plot(0:150,lx((0:150)*0.5,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5]),type="l",lty=1,main="Males surviving to age x",ylab="Proportion",xlab="Age",ylim=c(0,1))
lines(0:150,lx((0:150),pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5]))
for(i in 1:length(stretching)){
lines(0:150,lx((0:150)*stretching[i],pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5]),col="grey")
}
The corresponding survival probabilities are also shown below.
stretching <- seq(0.51,0.99,by=0.01)
#survival
#females
par(mfrow=c(1,2))
plot(0:150,px((0:150)*0.5,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5],scaling=0.5),type="l",lty=1,main="Female survival at age x",ylab="Proportion",xlab="Age",ylim=c(0,1))
lines(0:150,px((0:150),pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5]))
for(i in 1:length(stretching)){
lines(0:150,px((0:150)*stretching[i],pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5],scaling=stretching[i]),col="grey")
}
#males
plot(0:150,px((0:150)*0.5,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5],scaling=0.5),type="l",lty=1,main="Male survival at age x",ylab="Proportion",xlab="Age",ylim=c(0,1))
lines(0:150,px((0:150),pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5]))
for(i in 1:length(stretching)){
lines(0:150,px((0:150)*stretching[i],pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5],scaling=stretching[i]),col="grey")
}
When extending the function we refer to the scaling factor as the quotient between the original age and the target age. Therefore, considering the x axis, when using a scaling factor lower than one we stretch the function, and when we use a factor higher than one we shrink the function.
Next we evaluate by simulation (see below this is easier and faster to do by integrating lx
) the mean age at death as a function of this scaling factor. We are trying to stretch the function, hence we consider scaling factors lower than 1.
ls <- length(stretching)
#maximum age
maxage <- 150
#vector to hold mean age after stretching
mean.ages <- numeric(ls)
#range of ages
ages<-0:maxage
for(s in 1:ls){
#calculate stretched survivals
# I have this - I had here
# px(ages[-maxage]*stretching[s]...
# but removed it - not sure why - but things make no sence otherwise :()
survsF <- px(ages[-maxage]*stretching[s],pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5],scaling=stretching[s])
survsM <- px(ages[-maxage]*stretching[s],pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5],scaling=stretching[s])
#because the px function is a division, if both numerator and denominator are 0 we get a NaN - replace by 0's as that is the real survival probability in that case
survsF[is.nan(survsF)]<-0
survsM[is.nan(survsM)]<-0
#build Survival matrix
Tr <- matrix(0,nrow=maxage*2,ncol=maxage*2)
for(i in 1:maxage){
for (i in 1:(maxage-1)){
Tr[i+1,i] <- survsF[i]
Tr[maxage+i+1,maxage+i] <- survsM[i]
}
}
#define inital population
nstart <- 500000
Nstart=c(nstart/2,rep(0,maxage-1),nstart/2,rep(0,maxage-1))
#propagate it forward
years <- 200
Ns<-matrix(0,nrow=maxage*2,ncol=years)
Ns[,1] <- Nstart
for(i in 2:years){
Ns[,i] <- Tr %*% Ns[,i-1]
}
#see the pop size as a function of time
#plot(1:years,colSums(Ns)/nstart,ylim=c(0,1),xlab="year",ylab="P(survival up to age)")
#observed ages at death
ages.at.death <- rep(1:(years-1),times=abs(diff(colSums(Ns))))
#mean age at death
mean.ages[s] <- mean(ages.at.death)
}
#see the result
plot(stretching,mean.ages)
The effect is as expected.
In the next plot, just for comparison, we not only stretch the function horizontally, as we multiply the survival probabilities by a small number (1.005). That naturally leads to a much faster increase in the mean ages. (this was done to illustrate what might need to happen in practice, since just stretching the function might not be enough).
ls <- length(stretching)
#maximum age
maxage <- 150
#vector to hold mean age after stretching
mean.ages <- numeric(ls)
#range of ages
ages<-0:maxage
for(s in 1:ls){
#calculate stretched survivals
survsF <- px(ages[-maxage]*stretching[s],pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5],scaling=stretching[s])
survsM <- px(ages[-maxage]*stretching[s],pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5],scaling=stretching[s])
#because the px function is a division, if both numerator and denominator are 0 we get a NaN
survsF[is.nan(survsF)]<-0
survsM[is.nan(survsM)]<-0
survsF=1.005*survsF
survsM=1.005*survsM
#build Survival matrix
Tr <- matrix(0,nrow=maxage*2,ncol=maxage*2)
for(i in 1:maxage){
for (i in 1:(maxage-1)){
Tr[i+1,i] <- survsF[i]
Tr[maxage+i+1,maxage+i] <- survsM[i]
}
}
#define inital population
nstart <- 500000
Nstart=c(nstart/2,rep(0,maxage-1),nstart/2,rep(0,maxage-1))
#propagate it forward
years <- 200
Ns<-matrix(0,nrow=maxage*2,ncol=years)
Ns[,1] <- Nstart
for(i in 2:years){
Ns[,i] <- Tr %*% Ns[,i-1]
}
#see the pop size as a function of time
#plot(1:years,colSums(Ns)/nstart,ylim=c(0,1),xlab="year",ylab="P(survival up to age)")
#observed ages at death
ages.at.death <- rep(1:(years-1),times=abs(diff(colSums(Ns))))
#mean age at death
mean.ages[s] <- mean(ages.at.death)
}
#see the result
plot(stretching,mean.ages)
So the question now becomes, how to do it in practice.
Note the value of the mean age at death under no stretching is 15.3, a bit higher than the value originally found for the BND of around 14.
This might just be a consequence of using different survival probabilities, before used a random sample from the survivals resulting from the Siler model, here another random sample derived from a single iteration of the Siler model.
If we look into the distribution of possible mean ages, we see that the value 14 is well within the range expected. This reflects something that we will need to account for: we need to find out a good way to propagate the uncertainty in the Siler model, via this stretching procedure, to the resulting age structure model.
#maximum age
maxage <- 150
# number of repeats
repeats <- 1000
#vector to hold mean ages
mean.ages2 <- numeric(repeats)
#range of ages
ages<-0:maxage
for(s in 1:repeats){
#calculate stretched survivals
survsF <- px(ages[-maxage],pf[s,1],pf[s,2],pf[s,3],pf[s,4],pf[s,5])
survsM <- px(ages[-maxage],pm[s,1],pm[s,2],pm[s,3],pm[s,4],pm[s,5])
#because the px function is a division, if both numerator and denominator are 0 we get a NaN
survsF[is.nan(survsF)]<-0
survsM[is.nan(survsM)]<-0
#build Survival matrix
Tr <- matrix(0,nrow=maxage*2,ncol=maxage*2)
for(i in 1:maxage){
for (i in 1:(maxage-1)){
Tr[i+1,i] <- survsF[i]
Tr[maxage+i+1,maxage+i] <- survsM[i]
}
}
#define inital population
nstart <- 500000
Nstart=c(nstart/2,rep(0,maxage-1),nstart/2,rep(0,maxage-1))
#propagate it forward
years <- 200
Ns<-matrix(0,nrow=maxage*2,ncol=years)
Ns[,1] <- Nstart
for(i in 2:years){
Ns[,i] <- Tr %*% Ns[,i-1]
}
#see the pop size as a function of time
#plot(1:years,colSums(Ns)/nstart,ylim=c(0,1),xlab="year",ylab="P(survival up to age)")
#observed ages at death
ages.at.death <- rep(1:(years-1),times=abs(diff(colSums(Ns))))
#mean age at death
mean.ages[s] <- mean(ages.at.death)
}
#see the result
hist(mean.ages,xlab="Age",freq=FALSE,main="")
Do we, and how do we, account for variability in mean age in this exercise?
The paper
Lahdenpera et al 2018. Differences in age-specific mortality between wild-caught and captive-born Asian elephants. Nature Communications 9: 3023 DOI: http://doi.org/10.1038/s41467-018-05515-8
presents survival estimates from a sample of elephants, also considering a Siler model. The minor difference is that these authors included covariates on the Siler model parameters.
A way to see if the above procedure is somewhat reasonable (hard to know upfront what this means, as we are comparing apples and oranges- well literally, dolphins and elephants) is to check if we manage to extend the dolphins survivals to match the elephants key age statistics. The resulting survivorship functions for the elephants are shown in Figure . We can see that these long lived animals will outlive the dolphins.
Screenshot of figure 4 in Lahdenpera et al 2018. Survivorship functions resulting from the fitted Siler models.
However, at first sight, what seems obvious to me is that while the survivorship in the dolphins is concave, in the elephants it is convex. And this is something that our “expanding” method can’t cope with, because a concave function will not turn convex by “expansion”, or vice versa.
The baseline mortality rates are also presented here:
Screenshot of Figure 1 in Lahdenpera et al 2018.
Any way, lets proceed and see where this gets us.
To facilitate parameter matching, the parametrization used in the elephant model is shown here (note we are not considering the second line - that is only involved with additional covariates, not present for the dolphins - hence only the first line applies here)
Screenshot of Siler model parametrization in Lahdenpera et al 2018.
and that used in Lori’s paper is here:
Screenshot of Siler model parametrization in Schwacke et al 2017.
It is key to note that while \(p\) from the elephant paper (iEp) represents the hazard (probability of dying at a given age), in Lori’s paper (iLp) \(l_S(x)\) represents the probability of surviving up to age x, and so in fact
\[ 1-\frac{l(x+1)}{l(x)}=1-p(x)=q(x)= in~Lori's~paper~is~equivalent~to~the~Elephant's ~paper~=p \]
A source of potential confusion is worth bringing upfront: there are two functions denoted as \(p\), and the Elephant’s paper \(p\) corresponds to Lori code \(q(x)\), not \(p(x)\) which actually corresponds to \(1-q(x)\).
Note therefore that, irrespective of there being reparametrizations, we have the following connections:
First, we need to read in the model parameter estimates from the Elephants. From the supplementary material in the elephant paper we have estimates for the Siler model parameters for both males and females (there are actually a large number of parameters, model averaging intricacies, etc, but I ignore these for now). This is shown in , and I note explicitly we are considering the \(w_3\) for the year 2000 cohort.
Screenshot of Supplementary Table 1 in Lahdenpera et al 2018. Siler model parameter estimates.
Now we use these parameters and check that we can get results that look like those presented in the Elephant paper.
#first row females
#second row males
eleSp<-data.frame(w1=c(0.0807,0.10465),w2=c(0.00101,0.0027),w3=c(0.00426,0.00426),b1=c(0.34585,0.30425),b2=c(0.07656,0.06513),row.names = c("females","males"))
kable(eleSp)
w1 | w2 | w3 | b1 | b2 | |
---|---|---|---|---|---|
females | 0.08070 | 0.00101 | 0.00426 | 0.34585 | 0.07656 |
males | 0.10465 | 0.00270 | 0.00426 | 0.30425 | 0.06513 |
We can now plot the survivorship function:
ages <- 0:60
plot(ages,lx(ages,eleSp[1,1],eleSp[1,3],eleSp[1,2],eleSp[1,4],eleSp[1,5]),type="l",lty=1,main="Proportion surviving to age x",ylab="Proportion",xlab="Age",ylim=c(0,1))
lines(ages,lx(ages,eleSp[2,1],eleSp[2,3],eleSp[2,2],eleSp[2,4],eleSp[2,5]),lty=2)
legend("bottomleft",legend=c("Males","Females"),lty=2:1, inset=0.05)
Something went terribly wrong, as these are much larger survivorship values than what is presented in the image above.
Might I have misunderstood what parameter is what across parametrizations?
We can plot equation 1 in the elephant paper
p <- function(x,w1,w2,w3,b1,b2){
w1*exp(-b1*x)+w2*exp(b2*x)+w3}
ages <- 0:60
plot(ages,p(ages,eleSp[1,1],eleSp[1,2],eleSp[1,3],eleSp[1,4],eleSp[1,5]),type="l",lty=1,main="Equation 1 Elephant",ylab="Proportion",xlab="Age",ylim=c(0,0.15))
lines(ages,p(ages,eleSp[2,1],eleSp[2,2],eleSp[2,3],eleSp[2,4],eleSp[2,5]),lty=2)
legend("bottomright",legend=c("Males","Females"),lty=2:1, inset=0.02)
and that seems to add up well to what they have in the paper (again, note different baseline so not really the same necessarily).
Now we plot the age specific mortality rate, using Lori’s plotting code, with the Elephant parameter values.
ages <- 0:60
plot(ages,qx(ages,eleSp[1,1],eleSp[1,3],eleSp[1,2],eleSp[1,4],eleSp[1,5]),type="l",lty=1,main="Mortality rate at age x",ylab="Proportion",xlab="Age",ylim=c(0,0.15))
lines(ages,qx(ages,eleSp[2,1],eleSp[2,3],eleSp[2,2],eleSp[2,4],eleSp[2,5]),lty=2)
legend("topright",legend=c("Males","Females"),lty=2:1, inset=0.02)
This is much lower and does not add up - so this is likely an issue with the parametrization.
Back to the blackboard, considering the expressions in Siler’s paper,
Screenshot of Figure 1 in Lahdenpera et al 2018.
a slight reparametrization might be at stake, where \(a_1\) is divided by \(b_1\), and \(a_3\) is divided by \(b_3\), and so
\[a1=\frac{w1}{b1}\] \[a3=\frac{w2}{b2}\]
We can now check the Mortality rate, using the Elephant parameter values and Lori’s code with the new parametrization:
ages <- 0:60
plot(ages,qx(ages,eleSp[1,1]/eleSp[1,4],eleSp[1,3],eleSp[1,2]/eleSp[1,5],eleSp[1,4],eleSp[1,5]),type="l",lty=1,main="Mortality rate at age x",ylab="Proportion",xlab="Age",ylim=c(0,0.15))
lines(ages,qx(ages,eleSp[2,1]/eleSp[2,4],eleSp[2,3],eleSp[2,2]/eleSp[2,5],eleSp[2,4],eleSp[2,5]),lty=2)
legend("topright",legend=c("Males","Females"),lty=2:1, inset=0.02)
And now the lx
function, i.e. the proportion surviving to a given age, and that seems now to add up.
ages <- 0:60
plot(ages,lx(ages,eleSp[1,1]/eleSp[1,4],eleSp[1,3],eleSp[1,2]/eleSp[1,5],eleSp[1,4],eleSp[1,5]),type="l",lty=1,main="Proportion surviving to age x",ylab="Proportion",xlab="Age",ylim=c(0,1))
lines(ages,lx(ages,eleSp[2,1]/eleSp[2,4],eleSp[2,3],eleSp[2,2]/eleSp[2,5],eleSp[2,4],eleSp[2,5]),lty=2)
legend("topright",legend=c("Males","Females"),lty=2:1, inset=0.02)
It seems like all is good, and we probably are where we want to be, but it would be nice to demonstrate that if one starts at
\[ 1-\frac{l(x+1)}{l(x)}=1-p(x)=q(x)= in~Lori's~paper~is~equivalent~to~the~Elephant's ~paper~=p \]
and then check the parameters from p correspond to the reparametrizations found.
In any case, I can try to empirically show this is correct:
par(mfrow=c(1,1))
ages <- 0:60
#use Elephant formulation
plot(ages,p(ages,eleSp[1,1],eleSp[1,2],eleSp[1,3],eleSp[1,4],eleSp[1,5]),type="l",lty=1,main="Equation 1 Elephant (p) vs. Lori's q(x) ",ylab="Proportion",xlab="Age",ylim=c(0,0.15),col=3)
lines(ages,p(ages,eleSp[2,1],eleSp[2,2],eleSp[2,3],eleSp[2,4],eleSp[2,5]),lty=2,col=3)
#Use Lori's formulation with Elephants parameters
lines(ages,qx(ages,eleSp[1,1]/eleSp[1,4],eleSp[1,3],eleSp[1,2]/eleSp[1,5],eleSp[1,4],eleSp[1,5]),type="l",lty=1,main="Hazard at age x",ylab="Proportion",xlab="Age",ylim=c(0,0.15),col=4)
lines(ages,qx(ages,eleSp[2,1]/eleSp[2,4],eleSp[2,3],eleSp[2,2]/eleSp[2,5],eleSp[2,4],eleSp[2,5]),lty=2,col=4)
legend("top",legend=c("Males p","Females p","Males q(x)","Females q(x)"),lty=c(2,1,2,1),col=c(3,3,4,4), inset=0.02)
The result in green and blue is almost the same, but it should be exactly the same. Is this an issue with some numerical approximation gone wrong?
For the time being I proceed assuming the above is correct but there is something that is not really adding up in this. I suspect the consequences are negligible, but would be great to find out the slight discrepancy.
We can evaluate what would be the average life span of an Elephant. Given the probability of an elephant reaching a given age x is provided by \(l(x)\), as in the image below
ages <- 0:100
plot(ages,lx(ages,eleSp[1,1]/eleSp[1,4],eleSp[1,3],eleSp[1,2]/eleSp[1,5],eleSp[1,4],eleSp[1,5]),type="l",lty=1,main="Proportion surviving to age x",ylab="Proportion",xlab="Age",ylim=c(0,1))
lines(ages,lx(ages,eleSp[2,1]/eleSp[2,4],eleSp[2,3],eleSp[2,2]/eleSp[2,5],eleSp[2,4],eleSp[2,5]),lty=2)
legend("topright",legend=c("Males","Females"),lty=2:1, inset=0.02)
then the mean age of an elephant will be given by the integral under that function. Assuming that we have 1:1 sex ratio, the mean age is the mean age of males and females.
MeanAgeEmale <- integrate(lx,lower=0,upper=100,a1=eleSp[1,1]/eleSp[1,4],a2=eleSp[1,3],a3=eleSp[1,2]/eleSp[1,5],b1=eleSp[1,4],b3= eleSp[1,5])$value
MeanAgeEfemale <- integrate(lx,lower=0,upper=100,a1=eleSp[2,1]/eleSp[2,4],a2=eleSp[2,3],a3=eleSp[2,2]/eleSp[2,5],b1=eleSp[2,4],b3= eleSp[2,5])$value
MeanAgeE <- (MeanAgeEmale+MeanAgeEfemale)/2
This corresponds to a mean age of Male Elephants of 35.784 and of Female Elephants of 28.005, resulting in an overall mean age under a 1:1 sex ratio assumption of 31.895.
Just to check that we get a similar mean age using the simulation approach considered above
#maximum age
maxage <- 150
#range of ages
ages<-0:maxage
#calculate stretched survivals
survsF <- px(ages,eleSp[2,1]/eleSp[2,4],eleSp[2,3],eleSp[2,2]/eleSp[2,5],eleSp[2,4],eleSp[2,5])
survsM <- px(ages,eleSp[1,1]/eleSp[1,4],eleSp[1,3],eleSp[1,2]/eleSp[1,5],eleSp[1,4],eleSp[1,5])
#because the px function is a division, if both numerator and denominator are 0 we get a NaN
survsF[is.nan(survsF)]<-0
survsM[is.nan(survsM)]<-0
#build Survival matrix
Tr <- matrix(0,nrow=maxage*2,ncol=maxage*2)
for(i in 1:maxage){
for (i in 1:(maxage-1)){
Tr[i+1,i] <- survsF[i]
Tr[maxage+i+1,maxage+i] <- survsM[i]
}}
#define inital population
nstart <- 500000
Nstart=c(nstart/2,rep(0,maxage-1),nstart/2,rep(0,maxage-1))
#propagate it forward
years <- 200
Ns<-matrix(0,nrow=maxage*2,ncol=years)
Ns[,1] <- Nstart
for(i in 2:years){
Ns[,i] <- Tr %*% Ns[,i-1]
}
#observed ages at death
ages.at.death <- rep(1:(years-1),times=abs(diff(colSums(Ns))))
#mean age at death
mean.ages<- mean(ages.at.death)
mean.ages
## [1] 32.40203
This is a difference of 0.51. Half a year… is this a coincidence, or am I doing something wrong in the approximation by simulation? I actually think that it is the latter, because as currently implemented by simulation the minimum age one records an animal dying at is 1, i.e. you record the year the animal would have in it’s next birthday after dying, which will be on average about half a year later than its actual age, assuming uniform deaths during a year.
So now, let’s stretch dolphins to become elephants… (I always wanted to say that… it’s a remarkable Lamark remark :) - whether this is a direct consequence of social isolation for too long is unclear. But it’s not a big stretch to say so!
We can essentially redo what we attempted above, but now knowing we are aiming for a mean age at death of 31.89. That essentially means more “intense” stretching than before.
stretching <- seq(0.11,0.99,by=0.01)
maxage <- 300
#survival
#females
plot(0:maxage,lx((0:maxage)*0.5,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5]),type="l",lty=1,main="Females surviving to age x",ylab="Proportion",xlab="Age")
lines(0:maxage,lx((0:maxage),pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5]))
for(i in 1:length(stretching)){
lines(0:maxage,lx((0:maxage)*stretching[i],pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5]),col="grey")
}
ls <- length(stretching)
#vector to hold mean age after stretching
mean.ages <- numeric(ls)
#range of ages
ages<-0:maxage
for(s in 1:ls){
#calculate stretched survivals
survsF <- px(ages[-maxage]*stretching[s],pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5],scaling=stretching[s])
survsM <- px(ages[-maxage]*stretching[s],pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5],scaling=stretching[s])
#because the px function is a division, if both numerator and denominator are 0 we get a NaN
survsF[is.nan(survsF)]<-0
survsM[is.nan(survsM)]<-0
#build Survival matrix
Tr <- matrix(0,nrow=maxage*2,ncol=maxage*2)
for(i in 1:maxage){
for (i in 1:(maxage-1)){
Tr[i+1,i] <- survsF[i]
Tr[maxage+i+1,maxage+i] <- survsM[i]
}
}
#define inital population
nstart <- 500000
Nstart=c(nstart/2,rep(0,maxage-1),nstart/2,rep(0,maxage-1))
#propagate it forward
years <- 400
Ns<-matrix(0,nrow=maxage*2,ncol=years)
Ns[,1] <- Nstart
for(i in 2:years){
Ns[,i] <- Tr %*% Ns[,i-1]
}
#see the pop size as a function of time
#plot(1:years,colSums(Ns)/nstart,ylim=c(0,1),xlab="year",ylab="P(survival up to age)")
#observed ages at death
ages.at.death <- rep(1:(years-1),times=abs(diff(colSums(Ns))))
#mean age at death
mean.ages[s] <- mean(ages.at.death)
}
#see the result
plot(stretching,mean.ages)
We proceed with the comparison using the the direct integration method. This is more computationally efficient any way and seems to avoid a couple of issues (the age being incorrectly computed by half a year - see above - and this current unexpected result most likely due to some bug in the code).
We define a new lx
function, called lxstretch
, that can be numerically integrated by just adjusting the x’s.
#the required function
# Proportion surviving to age x
lxstretch <- function(x,stretch,a1,a2,a3,b1,b3) {
exp(-a2*x*stretch) * exp((-a1)*(1-exp(-b1*x*stretch))) * exp(a3*(1-exp(b3*x*stretch)))
}
stretching <- seq(0.11,0.99,by=0.01)
ls <- length(stretching)
#vector to hold mean age after stretching
mean.agesI <- numeric(ls)
for(s in 1:ls){
MeanAgeDmale <- integrate(lxstretch,lower=0,upper=500,stretch=stretching[s],a1=pf[1,1],a2=pf[1,2],a3=pf[1,3],b1=pf[1,4],b3= pf[1,5])$value
MeanAgeDfemale <- integrate(lxstretch,lower=0,upper=500,stretch=stretching[s],a1=pm[1,1],a2=pm[1,2],a3=pm[1,3],b1=pm[1,4],b3= pm[1,5])$value
MeanAgeD <- (MeanAgeDmale+MeanAgeDfemale)/2
mean.agesI[s] <- MeanAgeD}
Below we show the mean age of an animal from a population after stretching, with dashed lines representing the stretching required for:
#see the result
plot(stretching,mean.agesI,ylim=c(10,80))
abline(h=MeanAgeE,lty=2)
abline(v=0.46,lty=2)
abline(h=24.85,lty=3)
abline(v=0.6,lty=3)
legend("topright",legend=c("Elephant","Sperm whale"),lty=2:3,inset=0.05)
So the corresponding look and feel for the proportion of animals surviving to a given age, considering both
is shown in the figure below
ages <- 0:150
# Proportion surviving to age x
plot(ages,lx(ages*0.46,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5]),type="l",lty=1,main="Proportion surviving to age x",ylab="Proportion",xlab="Age",col=4)
lines(ages,lx(ages*0.46,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5]),lty=2,col=4)
lines(ages,lx(ages,eleSp[1,1]/eleSp[1,4],eleSp[1,3],eleSp[1,2]/eleSp[1,5],eleSp[1,4],eleSp[1,5]),col=3)
lines(ages,lx(ages,eleSp[2,1]/eleSp[2,4],eleSp[2,3],eleSp[2,2]/eleSp[2,5],eleSp[2,4],eleSp[2,5]),lty=2,col=3)
legend("top",legend=c("Males E","Females E","Males D2E","Females D2E"),lty=c(2,1,2,1),col=c(3,3,4,4), inset=0.02)
The conclusion of this exercise is somewhat expected. While we manage to extend the dolphins function, because the dolphins function had higher values for low age survival probabilities, and the convex shape, the only way we manage to extend the function to get the same mean age is by overestimating the ability to get really long lived animals, which is a consequence from a wider plateau in survival probabilities for adults. This is visible in the corresponding survival functions
ages <- 0:200
# Proportion surviving to age x
plot(ages,px(ages*0.46,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5],scaling=0.46),type="l",lty=1,main="Survival probability at age x",ylab="Proportion",xlab="Age",col=4,ylim=c(0,1))
lines(ages,px(ages*0.46,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5],scaling=0.46),lty=2,col=4)
lines(ages,px(ages,eleSp[1,1]/eleSp[1,4],eleSp[1,3],eleSp[1,2]/eleSp[1,5],eleSp[1,4],eleSp[1,5]),col=3)
lines(ages,px(ages,eleSp[2,1]/eleSp[2,4],eleSp[2,3],eleSp[2,2]/eleSp[2,5],eleSp[2,4],eleSp[2,5]),lty=2,col=3)
legend("bottomleft",legend=c("Males E","Females E","Males D2E","Females D2E"),lty=c(2,1,2,1),col=c(3,3,4,4), inset=0.02)
Just for an additional comparison, below is the comparison of the dolphin model and its extended version
ages <- 0:200
# Proportion surviving to age x
plot(ages,px(ages*0.46,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5],scaling=0.46),type="l",lty=1,main="Survival probability at age x",ylab="Proportion",xlab="Age",col=4,ylim=c(0,1))
lines(ages,px(ages*0.46,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5],scaling=0.46),lty=2,col=4)
lines(ages,px(ages,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5]),col=3)
lines(ages,px(ages,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5]),lty=2,col=3)
legend("bottomleft",legend=c("Males E","Females E","Males D","Females D"),lty=c(2,1,2,1),col=c(4,4,3,3), inset=0.02)
This seems to imply that simply extending the function in this way is not enough, since you can’t really go above the maximum survival probability from the original function. If the maximum survival probability, conditional on age, was say \(p_{max}\), then extending it this way, the new function cannot have any age with a survival probability higher than said \(p_{max}\).
And the reason I don’t like that is because I suspect that \(p_{max}\) will be correlated with the mean age at death, in that higher longevity species will have higher \(p_{max}\) than short lived ones. And hence, expanding sideways does not quite cut it, one would need to extend it upwards too…
The proportion surviving to age x corresponding for “dolphins extended to be sperm whales”, as well as survival probabilities at age x, using the mean age of 24.85 as described above, is shown below. For comparison, the values for the stage specific survival probabilities considered for the sperm whales are shown as dashed black lines
ages <- 0:100
# Proportion surviving to age x
plot(ages,lx(ages*0.6,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5]),type="l",lty=1,main="Proportion surviving to age x",ylab="Proportion",xlab="Age",col=4)
lines(ages,lx(ages*0.6,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5]),lty=2,col=4)
legend("top",legend=c("Males Pmac","Females Pmac"),lty=2,col=4, inset=0.02)
# Survival at age x
ages <- 0:200
plot(ages,px(ages*0.6,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5],scaling=0.6),type="l",lty=1,main="Survival probability at age x",ylab="Proportion",xlab="Age",col=4,ylim=c(0,1))
lines(ages,px(ages*0.6,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5],scaling=0.6),lty=2,col=4)
legend("bottomleft",legend=c("Males Pmac","Females Pmac","Stage structured values"),lty=c(1,2,2),col=c(4,4,1), inset=0.02,lwd=c(1,1,3))
#from Lance's files
PmacScalf <- 0.907
PmacSjuv <- 0.9424
PmacSFemale <- 0.9777
PmacSMale <- 0.918
#add lines
segments(x0=0,x1=2,y0=PmacScalf,y1=PmacScalf,lty=2,lwd=3)
segments(x0=2,x1=9,y0=PmacSjuv,y1=PmacSjuv,lty=2,lwd=3)
segments(x0=9,x1=100,y0=PmacSFemale,y1=PmacSFemale,lty=1,lwd=3)
segments(x0=9,x1=100,y0=PmacSMale,y1=PmacSMale,lty=2,lwd=3)
For an additional comparison, below is the comparison of the dolphin model and the sperm whale model
ages <- 0:200
# Survival probability at age x - Pmac
plot(ages,px(ages*0.6,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5],scaling=0.6),type="l",lty=1,main="Survival probability at age x",ylab="Proportion",xlab="Age",col=4,ylim=c(0,1))
lines(ages,px(ages*0.6,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5],scaling=0.6),lty=2,col=4)
# Survival probability at age x - dolphins
lines(ages,px(ages,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5]),col=3)
lines(ages,px(ages,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5]),lty=2,col=3)
legend("topright",legend=c("Males Pmac","Females Pmac","Males D","Females D"),lty=c(2,1,2,1),col=c(4,4,3,3), inset=0.02)
Motivated by the email that LS sent Sun 4/12/2020 4:21 AM, we plot below the corresponding lx
function taking the stage specific survival probabilities and associating them to each age - note there is a slight mismatch in doing it this way, as this is not a direct 1:1 translation from a stage to an age strutcured model.
PmacPaliveF=c(PmacScalf,PmacScalf,PmacSjuv,PmacSjuv,PmacSjuv,PmacSjuv,PmacSjuv,PmacSjuv,PmacSjuv,rep(PmacSFemale,141))
PmacPaliveM=c(PmacScalf,PmacScalf,PmacSjuv,PmacSjuv,PmacSjuv,PmacSjuv,PmacSjuv,PmacSjuv,PmacSjuv,rep(PmacSMale,141))
PmacPaliveF=cumprod(PmacPaliveF)
PmacPaliveM=cumprod(PmacPaliveM)
#the Pmac lx from the age structured survivals
plot(1:150,PmacPaliveF,col="black",type="l",ylim=c(0,1),main="Proportion surviving to age x",ylab="Proportion",xlab="Age")
lines(1:150,PmacPaliveM,lty=2,col="black")
#now the Pmac after streching
lines(ages,lx(ages*0.6,pf[1,1],pf[1,2],pf[1,3],pf[1,4],pf[1,5]),type="l",lty=1,col=4)
lines(ages,lx(ages*0.6,pm[1,1],pm[1,2],pm[1,3],pm[1,4],pm[1,5]),lty=2,col=4)
legend("topright",legend=c("Males Pmac age extended","Females Pmac age extended","Males stage structured","Females stage structured"),lty=c(2,1,2,1),col=c(4,4,1,1), inset=0.02)
What the above figure represents relates to what LS reported. The median age for Pmac (using an age structured model) is actually lower than that of dolphins (using the age structured model). When we consider the mean age, that is larger in the Pmac than in the BSE BND model. This comes at the cost of there bing lot’s of unrealistically old animals under the sperm whale stage structured model.
Not sure this approximation is good at all. Need to think further about this!
While reading Barlow & Boveng 1991 Modeling age-specific mortality for marine mammal populations. Marine Mammal Science 7: 50-65 I realized that actually, Jay has already thought about all of this. How am I not surprised that this was nothing new…! In a section below we actually reproduce some of Jay’s insights as they might be pertinent for what we are doing.
While reading Barlow & Boveng 1991 I was looking at this bit in particular:
Screenshot of expression in BBarlow and Boveng 1991.
and thinking this would provide a way to calculate the mean age of an animal sampled at random at the time of the oil spill, which is something that will be required to evaluate how long the recovery will take for animals that do recover from the oil after being exposed. To understand why we need this value the reader would have to read section 2.7 in CARMMHA’s internal document “How to land it - a road map to finishing CARMMHA”.
Given the above, that mean age would be given by
\[ A_m = E(x)= \int x~f_{pdf}(x) dx = \int x \frac{ l(x) }{\int l(x) dx}dx=\frac{ \int x l(x) dx}{\int l(x) dx}\]
We first define a function xlx
and then we use it
Amf <- integrate(xlx,lower=0,upper=500,a1=pf[1,1],a2=pf[1,2],a3=pf[1,3],b1=pf[1,4],b3= pf[1,5])$value/integrate(lx,lower=0,upper=500,a1=pf[1,1],a2=pf[1,2],a3=pf[1,3],b1=pf[1,4],b3= pf[1,5])$value
Amm<- integrate(xlx,lower=0,upper=500,a1=pm[1,1],a2=pm[1,2],a3=pm[1,3],b1=pm[1,4],b3= pm[1,5])$value/integrate(lx,lower=0,upper=500,a1=pm[1,1],a2=pm[1,2],a3=pm[1,3],b1=pm[1,4],b3= pm[1,5])$value
#50% sex ratio
Am <- (Amf+Amm)/2
The mean age of an animal taken from the population at any one time, hence at the time of the oil spill, would be 14.3 for females, 11.7 for males, and hence under a 1:1 sex ratio, 13 overall. Note that, given the number obtained above for the mean age at death (15.3), this poses a question. It seems like on average, at any one time, you only have 2.3 years to live…
This is a bit counter-intuitive…? What is my intuition failing to grasp?
A key thing that seems odd is that while the mean age of an animal at death does not depend on fecundity, the mean age of an animal alive at any one time should depend on fecundity too. Or at the very least, the above expression from Jay depends on it, because the stable age distribution will depend on fecundity, but in the above, I’ve not really used fecundity anywhere.
The Barlow and Boveng 1991 paper is interesting for another reason, that links to a topic discussed above. As per the screenshot below, the curve for human females, with the longest longevity (81 years, 34 for monkeys and 18 for seals), looks convex, as for the elephants.
Is that a “feature” of longer lived animals, and if so, is our extending the function approach not capable to capture that?
Screenshot of figure 1 in BBarlow and Boveng 1991.
These are pending questions collated from all the above that I feel would be great to answer:
How to incorporate the variability in the mean age when extending the function? Back to question 1
Why are the blue and green functions not exactly the same? Back to question 2
Why extending it the way we are doing it might not quite make it? Back to question 4
What am I missing here, how can the mean age at death and the age of a random individual from the population at any one time be so similar? Back to question 5
Did I ignored fecundity here? Back to question 6
Are longed lived species different from the others? Back to question 7
Is the age structured model sensible given the stage structured stage-conditional survival values? Back to question 8
Scale by a mean value or scale by an extreme quantile. I wonder what would Jay say? Back to question 9
Barlow, J. & Boveng, P. (1991). Modeling age-specific mortality for marine mammal populations. Marine Mammal Science, 7, 50–65.
Caswell, H., Brault, S., Read, A.J. & Smith, T.D. (1998). Harbor porpoise and fisheries: An uncertainty analysis of incidental mortality. Ecological Applications, 8, 1226–1238.
Eakin, T. (1994). Intrinsic time scaling in survival analysis: Application to biological populations. Bulletin of Mathematical Biology, 56, 1121–1141.
Krementz, D.G., Sauer, J.R. & Nichols, J.D. (1989). Model-based estimates of annual survival rate are preferable to observed maximum lifespan statistics for use in comparative life-history studies. Oikos, 56, 203.
Lahdenpera, M., Mar, K.U., Courtiol, A. & Lummaa, V. (2018). Differences in age-specific mortality between wild-caught and captive-born asian elephants. Nature Communications, 9, 3023.
Schwacke, L.H., Thomas, L., Wells, R.S., McFee, W.E., Hohn, A.A., Mullin, K.D., Zolman, E.S., Quigley, B.M., Rowles, T.K. & Schwacke, J.H. (2017). Quantifying injury to common bottlenose dolphins from the Deepwater Horizon oil spill using an age-, sex- and class-structured population model. Endangered Species Research, 33, 265–279.