R/mortality_weights.R
mortality_weights.RdMortality Model averaging via stacking of predictive distributions and Pseudo-BMA weighting. Based on Yao et al. (2018) but adapted by Barigou et al. (2021) for mortality forecasting.
mortality_weights(X)
| X | A list of stanfit objects. |
|---|
A matrix containing one weight for each model and each approach.
Mortality model averaging via stacking of predictive distributions or pseudo-BMA weighting. Both approaches were proposed in Yao et al. (2018) based leave-one-out cross-validation which is not suited for forecasting. Barigou et al. (2021) adapted both approaches based on leave-future-out validation which is more appropriate for mortality forecasting.
The stacking method combines all models by maximizing the leave-future-out predictive density of the combination distribution. That is, it finds the optimal linear combining weights for maximizing the leave-future-out log score.
The pseudo-BMA method finds the relative weights proportional to the expected log predictive density of each model.
Similar to Yao et al. (2018), we recommend stacking for averaging predictive distributions as pseudo-BMA tends to select only one model.
Yao, Y., Vehtari, A., Simpson, D., & Gelman, A. (2018). Using stacking to average Bayesian predictive distributions (with discussion). Bayesian Analysis, 13(3), 917-1007.
Barigou K., Goffard P-O., Loisel S., Salhi Y. (2021). Bayesian model averaging for mortality forecasting using leave-future-out validation. arXiv preprint arXiv:2103.15434.
# \donttest{ #10-year forecasts for French data for ages 50-90 and years 1970-2017 with a log-Poisson model #where the 10 last years are held out for validation. We search for the model weights between #the Lee-Carter model and the RH model (Lee-Carter with cohort effect). ages.fit<-50:90 years.fit<-1970:2017 deathFR<-FRMaleData$Dxt[formatC(ages.fit),formatC(years.fit)] exposureFR<-FRMaleData$Ext[formatC(ages.fit),formatC(years.fit)] iterations<-1000 # Toy example, consider at least 2000 iterations fitLC=lc_stan(death = deathFR,exposure=exposureFR, forecast = 10, validation=10, family = "poisson",iter=iterations,chains=1) #> #> SAMPLING FOR MODEL 'leecarter' NOW (CHAIN 1). #> Chain 1: #> Chain 1: Gradient evaluation took 0.001 seconds #> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10 seconds. #> Chain 1: Adjust your expectations accordingly! #> Chain 1: #> Chain 1: #> Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup) #> Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup) #> Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup) #> Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup) #> Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup) #> Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup) #> Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling) #> Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling) #> Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling) #> Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling) #> Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling) #> Chain 1: Iteration: 1000 / 1000 [100%] (Sampling) #> Chain 1: #> Chain 1: Elapsed Time: 11.236 seconds (Warm-up) #> Chain 1: 5.106 seconds (Sampling) #> Chain 1: 16.342 seconds (Total) #> Chain 1: fitRH=rh_stan(death = deathFR,exposure=exposureFR, forecast = 10, validation=10, family = "poisson",iter=iterations,chains=1) #> #> SAMPLING FOR MODEL 'RHmodel' NOW (CHAIN 1). #> Chain 1: #> Chain 1: Gradient evaluation took 0 seconds #> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds. #> Chain 1: Adjust your expectations accordingly! #> Chain 1: #> Chain 1: #> Chain 1: Iteration: 1 / 1000 [ 0%] (Warmup) #> Chain 1: Iteration: 100 / 1000 [ 10%] (Warmup) #> Chain 1: Iteration: 200 / 1000 [ 20%] (Warmup) #> Chain 1: Iteration: 300 / 1000 [ 30%] (Warmup) #> Chain 1: Iteration: 400 / 1000 [ 40%] (Warmup) #> Chain 1: Iteration: 500 / 1000 [ 50%] (Warmup) #> Chain 1: Iteration: 501 / 1000 [ 50%] (Sampling) #> Chain 1: Iteration: 600 / 1000 [ 60%] (Sampling) #> Chain 1: Iteration: 700 / 1000 [ 70%] (Sampling) #> Chain 1: Iteration: 800 / 1000 [ 80%] (Sampling) #> Chain 1: Iteration: 900 / 1000 [ 90%] (Sampling) #> Chain 1: Iteration: 1000 / 1000 [100%] (Sampling) #> Chain 1: #> Chain 1: Elapsed Time: 52.512 seconds (Warm-up) #> Chain 1: 21.197 seconds (Sampling) #> Chain 1: 73.709 seconds (Total) #> Chain 1: model_weights<-mortality_weights(list(fitLC,fitRH)) # }