An R function and code to estimate parameters of mortality models with maximum likelihood.

MortalityModels.R

#A likelihood function for GOMPERTZ, GOMPERTZ-MAKEHAM and SILER models.

likeLT <- function(lifetable,pars,type="GO"){

    # Extract data from life table
    Dx = lifetable$Dx
    Nx = lifetable$Nx
    StartInt = lifetable$StartAge
    EndInt = lifetable$EndAge
    LT.Type = as.character(lifetable$Type[1])

    # Read in parameters.
    if(type=="GM"){a = pars[1]; b = pars[2]; c = pars[3]}
    if(type=="GO"){a = pars[1]; b = pars[2]}
    # Siler has 5 free parameters: a1, b1, c, a2, b2.
    if(type=="SI"){a1 = pars[1]; b1 = pars[2]; c = pars[3]; a2 = pars[4]; b2 = pars[5]}

    # Calculated cumulative hazard
    if(type=="GO"){ H = (a / b) * (exp(b * EndInt) - exp(b * StartInt))}
    if(type=="GM"){ H = (a / b) * (exp(b * EndInt) - exp(b * StartInt)) + (c * (EndInt - StartInt))}
    if(type=="SI"){ H = -(a1 / b1) * (exp(-b1 * EndInt) - exp(-b1 * StartInt)) + c * (EndInt-StartInt) + (a2 / b2) * (exp(b2 * EndInt) - exp(b2 * StartInt))}

    # Calculate Likelihood
    # H is the cumulative hazard over the interval. The probability of dying in
    # the interval is q = 1 - exp(-H); this is what the binomial needs (passing
    # the raw cumulative hazard could give a probability > 1, hence NaN).
    if(LT.Type == "Cohort"){ LH = -sum(dbinom(Dx,Nx,1 - exp(-H),log=TRUE))}
    if(LT.Type == "Period"){ LH = -sum(dpois(Dx,H*Nx,log=TRUE))}

    return(LH)
    }

#You can use optim() to optimise to get a maximum likelihood fit.
#(Supply as many starting values as the model has parameters:
# GO = 2, GM = 3, SI = 5.)
out     = optim(par=c(0.1,0.1), likeLT, lifetable=lifetable, type="GO")
MLE     = out$par