#' 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.
#'
#' @seealso \code{\link{fit.StMoMo}}, \code{\link{lc}}, \code{\link{cbd}},
#' \code{\link{apc}}, \code{\link{rh}}, \code{\link{m6}}, \code{\link{m7}},
#' \code{\link{m8}}
#'
#' @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
#' that deaths follow a Binomial distribution and a logit link.
#'
#' @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{link}{a character string defining the link function of the model.}
#'
#' \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
#---------------------------------------------------------------------
#link
link <- match.arg(link)
if (link == "logit") {
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.")
}
out <- list(link = link, staticAgeFun = staticAgeFun,
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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.