I'm working on fitting a Poisson Lee-Carter model to the Danish male population data for calendar years 1990–2010 and ages 20–90 in R. After fitting the model, I want to plot the predicted mortality rates $\hat{\mu}_t^{L-C,(1)}(x)$ for the years 1990–2020 for ages 20, 40, 60, 75, and 90. However, I am struggling to plot the predictions for the calendar years 1990–2010, which were used for fitting the model.
The code below plots the predicted mortality rates for 2011–2020, but I need help on how to plot the predicted rates for the years 1990–2010 as well, since those years were used to train the model. How can I extend the plot to include the predicted values for the years 1990–2010, while keeping the raw mortality rates for comparison?
Here is the R code I have so far:
# Load necessary libraries
library("devtools")
library("gnm")
library("forecast")
library("StMoMo")
library("cobs")
library("ftsa")
library("demography")
# Read in the Swedish mortality data
MortalityDataSWE <- read.demogdata("SWE_Mx_1x1.txt", "SWE_Exposures_1x1.txt", type = "mortality", label = "Sweden")
# Prepare the data using the StMoMo package for female mortality series
StMoMoSWE <- StMoMoData(MortalityDataSWE, series = "female")
# Define the LC model with a log-link function
mortalityFunc <- lc(link = "log")
# Define the age and year ranges for fitting the model
fitAges <- 20:90
fitYrs <- 1990:2010
totYrs <- 1990:2020
# Calculate observed death rates
obsDeathRateSWE <- StMoMoSWE$Dxt[as.character(fitAges), as.character(totYrs)] / StMoMoSWE$Ext[as.character(fitAges), as.character(totYrs)]
# Fit the LC model
LCfitSWE <- fit(mortalityFunc, data = StMoMoSWE, ages.fit = fitAges, years.fit = fitYrs)
# Simulate the model to get forecasts for 2011-2020
forecastLC <- simulate(LCfitSWE, years = 2011:2020)
# Get residuals from the fitted model
residualsLC <- residuals(LCfitSWE)
# Calculate mean and percentiles (5% and 95%) for the simulated mortality rates
meanMortality <- apply(forecastLC$rates, c(1, 2), mean)
lowerBound <- apply(forecastLC$rates, c(1, 2), quantile, probs = 0.05)
upperBound <- apply(forecastLC$rates, c(1, 2), quantile, probs = 0.95)
# Define the ages to analyze
selectedAges <- c(20, 40, 60, 75, 90)
years <- 2011:2020
# Extract the predicted mortality rates and percentiles for the selected ages
meanMortalitySelected <- meanMortality[as.character(selectedAges), ]
lowerBoundSelected <- lowerBound[as.character(selectedAges), ]
upperBoundSelected <- upperBound[as.character(selectedAges), ]
# Subset the predicted data to match the years range (2011-2020)
meanMortalitySelected <- meanMortalitySelected[, as.character(years)]
lowerBoundSelected <- lowerBoundSelected[, as.character(years)]
upperBoundSelected <- upperBoundSelected[, as.character(years)]
# Get the raw observed death rates (if available)
rawMortalityRates <- obsDeathRateSWE[as.character(selectedAges), as.character(years)]
# Plot for Age 20
plot(years, meanMortalitySelected["20", ], type = "l",
ylim = range(c(meanMortalitySelected["20", ], lowerBoundSelected["20", ], upperBoundSelected["20", ], rawMortalityRates["20", ]), na.rm = TRUE),
col = "blue", lwd = 2, main = "Age: 20")
lines(years, lowerBoundSelected["20", ], col = "red", lty = 2)
lines(years, upperBoundSelected["20", ], col = "red", lty = 2)
if (!is.null(rawMortalityRates)) {
points(years, rawMortalityRates["20", ], col = "black", pch = 16)
}
legend("topright", legend = c("Mean Prediction", "5% and 95% Percentiles", "Raw Rates"),
col = c("blue", "red", "black"), lty = c(1, 2, NA), pch = c(NA, NA, 16), bty = "n")
Then I have one plot for each age (here I only show the one for age 20). The datasets are from mortality and are called SWE_Mx_1x1 and SWE_Exposures_1x1.
I would appreciate any suggestions or solutions for plotting the predicted mortality rates for the years 1990–2010, so I can compare them with the raw mortality rates for that period. Thanks!
Note: I have tried to simulate for years = 1990:2020 (forecastLC <- simulate(LCfitSWE, years = 2011:2020)) but i still only get the first forecast at year 2011, which i am guessing is since we are fitting the model with the years 1990-2010.