# R/StMoMo.R In StMoMo: Stochastic Mortality Modelling

#### Documented in StMoMoStMoMo

#' Create a new Stochastic Mortality Model
#'
#' Initialises a StMoMo object which represents a Generalised
#' Age-Period-Cohort Stochastic Mortality Model.
#'
#' @details
#'
#' R implementation of the family of Generalised Age-Period-Cohort
#' stochastic mortality models.  This family of models encompasses
#' many models proposed in the literature including the well-known
#' Lee-Carter model, CBD model and APC model.
#'
#' \code{StMoMo} defines an abstract representation of a Generalised
#' Age-Period-Cohort (GAPC) Stochastic model that fits within the
#' general class of generalised non-linear models defined as follows
#' \deqn{D_{xt} \sim Poisson(E_{xt}\mu_{xt}), D_{xt} \sim
#' Binomial(E_{xt},q_{xt})}
#' \deqn{\eta_{xt} = \log \mu_{xt}, \eta_{xt} = \mathrm{logit}\,
#' q_{xt}} \deqn{\eta_{xt} = \alpha_x + \sum_{i=1}^N \beta_x^{(i)}\kappa_t^{(i)}
#' + \beta_x^{(0)}\gamma_{t-x}} \deqn{v: \{\alpha_{x}, \beta_x^{(1)},...,
#' \beta_x^{(N)}, \kappa_t^{(1)},..., \kappa_t^{(N)}, \beta_x^{(0)}, \gamma_{t-x}\} \mapsto
#' \{\alpha_{x}, \beta_x^{(1)},..., \beta_x^{(N)}, \kappa_t^{(1)},..., \kappa_t^{(N)},
#' \beta_x^{(0)}, \gamma_{t-x}\},} where
#' \itemize{
#'  \item \eqn{\alpha_x} is a static age function;
#'  \item \eqn{\beta_x^{(i)}\kappa_t^{(i)}, i = 1,..N}, are age/period terms;
#'  \item \eqn{\beta_x^{(0)}\gamma_{t-x}} is the age/cohort term; and
#'  \item \eqn{v} is a  function defining the identifiability constraints of
#'  the model.
#' }
#' Most Stochastic mortality models proposed in the literature can be cast to
#' this representation (See Hunt and Blake (2015)).
#'
#' Parametric age functions should be scalar functions of the form
#' \code{f <- function(x, ages)} taking a scalar age \code{x} and a vector
#' of model fitting \code{ages} (see examples below).
#'
#' Do to limitation of functions \code{\link[gnm]{gnm}} within package
#' \pkg{gnm}, which is used for fitting \code{"StMoMo"} objects to data
#' (see \code{\link{fit.StMoMo}}), models combining parametric and
#' non-parametric age-modulating functions are not supported at the moment.
#'
#'
#' @param link defines the link function and random component associated with
#'   the mortality model. \code{"log"} would assume that deaths follow a
#'   Poisson distribution and use a log link while \code{"logit"} would assume
#'
#' @param staticAgeFun logical value indicating if a static age function
#'   \eqn{\alpha_x} is to be included.
#'
#' @param periodAgeFun  a list of length \eqn{N} with the definitions of the
#'   period age modulating parameters \eqn{\beta_x^{(i)}}. Each entry can take
#'   values: \code{"NP"} for non-parametric age terms, \code{"1"} for
#'   \eqn{\beta_x^{(i)}=1} or a predefined parametric function of
#'   age (see details). Set this to \code{NULL} if there are no period terms
#'   in the model.
#'
#' @param cohortAgeFun defines the cohort age modulating parameter
#'   \eqn{\beta_x^{(0)}}. It can take values: \code{"NP"} for non-parametric
#'   age terms, \code{"1"} for \eqn{\beta_x^{(0)}=1}, a predefined parametric
#'   function of age (see details) or \code{NULL} if there is no cohort effect.
#'
#' @param  constFun function defining the identifiability constraints of the
#'   model. It must be a function of the form
#'   \code{constFun <- function(ax, bx, kt, b0x, gc, wxt, ages)} taking a set
#'   of fitted model parameters and returning a list
#'   \code{list(ax = ax, bx = bx, kt = kt, b0x = b0x, gc = gc)}
#'   of the model parameters with the identifiability constraints applied. If
#'   omitted no identifiability constraints are applied to the model.
#'
#' @return A list with class \code{"StMoMo"} with components:
#'
#'
#'   \item{staticAgeFun}{a logical value indicating if the model has a static
#'   age function.}
#'
#'   \item{periodAgeFun}{a list defining the period age modulating parameters.}
#'
#'   \item{cohortAgeFun}{an object defining the cohort age modulating
#'   parameters.}
#'
#'   \item{constFun}{a function defining the identifiability constraints.}
#'
#'   \item{N}{an integer specifying The number of age-period terms in the
#'   model.}
#'
#'   \item{textFormula}{a character string of the model formula.}
#'
#'   \item{gnmFormula}{a formula that can be used for fitting the model with
#'   package \pkg{gnm}.}
#'
#' @references
#'
#' Plat, R. (2009). On stochastic mortality modeling. Insurance: Mathematics
#'   and Economics, 45(3), 393-404.
#'
#' Hunt, A., & Blake, D. (2015). On the Structure and Classification of
#' Mortality Models Mortality Models. Pension Institute Working Paper.
#' \url{http://www.pensions-institute.org/workingpapers/wp1506.pdf}.
#'
#'
#' @examples
#' #Lee-Carter model
#' constLC <- function(ax, bx, kt, b0x, gc, wxt, ages) {
#'  c1 <- mean(kt[1, ], na.rm = TRUE)
#'  c2 <- sum(bx[, 1], na.rm = TRUE)
#'  list(ax = ax + c1 * bx, bx = bx / c2, kt = c2 * (kt - c1))
#' }
#' LC <- StMoMo(link = "log", staticAgeFun = TRUE, periodAgeFun = "NP",
#'              constFun = constLC)
#'
#' plot(fit(LC, data = EWMaleData, ages.fit = 55:89))
#'
#'
#' #CBD model
#' f2 <- function(x, ages) x - mean(ages)
#' CBD <- StMoMo(link = "logit", staticAgeFun = FALSE,
#'               periodAgeFun = c("1", f2))
#' plot(fit(CBD, data = EWMaleData, ages.fit = 55:89))
#'
#' #Reduced Plat model (Plat, 2009)
#' f2 <- function(x, ages) mean(ages) - x
#' constPlat <- function(ax, bx, kt, b0x, gc, wxt, ages) {
#'  nYears <- dim(wxt)[2]
#'  x <- ages
#'  t <- 1:nYears
#'  c <- (1 - tail(ages, 1)):(nYears - ages[1])
#'  xbar <- mean(x)
#'  #nsum g(c)=0, nsum cg(c)=0, nsum c^2g(c)=0
#'  phiReg <- lm(gc ~ 1 + c + I(c^2), na.action = na.omit)
#'  phi <- coef(phiReg)
#'  gc <- gc - phi[1] - phi[2] * c - phi[3] * c^2
#'  kt[2, ] <- kt[2, ] + 2 * phi[3] * t
#'  kt[1, ] <- kt[1, ] + phi[2] * t + phi[3] * (t^2 - 2 * xbar * t)
#'  ax <- ax + phi[1] - phi[2] * x + phi[3] * x^2
#'  #nsum kt[i, ] = 0
#'  ci <- rowMeans(kt, na.rm = TRUE)
#'  ax <- ax + ci[1] + ci[2] * (xbar - x)
#'  kt[1, ] <- kt[1, ] - ci[1]
#'  kt[2, ] <- kt[2, ] - ci[2]
#'  list(ax = ax, bx = bx, kt = kt, b0x = b0x, gc = gc)
#' }
#' PLAT <- StMoMo(link = "log", staticAgeFun = TRUE,
#'                periodAgeFun = c("1", f2), cohortAgeFun = "1",
#'                constFun = constPlat)
#'
#' plot(fit(PLAT, data = EWMaleData, ages.fit = 55:89))
#'
#' #Models not supported
#' \dontrun{
#' MnotSup1 <- StMoMo(periodAgeFun = c(f2, "NP"))
#' MnotSup1 <- StMoMo(periodAgeFun = f2, cohortAgeFun = "NP")
#' }
#' @export
StMoMo  <- function(link = c("log", "logit"), staticAgeFun = TRUE,
periodAgeFun = 'NP', cohortAgeFun = NULL,
constFun = function(ax, bx, kt, b0x, gc, wxt, ages)
list(ax = ax, bx = bx, kt = kt, b0x = b0x, gc = gc)) {

#---------------------------------------------------------------------
# Construct the model formula
#---------------------------------------------------------------------

gnmFormula <- "D/E ~ -1 + offset(o)"
textFormula <- "logit q[x,t]"
} else if (link == "log") {
gnmFormula <- "D ~ -1 + offset(log(E)) + offset(o)"
textFormula <- "log m[x,t]"
}

#static life table
termSep <- " = "
if (staticAgeFun == TRUE) {
gnmFormula <- paste(gnmFormula, " + factor(x)", sep = "")
textFormula <- paste(textFormula, termSep, "a[x]", sep = "")
termSep <- " + "
}

#age-period terms
periodAgeFun <- c(periodAgeFun)
N <- length(periodAgeFun)
if (N>0) {
for (i in 1:N) {
if (is.function(periodAgeFun[[i]])) {
gnmFormula <- paste(gnmFormula, " + B", i, ":factor(t)", sep = "")
textFormula <- paste(textFormula, termSep, "f", i, "[x] k", i, "[t]",
sep = "")
} else if (periodAgeFun[[i]] == "1") {
gnmFormula <- paste(gnmFormula, " + factor(t)", sep = "")
textFormula <- paste(textFormula, termSep, "k", i, "[t]", sep = "")
} else if ( periodAgeFun[[i]] == "NP" ) {
ind <- i:N
inst <- 1 + sum(periodAgeFun[-ind] == periodAgeFun[[i]])
gnmFormula <- paste(gnmFormula,
" + Mult(factor(x), factor(t), inst = ", inst,
")", sep = "")
textFormula <- paste(textFormula, termSep, "b", i, "[x] k", i, "[t]",
sep = "")
} else {
stop("Not appropriate period age modulating function definition")
}
termSep <- " + "
}
totalP <- sum(sapply(periodAgeFun, is.function))
totalNP <- sum(periodAgeFun == "NP")
} else {
totalP <- 0
totalNP <- 0
}

#age-cohort term
if (!is.null(cohortAgeFun)) {
if (is.function(cohortAgeFun)) {
gnmFormula <- paste(gnmFormula, " + B0:factor(c)", sep = "")
textFormula <- paste(textFormula, termSep, "f0[x] g[t-x]", sep = "")
} else if (cohortAgeFun == "1") {
gnmFormula <- paste(gnmFormula, " + factor(c)", sep = "")
textFormula <- paste(textFormula, termSep, "g[t-x]", sep = "")
} else if (cohortAgeFun == "NP") {
gnmFormula <- paste(gnmFormula, " + Mult(factor(x), factor(c))",
sep = "")
textFormula <- paste(textFormula, termSep, "b0[x] g[t-x]", sep = "")
} else {
stop("Not appropriate cohort age modulating function definition")
}
totalP <- totalP + is.function(cohortAgeFun)
totalNP <-totalNP + sum(list(cohortAgeFun) == "NP")
}

#Check constraint on number of simulataneous parametric and non-parametric
#age functions
if (totalP >= 1 && totalNP >= 1) {
stop("Models combining parametric and non-parametric age-modulating functions are not supported at the moment.")
}

periodAgeFun = periodAgeFun, cohortAgeFun = cohortAgeFun,
constFun = constFun, N = N, textFormula = textFormula,
gnmFormula = gnmFormula)
class(out) <- "StMoMo"
out
}

#' @export
print.StMoMo <- function(x, ...) {
if (x$link == "logit") { cat("Binomial model with predictor: ") } else { cat("Poisson model with predictor: ") } cat("") cat(x$textFormula)
}


## Try the StMoMo package in your browser

Any scripts or data that you put into this service are public.

StMoMo documentation built on May 2, 2019, 11:42 a.m.