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

#### Documented in fit.rhrh

#' Create a Renshaw and Haberman (Lee-Carter with cohorts) mortality model
#'
#' Utility function to initialise a \code{StMoMo} object representing a
#' Renshaw and Haberman (Lee-Carter with cohorts) mortality model introduced
#' in Renshaw and Haberman (2006).
#'
#' The created model is either a log-Poisson or a  logit-Binomial version
#' of the Renshaw and Haberman model which has predictor structure
#' \deqn{\eta_{xt} = \alpha_x + \beta^{(1)}_x\kappa_t + \beta^{(0)} \gamma_{t-x}.}
#' or
#' \deqn{\eta_{xt} = \alpha_x + \beta^{(1)}_x\kappa_t + \gamma_{t-x}.}
#' depending on the value of argument \code{cohortAgeFun}.
#'
#' To ensure identifiability the following constraints are imposed
#' \deqn{\sum_t\kappa_t = 0, \sum_x\beta^{(1)}_x = 1, \sum_c\gamma_c = 0}
#' plus
#' \deqn{\sum_x\beta^{(0)}_x = 1}
#' if \code{cohortAgeFun = "NP"}
#'
#' In addition, if \code{approxConst=TRUE} then the approximate
#' identifiability constraint
#' \deqn{\sum_c (c-\bar{c})\gamma_c = 0}
#' is applied to improve the stability and robustness of the model
#' (see Hunt and Villegas (2015)).
#'
#' By default \eqn{\beta^{(0)}_x = 1} as this model has shown to be more
#' stable (see Haberman and Renshaw (2011) and Hunt and Villegas (2015)).
#'
#' @inheritParams StMoMo
#' @param cohortAgeFun defines the cohort age modulating parameter
#'   \eqn{\beta_x^{(0)}}. It can take values: \code{"NP"} for a non-parametric age
#'   term or \code{"1"} for \eqn{\beta_x^{(0)}=1} (the default).
#'
#' @param approxConst defines if the approximate identifiability constraint of
#' Hunt and Villegas (2015) is applied or not. If \code{TRUE}, the output object
#' is of class \code{rh} and subsequent model fitting is performed with
#' \code{\link{fit.rh}}. If \code{FALSE}, the output object is of class
#' \code{StMoMo} and subsequent model fitting is performed with
#'
#' @return An object of class \code{"StMoMo"} or \code{"rh"}.
#'
#'
#' @references
#'
#' Haberman, S., & Renshaw, A. (2011). A comparative study of parametric
#' mortality projection models. Insurance: Mathematics and Economics,
#' 48(1), 35-55.
#'
#' Hunt, A., & Villegas, A. M. (2015). Robustness and convergence in the
#' Lee-Carter model with cohorts. Insurance: Mathematics and Economics,
#' 64, 186-202.
#'
#' Renshaw, A. E., & Haberman, S. (2006). A cohort-based extension to the
#' Lee-Carter model for mortality reduction factors.
#' Insurance: Mathematics and Economics, 38(3), 556-570.
#'
#' @examples
#'
#' LCfit <-  fit(lc(), data = EWMaleData, ages.fit = 55:89)
#' wxt <- genWeightMat(55:89,  EWMaleData$years, clip = 3) #' RHfit <- fit(rh(), data = EWMaleData, ages.fit = 55:89, wxt = wxt, #' start.ax = LCfit$ax, start.bx = LCfit$bx, start.kt = LCfit$kt)
#' plot(RHfit)
#'
#' #Impose approximate constraint as in Hunt and Villegas (2015)
#' \dontrun{
#' RHapprox <- rh(approxConst = TRUE)
#' RHapproxfit <- fit(RHapprox, data = EWMaleData, ages.fit = 55:89,
#'                     wxt = wxt)
#' plot(RHapproxfit)
#' }
#'
#'
#' @export
rh <- function(link = c("log", "logit"), cohortAgeFun = c("1", "NP"),
approxConst = FALSE) {
cohortAgeFun <- match.arg(cohortAgeFun)
constRHgeneral <- function(ax, bx, kt, b0x, gc, wxt, ages, cohortAgeFun, approxConst) {
#\sum k[t] = 0
c1 <- mean(kt[1, ], na.rm = TRUE)
ax <- ax + c1 * bx[, 1]
kt[1, ] <- kt[1, ] - c1
#\sum b[x, 1] = 0
c2 <- sum(bx[, 1], na.rm = TRUE)
bx[, 1] <- bx[, 1] / c2
kt[1, ] <- kt[1, ] * c2
#\sum g[c] = 0
c3 <- mean(gc, na.rm = TRUE)
ax <- ax + c3 * b0x
gc <- gc - c3
#\sum b0[x] = 1
if (cohortAgeFun == "NP") {
c4 <- sum(b0x, na.rm = TRUE)
b0x <- b0x / c4
gc <- gc * c4
}

list(ax = ax, bx = bx, kt = kt, b0x = b0x, gc = gc)
}
constRH <- function(ax, bx, kt, b0x, gc, wxt, ages) {
constRHgeneral(ax, bx, kt, b0x, gc, wxt, ages, cohortAgeFun, approxConst)
}
cohortAgeFun = cohortAgeFun, constFun = constRH)
if (!is.logical(approxConst)) {
stop( "approxConst should be a logical variable")
}
out$approxConst <- approxConst if(approxConst) class(out) <- c("rh", "StMoMo") out } #' Fit a Renshaw and Haberman (Lee-Carter with cohorts) mortality model #' #' Fit a Renshaw and Haberman (Lee-Carter with cohorts) mortality model #' using the iterative Newton-Raphson procedure presented in Algorithm 1 #' of Hunt and Villegas (2015). This approach helps solve the #' well-known robustness and converges issues of the Lee-Carter model #' with cohort-effects. #' #' @inheritParams fit.StMoMo #' @param object an object of class \code{"rh"} created with function #' \code{\link{rh}}. #' @param tolerance a positive numeric value specifying the tolerance #' level for convergence. #' @param iterMax a positive integer specifying the maximum number of #' iterations to perform. #' @param ... arguments to be passed to or from other methods. #' #' @return #' #' \item{model}{ the object of class \code{"rh"} defining the fitted #' stochastic mortality model.} #' #' \item{ax}{ vector with the fitted values of the static age function #' \eqn{\alpha_x}. If the model does not have a static age function or #' failed to fit this is set to \code{NULL}.} #' #' \item{bx}{ matrix with the values of the period age-modulating functions #' \eqn{\beta_x^{(i)}, i=1, ..., N}. If the \eqn{i}-th age-modulating #' function is non-parametric (e.g. as in the Lee-Carter model) #' \code{bx[, i]} contains the estimated values. If the model does not have #' any age-period terms (i.e. \eqn{N=0}) or failed to fit this is set to #' \code{NULL}.} #' #' \item{kt}{ matrix with the values of the fitted period indexes #' \eqn{\kappa_t^{(i)}, i=1, ..., N}. \code{kt[i, ]} contains the estimated #' values of the \eqn{i}-th period index. If the model does not have any #' age-period terms (i.e. \eqn{N=0}) or failed to fit this is set to #' \code{NULL}.} #' \item{b0x}{ vector with the values of the cohort age-modulating function #' \eqn{\beta_x^{(0)}}. If the age-modulating function is non-parametric #' \code{b0x} contains the estimated values. If the model does not have a #' cohort effect or failed to fit this is set to \code{NULL}.} #' #' \item{gc}{ vector with the fitted cohort index \eqn{\gamma_{c}}. #' If the model does not have a cohort effect or failed to fit this is set #' to \code{NULL}.} #' #' \item{data}{ StMoMoData object provided for fitting the model.} #' #' \item{Dxt}{ matrix of deaths used in the fitting.} #' #' \item{Ext}{ matrix of exposures used in the fitting.} #' #' \item{oxt}{ matrix of known offset values used in the fitting.} #' #' \item{wxt}{ matrix of 0-1 weights used in the fitting.} #' #' \item{ages}{ vector of ages used in the fitting.} #' #' \item{years}{ vector of years used in the fitting.} #' #' \item{cohorts}{ vector of cohorts used in the fitting.} #' #' \item{fittingModel}{ output from the iterative fitting algorithm.} #' #' \item{loglik}{ log-likelihood of the model. If the fitting failed to #' converge this is set to \code{NULL}.} #' #' \item{deviance}{ deviance of the model. If the fitting failed to #' converge this is set to \code{NULL}.} #' #' \item{npar}{ effective number of parameters in the model. If the fitting #' failed to converge this is set to \code{NULL}.} #' #' \item{nobs}{ number of observations in the model fit. If the fitting #' failed to converge this is set to \code{NULL}.} #' #' \item{fail}{ \code{TRUE} if a model could not be fitted and #' \code{FALSE} otherwise.} #' #' \item{conv}{ \code{TRUE} if the model fitting converged and #' \code{FALSE} if it didn't.} #' #' @references #' #' Hunt, A., & Villegas, A. M. (2015). Robustness and convergence in the #' Lee-Carter model with cohorts. Insurance: Mathematics and Economics, #' 64, 186-202. #' #' @examples #' #' LCfit <- fit(lc(), data = EWMaleData, ages.fit = 55:89) #' wxt <- genWeightMat(55:89, EWMaleData$years, clip = 3)
#' RHfit <- fit(rh(), data = EWMaleData, ages.fit = 55:89,
#'              wxt = wxt, start.ax = LCfit$ax, #' start.bx = LCfit$bx, start.kt = LCfit$kt) #' plot(RHfit) #' #' #Impose approximate constraint as in Hunt and Villegas (2015) #' \dontrun{ #' RHapprox <- rh(approxConst = TRUE) #' RHapproxfit <- fit(RHapprox, data = EWMaleData, ages.fit = 55:89, #' wxt = wxt) #' plot(RHapproxfit) #' } #' #' @export fit.rh <- function(object, data = NULL, Dxt = NULL, Ext = NULL, ages = NULL, years = NULL, ages.fit = NULL, years.fit = NULL, oxt = NULL, wxt = NULL, start.ax = NULL, start.bx = NULL, start.kt = NULL, start.b0x = NULL, start.gc = NULL, verbose = TRUE, tolerance = 1e-04, iterMax = 10000, ...) { #Hack to remove notes in CRAN check x <- NULL # Select data from data or from Dxt, Ext, ages, years if(!is.null(data)) { if (class(data) != "StMoMoData") stop("Argument data needs to be of class StMoMoData.") Dxt <- data$Dxt
Ext <- data$Ext ages <- data$ages
years <- data$years } else { if (is.null(Dxt) || is.null(Ext)) stop("Either argument data or arguments Dxt and Ext need to be provided.") if (is.null(ages)) ages <- 1:nrow(Dxt) if (is.null(years)) years <- 1:ncol(Dxt) data <- structure(list(Dxt = Dxt, Ext = Ext, ages = ages, years = years, type = ifelse(object$link == "log", "central", "initial"),
series = "unknown", label = "unknown"), class = "StMoMoData")
}
if (is.null(ages.fit)) ages.fit <- ages
if (is.null(years.fit)) years.fit <- years

# Construct fitting data

nAges <- length(ages)
nYears <- length(years)
Dxt <- as.matrix(Dxt)
if (nrow(Dxt) != nAges ||  ncol(Dxt) != nYears) {
stop( "Mismatch between the dimension of Dxt and the
number of years or ages")
}
rownames(Dxt) <- ages
colnames(Dxt) <- years

Ext <- as.matrix(Ext)
if (nrow(Ext) != nAges ||  ncol(Ext) != nYears) {
stop( "Mismatch between the dimension of Ext and the
number of years or ages")
}
rownames(Ext) <- ages
colnames(Ext) <- years
if (data$type == "central" && object$link == "logit")
warning( "logit-Binomial model fitted to central exposure data\n")
if (data$type == "initial" && object$link == "log")
warning( "log-Poisson model fitted to initial exposure data\n")

#Extract the specific ages and years for fitting
if (length(ages.fit) != length(which(ages.fit %in% ages))) {
stop( "ages.fit is not a subset of ages")
}
if (length(years.fit) != length(which(years.fit %in% years))) {
stop( "years.fit is not a subset of years")
}
Dxt <- Dxt[which(ages %in% ages.fit), which(years %in% years.fit)]
Ext <- Ext[which(ages %in% ages.fit), which(years %in% years.fit)]
ages <- ages.fit
years <- years.fit

# Construct fitting data
nAges <- length(ages)
nYears <- length(years)
cohorts <- (years - ages[nAges]):(years[nYears] - ages)
nCohorts <- length(cohorts)

if (is.null(oxt)) {
oxt <- matrix(0, nrow = nAges, ncol = nYears)
rownames(oxt) <- ages
colnames(oxt) <- years
} else {
oxt <- matrix(oxt, nrow = nAges, ncol = nYears)
rownames(oxt) <- ages
colnames(oxt) <- years
}

if (is.null(wxt)) {
wxt <- matrix(1, nrow = nAges, ncol = nYears)
} else {
wxt <- as.matrix(wxt)
if ( nrow(wxt) != nAges ||  ncol(wxt) != nYears) {
stop( "Mismatch between the dimension of the weigth matrix wxt and the
number of fitting years or ages")
}
}
rownames(wxt) <- ages
colnames(wxt) <- years
if (any(Ext <= 0)) { #Non-positive exposures
indExt <- (Ext <= 0)
wxt[indExt] <- 0
warning(paste("StMoMo: ", sum(indExt), " data points have
non-positive exposures and have been zero weighted\n",
sep = ""))
}
if (any(is.na(Dxt / Ext))) { #Missing values
indqxt <- is.na(Dxt / Ext)
wxt[indqxt] <- 0
warning(paste("StMoMo: ", sum(indqxt),
" missing values which have been zero weighted\n", sep = ""))
}

# Identify the data ages, years or cohorts with 0 weight
wxtDF <- (reshape2::melt(wxt, value.name = "w", varnames = c("x", "t")))
wxtDF <- transform(wxtDF, c = t - x)
wxTemp <- aggregate(data = wxtDF, w ~ x, FUN = sum)
zeroWeigthAges <- as.character(wxTemp$x[which((wxTemp$w <= 0))])
indX <- which((!(as.character(ages) %in% zeroWeigthAges)))
wtTemp <- aggregate(data = wxtDF, w ~ t, FUN = sum)
zeroWeigthYears <- as.character(wtTemp$t[which((wtTemp$w <= 0))])
indT <- which((!(as.character(years) %in% zeroWeigthYears)))
wcTemp <- aggregate(data = wxtDF, w ~ c, FUN = sum)
zeroWeigthCohorts <- as.character(wcTemp$c[which((wcTemp$w <= 0))])
indC <- which((!(as.character(cohorts) %in% zeroWeigthCohorts)))
if (verbose) {
if (length(zeroWeigthAges) > 0) {
cat("StMoMo: The following ages have been zero weigthed:",
zeroWeigthAges, "\n")
}
if (length(zeroWeigthYears) > 0) {
cat("StMoMo: The following years have been zero weigthed:",
zeroWeigthYears, "\n")
}
if (length(zeroWeigthCohorts) > 0) {
cat("StMoMo: The following cohorts have been zero weigthed:",
zeroWeigthCohorts, "\n")
}
}

### Set starting values
if (is.null(start.ax)) {
if (object$link == "log") { logmxt <- log(Dxt / Ext) start.ax <- rowSums((logmxt-oxt) * wxt) / rowSums(wxt) }else{ logoddsxt <- log(Dxt / (Ext - Dxt)) start.ax <- rowSums((logoddsxt - oxt) * wxt) / rowSums(wxt) } } else { if (length(start.ax) != nAges) { stop( "Mismatch between the number of ages and start.ax") } } names(start.ax) <- ages start.ax[zeroWeigthAges] <- NA if (is.null(start.bx)) { start.bx <- array(1 / nAges, c(nAges, 1), dimnames = list(ages, 1)) } else { if (nrow(start.bx) != nAges) { stop( "Mismatch between the number of ages and start.bx.") } if (ncol(start.bx) != 1) { stop( "Mismatch between the number of age/period terms and start.bx.") } } dimnames(start.bx) <- list(ages, 1) start.bx[zeroWeigthAges, ] <- NA if (is.null(start.kt)) { start.kt <- array(0, c(1, nYears), dimnames = list(1, years)) } else { if (ncol(start.kt) != nYears) { stop( "Mismatch between the number of years and start.kt.") } if (nrow(start.kt) != 1) { stop( "Mismatch between the number of age/period terms and start.kt.") } } dimnames(start.kt) <- list(1, years) start.kt[, zeroWeigthYears] <- NA if (is.null(start.b0x) || object$cohortAgeFun == "1" ) {
start.b0x <- rep(1, nAges)
} else {
if (length(start.b0x) != nAges && object$cohortAgeFun == "NP") { stop( "Mismatch between the number of ages and start.b0x.") } } names(start.b0x) <- ages if (object$cohortAgeFun != "1" )  start.b0x[zeroWeigthAges] <- NA

if (is.null(start.gc)) {
start.gc <- rep(0, nCohorts)
} else {
if (length(start.gc) != nCohorts) {
stop( "Mismatch between the number of cohorts and start.gc.")
}
}
names(start.gc) <- cohorts
start.gc[zeroWeigthCohorts] <- NA

out <- list(model = object, ax = start.ax, bx = start.bx,
kt = start.kt, b0x = start.b0x, gc = start.gc,
data = data, Dxt = Dxt, Ext = Ext, oxt = oxt ,
wxt = wxt, ages = ages, years = years,
cohorts = cohorts, fittingModel = NULL,
loglik = NA, deviance = NA,
npar = (2 + (object$cohortAgeFun == "NP")) * (nAges - length(zeroWeigthAges)) + (nYears - length(zeroWeigthYears)) + (nCohorts - length(zeroWeigthCohorts)) - (3 + object$approxConst + (object$cohortAgeFun == "NP")), nobs = sum(wxt), conv = NA, fail = FALSE, call = match.call()) class(out) <- c("fitrh", "fitStMoMo") ### Do maximisation iterations #auxiliary variables for itertions matbx <- array(NA, dim = c(nAges, nYears)) matbx2 <- array(NA, dim = c(nAges, nYears)) matkt <- array(NA, dim = c(nYears, nAges)) matkt2 <- array(NA, dim = c(nYears, nAges)) matb0x <- array(NA, dim = c(nAges, nYears)) matb0x2 <- array(NA, dim = c(nAges, nYears)) diagonal <- col(Ext) - row(Ext) sumDiags <- function(A) sapply(split(A, diagonal), FUN = sum, na.rm = TRUE) matgc <- array(NA, dim = c(nAges, nYears)) matgc2 <- array(NA, dim = c(nAges, nYears)) gcIndex <- as.numeric(gl(nYears, nAges)) + seq(nAges - 1, 0) xbar <- 1:nAges names(xbar) <- ages xbar[zeroWeigthAges] <- NA xbar <- xbar - mean(xbar, na.rm = TRUE) tbar <- 1:nYears names(tbar) <- years tbar[zeroWeigthYears] <- NA tbar <- tbar - mean(tbar, na.rm = TRUE) sumtbar2 <- sum(tbar ^ 2, na.rm = TRUE); cbar <- 1:nCohorts names(cbar) <- cohorts cbar[zeroWeigthCohorts] <- NA cbar <- cbar - mean(cbar, na.rm = TRUE) sumcbar2 <- sum(cbar ^ 2, na.rm = TRUE); #initialisation mxt <- Dxt / Ext mhatxt <- fitted.fitStMoMo(out, type = "rates") Dhatxt <- Ext * mhatxt out$loglik <- ifelse(object$link == "log", computeLogLikPoisson(obs = Dxt, fit = Dhatxt, weight = wxt), computeLogLikBinomial(obs = mxt, fit = mhatxt, exposure = Ext, weight = wxt)) absErrL <- 10 iter <- 0 if (verbose) { cat("StMoMo: Start fitting with iterative Newton-Raphson\n") cat("Running iterations") } while(absErrL > tolerance & iter < iterMax){ #update ax sumNum <- rowSums(wxt * (Dxt - Dhatxt), na.rm = TRUE) sumDen <- switch(object$link,
log = rowSums(wxt * Dhatxt, na.rm = TRUE),
logit = rowSums(wxt * Dhatxt * (1 - mhatxt),
na.rm = TRUE))
out$ax[indX] <- out$ax[indX] + sumNum[indX] / sumDen[indX]

#update kt
mhatxt <- fitted.fitStMoMo(out, type = "rates")
Dhatxt <- Ext * mhatxt
matbx[] <- out$bx[, 1] matbx2[] <- out$bx[, 1] ^ 2
sumNum <- colSums(wxt * (Dxt - Dhatxt) * matbx, na.rm = TRUE)
sumDen <- switch(object$link, log = colSums(wxt * Dhatxt * matbx2, na.rm = TRUE), logit = colSums(wxt * Dhatxt * matbx2 * (1 - mhatxt), na.rm = TRUE)) out$kt[1, indT] <- out$kt[1, indT] + sumNum[indT] / sumDen[indT] #update bx mhatxt <- fitted.fitStMoMo(out, type = "rates") Dhatxt <- Ext * mhatxt matkt[] <- out$kt[1, ]
matkt2[] <- out$kt[1, ] ^ 2 sumNum <- rowSums(wxt * (Dxt - Dhatxt) * t(matkt), na.rm = TRUE) sumDen <- switch(object$link,
log = rowSums(wxt * Dhatxt * t(matkt2), na.rm = TRUE),
logit = rowSums(wxt * Dhatxt * t(matkt2) * (1 - mhatxt),
na.rm = TRUE))
out$bx[indX, 1] <- out$bx[indX, 1] + sumNum[indX] / sumDen[indX]

#update gc
mhatxt <- fitted.fitStMoMo(out, type = "rates")
Dhatxt <- Ext * mhatxt
if (object$cohortAgeFun == "NP"){ matb0x[] <- out$b0x
matb0x2[] <- out$b0x ^ 2 sumNum <- sumDiags(wxt * (Dxt - Dhatxt) * matb0x) sumDen <- switch(object$link,
log = sumDiags(wxt * Dhatxt * matb0x2),
logit = sumDiags(wxt * Dhatxt * matb0x2 * (1 - mhatxt)))
} else {
sumNum <- sumDiags(wxt * (Dxt - Dhatxt))
sumDen <- switch(object$link, log = sumDiags(wxt * Dhatxt), logit = sumDiags(wxt * Dhatxt * (1 - mhatxt))) } out$gc[indC] <- out$gc[indC] + sumNum[indC] / sumDen[indC] #update b0x if (object$cohortAgeFun == "NP"){
mhatxt <- fitted.fitStMoMo(out, type = "rates")
Dhatxt <- Ext * mhatxt
matgc[] <- out$gc[gcIndex] matgc2[] <- (out$gc ^ 2)[gcIndex]
sumNum <- rowSums(wxt * (Dxt - Dhatxt) * matgc, na.rm = TRUE)
sumDen <- switch(object$link, log = rowSums(wxt * Dhatxt * matgc2, na.rm = TRUE), logit = colSums(wxt * Dhatxt * matgc2 * (1 - mhatxt), na.rm = TRUE)) out$b0x[indX] <- out$b0x[indX] + sumNum[indX] / sumDen[indX] mhatxt <- fitted.fitStMoMo(out, type = "rates") Dhatxt <- Ext * mhatxt } #apply constraints constPar<- object$constFun(out$ax, out$bx, out$kt, out$b0x, out$gc, wxt, ages) out$ax <- constPar$ax out$bx <- constPar$bx out$kt <- constPar$kt out$b0x <- constPar$b0x out$gc <- constPar$gc #Apply approx constraint if (object$approxConst) {
K <- sum(tbar[indT] * out$kt[1, indT]) / sumtbar2 e <- -sum(cbar[indC] * out$gc[indC]) / sumcbar2
out$ax[indX] <- out$ax[indX] + e * out$b0x[indX] * xbar[indX] out$bx[indX, 1] <- K / (K - e) * out$bx[indX, 1] - e / (K - e) * out$b0x[indX]
out$kt[, indT] <- (K - e) / K * out$kt[, indT]
out$gc[indC] <- out$gc[indC] + e * cbar[indC]
}

#update iteration
oldLoglik <- out$loglik mhatxt <- fitted.fitStMoMo(out, type = "rates") Dhatxt <- Ext * mhatxt out$loglik <- ifelse(object$link == "log", computeLogLikPoisson(obs = Dxt, fit = Dhatxt, weight = wxt), computeLogLikBinomial(obs = mxt, fit = mhatxt, exposure = Ext, weight = wxt)) absErrL <- abs(out$loglik - oldLoglik)
iter <- iter + 1
if (verbose ) cat(".")

}

if (verbose ){
cat('\nDone\n')
cat('StMoMo: Finish fitting with iterative Newton-Raphson\n')
}

### Prepare output
out$deviance <- ifelse(object$link == "log",
computeDeviancePoisson(obs = Dxt, fit = Dhatxt,
weight = wxt),
computeDevianceBinomial(obs = mxt, fit = mhatxt,
exposure = Ext, weight = wxt))
out$conv <- (iter < iterMax) if (!out$conv) warning( "Iterative Newton-Raphson procedure did not converge.
See fittingModel for information about the last iteration.\n")
if (is.na(out$loglik)) { warning( "Iterative Newton-Raphson procedure failed. See fittingModel for information about the last iteration.\n") out$fail <- TRUE
out$conv <- FALSE } out$fittingModel <- list(tolerance = tolerance, iterMax = iterMax, iter = iter, absErr = absErrL)

out

}


## 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.