StMoMo 死亡率包中的完整 PLAT 模型
full PLAT model in StMoMo mortality package
我正在使用 R 包 StMoMo 进行随机死亡率建模。描述符号的论文可以在这里找到:https://openaccess.city.ac.uk/id/eprint/17378/7/StMoMoVignette.pdf
论文描述了完整的 PLAT 模型:
以及简化的 PLAT 模型:
然后它提供(见p13-14)简化的PLAT模型的代码。此代码工作正常。
#to get data
ages.fit = 12:84
years.fit = 2008:2017
gender = "male"
JPNdata = hmd.mx(country="JPN",username=username,password=password,label="Japan")
JPNStMoMo = StMoMoData(JPNdata, series = gender,type="initial")
#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)
reducedPlat %>% fit(data=JPNStMoMo,ages.fit = ages.fit,years.fit=years.fit)
但是,当我尝试稍微修改代码以获得完整的 Plat 模型时出现以下错误:
The parameter transformation function does not preserve the fitted rates.
Check the 'constFun' argument of StMoMo."
修改后的代码如下:
f2 <- function(x, ages) mean(ages) - x
f3 <- function(x, ages) max(f2(x,ages),0) #added
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) #modified
kt[1, ] <- kt[1, ] - ci[1]
kt[2, ] <- kt[2, ] - ci[2]
kt[3, ] <- kt[3, ] - ci[3] #added
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) #modified
fullPlat %>% fit(data=JPNStMoMo,ages.fit = ages.fit,years.fit=years.fit)
虽然我的改动真的很小,但我并没有发现自己的错误。如果有人发现了什么,请提前致谢!
在我的代码中,必须将 max 更改为 pmax。
此外,包的作者提供了完整模型的代码:
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)
我正在使用 R 包 StMoMo 进行随机死亡率建模。描述符号的论文可以在这里找到:https://openaccess.city.ac.uk/id/eprint/17378/7/StMoMoVignette.pdf
论文描述了完整的 PLAT 模型:
以及简化的 PLAT 模型:
然后它提供(见p13-14)简化的PLAT模型的代码。此代码工作正常。
#to get data
ages.fit = 12:84
years.fit = 2008:2017
gender = "male"
JPNdata = hmd.mx(country="JPN",username=username,password=password,label="Japan")
JPNStMoMo = StMoMoData(JPNdata, series = gender,type="initial")
#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)
reducedPlat %>% fit(data=JPNStMoMo,ages.fit = ages.fit,years.fit=years.fit)
但是,当我尝试稍微修改代码以获得完整的 Plat 模型时出现以下错误:
The parameter transformation function does not preserve the fitted rates.
Check the 'constFun' argument of StMoMo."
修改后的代码如下:
f2 <- function(x, ages) mean(ages) - x
f3 <- function(x, ages) max(f2(x,ages),0) #added
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) #modified
kt[1, ] <- kt[1, ] - ci[1]
kt[2, ] <- kt[2, ] - ci[2]
kt[3, ] <- kt[3, ] - ci[3] #added
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) #modified
fullPlat %>% fit(data=JPNStMoMo,ages.fit = ages.fit,years.fit=years.fit)
虽然我的改动真的很小,但我并没有发现自己的错误。如果有人发现了什么,请提前致谢!
在我的代码中,必须将 max 更改为 pmax。 此外,包的作者提供了完整模型的代码:
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)