[help wanted] full PLAT model with StMoMo

[ACL90]

First, thank you for this incredible StMoMo package!

When I apply the code in your companion paper for the reduced Plat model, everything is fine.

the reduced Plat model

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)
   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
    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)
}
 reducedPlat <- StMoMo(link = "logit", staticAgeFun = TRUE,
                       periodAgeFun = c("1", f2), cohortAgeFun = "1", constFun = constPlat)

However, I get the following error when I try to slightly modify the code to get the full Plat model, as defined in your companion paper (https://openaccess.city.ac.uk/id/eprint/17378/7/StMoMoVignette.pdf, p10): "Error in fit.StMoMo(., data = JPNStMoMo, ages.fit = ages.fit, years.fit = years.fit) : The parameter transformation function does not preserve the fitted rates. Check the 'constFun' argument of StMoMo."

Here is the modified code: f2 <- function(x, ages) mean(ages) - x f3 <- function(x, ages) max(f2(x,ages),0) 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) 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 ci <- rowMeans(kt, na.rm = TRUE) ax <- ax + ci[1] + ci[2] (xbar - x) + ci[3] max(xbar - x,0) kt[1, ] <- kt[1, ] - ci[1] kt[2, ] <- kt[2, ] - ci[2] kt[3, ] <- kt[3, ] - ci[3] list(ax = ax, bx = bx, kt = kt, b0x = b0x, gc = gc) } fullPlat <- StMoMo(link = "logit", staticAgeFun = TRUE, periodAgeFun = c("1", f2, f3), cohortAgeFun = "1", constFun = constPlat)

Although my changes are really small, I do not spot my mistake. Any tip will be welcome! Thank you in advance.

[ACL90]

duplicated here: https://stackoverflow.com/questions/60431889/full-plat-model-in-stmomo-mortality-package

[amvillegas]

Hello, try replacing max to pmax so that you get the argument by argument maximum.

Below is some code for implementing the Plat 2009 model. I hope it helps

library(StMoMo)
f2 <- function(x, ages) mean(ages) - x
f3 <- function(x, ages) pmax(mean(ages)-x,0)
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)
  #\sum g(c)=0, \sum cg(c)=0, \sum 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
  #\sum kt[i, ] = 0
  ci <- rowMeans(kt, na.rm = TRUE)
  ax <- ax + ci[1] + ci[2] * (xbar - x) + ci[3] * pmax(xbar - x, 0)
  kt[1, ] <- kt[1, ] - ci[1]
  kt[2, ] <- kt[2, ] - ci[2]
  kt[3, ] <- kt[3, ] - ci[3]
  list(ax = ax, bx = bx, kt = kt, b0x = b0x, gc = gc)
}
PLAT <- StMoMo(link = "log", staticAgeFun = TRUE,
               periodAgeFun = c("1", f2, f3), cohortAgeFun = "1",
               constFun = constPlat)

ages.fit <- 0:100
wxt <- genWeightMat(ages = ages.fit, years = EWMaleData$years, clip = 3)

PLATfit <- fit(PLAT, data = EWMaleData, ages.fit = ages.fit, wxt = wxt)
plot(PLATfit, parametricbx = FALSE)