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

MortalityModels.R

#A likelihood function for GOMPERTZ, GOMPERTZ-MAKEHANM 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]}
    if(type=="SI"){a1 = pars[1]; b1 = pars[2]; c = pars[3]; a2 = pars[1]; b2 = pars[2]}

    # 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
    if(LT.Type == "Cohort"){ LH = -sum(dbinom(Dx,Nx,H,log=TRUE))}
    if(LT.Type == "Period"){ LH = -sum(dpois(Dx,H*Nx,log=TRUE))}
   
    return(LH)
    }    
    
#You can use optim() to optimes to get a maximum likelihood fit.
out     = optim(par=c(0.1,0.1,0.05), likeLT, lifetable=lifetable, type="GO")
MLE     = out$par