Точность прогноза в R

Question

Я следовал инструкциям в этом пакетном документе StMoMo, чтобы сопоставить Ли Картер с данными о смертности в Канаде.

Следующим шагом в моем проекте является измерение точности прогноза модели Ли Картера, когда он соответствует этим канадским данным.

Чтобы сделать это, я попытался использовать точность (), но столкнулся с ошибкой, поскольку моя подгонка Ли Картера относится к классу "fitStMoMo", а не к классу "прогноз" или временному ряду.

Существует ли альтернативная функция точности прогноза, которую я могу использовать для объектов "fitStMoMo", которая будет вычислять среднюю ошибку, среднеквадратичную ошибку, среднюю абсолютную ошибку, среднюю процентную ошибку, среднюю абсолютную процентную ошибку и среднюю абсолютную масштабированную ошибкудля меня?

Представляет

Представляет, созданный с использованием EWMaleData, который используется в документе StMoMo для специальной отметки ошибки:

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

Я не совсем уверен в том, как работает accuracy() из forecast, но в некотором смысле он должен сравнивать реальные и прогнозируемые значения и возвращать метрики о том, насколько они различаются.В широком смысле это можно рассматривать как форму перекрестной проверки.Поскольку accuracy() не работает для StMoMo объектов, мы могли бы также самостоятельно разработать процедуру перекрестной проверки.
Для краткого учебного пособия по этой форме перекрестной проверки я бы порекомендовал Роба Хиндмананоты на tsCV() от forecast.Было бы неплохо, если бы мы могли использовать tsCV() здесь, но это работает только для одномерных временных рядов, а данные о смертности по сути являются многомерными временными рядами.
Я должен также упомянуть, что до сегодняшнего дня я никогда не слышал о моделировании смертности,так что это часть теории моделей, на которой я очень размыт.

Этот первый бит идентичен тому, что уже был опубликован

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)

Как видно из заметок Хиндмана, чтобы сделать это должным образом, мы должны провести это сравнение, используя прогнозы в нескольких точках нашего временного ряда и усредненные оценки.