Bayesian Markov-Chain-Monte-Carlo and Claims Data

using MortalityTables
using Turing
using UUIDs
using Random
Random.seed!(20261012)
using DataFramesMeta
using MCMCChains
using CairoMakie
using LinearAlgebra
using Pipe
using StatsFuns

# Apply the site's custom Makie plot theme (cassette futurism).
include(joinpath(@__DIR__, "..", "..", "assets", "themes", "cassette_futurism.jl"))
set_theme!(cassette_futurism_theme())

# Small Makie equivalent of StatsPlots' plot(chain): trace + density per
# parameter, one row per parameter, chains overlaid in palette colors.
function chainplot(ch::Chains)
    params = names(ch, :parameters)
    nchains = size(ch, 3)
    fig = Figure(size=(900, 160 * length(params)))
    for (i, p) in enumerate(params)
        ax_t = Axis(fig[i, 1]; title=string(p), xlabel="iteration")
        ax_d = Axis(fig[i, 2]; title="density")
        for c in 1:nchains
            col = CASSETTE_PALETTE[mod1(c, length(CASSETTE_PALETTE))]
            v = vec(Array(ch[:, p, c]))
            lines!(ax_t, v; color=(col, 0.7), linewidth=0.8)
            density!(ax_d, v; color=(col, 0.0),
                strokecolor=col, strokewidth=1.5)
        end
    end
    return fig
end

chainplot (generic function with 1 method)

Generating fake data

The problem of interest is to look at mortality rates, which are given in terms of exposures (whether or not a life experienced a death in a given year).

We’ll grab some example rates from an insurance table, which has a “selection” component: When someone enters observation, say at age 50, their mortality is path dependent (so for someone who started being observed at 50 will have a different risk/mortality rate at age 55 than someone who started being observed at 45).

Addtionally, there may be additional groups of interest, such as: - high/medium/low risk classification - sex - group (e.g. company, data source, etc.) - type of insurance product offered

The example data will start with only the risk classification above

src = MortalityTables.table("2001 VBT Residual Standard Select and Ultimate - Male Nonsmoker, ANB")

MortalityTable (Insured Lives Mortality):
   Name:
       2001 VBT Residual Standard Select and Ultimate - Male Nonsmoker, ANB
   Fields:
       (:select, :ultimate, :metadata)
   Provider:
       Society of Actuaries
   mort.SOA.org ID:
       1118
   mort.SOA.org link:
       https://mort.soa.org/ViewTable.aspx?&TableIdentity=1118
   Description:
       2001 Valuation Basic Table (VBT) Residual Standard Select and Ultimate Table -  Male Nonsmoker. Basis: Age Nearest Birthday. Minimum Select Age: 0. Maximum Select Age: 99. Minimum Ultimate Age: 25. Maximum Ultimate Age: 120

n = 10_000

function generate_data_individual(tbl, issue_age=rand(50:55), inforce_years=rand(1:30), risklevel=rand(1:3))
    # risk_factors will scale the "true" parameter up or down
    # we observe the assigned risklevel, but not risk_factor
    risk_factors = [0.7, 1.0, 1.5]
    rf = risk_factors[risklevel]
    deaths = rand(inforce_years) .< (tbl.select[issue_age][issue_age .+ (1:inforce_years) .- 1] .* rf)

    endpoint = if sum(deaths) == 0
        inforce_years
    else
        findfirst(deaths)
    end
    id = uuid1()
    map(1:endpoint) do i
        (
            issue_age=issue_age,
            risklevel=risklevel,
            att_age=issue_age + i - 1,
            death=deaths[i],
            id=id,
        )
    end

end

exposures = vcat([generate_data_individual(src) for _ in 1:n]...) |> DataFrame

144068×5 DataFrame

144043 rows omitted

Row	issue_age	risklevel	att_age	death	id
	Int64	Int64	Int64	Bool	Base.UUID
1	54	3	54	false	1c3a8ef0-c66b-11f1-beec-7b302928fd1a
2	54	3	55	false	1c3a8ef0-c66b-11f1-beec-7b302928fd1a
3	54	3	56	false	1c3a8ef0-c66b-11f1-beec-7b302928fd1a
4	54	3	57	false	1c3a8ef0-c66b-11f1-beec-7b302928fd1a
5	54	3	58	false	1c3a8ef0-c66b-11f1-beec-7b302928fd1a
6	54	3	59	false	1c3a8ef0-c66b-11f1-beec-7b302928fd1a
7	54	3	60	false	1c3a8ef0-c66b-11f1-beec-7b302928fd1a
8	54	3	61	false	1c3a8ef0-c66b-11f1-beec-7b302928fd1a
9	54	3	62	false	1c3a8ef0-c66b-11f1-beec-7b302928fd1a
10	54	3	63	false	1c3a8ef0-c66b-11f1-beec-7b302928fd1a
11	54	3	64	false	1c3a8ef0-c66b-11f1-beec-7b302928fd1a
12	54	3	65	false	1c3a8ef0-c66b-11f1-beec-7b302928fd1a
13	54	3	66	false	1c3a8ef0-c66b-11f1-beec-7b302928fd1a
⋮	⋮	⋮	⋮	⋮	⋮
144057	55	1	68	false	1c4111da-c66b-11f1-b0c7-2b66a9cc9947
144058	55	1	69	false	1c4111da-c66b-11f1-b0c7-2b66a9cc9947
144059	55	1	70	false	1c4111da-c66b-11f1-b0c7-2b66a9cc9947
144060	55	1	71	false	1c4111da-c66b-11f1-b0c7-2b66a9cc9947
144061	55	1	72	false	1c4111da-c66b-11f1-b0c7-2b66a9cc9947
144062	55	1	73	false	1c4111da-c66b-11f1-b0c7-2b66a9cc9947
144063	55	1	74	false	1c4111da-c66b-11f1-b0c7-2b66a9cc9947
144064	55	1	75	false	1c4111da-c66b-11f1-b0c7-2b66a9cc9947
144065	55	1	76	false	1c4111da-c66b-11f1-b0c7-2b66a9cc9947
144066	55	1	77	false	1c4111da-c66b-11f1-b0c7-2b66a9cc9947
144067	55	1	78	false	1c4111da-c66b-11f1-b0c7-2b66a9cc9947
144068	55	1	79	false	1c4111da-c66b-11f1-b0c7-2b66a9cc9947

Two groupings, one with and without the risk level:

data = combine(groupby(exposures, [:issue_age, :att_age])) do subdf
    (exposures=nrow(subdf),
        deaths=sum(subdf.death),
        fraction=sum(subdf.death) / nrow(subdf))
end

180×5 DataFrame

155 rows omitted

Row	issue_age	att_age	exposures	deaths	fraction
	Int64	Int64	Int64	Int64	Float64
1	50	50	1653	2	0.00120992
2	50	51	1588	1	0.000629723
3	50	52	1533	5	0.00326158
4	50	53	1477	5	0.00338524
5	50	54	1412	3	0.00212465
6	50	55	1359	4	0.00294334
7	50	56	1292	10	0.00773994
8	50	57	1242	3	0.00241546
9	50	58	1173	6	0.00511509
10	50	59	1107	8	0.00722674
11	50	60	1048	8	0.00763359
12	50	61	989	12	0.0121335
13	50	62	925	9	0.00972973
⋮	⋮	⋮	⋮	⋮	⋮
169	55	73	534	25	0.0468165
170	55	74	471	7	0.014862
171	55	75	423	13	0.0307329
172	55	76	368	16	0.0434783
173	55	77	302	17	0.0562914
174	55	78	248	17	0.0685484
175	55	79	202	12	0.0594059
176	55	80	143	8	0.0559441
177	55	81	108	5	0.0462963
178	55	82	80	4	0.05
179	55	83	49	5	0.102041
180	55	84	22	0	0.0

data2 = combine(groupby(exposures, [:issue_age, :att_age, :risklevel])) do subdf
    (exposures=nrow(subdf),
        deaths=sum(subdf.death),
        fraction=sum(subdf.death) / nrow(subdf))
end

540×6 DataFrame

515 rows omitted

Row	issue_age	att_age	risklevel	exposures	deaths	fraction
	Int64	Int64	Int64	Int64	Int64	Float64
1	50	50	1	565	0	0.0
2	50	50	2	557	0	0.0
3	50	50	3	531	2	0.00376648
4	50	51	1	539	0	0.0
5	50	51	2	536	1	0.00186567
6	50	51	3	513	0	0.0
7	50	52	1	522	0	0.0
8	50	52	2	515	2	0.0038835
9	50	52	3	496	3	0.00604839
10	50	53	1	504	2	0.00396825
11	50	53	2	500	1	0.002
12	50	53	3	473	2	0.00422833
13	50	54	1	478	1	0.00209205
⋮	⋮	⋮	⋮	⋮	⋮	⋮
529	55	81	1	42	1	0.0238095
530	55	81	2	33	2	0.0606061
531	55	81	3	33	2	0.0606061
532	55	82	1	32	2	0.0625
533	55	82	2	22	1	0.0454545
534	55	82	3	26	1	0.0384615
535	55	83	1	21	2	0.0952381
536	55	83	2	11	0	0.0
537	55	83	3	17	3	0.176471
538	55	84	1	13	0	0.0
539	55	84	2	4	0	0.0
540	55	84	3	5	0	0.0

1: A single binomial parameter model

Estimate \(p\), the average mortality rate, not accounting for any variation within the population/sample:

@model function mortality(data, deaths)
    p ~ Beta(1, 1)
    for i = 1:nrow(data)
        deaths[i] ~ Binomial(data.exposures[i], p)
    end
end


m1 = mortality(data, data.deaths)

num_chains = 4

4

Sampling from the posterior

We use a No-U-Turn-Sampler (NUTS) technique to sample multiple chains at once:

chain = MCMCChains.Chains(sample(m1, NUTS(), MCMCThreads(), 1000, num_chains))

chainplot(chain)

Sampling (4 threads)   0%|                              |  ETA: N/A
┌ Info: Found initial step size
└   ϵ = 0.00625
┌ Info: Found initial step size
└   ϵ = 0.00625
┌ Info: Found initial step size
└   ϵ = 0.0125
┌ Info: Found initial step size
└   ϵ = 0.0125
Sampling (4 threads)   0%|▏                             |  ETA: 0:06:19
Sampling (4 threads)   1%|▎                             |  ETA: 0:03:10
Sampling (4 threads)   2%|▌                             |  ETA: 0:02:06
Sampling (4 threads)   2%|▋                             |  ETA: 0:01:34
Sampling (4 threads)   2%|▊                             |  ETA: 0:01:15
Sampling (4 threads)   3%|▉                             |  ETA: 0:01:02
Sampling (4 threads)   4%|█                             |  ETA: 0:00:53
Sampling (4 threads)   4%|█▎                            |  ETA: 0:00:46
Sampling (4 threads)   4%|█▍                            |  ETA: 0:00:41
Sampling (4 threads)   5%|█▌                            |  ETA: 0:00:37
Sampling (4 threads)   6%|█▋                            |  ETA: 0:00:33
Sampling (4 threads)   6%|█▊                            |  ETA: 0:00:30
Sampling (4 threads)   6%|██                            |  ETA: 0:00:28
Sampling (4 threads)   7%|██▏                           |  ETA: 0:00:26
Sampling (4 threads)   8%|██▎                           |  ETA: 0:00:24
Sampling (4 threads)   8%|██▍                           |  ETA: 0:00:23
Sampling (4 threads)   8%|██▌                           |  ETA: 0:00:21
Sampling (4 threads)   9%|██▊                           |  ETA: 0:00:20
Sampling (4 threads)  10%|██▉                           |  ETA: 0:00:19
Sampling (4 threads)  10%|███                           |  ETA: 0:00:18
Sampling (4 threads)  10%|███▏                          |  ETA: 0:00:17
Sampling (4 threads)  11%|███▎                          |  ETA: 0:00:16
Sampling (4 threads)  12%|███▌                          |  ETA: 0:00:15
Sampling (4 threads)  12%|███▋                          |  ETA: 0:00:15
Sampling (4 threads)  12%|███▊                          |  ETA: 0:00:14
Sampling (4 threads)  13%|███▉                          |  ETA: 0:00:13
Sampling (4 threads)  14%|████                          |  ETA: 0:00:13
Sampling (4 threads)  14%|████▎                         |  ETA: 0:00:12
Sampling (4 threads)  14%|████▍                         |  ETA: 0:00:12
Sampling (4 threads)  15%|████▌                         |  ETA: 0:00:11
Sampling (4 threads)  16%|████▋                         |  ETA: 0:00:11
Sampling (4 threads)  16%|████▊                         |  ETA: 0:00:11
Sampling (4 threads)  16%|█████                         |  ETA: 0:00:10
Sampling (4 threads)  17%|█████▏                        |  ETA: 0:00:10
Sampling (4 threads)  18%|█████▎                        |  ETA: 0:00:09
Sampling (4 threads)  18%|█████▍                        |  ETA: 0:00:09
Sampling (4 threads)  18%|█████▌                        |  ETA: 0:00:09
Sampling (4 threads)  19%|█████▊                        |  ETA: 0:00:09
Sampling (4 threads)  20%|█████▉                        |  ETA: 0:00:08
Sampling (4 threads)  20%|██████                        |  ETA: 0:00:08
Sampling (4 threads)  20%|██████▏                       |  ETA: 0:00:08
Sampling (4 threads)  21%|██████▎                       |  ETA: 0:00:08
Sampling (4 threads)  22%|██████▌                       |  ETA: 0:00:07
Sampling (4 threads)  22%|██████▋                       |  ETA: 0:00:07
Sampling (4 threads)  22%|██████▊                       |  ETA: 0:00:07
Sampling (4 threads)  23%|██████▉                       |  ETA: 0:00:07
Sampling (4 threads)  24%|███████                       |  ETA: 0:00:07
Sampling (4 threads)  24%|███████▎                      |  ETA: 0:00:06
Sampling (4 threads)  24%|███████▍                      |  ETA: 0:00:06
Sampling (4 threads)  25%|███████▌                      |  ETA: 0:00:06
Sampling (4 threads)  26%|███████▋                      |  ETA: 0:00:06
Sampling (4 threads)  26%|███████▊                      |  ETA: 0:00:06
Sampling (4 threads)  26%|████████                      |  ETA: 0:00:06
Sampling (4 threads)  27%|████████▏                     |  ETA: 0:00:05
Sampling (4 threads)  28%|████████▎                     |  ETA: 0:00:05
Sampling (4 threads)  28%|████████▍                     |  ETA: 0:00:05
Sampling (4 threads)  28%|████████▌                     |  ETA: 0:00:05
Sampling (4 threads)  29%|████████▊                     |  ETA: 0:00:05
Sampling (4 threads)  30%|████████▉                     |  ETA: 0:00:05
Sampling (4 threads)  30%|█████████                     |  ETA: 0:00:05
Sampling (4 threads)  30%|█████████▏                    |  ETA: 0:00:05
Sampling (4 threads)  31%|█████████▎                    |  ETA: 0:00:05
Sampling (4 threads)  32%|█████████▌                    |  ETA: 0:00:04
Sampling (4 threads)  32%|█████████▋                    |  ETA: 0:00:04
Sampling (4 threads)  32%|█████████▊                    |  ETA: 0:00:04
Sampling (4 threads)  33%|█████████▉                    |  ETA: 0:00:04
Sampling (4 threads)  34%|██████████                    |  ETA: 0:00:04
Sampling (4 threads)  34%|██████████▎                   |  ETA: 0:00:04
Sampling (4 threads)  34%|██████████▍                   |  ETA: 0:00:04
Sampling (4 threads)  35%|██████████▌                   |  ETA: 0:00:04
Sampling (4 threads)  36%|██████████▋                   |  ETA: 0:00:04
Sampling (4 threads)  36%|██████████▊                   |  ETA: 0:00:04
Sampling (4 threads)  36%|███████████                   |  ETA: 0:00:04
Sampling (4 threads)  37%|███████████▏                  |  ETA: 0:00:04
Sampling (4 threads)  38%|███████████▎                  |  ETA: 0:00:03
Sampling (4 threads)  38%|███████████▍                  |  ETA: 0:00:03
Sampling (4 threads)  38%|███████████▌                  |  ETA: 0:00:03
Sampling (4 threads)  39%|███████████▊                  |  ETA: 0:00:03
Sampling (4 threads)  40%|███████████▉                  |  ETA: 0:00:03
Sampling (4 threads)  40%|████████████                  |  ETA: 0:00:03
Sampling (4 threads)  40%|████████████▏                 |  ETA: 0:00:03
Sampling (4 threads)  41%|████████████▎                 |  ETA: 0:00:03
Sampling (4 threads)  42%|████████████▌                 |  ETA: 0:00:03
Sampling (4 threads)  42%|████████████▋                 |  ETA: 0:00:03
Sampling (4 threads)  42%|████████████▊                 |  ETA: 0:00:03
Sampling (4 threads)  43%|████████████▉                 |  ETA: 0:00:03
Sampling (4 threads)  44%|█████████████                 |  ETA: 0:00:03
Sampling (4 threads)  44%|█████████████▎                |  ETA: 0:00:03
Sampling (4 threads)  44%|█████████████▍                |  ETA: 0:00:03
Sampling (4 threads)  45%|█████████████▌                |  ETA: 0:00:03
Sampling (4 threads)  46%|█████████████▋                |  ETA: 0:00:03
Sampling (4 threads)  46%|█████████████▊                |  ETA: 0:00:02
Sampling (4 threads)  46%|██████████████                |  ETA: 0:00:02
Sampling (4 threads)  47%|██████████████▏               |  ETA: 0:00:02
Sampling (4 threads)  48%|██████████████▎               |  ETA: 0:00:02
Sampling (4 threads)  48%|██████████████▍               |  ETA: 0:00:02
Sampling (4 threads)  48%|██████████████▌               |  ETA: 0:00:02
Sampling (4 threads)  49%|██████████████▊               |  ETA: 0:00:02
Sampling (4 threads)  50%|██████████████▉               |  ETA: 0:00:02
Sampling (4 threads)  50%|███████████████               |  ETA: 0:00:02
Sampling (4 threads)  50%|███████████████▏              |  ETA: 0:00:02
Sampling (4 threads)  51%|███████████████▎              |  ETA: 0:00:02
Sampling (4 threads)  52%|███████████████▌              |  ETA: 0:00:02
Sampling (4 threads)  52%|███████████████▋              |  ETA: 0:00:02
Sampling (4 threads)  52%|███████████████▊              |  ETA: 0:00:02
Sampling (4 threads)  53%|███████████████▉              |  ETA: 0:00:02
Sampling (4 threads)  54%|████████████████              |  ETA: 0:00:02
Sampling (4 threads)  54%|████████████████▎             |  ETA: 0:00:02
Sampling (4 threads)  55%|████████████████▍             |  ETA: 0:00:02
Sampling (4 threads)  55%|████████████████▌             |  ETA: 0:00:02
Sampling (4 threads)  56%|████████████████▋             |  ETA: 0:00:02
Sampling (4 threads)  56%|████████████████▊             |  ETA: 0:00:02
Sampling (4 threads)  56%|█████████████████             |  ETA: 0:00:02
Sampling (4 threads)  57%|█████████████████▏            |  ETA: 0:00:02
Sampling (4 threads)  57%|█████████████████▎            |  ETA: 0:00:02
Sampling (4 threads)  58%|█████████████████▍            |  ETA: 0:00:02
Sampling (4 threads)  58%|█████████████████▌            |  ETA: 0:00:02
Sampling (4 threads)  59%|█████████████████▊            |  ETA: 0:00:01
Sampling (4 threads)  60%|█████████████████▉            |  ETA: 0:00:01
Sampling (4 threads)  60%|██████████████████            |  ETA: 0:00:01
Sampling (4 threads)  60%|██████████████████▏           |  ETA: 0:00:01
Sampling (4 threads)  61%|██████████████████▎           |  ETA: 0:00:01
Sampling (4 threads)  62%|██████████████████▌           |  ETA: 0:00:01
Sampling (4 threads)  62%|██████████████████▋           |  ETA: 0:00:01
Sampling (4 threads)  62%|██████████████████▊           |  ETA: 0:00:01
Sampling (4 threads)  63%|██████████████████▉           |  ETA: 0:00:01
Sampling (4 threads)  64%|███████████████████           |  ETA: 0:00:01
Sampling (4 threads)  64%|███████████████████▎          |  ETA: 0:00:01
Sampling (4 threads)  64%|███████████████████▍          |  ETA: 0:00:01
Sampling (4 threads)  65%|███████████████████▌          |  ETA: 0:00:01
Sampling (4 threads)  66%|███████████████████▋          |  ETA: 0:00:01
Sampling (4 threads)  66%|███████████████████▊          |  ETA: 0:00:01
Sampling (4 threads)  66%|████████████████████          |  ETA: 0:00:01
Sampling (4 threads)  67%|████████████████████▏         |  ETA: 0:00:01
Sampling (4 threads)  68%|████████████████████▎         |  ETA: 0:00:01
Sampling (4 threads)  68%|████████████████████▍         |  ETA: 0:00:01
Sampling (4 threads)  68%|████████████████████▌         |  ETA: 0:00:01
Sampling (4 threads)  69%|████████████████████▊         |  ETA: 0:00:01
Sampling (4 threads)  70%|████████████████████▉         |  ETA: 0:00:01
Sampling (4 threads)  70%|█████████████████████         |  ETA: 0:00:01
Sampling (4 threads)  70%|█████████████████████▏        |  ETA: 0:00:01
Sampling (4 threads)  71%|█████████████████████▎        |  ETA: 0:00:01
Sampling (4 threads)  72%|█████████████████████▌        |  ETA: 0:00:01
Sampling (4 threads)  72%|█████████████████████▋        |  ETA: 0:00:01
Sampling (4 threads)  72%|█████████████████████▊        |  ETA: 0:00:01
Sampling (4 threads)  73%|█████████████████████▉        |  ETA: 0:00:01
Sampling (4 threads)  74%|██████████████████████        |  ETA: 0:00:01
Sampling (4 threads)  74%|██████████████████████▎       |  ETA: 0:00:01
Sampling (4 threads)  74%|██████████████████████▍       |  ETA: 0:00:01
Sampling (4 threads)  75%|██████████████████████▌       |  ETA: 0:00:01
Sampling (4 threads)  76%|██████████████████████▋       |  ETA: 0:00:01
Sampling (4 threads)  76%|██████████████████████▊       |  ETA: 0:00:01
Sampling (4 threads)  76%|███████████████████████       |  ETA: 0:00:01
Sampling (4 threads)  77%|███████████████████████▏      |  ETA: 0:00:01
Sampling (4 threads)  78%|███████████████████████▎      |  ETA: 0:00:01
Sampling (4 threads)  78%|███████████████████████▍      |  ETA: 0:00:01
Sampling (4 threads)  78%|███████████████████████▌      |  ETA: 0:00:01
Sampling (4 threads)  79%|███████████████████████▊      |  ETA: 0:00:01
Sampling (4 threads)  80%|███████████████████████▉      |  ETA: 0:00:01
Sampling (4 threads)  80%|████████████████████████      |  ETA: 0:00:01
Sampling (4 threads)  80%|████████████████████████▏     |  ETA: 0:00:01
Sampling (4 threads)  81%|████████████████████████▎     |  ETA: 0:00:01
Sampling (4 threads)  82%|████████████████████████▌     |  ETA: 0:00:01
Sampling (4 threads)  82%|████████████████████████▋     |  ETA: 0:00:00
Sampling (4 threads)  82%|████████████████████████▊     |  ETA: 0:00:00
Sampling (4 threads)  83%|████████████████████████▉     |  ETA: 0:00:00
Sampling (4 threads)  84%|█████████████████████████     |  ETA: 0:00:00
Sampling (4 threads)  84%|█████████████████████████▎    |  ETA: 0:00:00
Sampling (4 threads)  84%|█████████████████████████▍    |  ETA: 0:00:00
Sampling (4 threads)  85%|█████████████████████████▌    |  ETA: 0:00:00
Sampling (4 threads)  86%|█████████████████████████▋    |  ETA: 0:00:00
Sampling (4 threads)  86%|█████████████████████████▊    |  ETA: 0:00:00
Sampling (4 threads)  86%|██████████████████████████    |  ETA: 0:00:00
Sampling (4 threads)  87%|██████████████████████████▏   |  ETA: 0:00:00
Sampling (4 threads)  88%|██████████████████████████▎   |  ETA: 0:00:00
Sampling (4 threads)  88%|██████████████████████████▍   |  ETA: 0:00:00
Sampling (4 threads)  88%|██████████████████████████▌   |  ETA: 0:00:00
Sampling (4 threads)  89%|██████████████████████████▊   |  ETA: 0:00:00
Sampling (4 threads)  90%|██████████████████████████▉   |  ETA: 0:00:00
Sampling (4 threads)  90%|███████████████████████████   |  ETA: 0:00:00
Sampling (4 threads)  90%|███████████████████████████▏  |  ETA: 0:00:00
Sampling (4 threads)  91%|███████████████████████████▎  |  ETA: 0:00:00
Sampling (4 threads)  92%|███████████████████████████▌  |  ETA: 0:00:00
Sampling (4 threads)  92%|███████████████████████████▋  |  ETA: 0:00:00
Sampling (4 threads)  92%|███████████████████████████▊  |  ETA: 0:00:00
Sampling (4 threads)  93%|███████████████████████████▉  |  ETA: 0:00:00
Sampling (4 threads)  94%|████████████████████████████  |  ETA: 0:00:00
Sampling (4 threads)  94%|████████████████████████████▎ |  ETA: 0:00:00
Sampling (4 threads)  94%|████████████████████████████▍ |  ETA: 0:00:00
Sampling (4 threads)  95%|████████████████████████████▌ |  ETA: 0:00:00
Sampling (4 threads)  96%|████████████████████████████▋ |  ETA: 0:00:00
Sampling (4 threads)  96%|████████████████████████████▊ |  ETA: 0:00:00
Sampling (4 threads)  96%|█████████████████████████████ |  ETA: 0:00:00
Sampling (4 threads)  97%|█████████████████████████████▏|  ETA: 0:00:00
Sampling (4 threads)  98%|█████████████████████████████▎|  ETA: 0:00:00
Sampling (4 threads)  98%|█████████████████████████████▍|  ETA: 0:00:00
Sampling (4 threads)  98%|█████████████████████████████▌|  ETA: 0:00:00
Sampling (4 threads)  99%|█████████████████████████████▊|  ETA: 0:00:00
Sampling (4 threads) 100%|█████████████████████████████▉|  ETA: 0:00:00
Sampling (4 threads) 100%|██████████████████████████████| Time: 0:00:02
Sampling (4 threads) 100%|██████████████████████████████| Time: 0:00:02

Plotting samples from the posterior

We can see that the sampling of possible posterior parameters doesn’t really fit the data very well since our model was so simplified. The lines represent the posterior binomial probability.

This is saying that for the observed data, if there really is just a single probability p that governs the true process that came up with the data, there’s a pretty narrow range of values it could possibly be:

let
    data_weight = data.exposures ./ sum(data.exposures)
    data_weight = .√(data_weight ./ maximum(data_weight) .* 20) .* 4  # Makie markersize ≈ 4× Plots

    fig = Figure()
    ax = Axis(fig[1,1]; xlabel="age", ylabel="mortality rate",
        title="Parametric Bayesian Mortality")

    # show n samples from the posterior plotted on the graph
    n = 300
    for _ in 1:n
        p_posterior = sample(chain, 1)[:p][1]
        hlines!(ax, [p_posterior]; color=(CASSETTE_PALETTE[3], 0.08))
    end

    scatter!(ax, data.att_age, data.fraction;
        markersize=data_weight,
        color=(CASSETTE_PALETTE[1], 0.5),
        strokewidth=0,
        label="Experience data point (size = relative exposure)")

    ylims!(ax, 0.0, 0.25)
    axislegend(ax, position=:rt)
    fig
end

The posterior mean of p is of course very close to the simple proportoin of claims to exposures:

let
    a = mean(chain, :p)
    b = sum(data.deaths) / sum(data.exposures)
    a, b
end

(0.01191242548308005, 0.0119041008412694)

2. Parametric model

We use the Makeham–Beard mortality law to describe how the force of mortality varies with age. hazard(model, age) is an instantaneous rate, not a death probability. The binomial likelihood needs the annual probability decrement(model, age, age + 1), obtained by integrating the hazard over that year. This probability is between zero and one for a valid mortality law.

The priors keep a, b and k positive, with c = 0. The binomial probability comes from decrement; the instantaneous hazard itself need not be bounded by one.

This isn’t an ideal parameterization (e.g. we aren’t including information about the select underwriting period), but is an example of utilizing Bayesian techniques on life experience data.

@model function mortality2(data, deaths)
    a ~ Exponential(0.1)
    b ~ Exponential(0.1)
    c = 0.0
    k ~ truncated(Exponential(1), 1, Inf)

    # use the variables to create a parametric mortality model
    m = MortalityTables.MakehamBeard(; a, b, c, k)

    # loop through the rows of the dataframe to let Turing observe the data
    # and how consistent the parameters are with the data
    for i = 1:nrow(data)
        age = data.att_age[i]
        q = clamp(MortalityTables.decrement(m, age, age + 1), 1e-10, 1 - 1e-10)
        deaths[i] ~ Binomial(data.exposures[i], q)
    end
end

mortality2 (generic function with 2 methods)

Sampling from the posterior

We combine the model with the data and a use use a No-U-Turn-Sampler (NUTS) technique to sample:

m2 = mortality2(data, data.deaths)
chain2 = MCMCChains.Chains(sample(m2, NUTS(), 1000))
summarize(chain2)

Sampling   0%|                                          |  ETA: N/A
┌ Info: Found initial step size
└   ϵ = 0.025
Sampling   1%|▎                                         |  ETA: 0:02:28
Sampling   1%|▍                                         |  ETA: 0:01:26
Sampling   2%|▋                                         |  ETA: 0:00:56
Sampling   2%|▉                                         |  ETA: 0:00:43
Sampling   3%|█▏                                        |  ETA: 0:00:34
Sampling   3%|█▎                                        |  ETA: 0:00:29
Sampling   4%|█▌                                        |  ETA: 0:00:25
Sampling   4%|█▋                                        |  ETA: 0:00:23
Sampling   5%|█▉                                        |  ETA: 0:00:22
Sampling   5%|██▏                                       |  ETA: 0:00:21
Sampling   6%|██▍                                       |  ETA: 0:00:21
Sampling   6%|██▌                                       |  ETA: 0:00:24
Sampling   7%|██▊                                       |  ETA: 0:00:32
Sampling   7%|███                                       |  ETA: 0:00:32
Sampling   8%|███▏                                      |  ETA: 0:00:33
Sampling   8%|███▍                                      |  ETA: 0:00:35
Sampling   9%|███▋                                      |  ETA: 0:00:37
Sampling   9%|███▊                                      |  ETA: 0:00:38
Sampling  10%|████                                      |  ETA: 0:00:36
Sampling  10%|████▎                                     |  ETA: 0:00:35
Sampling  11%|████▍                                     |  ETA: 0:00:34
Sampling  11%|████▋                                     |  ETA: 0:00:34
Sampling  12%|████▉                                     |  ETA: 0:00:33
Sampling  12%|█████                                     |  ETA: 0:00:34
Sampling  13%|█████▎                                    |  ETA: 0:00:33
Sampling  13%|█████▌                                    |  ETA: 0:00:33
Sampling  14%|█████▋                                    |  ETA: 0:00:33
Sampling  14%|█████▉                                    |  ETA: 0:00:32
Sampling  15%|██████▏                                   |  ETA: 0:00:31
Sampling  15%|██████▎                                   |  ETA: 0:00:31
Sampling  16%|██████▌                                   |  ETA: 0:00:31
Sampling  16%|██████▊                                   |  ETA: 0:00:30
Sampling  17%|███████                                   |  ETA: 0:00:30
Sampling  17%|███████▏                                  |  ETA: 0:00:30
Sampling  18%|███████▍                                  |  ETA: 0:00:30
Sampling  18%|███████▌                                  |  ETA: 0:00:30
Sampling  19%|███████▊                                  |  ETA: 0:00:29
Sampling  19%|████████                                  |  ETA: 0:00:29
Sampling  20%|████████▎                                 |  ETA: 0:00:28
Sampling  20%|████████▍                                 |  ETA: 0:00:28
Sampling  21%|████████▋                                 |  ETA: 0:00:27
Sampling  21%|████████▉                                 |  ETA: 0:00:27
Sampling  22%|█████████                                 |  ETA: 0:00:27
Sampling  22%|█████████▎                                |  ETA: 0:00:27
Sampling  23%|█████████▌                                |  ETA: 0:00:26
Sampling  23%|█████████▋                                |  ETA: 0:00:26
Sampling  24%|█████████▉                                |  ETA: 0:00:26
Sampling  24%|██████████▏                               |  ETA: 0:00:25
Sampling  25%|██████████▎                               |  ETA: 0:00:25
Sampling  25%|██████████▌                               |  ETA: 0:00:25
Sampling  26%|██████████▊                               |  ETA: 0:00:24
Sampling  26%|██████████▉                               |  ETA: 0:00:24
Sampling  27%|███████████▏                              |  ETA: 0:00:23
Sampling  27%|███████████▍                              |  ETA: 0:00:23
Sampling  28%|███████████▋                              |  ETA: 0:00:23
Sampling  28%|███████████▊                              |  ETA: 0:00:23
Sampling  29%|████████████                              |  ETA: 0:00:23
Sampling  29%|████████████▏                             |  ETA: 0:00:22
Sampling  30%|████████████▍                             |  ETA: 0:00:22
Sampling  30%|████████████▋                             |  ETA: 0:00:22
Sampling  31%|████████████▉                             |  ETA: 0:00:22
Sampling  31%|█████████████                             |  ETA: 0:00:22
Sampling  32%|█████████████▎                            |  ETA: 0:00:21
Sampling  32%|█████████████▌                            |  ETA: 0:00:21
Sampling  33%|█████████████▋                            |  ETA: 0:00:21
Sampling  33%|█████████████▉                            |  ETA: 0:00:20
Sampling  34%|██████████████▏                           |  ETA: 0:00:20
Sampling  34%|██████████████▎                           |  ETA: 0:00:20
Sampling  35%|██████████████▌                           |  ETA: 0:00:20
Sampling  35%|██████████████▊                           |  ETA: 0:00:20
Sampling  36%|██████████████▉                           |  ETA: 0:00:20
Sampling  36%|███████████████▏                          |  ETA: 0:00:20
Sampling  37%|███████████████▍                          |  ETA: 0:00:19
Sampling  37%|███████████████▌                          |  ETA: 0:00:19
Sampling  38%|███████████████▊                          |  ETA: 0:00:19
Sampling  38%|████████████████                          |  ETA: 0:00:19
Sampling  39%|████████████████▏                         |  ETA: 0:00:19
Sampling  39%|████████████████▍                         |  ETA: 0:00:19
Sampling  40%|████████████████▋                         |  ETA: 0:00:18
Sampling  40%|████████████████▊                         |  ETA: 0:00:18
Sampling  41%|█████████████████                         |  ETA: 0:00:18
Sampling  41%|█████████████████▎                        |  ETA: 0:00:18
Sampling  42%|█████████████████▌                        |  ETA: 0:00:18
Sampling  42%|█████████████████▋                        |  ETA: 0:00:18
Sampling  43%|█████████████████▉                        |  ETA: 0:00:17
Sampling  43%|██████████████████                        |  ETA: 0:00:17
Sampling  44%|██████████████████▎                       |  ETA: 0:00:17
Sampling  44%|██████████████████▌                       |  ETA: 0:00:17
Sampling  45%|██████████████████▊                       |  ETA: 0:00:17
Sampling  45%|██████████████████▉                       |  ETA: 0:00:17
Sampling  46%|███████████████████▏                      |  ETA: 0:00:16
Sampling  46%|███████████████████▍                      |  ETA: 0:00:16
Sampling  47%|███████████████████▌                      |  ETA: 0:00:16
Sampling  47%|███████████████████▊                      |  ETA: 0:00:16
Sampling  48%|████████████████████                      |  ETA: 0:00:16
Sampling  48%|████████████████████▏                     |  ETA: 0:00:15
Sampling  49%|████████████████████▍                     |  ETA: 0:00:15
Sampling  49%|████████████████████▋                     |  ETA: 0:00:15
Sampling  50%|████████████████████▊                     |  ETA: 0:00:15
Sampling  50%|█████████████████████                     |  ETA: 0:00:15
Sampling  51%|█████████████████████▎                    |  ETA: 0:00:15
Sampling  51%|█████████████████████▍                    |  ETA: 0:00:15
Sampling  52%|█████████████████████▋                    |  ETA: 0:00:15
Sampling  52%|█████████████████████▉                    |  ETA: 0:00:15
Sampling  53%|██████████████████████▏                   |  ETA: 0:00:15
Sampling  53%|██████████████████████▎                   |  ETA: 0:00:15
Sampling  54%|██████████████████████▌                   |  ETA: 0:00:14
Sampling  54%|██████████████████████▋                   |  ETA: 0:00:14
Sampling  55%|██████████████████████▉                   |  ETA: 0:00:14
Sampling  55%|███████████████████████▏                  |  ETA: 0:00:14
Sampling  56%|███████████████████████▍                  |  ETA: 0:00:14
Sampling  56%|███████████████████████▌                  |  ETA: 0:00:14
Sampling  57%|███████████████████████▊                  |  ETA: 0:00:13
Sampling  57%|████████████████████████                  |  ETA: 0:00:13
Sampling  58%|████████████████████████▏                 |  ETA: 0:00:13
Sampling  58%|████████████████████████▍                 |  ETA: 0:00:13
Sampling  59%|████████████████████████▋                 |  ETA: 0:00:13
Sampling  59%|████████████████████████▊                 |  ETA: 0:00:13
Sampling  60%|█████████████████████████                 |  ETA: 0:00:12
Sampling  60%|█████████████████████████▎                |  ETA: 0:00:12
Sampling  61%|█████████████████████████▍                |  ETA: 0:00:12
Sampling  61%|█████████████████████████▋                |  ETA: 0:00:12
Sampling  62%|█████████████████████████▉                |  ETA: 0:00:12
Sampling  62%|██████████████████████████                |  ETA: 0:00:12
Sampling  63%|██████████████████████████▎               |  ETA: 0:00:12
Sampling  63%|██████████████████████████▌               |  ETA: 0:00:11
Sampling  64%|██████████████████████████▋               |  ETA: 0:00:11
Sampling  64%|██████████████████████████▉               |  ETA: 0:00:11
Sampling  65%|███████████████████████████▏              |  ETA: 0:00:11
Sampling  65%|███████████████████████████▎              |  ETA: 0:00:11
Sampling  66%|███████████████████████████▌              |  ETA: 0:00:11
Sampling  66%|███████████████████████████▊              |  ETA: 0:00:11
Sampling  67%|████████████████████████████              |  ETA: 0:00:10
Sampling  67%|████████████████████████████▏             |  ETA: 0:00:10
Sampling  68%|████████████████████████████▍             |  ETA: 0:00:10
Sampling  68%|████████████████████████████▌             |  ETA: 0:00:10
Sampling  69%|████████████████████████████▊             |  ETA: 0:00:10
Sampling  69%|█████████████████████████████             |  ETA: 0:00:10
Sampling  70%|█████████████████████████████▎            |  ETA: 0:00:09
Sampling  70%|█████████████████████████████▍            |  ETA: 0:00:09
Sampling  71%|█████████████████████████████▋            |  ETA: 0:00:09
Sampling  71%|█████████████████████████████▉            |  ETA: 0:00:09
Sampling  72%|██████████████████████████████            |  ETA: 0:00:09
Sampling  72%|██████████████████████████████▎           |  ETA: 0:00:09
Sampling  73%|██████████████████████████████▌           |  ETA: 0:00:08
Sampling  73%|██████████████████████████████▋           |  ETA: 0:00:08
Sampling  74%|██████████████████████████████▉           |  ETA: 0:00:08
Sampling  74%|███████████████████████████████▏          |  ETA: 0:00:08
Sampling  75%|███████████████████████████████▎          |  ETA: 0:00:08
Sampling  75%|███████████████████████████████▌          |  ETA: 0:00:08
Sampling  76%|███████████████████████████████▊          |  ETA: 0:00:08
Sampling  76%|███████████████████████████████▉          |  ETA: 0:00:07
Sampling  77%|████████████████████████████████▏         |  ETA: 0:00:07
Sampling  77%|████████████████████████████████▍         |  ETA: 0:00:07
Sampling  78%|████████████████████████████████▋         |  ETA: 0:00:07
Sampling  78%|████████████████████████████████▊         |  ETA: 0:00:07
Sampling  79%|█████████████████████████████████         |  ETA: 0:00:07
Sampling  79%|█████████████████████████████████▏        |  ETA: 0:00:07
Sampling  80%|█████████████████████████████████▍        |  ETA: 0:00:06
Sampling  80%|█████████████████████████████████▋        |  ETA: 0:00:06
Sampling  81%|█████████████████████████████████▉        |  ETA: 0:00:06
Sampling  81%|██████████████████████████████████        |  ETA: 0:00:06
Sampling  82%|██████████████████████████████████▎       |  ETA: 0:00:06
Sampling  82%|██████████████████████████████████▌       |  ETA: 0:00:06
Sampling  83%|██████████████████████████████████▋       |  ETA: 0:00:05
Sampling  83%|██████████████████████████████████▉       |  ETA: 0:00:05
Sampling  84%|███████████████████████████████████▏      |  ETA: 0:00:05
Sampling  84%|███████████████████████████████████▎      |  ETA: 0:00:05
Sampling  85%|███████████████████████████████████▌      |  ETA: 0:00:05
Sampling  85%|███████████████████████████████████▊      |  ETA: 0:00:05
Sampling  86%|███████████████████████████████████▉      |  ETA: 0:00:05
Sampling  86%|████████████████████████████████████▏     |  ETA: 0:00:04
Sampling  87%|████████████████████████████████████▍     |  ETA: 0:00:04
Sampling  87%|████████████████████████████████████▌     |  ETA: 0:00:04
Sampling  88%|████████████████████████████████████▊     |  ETA: 0:00:04
Sampling  88%|█████████████████████████████████████     |  ETA: 0:00:04
Sampling  89%|█████████████████████████████████████▏    |  ETA: 0:00:04
Sampling  89%|█████████████████████████████████████▍    |  ETA: 0:00:03
Sampling  90%|█████████████████████████████████████▋    |  ETA: 0:00:03
Sampling  90%|█████████████████████████████████████▊    |  ETA: 0:00:03
Sampling  91%|██████████████████████████████████████    |  ETA: 0:00:03
Sampling  91%|██████████████████████████████████████▎   |  ETA: 0:00:03
Sampling  92%|██████████████████████████████████████▌   |  ETA: 0:00:03
Sampling  92%|██████████████████████████████████████▋   |  ETA: 0:00:02
Sampling  93%|██████████████████████████████████████▉   |  ETA: 0:00:02
Sampling  93%|███████████████████████████████████████   |  ETA: 0:00:02
Sampling  94%|███████████████████████████████████████▎  |  ETA: 0:00:02
Sampling  94%|███████████████████████████████████████▌  |  ETA: 0:00:02
Sampling  95%|███████████████████████████████████████▊  |  ETA: 0:00:02
Sampling  95%|███████████████████████████████████████▉  |  ETA: 0:00:02
Sampling  96%|████████████████████████████████████████▏ |  ETA: 0:00:01
Sampling  96%|████████████████████████████████████████▍ |  ETA: 0:00:01
Sampling  97%|████████████████████████████████████████▌ |  ETA: 0:00:01
Sampling  97%|████████████████████████████████████████▊ |  ETA: 0:00:01
Sampling  98%|█████████████████████████████████████████ |  ETA: 0:00:01
Sampling  98%|█████████████████████████████████████████▏|  ETA: 0:00:01
Sampling  99%|█████████████████████████████████████████▍|  ETA: 0:00:00
Sampling  99%|█████████████████████████████████████████▋|  ETA: 0:00:00
Sampling 100%|█████████████████████████████████████████▊|  ETA: 0:00:00
Sampling 100%|██████████████████████████████████████████| Time: 0:00:31
Sampling 100%|██████████████████████████████████████████| Time: 0:00:31

  parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   e ⋯
      Symbol   Float64   Float64   Float64    Float64    Float64   Float64     ⋯

           a    0.0000    0.0000    0.0000   159.8696   159.3169    0.9999     ⋯
           b    0.1460    0.0096    0.0008   156.2018   163.3244    1.0018     ⋯
           k    5.7711    2.1983    0.1737   162.7053   240.7273    1.0048     ⋯

                                                                1 column omitted

chainplot(chain2)

Plotting samples from the posterior

We can see that the sampling of possible posterior parameters fits the data well:

let
    data_weight = data.exposures ./ sum(data.exposures)
    data_weight = .√(data_weight ./ maximum(data_weight) .* 20) .* 4

    fig = Figure()
    ax = Axis(fig[1,1]; xlabel="age", ylabel="mortality rate",
        title="Parametric Bayesian Mortality")

    n = 300
    ages = sort!(unique(data.att_age))
    for _ in 1:n
        s = sample(chain2, 1)
        a = only(s[:a]); b = only(s[:b]); k = only(s[:k]); c = 0
        m = MortalityTables.MakehamBeard(; a, b, c, k)
        lines!(ax, ages, age -> MortalityTables.decrement(m, age, age + 1);
            color=(CASSETTE_PALETTE[3], 0.08), linewidth=0.8)
    end

    scatter!(ax, data.att_age, data.fraction;
        markersize=data_weight, color=(CASSETTE_PALETTE[1], 0.5),
        strokewidth=0,
        label="Experience data point (size = relative exposure)")

    ylims!(ax, 0.0, 0.25)
    axislegend(ax, position=:rt)
    fig
end

3. Parametric model

This model extends the prior to create a multi-level model. Each risk class (risklevel) gets its own \(a\) paramater in the MakhamBeard model. The prior for \(a_i\) is determined by the hyperparameter \(\bar{a}\).

@model function mortality3(data, deaths)
    risk_levels = length(levels(data.risklevel))
    b ~ Exponential(0.1)
    ā ~ Exponential(0.1)
    a ~ filldist(Exponential(ā), risk_levels)
    c = 0
    k ~ truncated(Exponential(1), 1, Inf)

    # use the variables to create a parametric mortality model

    # loop through the rows of the dataframe to let Turing observe the data 
    # and how consistent the parameters are with the data
    for i = 1:nrow(data)
        risk = data.risklevel[i]

        m = MortalityTables.MakehamBeard(; a=a[risk], b, c, k)
        age = data.att_age[i]
        q = clamp(MortalityTables.decrement(m, age, age + 1), 1e-10, 1 - 1e-10)
        deaths[i] ~ Binomial(data.exposures[i], q)
    end
end

mortality3 (generic function with 2 methods)

Instantiate model with the data, sample, and summarize:

m3 = mortality3(data2, data2.deaths)
chain3 = MCMCChains.Chains(sample(m3, NUTS(), 1000))
summarize(chain3)

Sampling   0%|                                          |  ETA: N/A
┌ Info: Found initial step size
└   ϵ = 0.025
Sampling   1%|▎                                         |  ETA: 0:04:25
Sampling   1%|▍                                         |  ETA: 0:02:51
Sampling   2%|▋                                         |  ETA: 0:02:10
Sampling   2%|▉                                         |  ETA: 0:01:46
Sampling   3%|█▏                                        |  ETA: 0:01:33
Sampling   3%|█▎                                        |  ETA: 0:01:32
Sampling   4%|█▌                                        |  ETA: 0:01:31
Sampling   4%|█▋                                        |  ETA: 0:02:01
Sampling   5%|█▉                                        |  ETA: 0:04:19
Sampling   5%|██▏                                       |  ETA: 0:06:12
Sampling   6%|██▍                                       |  ETA: 0:07:33
Sampling   6%|██▌                                       |  ETA: 0:08:42
Sampling   7%|██▊                                       |  ETA: 0:09:41
Sampling   7%|███                                       |  ETA: 0:09:32
Sampling   8%|███▏                                      |  ETA: 0:08:54
Sampling   8%|███▍                                      |  ETA: 0:08:29
Sampling   9%|███▋                                      |  ETA: 0:08:04
Sampling   9%|███▊                                      |  ETA: 0:07:44
Sampling  10%|████                                      |  ETA: 0:07:23
Sampling  10%|████▎                                     |  ETA: 0:07:04
Sampling  11%|████▍                                     |  ETA: 0:06:45
Sampling  11%|████▋                                     |  ETA: 0:06:32
Sampling  12%|████▉                                     |  ETA: 0:06:27
Sampling  12%|█████                                     |  ETA: 0:06:14
Sampling  13%|█████▎                                    |  ETA: 0:06:00
Sampling  13%|█████▌                                    |  ETA: 0:05:48
Sampling  14%|█████▋                                    |  ETA: 0:05:38
Sampling  14%|█████▉                                    |  ETA: 0:05:30
Sampling  15%|██████▏                                   |  ETA: 0:05:18
Sampling  15%|██████▎                                   |  ETA: 0:05:08
Sampling  16%|██████▌                                   |  ETA: 0:04:59
Sampling  16%|██████▊                                   |  ETA: 0:04:52
Sampling  17%|███████                                   |  ETA: 0:04:43
Sampling  17%|███████▏                                  |  ETA: 0:04:38
Sampling  18%|███████▍                                  |  ETA: 0:04:30
Sampling  18%|███████▌                                  |  ETA: 0:04:27
Sampling  19%|███████▊                                  |  ETA: 0:04:23
Sampling  19%|████████                                  |  ETA: 0:04:16
Sampling  20%|████████▎                                 |  ETA: 0:04:11
Sampling  20%|████████▍                                 |  ETA: 0:04:07
Sampling  21%|████████▋                                 |  ETA: 0:04:04
Sampling  21%|████████▉                                 |  ETA: 0:03:58
Sampling  22%|█████████                                 |  ETA: 0:03:52
Sampling  22%|█████████▎                                |  ETA: 0:03:49
Sampling  23%|█████████▌                                |  ETA: 0:03:45
Sampling  23%|█████████▋                                |  ETA: 0:03:42
Sampling  24%|█████████▉                                |  ETA: 0:03:38
Sampling  24%|██████████▏                               |  ETA: 0:03:33
Sampling  25%|██████████▎                               |  ETA: 0:03:29
Sampling  25%|██████████▌                               |  ETA: 0:03:26
Sampling  26%|██████████▊                               |  ETA: 0:03:22
Sampling  26%|██████████▉                               |  ETA: 0:03:19
Sampling  27%|███████████▏                              |  ETA: 0:03:16
Sampling  27%|███████████▍                              |  ETA: 0:03:13
Sampling  28%|███████████▋                              |  ETA: 0:03:10
Sampling  28%|███████████▊                              |  ETA: 0:03:07
Sampling  29%|████████████                              |  ETA: 0:03:03
Sampling  29%|████████████▏                             |  ETA: 0:03:00
Sampling  30%|████████████▍                             |  ETA: 0:02:57
Sampling  30%|████████████▋                             |  ETA: 0:02:55
Sampling  31%|████████████▉                             |  ETA: 0:02:52
Sampling  31%|█████████████                             |  ETA: 0:02:50
Sampling  32%|█████████████▎                            |  ETA: 0:02:48
Sampling  32%|█████████████▌                            |  ETA: 0:02:46
Sampling  33%|█████████████▋                            |  ETA: 0:02:44
Sampling  33%|█████████████▉                            |  ETA: 0:02:42
Sampling  34%|██████████████▏                           |  ETA: 0:02:40
Sampling  34%|██████████████▎                           |  ETA: 0:02:39
Sampling  35%|██████████████▌                           |  ETA: 0:02:37
Sampling  35%|██████████████▊                           |  ETA: 0:02:35
Sampling  36%|██████████████▉                           |  ETA: 0:02:32
Sampling  36%|███████████████▏                          |  ETA: 0:02:30
Sampling  37%|███████████████▍                          |  ETA: 0:02:29
Sampling  37%|███████████████▌                          |  ETA: 0:02:27
Sampling  38%|███████████████▊                          |  ETA: 0:02:24
Sampling  38%|████████████████                          |  ETA: 0:02:23
Sampling  39%|████████████████▏                         |  ETA: 0:02:21
Sampling  39%|████████████████▍                         |  ETA: 0:02:19
Sampling  40%|████████████████▋                         |  ETA: 0:02:17
Sampling  40%|████████████████▊                         |  ETA: 0:02:15
Sampling  41%|█████████████████                         |  ETA: 0:02:14
Sampling  41%|█████████████████▎                        |  ETA: 0:02:12
Sampling  42%|█████████████████▌                        |  ETA: 0:02:11
Sampling  42%|█████████████████▋                        |  ETA: 0:02:09
Sampling  43%|█████████████████▉                        |  ETA: 0:02:07
Sampling  43%|██████████████████                        |  ETA: 0:02:06
Sampling  44%|██████████████████▎                       |  ETA: 0:02:04
Sampling  44%|██████████████████▌                       |  ETA: 0:02:03
Sampling  45%|██████████████████▊                       |  ETA: 0:02:02
Sampling  45%|██████████████████▉                       |  ETA: 0:02:00
Sampling  46%|███████████████████▏                      |  ETA: 0:01:59
Sampling  46%|███████████████████▍                      |  ETA: 0:01:57
Sampling  47%|███████████████████▌                      |  ETA: 0:01:55
Sampling  47%|███████████████████▊                      |  ETA: 0:01:54
Sampling  48%|████████████████████                      |  ETA: 0:01:53
Sampling  48%|████████████████████▏                     |  ETA: 0:01:51
Sampling  49%|████████████████████▍                     |  ETA: 0:01:50
Sampling  49%|████████████████████▋                     |  ETA: 0:01:48
Sampling  50%|████████████████████▊                     |  ETA: 0:01:47
Sampling  50%|█████████████████████                     |  ETA: 0:01:45
Sampling  51%|█████████████████████▎                    |  ETA: 0:01:44
Sampling  51%|█████████████████████▍                    |  ETA: 0:01:43
Sampling  52%|█████████████████████▋                    |  ETA: 0:01:41
Sampling  52%|█████████████████████▉                    |  ETA: 0:01:40
Sampling  53%|██████████████████████▏                   |  ETA: 0:01:39
Sampling  53%|██████████████████████▎                   |  ETA: 0:01:38
Sampling  54%|██████████████████████▌                   |  ETA: 0:01:36
Sampling  54%|██████████████████████▋                   |  ETA: 0:01:36
Sampling  55%|██████████████████████▉                   |  ETA: 0:01:34
Sampling  55%|███████████████████████▏                  |  ETA: 0:01:33
Sampling  56%|███████████████████████▍                  |  ETA: 0:01:31
Sampling  56%|███████████████████████▌                  |  ETA: 0:01:30
Sampling  57%|███████████████████████▊                  |  ETA: 0:01:29
Sampling  57%|████████████████████████                  |  ETA: 0:01:28
Sampling  58%|████████████████████████▏                 |  ETA: 0:01:26
Sampling  58%|████████████████████████▍                 |  ETA: 0:01:25
Sampling  59%|████████████████████████▋                 |  ETA: 0:01:24
Sampling  59%|████████████████████████▊                 |  ETA: 0:01:23
Sampling  60%|█████████████████████████                 |  ETA: 0:01:22
Sampling  60%|█████████████████████████▎                |  ETA: 0:01:21
Sampling  61%|█████████████████████████▍                |  ETA: 0:01:20
Sampling  61%|█████████████████████████▋                |  ETA: 0:01:19
Sampling  62%|█████████████████████████▉                |  ETA: 0:01:17
Sampling  62%|██████████████████████████                |  ETA: 0:01:16
Sampling  63%|██████████████████████████▎               |  ETA: 0:01:15
Sampling  63%|██████████████████████████▌               |  ETA: 0:01:14
Sampling  64%|██████████████████████████▋               |  ETA: 0:01:12
Sampling  64%|██████████████████████████▉               |  ETA: 0:01:12
Sampling  65%|███████████████████████████▏              |  ETA: 0:01:10
Sampling  65%|███████████████████████████▎              |  ETA: 0:01:10
Sampling  66%|███████████████████████████▌              |  ETA: 0:01:08
Sampling  66%|███████████████████████████▊              |  ETA: 0:01:07
Sampling  67%|████████████████████████████              |  ETA: 0:01:06
Sampling  67%|████████████████████████████▏             |  ETA: 0:01:05
Sampling  68%|████████████████████████████▍             |  ETA: 0:01:04
Sampling  68%|████████████████████████████▌             |  ETA: 0:01:03
Sampling  69%|████████████████████████████▊             |  ETA: 0:01:02
Sampling  69%|█████████████████████████████             |  ETA: 0:01:01
Sampling  70%|█████████████████████████████▎            |  ETA: 0:01:00
Sampling  70%|█████████████████████████████▍            |  ETA: 0:00:59
Sampling  71%|█████████████████████████████▋            |  ETA: 0:00:58
Sampling  71%|█████████████████████████████▉            |  ETA: 0:00:57
Sampling  72%|██████████████████████████████            |  ETA: 0:00:56
Sampling  72%|██████████████████████████████▎           |  ETA: 0:00:55
Sampling  73%|██████████████████████████████▌           |  ETA: 0:00:54
Sampling  73%|██████████████████████████████▋           |  ETA: 0:00:53
Sampling  74%|██████████████████████████████▉           |  ETA: 0:00:51
Sampling  74%|███████████████████████████████▏          |  ETA: 0:00:50
Sampling  75%|███████████████████████████████▎          |  ETA: 0:00:49
Sampling  75%|███████████████████████████████▌          |  ETA: 0:00:48
Sampling  76%|███████████████████████████████▊          |  ETA: 0:00:47
Sampling  76%|███████████████████████████████▉          |  ETA: 0:00:46
Sampling  77%|████████████████████████████████▏         |  ETA: 0:00:45
Sampling  77%|████████████████████████████████▍         |  ETA: 0:00:44
Sampling  78%|████████████████████████████████▋         |  ETA: 0:00:43
Sampling  78%|████████████████████████████████▊         |  ETA: 0:00:42
Sampling  79%|█████████████████████████████████         |  ETA: 0:00:41
Sampling  79%|█████████████████████████████████▏        |  ETA: 0:00:40
Sampling  80%|█████████████████████████████████▍        |  ETA: 0:00:39
Sampling  80%|█████████████████████████████████▋        |  ETA: 0:00:38
Sampling  81%|█████████████████████████████████▉        |  ETA: 0:00:37
Sampling  81%|██████████████████████████████████        |  ETA: 0:00:36
Sampling  82%|██████████████████████████████████▎       |  ETA: 0:00:35
Sampling  82%|██████████████████████████████████▌       |  ETA: 0:00:34
Sampling  83%|██████████████████████████████████▋       |  ETA: 0:00:33
Sampling  83%|██████████████████████████████████▉       |  ETA: 0:00:32
Sampling  84%|███████████████████████████████████▏      |  ETA: 0:00:31
Sampling  84%|███████████████████████████████████▎      |  ETA: 0:00:30
Sampling  85%|███████████████████████████████████▌      |  ETA: 0:00:29
Sampling  85%|███████████████████████████████████▊      |  ETA: 0:00:28
Sampling  86%|███████████████████████████████████▉      |  ETA: 0:00:27
Sampling  86%|████████████████████████████████████▏     |  ETA: 0:00:26
Sampling  87%|████████████████████████████████████▍     |  ETA: 0:00:25
Sampling  87%|████████████████████████████████████▌     |  ETA: 0:00:24
Sampling  88%|████████████████████████████████████▊     |  ETA: 0:00:23
Sampling  88%|█████████████████████████████████████     |  ETA: 0:00:22
Sampling  89%|█████████████████████████████████████▏    |  ETA: 0:00:21
Sampling  89%|█████████████████████████████████████▍    |  ETA: 0:00:20
Sampling  90%|█████████████████████████████████████▋    |  ETA: 0:00:19
Sampling  90%|█████████████████████████████████████▊    |  ETA: 0:00:19
Sampling  91%|██████████████████████████████████████    |  ETA: 0:00:18
Sampling  91%|██████████████████████████████████████▎   |  ETA: 0:00:17
Sampling  92%|██████████████████████████████████████▌   |  ETA: 0:00:16
Sampling  92%|██████████████████████████████████████▋   |  ETA: 0:00:15
Sampling  93%|██████████████████████████████████████▉   |  ETA: 0:00:14
Sampling  93%|███████████████████████████████████████   |  ETA: 0:00:13
Sampling  94%|███████████████████████████████████████▎  |  ETA: 0:00:12
Sampling  94%|███████████████████████████████████████▌  |  ETA: 0:00:11
Sampling  95%|███████████████████████████████████████▊  |  ETA: 0:00:10
Sampling  95%|███████████████████████████████████████▉  |  ETA: 0:00:09
Sampling  96%|████████████████████████████████████████▏ |  ETA: 0:00:08
Sampling  96%|████████████████████████████████████████▍ |  ETA: 0:00:07
Sampling  97%|████████████████████████████████████████▌ |  ETA: 0:00:06
Sampling  97%|████████████████████████████████████████▊ |  ETA: 0:00:05
Sampling  98%|█████████████████████████████████████████ |  ETA: 0:00:04
Sampling  98%|█████████████████████████████████████████▏|  ETA: 0:00:04
Sampling  99%|█████████████████████████████████████████▍|  ETA: 0:00:03
Sampling  99%|█████████████████████████████████████████▋|  ETA: 0:00:02
Sampling 100%|█████████████████████████████████████████▊|  ETA: 0:00:01
Sampling 100%|██████████████████████████████████████████| Time: 0:03:00
Sampling 100%|██████████████████████████████████████████| Time: 0:03:00

  parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   e ⋯
      Symbol   Float64   Float64   Float64    Float64    Float64   Float64     ⋯

           b    0.1385    0.0059    0.0004   201.3491   316.9654    1.0042     ⋯
           ā    0.0000    0.0000    0.0000   253.4258   301.7022    0.9994     ⋯
        a[1]    0.0000    0.0000    0.0000   208.3459   349.5525    1.0043     ⋯
        a[2]    0.0000    0.0000    0.0000   209.3108   316.1254    1.0040     ⋯
        a[3]    0.0000    0.0000    0.0000   204.9295   334.8008    1.0042     ⋯
           k    3.2459    1.3230    0.0896   193.2515   116.3892    1.0033     ⋯

                                                                1 column omitted

let data = data2
    data_weight = data.exposures ./ sum(data.exposures)
    data_weight = .√(data_weight ./ maximum(data_weight) .* 20) .* 4

    risk_colors = CASSETTE_PALETTE[1:3]
    point_colors = [risk_colors[r] for r in Int.(data.risklevel)]

    fig = Figure()
    ax = Axis(fig[1,1]; xlabel="age", ylabel="mortality rate",
        title="Parametric Bayesian Mortality")

    n = 100
    ages = sort!(unique(data.att_age))
    for r in 1:3, i in 1:n
        s = sample(chain3, 1)
        a = only(s[Symbol("a[$r]")])
        b = only(s[:b]); k = only(s[:k]); c = 0
        m = MortalityTables.MakehamBeard(; a, b, c, k)
        lines!(ax, ages, age -> MortalityTables.decrement(m, age, age + 1);
            color=(risk_colors[r], 0.12), linewidth=0.8,
            label = (i == 1 ? "risk level $r" : nothing))
    end

    scatter!(ax, data.att_age, data.fraction;
        markersize=data_weight, color=point_colors, strokewidth=0)

    ylims!(ax, 0.0, 0.25)
    axislegend(ax, position=:rt, unique=true)
    fig
end

Handling non-unit exposures

The key is to use the Poisson distribution:

@model function mortality4(data, deaths)
    risk_levels = length(levels(data.risklevel))
    b ~ Exponential(0.1)
    ā ~ Exponential(0.1)
    a ~ filldist(Exponential(ā), risk_levels)
    c ~ Beta(4, 18)
    k ~ truncated(Exponential(1), 1, Inf)

    # use the variables to create a parametric mortality model

    # loop through the rows of the dataframe to let Turing observe the data 
    # and how consistent the parameters are with the data
    for i = 1:nrow(data)
        risk = data.risklevel[i]

        m = MortalityTables.MakehamBeard(; a=a[risk], b, c, k)
        age = data.att_age[i]
        q = MortalityTables.decrement(m, age, age + 1)
        deaths[i] ~ Poisson(data.exposures[i] * q)
    end
end

mortality4 (generic function with 2 methods)

Instantiate model with the data, sample, and summarize:

m4 = mortality4(data2, data2.deaths)
chain4 = MCMCChains.Chains(sample(m4, NUTS(), 1000))
summarize(chain4)

Sampling   0%|                                          |  ETA: N/A
┌ Info: Found initial step size
└   ϵ = 0.0013671875000000001
Sampling   1%|▎                                         |  ETA: 0:04:15
Sampling   1%|▍                                         |  ETA: 0:02:47
Sampling   2%|▋                                         |  ETA: 0:02:10
Sampling   2%|▉                                         |  ETA: 0:06:22
Sampling   3%|█▏                                        |  ETA: 0:13:07
Sampling   3%|█▎                                        |  ETA: 0:15:54
Sampling   4%|█▌                                        |  ETA: 0:18:39
Sampling   4%|█▋                                        |  ETA: 0:16:55
Sampling   5%|█▉                                        |  ETA: 0:15:02
Sampling   5%|██▏                                       |  ETA: 0:13:41
Sampling   6%|██▍                                       |  ETA: 0:12:28
Sampling   6%|██▌                                       |  ETA: 0:12:02
Sampling   7%|██▊                                       |  ETA: 0:11:55
Sampling   7%|███                                       |  ETA: 0:12:09
Sampling   8%|███▏                                      |  ETA: 0:11:27
Sampling   8%|███▍                                      |  ETA: 0:11:17
Sampling   9%|███▋                                      |  ETA: 0:11:21
Sampling   9%|███▊                                      |  ETA: 0:11:11
Sampling  10%|████                                      |  ETA: 0:10:46
Sampling  10%|████▎                                     |  ETA: 0:10:31
Sampling  11%|████▍                                     |  ETA: 0:10:00
Sampling  11%|████▋                                     |  ETA: 0:09:44
Sampling  12%|████▉                                     |  ETA: 0:09:29
Sampling  12%|█████                                     |  ETA: 0:09:15
Sampling  13%|█████▎                                    |  ETA: 0:08:51
Sampling  13%|█████▌                                    |  ETA: 0:08:34
Sampling  14%|█████▋                                    |  ETA: 0:08:17
Sampling  14%|█████▉                                    |  ETA: 0:08:03
Sampling  15%|██████▏                                   |  ETA: 0:07:53
Sampling  15%|██████▎                                   |  ETA: 0:07:42
Sampling  16%|██████▌                                   |  ETA: 0:07:30
Sampling  16%|██████▊                                   |  ETA: 0:07:20
Sampling  17%|███████                                   |  ETA: 0:07:10
Sampling  17%|███████▏                                  |  ETA: 0:06:59
Sampling  18%|███████▍                                  |  ETA: 0:06:49
Sampling  18%|███████▌                                  |  ETA: 0:06:40
Sampling  19%|███████▊                                  |  ETA: 0:06:30
Sampling  19%|████████                                  |  ETA: 0:06:21
Sampling  20%|████████▎                                 |  ETA: 0:06:13
Sampling  20%|████████▍                                 |  ETA: 0:06:05
Sampling  21%|████████▋                                 |  ETA: 0:05:57
Sampling  21%|████████▉                                 |  ETA: 0:05:49
Sampling  22%|█████████                                 |  ETA: 0:05:41
Sampling  22%|█████████▎                                |  ETA: 0:05:33
Sampling  23%|█████████▌                                |  ETA: 0:05:25
Sampling  23%|█████████▋                                |  ETA: 0:05:21
Sampling  24%|█████████▉                                |  ETA: 0:05:14
Sampling  24%|██████████▏                               |  ETA: 0:05:09
Sampling  25%|██████████▎                               |  ETA: 0:05:01
Sampling  25%|██████████▌                               |  ETA: 0:04:57
Sampling  26%|██████████▊                               |  ETA: 0:04:51
Sampling  26%|██████████▉                               |  ETA: 0:04:47
Sampling  27%|███████████▏                              |  ETA: 0:04:43
Sampling  27%|███████████▍                              |  ETA: 0:04:40
Sampling  28%|███████████▋                              |  ETA: 0:04:34
Sampling  28%|███████████▊                              |  ETA: 0:04:30
Sampling  29%|████████████                              |  ETA: 0:04:25
Sampling  29%|████████████▏                             |  ETA: 0:04:23
Sampling  30%|████████████▍                             |  ETA: 0:04:18
Sampling  30%|████████████▋                             |  ETA: 0:04:14
Sampling  31%|████████████▉                             |  ETA: 0:04:11
Sampling  31%|█████████████                             |  ETA: 0:04:07
Sampling  32%|█████████████▎                            |  ETA: 0:04:04
Sampling  32%|█████████████▌                            |  ETA: 0:04:01
Sampling  33%|█████████████▋                            |  ETA: 0:03:57
Sampling  33%|█████████████▉                            |  ETA: 0:03:53
Sampling  34%|██████████████▏                           |  ETA: 0:03:50
Sampling  34%|██████████████▎                           |  ETA: 0:03:47
Sampling  35%|██████████████▌                           |  ETA: 0:03:44
Sampling  35%|██████████████▊                           |  ETA: 0:03:41
Sampling  36%|██████████████▉                           |  ETA: 0:03:39
Sampling  36%|███████████████▏                          |  ETA: 0:03:37
Sampling  37%|███████████████▍                          |  ETA: 0:03:34
Sampling  37%|███████████████▌                          |  ETA: 0:03:32
Sampling  38%|███████████████▊                          |  ETA: 0:03:29
Sampling  38%|████████████████                          |  ETA: 0:03:27
Sampling  39%|████████████████▏                         |  ETA: 0:03:24
Sampling  39%|████████████████▍                         |  ETA: 0:03:21
Sampling  40%|████████████████▋                         |  ETA: 0:03:19
Sampling  40%|████████████████▊                         |  ETA: 0:03:16
Sampling  41%|█████████████████                         |  ETA: 0:03:14
Sampling  41%|█████████████████▎                        |  ETA: 0:03:12
Sampling  42%|█████████████████▌                        |  ETA: 0:03:09
Sampling  42%|█████████████████▋                        |  ETA: 0:03:07
Sampling  43%|█████████████████▉                        |  ETA: 0:03:05
Sampling  43%|██████████████████                        |  ETA: 0:03:02
Sampling  44%|██████████████████▎                       |  ETA: 0:03:00
Sampling  44%|██████████████████▌                       |  ETA: 0:02:58
Sampling  45%|██████████████████▊                       |  ETA: 0:02:56
Sampling  45%|██████████████████▉                       |  ETA: 0:02:54
Sampling  46%|███████████████████▏                      |  ETA: 0:02:52
Sampling  46%|███████████████████▍                      |  ETA: 0:02:50
Sampling  47%|███████████████████▌                      |  ETA: 0:02:47
Sampling  47%|███████████████████▊                      |  ETA: 0:02:45
Sampling  48%|████████████████████                      |  ETA: 0:02:43
Sampling  48%|████████████████████▏                     |  ETA: 0:02:42
Sampling  49%|████████████████████▍                     |  ETA: 0:02:39
Sampling  49%|████████████████████▋                     |  ETA: 0:02:37
Sampling  50%|████████████████████▊                     |  ETA: 0:02:35
Sampling  50%|█████████████████████                     |  ETA: 0:02:33
Sampling  51%|█████████████████████▎                    |  ETA: 0:02:31
Sampling  51%|█████████████████████▍                    |  ETA: 0:02:29
Sampling  52%|█████████████████████▋                    |  ETA: 0:02:27
Sampling  52%|█████████████████████▉                    |  ETA: 0:02:25
Sampling  53%|██████████████████████▏                   |  ETA: 0:02:23
Sampling  53%|██████████████████████▎                   |  ETA: 0:02:21
Sampling  54%|██████████████████████▌                   |  ETA: 0:02:19
Sampling  54%|██████████████████████▋                   |  ETA: 0:02:18
Sampling  55%|██████████████████████▉                   |  ETA: 0:02:16
Sampling  55%|███████████████████████▏                  |  ETA: 0:02:14
Sampling  56%|███████████████████████▍                  |  ETA: 0:02:13
Sampling  56%|███████████████████████▌                  |  ETA: 0:02:11
Sampling  57%|███████████████████████▊                  |  ETA: 0:02:09
Sampling  57%|████████████████████████                  |  ETA: 0:02:08
Sampling  58%|████████████████████████▏                 |  ETA: 0:02:06
Sampling  58%|████████████████████████▍                 |  ETA: 0:02:04
Sampling  59%|████████████████████████▋                 |  ETA: 0:02:02
Sampling  59%|████████████████████████▊                 |  ETA: 0:02:01
Sampling  60%|█████████████████████████                 |  ETA: 0:01:59
Sampling  60%|█████████████████████████▎                |  ETA: 0:01:58
Sampling  61%|█████████████████████████▍                |  ETA: 0:01:56
Sampling  61%|█████████████████████████▋                |  ETA: 0:01:54
Sampling  62%|█████████████████████████▉                |  ETA: 0:01:52
Sampling  62%|██████████████████████████                |  ETA: 0:01:51
Sampling  63%|██████████████████████████▎               |  ETA: 0:01:49
Sampling  63%|██████████████████████████▌               |  ETA: 0:01:48
Sampling  64%|██████████████████████████▋               |  ETA: 0:01:46
Sampling  64%|██████████████████████████▉               |  ETA: 0:01:45
Sampling  65%|███████████████████████████▏              |  ETA: 0:01:43
Sampling  65%|███████████████████████████▎              |  ETA: 0:01:41
Sampling  66%|███████████████████████████▌              |  ETA: 0:01:40
Sampling  66%|███████████████████████████▊              |  ETA: 0:01:38
Sampling  67%|████████████████████████████              |  ETA: 0:01:36
Sampling  67%|████████████████████████████▏             |  ETA: 0:01:35
Sampling  68%|████████████████████████████▍             |  ETA: 0:01:33
Sampling  68%|████████████████████████████▌             |  ETA: 0:01:32
Sampling  69%|████████████████████████████▊             |  ETA: 0:01:31
Sampling  69%|█████████████████████████████             |  ETA: 0:01:29
Sampling  70%|█████████████████████████████▎            |  ETA: 0:01:27
Sampling  70%|█████████████████████████████▍            |  ETA: 0:01:26
Sampling  71%|█████████████████████████████▋            |  ETA: 0:01:24
Sampling  71%|█████████████████████████████▉            |  ETA: 0:01:23
Sampling  72%|██████████████████████████████            |  ETA: 0:01:21
Sampling  72%|██████████████████████████████▎           |  ETA: 0:01:19
Sampling  73%|██████████████████████████████▌           |  ETA: 0:01:18
Sampling  73%|██████████████████████████████▋           |  ETA: 0:01:16
Sampling  74%|██████████████████████████████▉           |  ETA: 0:01:15
Sampling  74%|███████████████████████████████▏          |  ETA: 0:01:13
Sampling  75%|███████████████████████████████▎          |  ETA: 0:01:12
Sampling  75%|███████████████████████████████▌          |  ETA: 0:01:11
Sampling  76%|███████████████████████████████▊          |  ETA: 0:01:09
Sampling  76%|███████████████████████████████▉          |  ETA: 0:01:08
Sampling  77%|████████████████████████████████▏         |  ETA: 0:01:06
Sampling  77%|████████████████████████████████▍         |  ETA: 0:01:05
Sampling  78%|████████████████████████████████▋         |  ETA: 0:01:03
Sampling  78%|████████████████████████████████▊         |  ETA: 0:01:02
Sampling  79%|█████████████████████████████████         |  ETA: 0:01:00
Sampling  79%|█████████████████████████████████▏        |  ETA: 0:00:59
Sampling  80%|█████████████████████████████████▍        |  ETA: 0:00:57
Sampling  80%|█████████████████████████████████▋        |  ETA: 0:00:56
Sampling  81%|█████████████████████████████████▉        |  ETA: 0:00:54
Sampling  81%|██████████████████████████████████        |  ETA: 0:00:53
Sampling  82%|██████████████████████████████████▎       |  ETA: 0:00:51
Sampling  82%|██████████████████████████████████▌       |  ETA: 0:00:50
Sampling  83%|██████████████████████████████████▋       |  ETA: 0:00:48
Sampling  83%|██████████████████████████████████▉       |  ETA: 0:00:47
Sampling  84%|███████████████████████████████████▏      |  ETA: 0:00:45
Sampling  84%|███████████████████████████████████▎      |  ETA: 0:00:44
Sampling  85%|███████████████████████████████████▌      |  ETA: 0:00:43
Sampling  85%|███████████████████████████████████▊      |  ETA: 0:00:41
Sampling  86%|███████████████████████████████████▉      |  ETA: 0:00:40
Sampling  86%|████████████████████████████████████▏     |  ETA: 0:00:38
Sampling  87%|████████████████████████████████████▍     |  ETA: 0:00:37
Sampling  87%|████████████████████████████████████▌     |  ETA: 0:00:36
Sampling  88%|████████████████████████████████████▊     |  ETA: 0:00:34
Sampling  88%|█████████████████████████████████████     |  ETA: 0:00:33
Sampling  89%|█████████████████████████████████████▏    |  ETA: 0:00:31
Sampling  89%|█████████████████████████████████████▍    |  ETA: 0:00:30
Sampling  90%|█████████████████████████████████████▋    |  ETA: 0:00:29
Sampling  90%|█████████████████████████████████████▊    |  ETA: 0:00:27
Sampling  91%|██████████████████████████████████████    |  ETA: 0:00:26
Sampling  91%|██████████████████████████████████████▎   |  ETA: 0:00:24
Sampling  92%|██████████████████████████████████████▌   |  ETA: 0:00:23
Sampling  92%|██████████████████████████████████████▋   |  ETA: 0:00:22
Sampling  93%|██████████████████████████████████████▉   |  ETA: 0:00:20
Sampling  93%|███████████████████████████████████████   |  ETA: 0:00:19
Sampling  94%|███████████████████████████████████████▎  |  ETA: 0:00:17
Sampling  94%|███████████████████████████████████████▌  |  ETA: 0:00:16
Sampling  95%|███████████████████████████████████████▊  |  ETA: 0:00:15
Sampling  95%|███████████████████████████████████████▉  |  ETA: 0:00:13
Sampling  96%|████████████████████████████████████████▏ |  ETA: 0:00:12
Sampling  96%|████████████████████████████████████████▍ |  ETA: 0:00:11
Sampling  97%|████████████████████████████████████████▌ |  ETA: 0:00:09
Sampling  97%|████████████████████████████████████████▊ |  ETA: 0:00:08
Sampling  98%|█████████████████████████████████████████ |  ETA: 0:00:07
Sampling  98%|█████████████████████████████████████████▏|  ETA: 0:00:05
Sampling  99%|█████████████████████████████████████████▍|  ETA: 0:00:04
Sampling  99%|█████████████████████████████████████████▋|  ETA: 0:00:03
Sampling 100%|█████████████████████████████████████████▊|  ETA: 0:00:01
Sampling 100%|██████████████████████████████████████████| Time: 0:04:28
Sampling 100%|██████████████████████████████████████████| Time: 0:04:28

  parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   e ⋯
      Symbol   Float64   Float64   Float64    Float64    Float64   Float64     ⋯

           b    0.1478    0.0088    0.0007   140.9338   270.9180    1.0103     ⋯
           ā    0.0000    0.0000    0.0000   274.2959   357.8357    1.0102     ⋯
        a[1]    0.0000    0.0000    0.0000   144.9663   288.4517    1.0091     ⋯
        a[2]    0.0000    0.0000    0.0000   143.7253   280.2862    1.0101     ⋯
        a[3]    0.0000    0.0000    0.0000   144.6586   259.1196    1.0099     ⋯
           c    0.0006    0.0002    0.0000   275.1041   326.8779    1.0005     ⋯
           k    3.9083    1.6700    0.1320   134.0132    97.2392    1.0111     ⋯

                                                                1 column omitted

risk_factors4 = [mean(chain4[Symbol("a[$f]")]) for f in 1:3]

3-element Vector{Float64}:
 5.777698173226399e-7
 8.685330222156446e-7
 1.2745556502583303e-6

let data = data2
    data_weight = data.exposures ./ sum(data.exposures)
    data_weight = .√(data_weight ./ maximum(data_weight) .* 20) .* 4

    risk_colors = CASSETTE_PALETTE[1:3]
    point_colors = [risk_colors[r] for r in Int.(data.risklevel)]

    fig = Figure()
    ax = Axis(fig[1,1]; xlabel="age", ylabel="mortality rate",
        title="Parametric Bayesian Mortality")

    n = 100
    ages = sort!(unique(data.att_age))
    for r in 1:3, i in 1:n
        s = sample(chain4, 1)
        a = only(s[Symbol("a[$r]")])
        b = only(s[:b]); k = only(s[:k]); c = only(s[:c])
        m = MortalityTables.MakehamBeard(; a, b, c, k)
        lines!(ax, ages, age -> MortalityTables.decrement(m, age, age + 1);
            color=(risk_colors[r], 0.12), linewidth=0.8,
            label = (i == 1 ? "risk level $r" : nothing))
    end

    scatter!(ax, data.att_age, data.fraction;
        markersize=data_weight, color=point_colors, strokewidth=0)

    ylims!(ax, 0.0, 0.25)
    axislegend(ax, position=:rt, unique=true)
    fig
end

Predictions

We can generate predictive estimates by passing a vector of missing in place of the outcome variables and then calling predict.

We get a table of values where each row is the the prediction implied by the corresponding chain sample, and the columns are the predicted value for each of the outcomes in our original dataset.

preds = predict(mortality4(data2, fill(missing, length(data2.deaths))), chain4)

Chains MCMC chain (1000×540×1 Array{Float64, 3}):

Iterations        = 1:1:1000
Number of chains  = 1
Samples per chain = 1000
parameters        = deaths[1], deaths[2], deaths[3], deaths[4], deaths[5], deaths[6], deaths[7], deaths[8], deaths[9], deaths[10], deaths[11], deaths[12], deaths[13], deaths[14], deaths[15], deaths[16], deaths[17], deaths[18], deaths[19], deaths[20], deaths[21], deaths[22], deaths[23], deaths[24], deaths[25], deaths[26], deaths[27], deaths[28], deaths[29], deaths[30], deaths[31], deaths[32], deaths[33], deaths[34], deaths[35], deaths[36], deaths[37], deaths[38], deaths[39], deaths[40], deaths[41], deaths[42], deaths[43], deaths[44], deaths[45], deaths[46], deaths[47], deaths[48], deaths[49], deaths[50], deaths[51], deaths[52], deaths[53], deaths[54], deaths[55], deaths[56], deaths[57], deaths[58], deaths[59], deaths[60], deaths[61], deaths[62], deaths[63], deaths[64], deaths[65], deaths[66], deaths[67], deaths[68], deaths[69], deaths[70], deaths[71], deaths[72], deaths[73], deaths[74], deaths[75], deaths[76], deaths[77], deaths[78], deaths[79], deaths[80], deaths[81], deaths[82], deaths[83], deaths[84], deaths[85], deaths[86], deaths[87], deaths[88], deaths[89], deaths[90], deaths[91], deaths[92], deaths[93], deaths[94], deaths[95], deaths[96], deaths[97], deaths[98], deaths[99], deaths[100], deaths[101], deaths[102], deaths[103], deaths[104], deaths[105], deaths[106], deaths[107], deaths[108], deaths[109], deaths[110], deaths[111], deaths[112], deaths[113], deaths[114], deaths[115], deaths[116], deaths[117], deaths[118], deaths[119], deaths[120], deaths[121], deaths[122], deaths[123], deaths[124], deaths[125], deaths[126], deaths[127], deaths[128], deaths[129], deaths[130], deaths[131], deaths[132], deaths[133], deaths[134], deaths[135], deaths[136], deaths[137], deaths[138], deaths[139], deaths[140], deaths[141], deaths[142], deaths[143], deaths[144], deaths[145], deaths[146], deaths[147], deaths[148], deaths[149], deaths[150], deaths[151], deaths[152], deaths[153], deaths[154], deaths[155], deaths[156], deaths[157], deaths[158], deaths[159], deaths[160], deaths[161], deaths[162], deaths[163], deaths[164], deaths[165], deaths[166], deaths[167], deaths[168], deaths[169], deaths[170], deaths[171], deaths[172], deaths[173], deaths[174], deaths[175], deaths[176], deaths[177], deaths[178], deaths[179], deaths[180], deaths[181], deaths[182], deaths[183], deaths[184], deaths[185], deaths[186], deaths[187], deaths[188], deaths[189], deaths[190], deaths[191], deaths[192], deaths[193], deaths[194], deaths[195], deaths[196], deaths[197], deaths[198], deaths[199], deaths[200], deaths[201], deaths[202], deaths[203], deaths[204], deaths[205], deaths[206], deaths[207], deaths[208], deaths[209], deaths[210], deaths[211], deaths[212], deaths[213], deaths[214], deaths[215], deaths[216], deaths[217], deaths[218], deaths[219], deaths[220], deaths[221], deaths[222], deaths[223], deaths[224], deaths[225], deaths[226], deaths[227], deaths[228], deaths[229], deaths[230], deaths[231], deaths[232], deaths[233], deaths[234], deaths[235], deaths[236], deaths[237], deaths[238], deaths[239], deaths[240], deaths[241], deaths[242], deaths[243], deaths[244], deaths[245], deaths[246], deaths[247], deaths[248], deaths[249], deaths[250], deaths[251], deaths[252], deaths[253], deaths[254], deaths[255], deaths[256], deaths[257], deaths[258], deaths[259], deaths[260], deaths[261], deaths[262], deaths[263], deaths[264], deaths[265], deaths[266], deaths[267], deaths[268], deaths[269], deaths[270], deaths[271], deaths[272], deaths[273], deaths[274], deaths[275], deaths[276], deaths[277], deaths[278], deaths[279], deaths[280], deaths[281], deaths[282], deaths[283], deaths[284], deaths[285], deaths[286], deaths[287], deaths[288], deaths[289], deaths[290], deaths[291], deaths[292], deaths[293], deaths[294], deaths[295], deaths[296], deaths[297], deaths[298], deaths[299], deaths[300], deaths[301], deaths[302], deaths[303], deaths[304], deaths[305], deaths[306], deaths[307], deaths[308], deaths[309], deaths[310], deaths[311], deaths[312], deaths[313], deaths[314], deaths[315], deaths[316], deaths[317], deaths[318], deaths[319], deaths[320], deaths[321], deaths[322], deaths[323], deaths[324], deaths[325], deaths[326], deaths[327], deaths[328], deaths[329], deaths[330], deaths[331], deaths[332], deaths[333], deaths[334], deaths[335], deaths[336], deaths[337], deaths[338], deaths[339], deaths[340], deaths[341], deaths[342], deaths[343], deaths[344], deaths[345], deaths[346], deaths[347], deaths[348], deaths[349], deaths[350], deaths[351], deaths[352], deaths[353], deaths[354], deaths[355], deaths[356], deaths[357], deaths[358], deaths[359], deaths[360], deaths[361], deaths[362], deaths[363], deaths[364], deaths[365], deaths[366], deaths[367], deaths[368], deaths[369], deaths[370], deaths[371], deaths[372], deaths[373], deaths[374], deaths[375], deaths[376], deaths[377], deaths[378], deaths[379], deaths[380], deaths[381], deaths[382], deaths[383], deaths[384], deaths[385], deaths[386], deaths[387], deaths[388], deaths[389], deaths[390], deaths[391], deaths[392], deaths[393], deaths[394], deaths[395], deaths[396], deaths[397], deaths[398], deaths[399], deaths[400], deaths[401], deaths[402], deaths[403], deaths[404], deaths[405], deaths[406], deaths[407], deaths[408], deaths[409], deaths[410], deaths[411], deaths[412], deaths[413], deaths[414], deaths[415], deaths[416], deaths[417], deaths[418], deaths[419], deaths[420], deaths[421], deaths[422], deaths[423], deaths[424], deaths[425], deaths[426], deaths[427], deaths[428], deaths[429], deaths[430], deaths[431], deaths[432], deaths[433], deaths[434], deaths[435], deaths[436], deaths[437], deaths[438], deaths[439], deaths[440], deaths[441], deaths[442], deaths[443], deaths[444], deaths[445], deaths[446], deaths[447], deaths[448], deaths[449], deaths[450], deaths[451], deaths[452], deaths[453], deaths[454], deaths[455], deaths[456], deaths[457], deaths[458], deaths[459], deaths[460], deaths[461], deaths[462], deaths[463], deaths[464], deaths[465], deaths[466], deaths[467], deaths[468], deaths[469], deaths[470], deaths[471], deaths[472], deaths[473], deaths[474], deaths[475], deaths[476], deaths[477], deaths[478], deaths[479], deaths[480], deaths[481], deaths[482], deaths[483], deaths[484], deaths[485], deaths[486], deaths[487], deaths[488], deaths[489], deaths[490], deaths[491], deaths[492], deaths[493], deaths[494], deaths[495], deaths[496], deaths[497], deaths[498], deaths[499], deaths[500], deaths[501], deaths[502], deaths[503], deaths[504], deaths[505], deaths[506], deaths[507], deaths[508], deaths[509], deaths[510], deaths[511], deaths[512], deaths[513], deaths[514], deaths[515], deaths[516], deaths[517], deaths[518], deaths[519], deaths[520], deaths[521], deaths[522], deaths[523], deaths[524], deaths[525], deaths[526], deaths[527], deaths[528], deaths[529], deaths[530], deaths[531], deaths[532], deaths[533], deaths[534], deaths[535], deaths[536], deaths[537], deaths[538], deaths[539], deaths[540]
internals         = 

Use `describe(chains)` for summary statistics and quantiles.