1
$\begingroup$

I am fitting USA mortality data from Human Mortality Database (with data downloaded from here after registering for a free account https://www.mortality.org/cgi-bin/hmd/hmd_download.php) with the StMoMo package in R. The model fit has not been great. I have looked into the source code of the package, and did not really find any useful information that could help me improve the model accuracy. (I noticed that it is using generalized nonlinear model in fitting the model parameters).

Thus, I am asking for any advise that could help me improve accuracy.

My code is below:

library("StMoMo")

label="U.S.A."

Variable "mx" contains the data on CENTRAL DEATH RATE; the data is divided by gender, age, and year. Cohort effect can be extracted.

mx <- try(utils::read.table("C:/Users/myself/Downloads/USA/STATS/Mx_1x1.txt", skip = 2, header = TRUE, 
    na.strings = "."), TRUE)
mx = mx[mx$Year<=2014 & mx$Year>=1977,]

Variable "pop" contains population data (number of exposure at risk) divided by gender, age, and year. Please note that the path should point to the USA folder downloaded from HMD for USA mortality data. Please also note that the year has been restricted to the recent years of 1977-2014. Although HMD has data from 1933-2016.

pop <- try(utils::read.table("C:/Users/myself/Downloads/USA/STATS/Exposures_1x1.txt", skip = 2, header = TRUE, 
    na.strings = "."), TRUE)
pop = pop[pop$Year<=2014 & pop$Year>=1977,]
obj <- list(type = "mortality", label = label, lambda = 0)
obj$year <- sort(unique(mx[, 1]))
n <- length(obj$year)
m <- length(unique(mx[, 2]))
obj$age <- mx[1:m, 2]
mnames <- names(mx)[-c(1, 2)]
n.mort <- length(mnames)
obj$rate <- obj$pop <- list()
for (i in 1:n.mort) {
obj$rate[[i]] <- matrix(mx[, i + 2], nrow = m, ncol = n)
obj$rate[[i]][obj$rate[[i]] < 0] <- NA
obj$pop[[i]] <- matrix(pop[, i + 2], nrow = m, ncol = n)
obj$pop[[i]][obj$pop[[i]] < 0] <- NA
dimnames(obj$rate[[i]]) <- dimnames(obj$pop[[i]]) <- list(obj$age, 
    obj$year)
}
names(obj$pop) = names(obj$rate) <- tolower(mnames)
obj$age <- as.numeric(as.character(obj$age))
if (is.na(obj$age[m])) 
    obj$age[m] <- 2 * obj$age[m - 1] - obj$age[m - 2]
USATrain = (structure(obj, class = "demogdata"))

Here we read in full data to compare against projections.

label="U.S.A."

mx <- try(utils::read.table("C:/Users/myself/Downloads/USA/STATS/Mx_1x1.txt", skip = 2, header = TRUE, 
    na.strings = "."), TRUE)

pop <- try(utils::read.table("C:/Users/myself/Downloads/USA/STATS/Exposures_1x1.txt", skip = 2, header = TRUE, 
    na.strings = "."), TRUE)
obj <- list(type = "mortality", label = label, lambda = 0)
obj$year <- sort(unique(mx[, 1]))
n <- length(obj$year)
m <- length(unique(mx[, 2]))
obj$age <- mx[1:m, 2]
mnames <- names(mx)[-c(1, 2)]
n.mort <- length(mnames)
obj$rate <- obj$pop <- list()
for (i in 1:n.mort) {
obj$rate[[i]] <- matrix(mx[, i + 2], nrow = m, ncol = n)
obj$rate[[i]][obj$rate[[i]] < 0] <- NA
obj$pop[[i]] <- matrix(pop[, i + 2], nrow = m, ncol = n)
obj$pop[[i]][obj$pop[[i]] < 0] <- NA
dimnames(obj$rate[[i]]) <- dimnames(obj$pop[[i]]) <- list(obj$age, 
    obj$year)
}
names(obj$pop) = names(obj$rate) <- tolower(mnames)
obj$age <- as.numeric(as.character(obj$age))
if (is.na(obj$age[m])) 
    obj$age[m] <- 2 * obj$age[m - 1] - obj$age[m - 2]
FullData = (structure(obj, class = "demogdata"))

Setting up constraint functions. I am fitting age 30-95.

usaMale = StMoMoData(USATrain, series = "male")
usaMaleStMoMo <- central2initial(usaMale)
ages.fit <- 30:95
LC = lc(link = "logit")
APC = apc(link = "logit")
M7 = m7()
M6 = m6()
M8 = m8(xc = 89)
CBD = cbd()
RH <- rh(link = "logit", cohortAgeFun = "1")

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)
 }



PLAT <- StMoMo(link = "logit", staticAgeFun = TRUE, periodAgeFun = c("1", f2), cohortAgeFun = "1", constFun = constPlat)

wxt <- genWeightMat(ages = ages.fit, years = usaMaleStMoMo$years, clip = 3)

Fitting the models.

LCfit <- fit(LC, data = usaMaleStMoMo, ages.fit = ages.fit, wxt = wxt)
APCfit <- fit(APC, data = usaMaleStMoMo, ages.fit = ages.fit, wxt = wxt)
CBDfit <- fit(CBD, data = usaMaleStMoMo, ages.fit = ages.fit, wxt = wxt)
M7fit <- fit(M7, data = usaMaleStMoMo, ages.fit = ages.fit, wxt = wxt)
M6fit <- fit(M6, data = usaMaleStMoMo, ages.fit = ages.fit, wxt = wxt)
M8fit <- fit(M8, data = usaMaleStMoMo, ages.fit = ages.fit, wxt = wxt)
PLATfit <- fit(PLAT, data = usaMaleStMoMo, ages.fit = ages.fit, wxt = wxt)
RHfit <- fit(RH, data = usaMaleStMoMo, ages.fit = ages.fit, wxt = wxt, start.ax = LCfit$ax, start.bx = LCfit$bx, start.kt = LCfit$kt)

Forecast for 50 years

LCfor <- forecast(LCfit, h = 50)
CBDfor <- forecast(CBDfit, h = 50)
APCfor <- forecast(APCfit, h = 50, gc.order = c(1, 1, 0))
RHfor <- forecast(RHfit, h = 50, gc.order = c(1, 1, 0))
M7for <- forecast(M7fit, h = 50, gc.order = c(2, 0, 0))
M6for <- forecast(M6fit, h = 50, gc.order = c(2, 0, 0))
M8for <- forecast(M8fit, h = 50, gc.order = c(2, 0, 0))
PLATfor <- forecast(PLATfit, h = 50, gc.order = c(2, 0, 0))

Grabbing the actual mortality rate for 2015 USA male and the forecasted rates from all models. I am comparing the 2015 projected rate against the actual 2015 rates.

actual2015 = FullData$rate$male[ages.fit+1,c("2015")]
LC2015 = LCfor$rates[,1]
CBD2015 = CBDfor$rates[,1]
APC2015 = APCfor$rates[,1]
RH2015 =  RHfor$rates[,1]
M72015 = M7for$rates[,1]
M62015 = M6for$rates[,1]
M82015 = M8for$rates[,1]
PLAT2015 = PLATfor$rates[,1]

Consolidate results. I admit that the code below is quite ugly, and mean squared error is probably better than average percent error. However, I think the fact that the percent error is off by the amount shown in the picture below indicates that the model is not fitting very well. Moreover, the residual plot also shows that there are unexplained structure in the data.

aicbic = matrix(ncol = 8, nrow = 4)
colnames(aicbic) = c("LC", "CBD", "APC", "RH", "M7", "PLAT", "M6", "M8")
rownames(aicbic) = c("Number Parameters", "AIC", "BIC", "Average Pct Error for t+1")

aicbic[1,1]=LCfit$npar
aicbic[1,2]=CBDfit$npar
aicbic[1,3]=APCfit$npar
aicbic[1,4]=RHfit$npar
aicbic[1,5]=M7fit$npar
aicbic[1,6]=PLATfit$npar
aicbic[1,7]=M6fit$npar
aicbic[1,8]=M8fit$npar
aicbic[2,1]=AIC(LCfit)
aicbic[2,2]=AIC(CBDfit)
aicbic[2,3]=AIC(APCfit)
aicbic[2,4]=AIC(RHfit)
aicbic[2,5]=AIC(M7fit)
aicbic[2,6]=AIC(PLATfit)
aicbic[2,7]=AIC(M6fit)
aicbic[2,8]=AIC(M8fit)
aicbic[3,1]=BIC(LCfit)
aicbic[3,2]=BIC(CBDfit)
aicbic[3,3]=BIC(APCfit)
aicbic[3,4]=BIC(RHfit)
aicbic[3,5]=BIC(M7fit)
aicbic[3,6]=BIC(PLATfit)
aicbic[3,7]=BIC(M6fit)
aicbic[3,8]=BIC(M8fit)
aicbic[4,1]=mean((LC2015 - actual2015)/actual2015)*100
aicbic[4,2]=mean((CBD2015 - actual2015)/actual2015)*100
aicbic[4,3]=mean((APC2015 - actual2015)/actual2015)*100
aicbic[4,4]=mean((RH2015 - actual2015)/actual2015)*100
aicbic[4,5]=mean((M72015 - actual2015)/actual2015)*100
aicbic[4,6]=mean((PLAT2015 - actual2015)/actual2015)*100
aicbic[4,7]=mean((M62015 - actual2015)/actual2015)*100
aicbic[4,8]=mean((M82015 - actual2015)/actual2015)*100

View(aicbic)

Here is the result

Residual plot.

For the best model in terms of AIC and BIC, the Plat model, exhibits the following residuals.

Code for residual plot. CBDRes = residuals(CBDfit) plot(CBDRes, type = "colourmap", main = "CBD Residuals")

LCRes = residuals(LCfit)
plot(LCRes, type = "colourmap", main = "Lee-Carter Residuals")

APCRes = residuals(APCfit)
plot(APCRes, type = "colourmap", main = "Age-Period-Cohort Residuals")

RHRes = residuals(RHfit)
plot(RHRes, type = "colourmap", main = "Renshaw and Haberman Residuals")

M7Res = residuals(M7fit)
plot(M7Res, type = "colourmap", main = "M7 Residuals")

PLATRes = residuals(PLATfit)
plot(PLATRes, type = "colourmap", main = "Plat Residuals")

M6Res = residuals(M6fit)
plot(M6Res, type = "colourmap", main = "M6 Residuals")

M8Res = residuals(M8fit)
plot(M8Res, type = "colourmap", main = "M8 Residuals")
$\endgroup$
  • $\begingroup$ Way too long, you should aim for short and to the point. $\endgroup$ – user2974951 Dec 13 '18 at 7:25
  • $\begingroup$ @user2974951 I think a lot of it is just reading in the data and setting up constraint functions. Code chunk 3 and 4 are long but they're necessary to make the object. I included many models so it's more complete. Should I include only one model? $\endgroup$ – user101998 Dec 13 '18 at 16:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.