## load package: MortalitySmooth ## construction of B-splines basis ## Section 2, page 4 library("MortalitySmooth") x <- runif(100) B3 <- MortSmooth_bbase(x=x, xl=min(x), xr=max(x), ndx=10, deg=3) ## construction of matrices of differences ## Section 2.1, page 4 D1 <- diff(diag(5), diff=1) D2 <- diff(diag(5), diff=2) ## comparison direct and GLAM computation ## of the inner product in Equation 8 ## Section 2.3, page 6-7 a <- 1:50 m <- length(a) y <- 1:100 n <- length(y) Ba <- MortSmooth_bbase(x=a,xl=min(a),xr=max(a),ndx=7,deg=3) ka <- ncol(Ba) By <- MortSmooth_bbase(x=y,xl=min(y),xr=max(y),ndx=17,deg=3) ky <- ncol(By) w <- runif(m*n) B <- kronecker(By, Ba) system.time(tBWB1 <- t(B) %*% (w*B)) Ba1 <- kronecker(matrix(1, ncol = ka, nrow = 1), Ba) Ba2 <- kronecker(Ba, matrix(1, ncol = ka, nrow = 1)) Ga <- Ba1 * Ba2 By1 <- kronecker(matrix(1, ncol = ky, nrow = 1), By) By2 <- kronecker(By, matrix(1, ncol = ky, nrow = 1)) Gy <- By1 * By2 Wbreve <- matrix(w, m, n) system.time(tBWB2 <- MortSmooth_BWB(Ga, Gy, ka, ky, Wbreve)) size1 <- object.size(list(B, w)) size2 <- object.size(list(Ga, Gy, ka, ky, Wbreve)) size1/size2 ## load mortality data, ## extract data with selectHMDdata: ## example of Danish death rates ## Section 3, page 10-11 x <- 50:100 y <- 1950:2009 mydata <- selectHMDdata(country="Denmark", data="Rates", sex="Females", ages=x, years=y) ## print the selected data ## Section 3, page 11 mydata ## !! note that in the article I shorthen the outcome ## plot selected data ## Section 3, page 11 (Figure 1, page 11) plot(mydata, cut=10, col.regions=terrain.colors(11)) ## select data, ## smooth one-dimensional mortality data (default arguments) ## and print fitted object ## Section 4.1, page 12 x <- 1945:2000 y <- selectHMDdata("Denmark", "Deaths", "Males", ages = 60, years = x) e <- selectHMDdata("Denmark", "Exposures", "Males", ages = 60, years = x) fit1D <- Mort1Dsmooth(x=x, y=y, offset=log(e)) fit1D ## extract and print deviance from the fitted object ## Section 4.1, page 12 fit1D$deviance ## print the summary of the fitted object ## Section 4.1, page 12-13 summary(fit1D) ## smooth the previous data ## using different model selection: ## AIC, subjective choice and given ED ## Section 4.1, page 13 fit1Daic <- Mort1Dsmooth(x=x, y=y, offset=log(e), method=2) fit1Dsub <- Mort1Dsmooth(x=x, y=y, offset=log(e), method=3, lambda=10000) fit1Dfix <- Mort1Dsmooth(x=x, y=y, offset=log(e), method=4, df=8) ## plot fitted object (with default arguments) ## Section 4.1, page 13 (Figure 2, page 14) plot(fit1D) ## add to the previous plot ## the others fitted objects ## Section 4.1, page 13 (Figure 2, page 14) lines(x, fit1Daic$logmortality, col=3, lwd=2) lines(x, fit1Dsub$logmortality, col=4, lwd=2) lines(x, fit1Dfix$logmortality, col=5, lwd=2) legend("bottom", c("Actual", "BIC", "AIC", "Subjective", "Fix df=8"), col=1:5, pch=c(1,-1,-1,-1,-1), lwd=c(1,2,2,2,2), lty=c(0,1,1,1,1), bty="n") ## the analogous models fitted with mgcv ## Section 4.1, page 14 library("mgcv") kk <- fit1D$ndx+fit1D$deg mgam <- gam(y ~ offset(log(e)) + s(x, bs="ps", k=kk), family="poisson") mgam2 <- gam(y ~ offset(log(e)) + s(x, bs="ps", k=kk), family="poisson", method="REML") ## outcomes from mgcv ## Section 4.1, page 14 ## !! these plotting commands are not shown in the paper ## !! but the outcomes are given in Figure 2, page 14 lines(x, log(fitted(mgam)/e), lwd=2, lty=2, col="green4") lines(x, log(fitted(mgam2)/e), lty=2, lwd=2, col="red4") ## extrapolating mortality data ## Section 4.1, page 14 x.new <- x[1]:2020 fit1Dfor <- predict(fit1D, newdata=x.new, se.fit=TRUE) ## smooth and extrapolate using different ## order of difference for the penalty term ## Section 4.1, page 15 fit1D3d <- Mort1Dsmooth(x=x, y=y, offset=log(e), pord=3) fit1D1d <- Mort1Dsmooth(x=x, y=y, offset=log(e), pord=1) fit1Dfor3d <- predict(fit1D3d, newdata=x.new) fit1Dfor1d <- predict(fit1D1d, newdata=x.new) ## plotting actual, fitted and forecast ## and associate 95% confidence interval ## for the death rates on a log-scale ## Section 4.1, page 16 (Figure 3, page 15) plot(x, log(y/e), xlim=range(x.new), ylim=range(fit1Dfor$fit, log(y/e))) lines(x.new, fit1Dfor$fit, lwd=2, col=2) lines(x.new, fit1Dfor$fit + 2*fit1Dfor$se.fit, lty=2, col=2) lines(x.new, fit1Dfor$fit - 2*fit1Dfor$se.fit, lty=2, col=2) lines(x.new, fit1Dfor3d, lwd=2, col=3) lines(x.new, fit1Dfor1d, lwd=2, col=4) ## in Figure 3 at page 15, we also plot the actual ## death rates for 2001:2009 without showing ## the following lines in the paper y2 <- selectHMDdata("Denmark", "Deaths", "Males", ages = 60, years = 2001:2009) e2 <- selectHMDdata("Denmark", "Exposures", "Males", ages = 60, years = 2001:2009) points(2001:2009, log(y2/e2), pch=16, col=1) ## also the R-code for the legend in skiped in the paper legend("bottom", c("Actual", "d=2", "d=3", "d=1"), col=1:4, pch=c(1,-1,-1,-1), lwd=c(1,2,2,2), lty=c(0,1,1,1)) ## smooth mortality data accommodating for ## the presence of overdispersion in ## the smoothing parameter selection ## Section 4.1, page 16 fit1Dover <- Mort1Dsmooth(x=x, y=y, offset=log(e), overdispersion=TRUE) ## compare optimal smoothing parameters and ## effective dimension ## with and without accommodation for overdispersion ## Section 4.1, page 16 c(fit1D$lambda, fit1D$df) c(fit1Dover$lambda, fit1Dover$df) ## select data, ## smooth two-dimensional mortality data (default arguments) ## and print fitted object ## Section 4.2, page 17 x <- 10:100 y <- 1930:2010 Y <- selectHMDdata("Sweden", "Deaths", "Females", ages = x, years = y) E <- selectHMDdata("Sweden", "Exposures", "Females", ages = x, years = y) fit2D <- Mort2Dsmooth(x=x, y=y, Z=Y, offset=log(E)) fit2D ## plot fitted object ## Section 4.2, page 17 (Figure 4, page 18) plot(fit2D, palette="terrain.colors") ## extract (pearson) residuals from ## the previous fitted object ## Section 4.2, page 17 res2D <- residuals(fit2D, type="pearson") ## two-dimensional plot ## of the person residuals ## Section 4.2, page 18 (Figure 5, page 19) res.grid <- expand.grid(list(x=x,y=y)) res.grid$res <- c(res2D) levelplot(res ~ y*x, data=res.grid, at=c(min(res2D), -2, -1, 1, 2, max(res2D))) ## interpolating mortality data ## Section 4.2, page 18-19 m <- length(x) n <- length(y) set.seed(1) whi0 <- sample(x=1:(m*n), size=5900) W <- matrix(1,m,n) W[whi0] <- 0 fit2Dint <- Mort2Dsmooth(x=x, y=y, Z=Y, offset=log(E), W=W) ## plotting the interpolated mortality data ## over ages for selected years ## Section 4.2, page 19 (Figure 6, page 20) lmx.hat <- fit2Dint$logmortality[,c(1,27,54,81)] lmx.act0 <- log(Y/E) lmx.act0[whi0] <- NA lmx.act <- lmx.act0[,c(1,27,54,81)] matplot(x, lmx.hat, t="l", lty=1) matpoints(x, lmx.act, pch=1) legend("topleft", legend=y[c(1,27,54,81)], col=1:4, pch=1, lty=1)