R中的预测准确性

Question

我已按照 thisStMoMo 软件包文档中的说明将 Lee Carter 与加拿大的死亡率数据相匹配。

我项目的下一步是测量 Lee Carter 模型在拟合加拿大数据时的预测准确性。

为此，我尝试使用 accuracy() 但遇到了错误，因为我的 Lee Carter 拟合属于“fitStMoMo”类而不是“预测”类或时间序列。

我是否可以在“fitStMoMo”对象上使用替代的预测精度函数来计算平均误差、均方根误差、平均绝对误差、平均百分比误差、平均绝对百分比误差和平均绝对比例误差？

代表

使用 StMoMo 文档中使用的 EWMaleData 创建的 Reprex 专门标记错误：

library("StMoMo")
library("demography")
library("forecast")

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 = "logit", staticAgeFun = TRUE, periodAgeFun = "NP",
constFun = constLC)
LC <- lc(link = "logit")
LC$gnmFormula
#> [1] "D/E ~ -1 + offset(o) + factor(x) + Mult(factor(x), factor(t), inst = 1)"

EWMaleData
#> Mortality data for England and Wales
#>     Series:  male
#>     Years: 1961 - 2011
#>     Ages:  0 - 100
#>     Exposure:  central

EWMaleIniData <- central2initial(EWMaleData)
ages.fit <- 55:89
wxt <- genWeightMat(ages = ages.fit, years = EWMaleIniData$years,
clip = 3)
LCfit <- fit(LC, data = EWMaleIniData, ages.fit = ages.fit, wxt = wxt)
#> StMoMo: The following cohorts have been zero weigthed: 
#>   1872 1873 1874 1954 1955 1956 
#> StMoMo: Start fitting with gnm
#> Initialising
#> Running start-up iterations..
#> Running main iterations.....
#> Done
#> StMoMo: Finish fitting with gnm

LCfor <- forecast(LCfit, h = 50)
class(LCfit)
#> [1] "fitStMoMo"
class(LCfor)
#> [1] "forStMoMo"
accuracy(LCfit)
#> Error in accuracy.default(LCfit): First argument should be a forecast object 
#>   or a time series.
accuracy(LCfor)
#> Error in accuracy.default(LCfor): First argument should be a forecast object
#>   or a time series.

Answer 1

我不完全确定来自forecast 的accuracy() 是如何工作的，但在某种程度上它必须比较真实值和预测值，并返回关于它们差异多少的指标。从广义上讲，这可以被视为一种交叉验证的形式。由于accuracy() 不适用于StMoMo 对象，我们不妨自己开发一个交叉验证例程。
对于这种交叉验证形式的简短入门，我建议 Rob Hyndman's notes 上的 tsCV() 来自 forecast。如果我们可以在这里使用tsCV() 那就太好了，但它只适用于单变量时间序列，而死亡率数据本质上是多变量时间序列。
我还应该提一下，在今天之前我从未听说过 Mortality Modeling，所以我对此的模型理论部分非常模糊。

第一个部分与已经发布的内容相同

library(StMoMo)
library(demography)
library(forecast)

data(EWMaleData)

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="logit", staticAgeFun=TRUE, periodAgeFun="NP", constFun=constLC)
LC <- lc(link="logit")

EWMaleIniData <- central2initial(EWMaleData)

然后事情变得有点不同。在时间序列上执行 CV 的中心点是对我们实际拥有的数据进行预测，但我们假装没有。因此，我们必须对我们的数据进行子集化，以便我们想要预测的块不是模型的一部分。在这个具体示例中，我们将使用前 30 年，然后预测接下来的 10 年

ages.fit <- 55:89
years.fit <- EWMaleIniData$years[1]:(EWMaleIniData$years[1] + 30)
years.for <- 10

wxt <- genWeightMat(ages=ages.fit, years=years.fit, clip=3)
LCfit <- fit(LC, data=EWMaleIniData, ages.fit=ages.fit,
  years.fit=years.fit, wxt=wxt)
LCfor <- forecast(LCfit, h=years.for)

现在我们有了一个十年的预测，我们可以将这些年与我们的实际数据进行比较，并使用我们想要的任何错误度量来查看预测的准确度。

cvy <- LCfor$years  # years used in forecast
cva <- LCfor$ages   # ages used in forecast

pred <- LCfor$rates # predicted mortality rates

# actual mortality rates subset to the same ages and years as forecast
actual <- EWMaleIniData$Dxt/EWMaleIniData$Ext
actual <- actual[rownames(actual) %in% cva,
                 colnames(actual) %in% cvy]

# A collection of error measures. plenty of others can be devised
err <- pred - actual
Q <- pred/actual
rmse <- sqrt(rowMeans(err^2))
mae <- rowMeans(abs(err))
smape <- 100 * (rowMeans(exp(abs(log(Q)))) - 1)

这个位纯粹是为了显示结果

par(mfrow=c(3, 1), mar=c(3, 3, 1, 1), mgp=c(2, 0.8, 0), oma=c(0, 0, 2, 0))
plot(as.numeric(names(rmse)), rmse, type="h", xlab="")
plot(as.numeric(names(mae)), mae, type="h", xlab="")
plot(as.numeric(names(smape)), smape, type="h", xlab="Ages")
mtext(paste("Forecast accuracy for the years", 
  paste(cvy[c(1, years.for)], collapse=" - ")), 
  3, outer=TRUE)

正如在 Hyndman 的笔记中所看到的，要正确地做到这一点，我们必须使用时间序列中几个点的预测以及平均得分来进行比较。