using MortalityTables
using Turing
using UUIDs
using DataFramesMeta
using MCMCChains, Plots, StatsPlots
using LinearAlgebra
using Pipe
using StatsFuns
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.+inforce_years.-1] .* rf)
endpoint = if sum(deaths) == 0
last(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
Row | issue_age | risklevel | att_age | death | id |
---|---|---|---|---|---|
Int64 | Int64 | Int64 | Bool | Base.UUID | |
1 | 50 | 2 | 50 | false | 95a1a3d0-4e5f-11f0-182e-97959ef3003d |
2 | 50 | 2 | 51 | false | 95a1a3d0-4e5f-11f0-182e-97959ef3003d |
3 | 50 | 2 | 52 | false | 95a1a3d0-4e5f-11f0-182e-97959ef3003d |
4 | 50 | 2 | 53 | false | 95a1a3d0-4e5f-11f0-182e-97959ef3003d |
5 | 50 | 2 | 54 | false | 95a1a3d0-4e5f-11f0-182e-97959ef3003d |
6 | 50 | 2 | 55 | false | 95a1a3d0-4e5f-11f0-182e-97959ef3003d |
7 | 50 | 2 | 56 | false | 95a1a3d0-4e5f-11f0-182e-97959ef3003d |
8 | 50 | 2 | 57 | false | 95a1a3d0-4e5f-11f0-182e-97959ef3003d |
9 | 50 | 2 | 58 | false | 95a1a3d0-4e5f-11f0-182e-97959ef3003d |
10 | 50 | 2 | 59 | false | 95a1a3d0-4e5f-11f0-182e-97959ef3003d |
11 | 50 | 2 | 60 | false | 95a1a3d0-4e5f-11f0-182e-97959ef3003d |
12 | 50 | 3 | 50 | false | 95a575b4-4e5f-11f0-032d-effd0f5f8a93 |
13 | 50 | 3 | 51 | false | 95a575b4-4e5f-11f0-032d-effd0f5f8a93 |
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
110404 | 51 | 2 | 66 | false | 95a81b8e-4e5f-11f0-2d98-bf8b7a139ca5 |
110405 | 51 | 2 | 67 | false | 95a81b8e-4e5f-11f0-2d98-bf8b7a139ca5 |
110406 | 51 | 2 | 68 | false | 95a81b8e-4e5f-11f0-2d98-bf8b7a139ca5 |
110407 | 51 | 2 | 69 | false | 95a81b8e-4e5f-11f0-2d98-bf8b7a139ca5 |
110408 | 51 | 2 | 70 | false | 95a81b8e-4e5f-11f0-2d98-bf8b7a139ca5 |
110409 | 54 | 1 | 54 | false | 95a81ba0-4e5f-11f0-1a2f-2f3072a3c1da |
110410 | 54 | 1 | 55 | false | 95a81ba0-4e5f-11f0-1a2f-2f3072a3c1da |
110411 | 54 | 1 | 56 | false | 95a81ba0-4e5f-11f0-1a2f-2f3072a3c1da |
110412 | 54 | 1 | 57 | false | 95a81ba0-4e5f-11f0-1a2f-2f3072a3c1da |
110413 | 54 | 1 | 58 | false | 95a81ba0-4e5f-11f0-1a2f-2f3072a3c1da |
110414 | 54 | 1 | 59 | false | 95a81ba0-4e5f-11f0-1a2f-2f3072a3c1da |
110415 | 54 | 1 | 60 | false | 95a81ba0-4e5f-11f0-1a2f-2f3072a3c1da |
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
Row | issue_age | att_age | exposures | deaths | fraction |
---|---|---|---|---|---|
Int64 | Int64 | Int64 | Int64 | Float64 | |
1 | 50 | 50 | 1593 | 29 | 0.0182046 |
2 | 50 | 51 | 1511 | 21 | 0.0138981 |
3 | 50 | 52 | 1438 | 23 | 0.0159944 |
4 | 50 | 53 | 1357 | 33 | 0.0243183 |
5 | 50 | 54 | 1280 | 25 | 0.0195312 |
6 | 50 | 55 | 1213 | 28 | 0.0230833 |
7 | 50 | 56 | 1132 | 15 | 0.0132509 |
8 | 50 | 57 | 1070 | 17 | 0.0158879 |
9 | 50 | 58 | 997 | 22 | 0.0220662 |
10 | 50 | 59 | 923 | 31 | 0.0335861 |
11 | 50 | 60 | 836 | 21 | 0.0251196 |
12 | 50 | 61 | 762 | 22 | 0.0288714 |
13 | 50 | 62 | 695 | 15 | 0.0215827 |
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
169 | 55 | 73 | 226 | 17 | 0.0752212 |
170 | 55 | 74 | 174 | 6 | 0.0344828 |
171 | 55 | 75 | 137 | 7 | 0.0510949 |
172 | 55 | 76 | 107 | 7 | 0.0654206 |
173 | 55 | 77 | 75 | 4 | 0.0533333 |
174 | 55 | 78 | 57 | 1 | 0.0175439 |
175 | 55 | 79 | 42 | 2 | 0.047619 |
176 | 55 | 80 | 34 | 2 | 0.0588235 |
177 | 55 | 81 | 17 | 2 | 0.117647 |
178 | 55 | 82 | 11 | 1 | 0.0909091 |
179 | 55 | 83 | 7 | 1 | 0.142857 |
180 | 55 | 84 | 3 | 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
Row | issue_age | att_age | risklevel | exposures | deaths | fraction |
---|---|---|---|---|---|---|
Int64 | Int64 | Int64 | Int64 | Int64 | Float64 | |
1 | 50 | 50 | 1 | 511 | 4 | 0.00782779 |
2 | 50 | 50 | 2 | 537 | 8 | 0.0148976 |
3 | 50 | 50 | 3 | 545 | 17 | 0.0311927 |
4 | 50 | 51 | 1 | 489 | 5 | 0.0102249 |
5 | 50 | 51 | 2 | 513 | 7 | 0.0136452 |
6 | 50 | 51 | 3 | 509 | 9 | 0.0176817 |
7 | 50 | 52 | 1 | 472 | 2 | 0.00423729 |
8 | 50 | 52 | 2 | 487 | 10 | 0.0205339 |
9 | 50 | 52 | 3 | 479 | 11 | 0.0229645 |
10 | 50 | 53 | 1 | 454 | 6 | 0.0132159 |
11 | 50 | 53 | 2 | 454 | 12 | 0.0264317 |
12 | 50 | 53 | 3 | 449 | 15 | 0.0334076 |
13 | 50 | 54 | 1 | 438 | 3 | 0.00684932 |
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
522 | 55 | 80 | 2 | 9 | 0 | 0.0 |
523 | 55 | 80 | 3 | 7 | 1 | 0.142857 |
524 | 55 | 81 | 1 | 10 | 1 | 0.1 |
525 | 55 | 81 | 2 | 4 | 1 | 0.25 |
526 | 55 | 81 | 3 | 3 | 0 | 0.0 |
527 | 55 | 82 | 1 | 7 | 0 | 0.0 |
528 | 55 | 82 | 2 | 2 | 1 | 0.5 |
529 | 55 | 82 | 3 | 2 | 0 | 0.0 |
530 | 55 | 83 | 1 | 6 | 1 | 0.166667 |
531 | 55 | 83 | 3 | 1 | 0 | 0.0 |
532 | 55 | 84 | 1 | 2 | 0 | 0.0 |
533 | 55 | 84 | 3 | 1 | 0 | 0.0 |
1: A single binomial parameter model
Estimate
@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 = sample(m1, NUTS(), 1000)
plot(chain)
Sampling 0%| | ETA: N/A ┌ Info: Found initial step size └ ϵ = 0.0125 Sampling 0%|▎ | ETA: 0:09:18 Sampling 1%|▍ | ETA: 0:04:39 Sampling 1%|▋ | ETA: 0:03:05 Sampling 2%|▊ | ETA: 0:02:18 Sampling 2%|█ | ETA: 0:01:50 Sampling 3%|█▏ | ETA: 0:01:31 Sampling 3%|█▍ | ETA: 0:01:18 Sampling 4%|█▋ | ETA: 0:01:08 Sampling 4%|█▊ | ETA: 0:01:00 Sampling 5%|██ | ETA: 0:00:54 Sampling 5%|██▏ | ETA: 0:00:49 Sampling 6%|██▍ | ETA: 0:00:44 Sampling 6%|██▌ | ETA: 0:00:41 Sampling 7%|██▊ | ETA: 0:00:38 Sampling 7%|███ | ETA: 0:00:35 Sampling 7%|███▏ | ETA: 0:00:33 Sampling 8%|███▍ | ETA: 0:00:31 Sampling 8%|███▌ | ETA: 0:00:29 Sampling 9%|███▊ | ETA: 0:00:27 Sampling 9%|███▉ | ETA: 0:00:26 Sampling 10%|████▏ | ETA: 0:00:24 Sampling 10%|████▎ | ETA: 0:00:23 Sampling 11%|████▌ | ETA: 0:00:22 Sampling 11%|████▊ | ETA: 0:00:21 Sampling 12%|████▉ | ETA: 0:00:20 Sampling 12%|█████▏ | ETA: 0:00:19 Sampling 13%|█████▎ | ETA: 0:00:18 Sampling 13%|█████▌ | ETA: 0:00:18 Sampling 14%|█████▋ | ETA: 0:00:17 Sampling 14%|█████▉ | ETA: 0:00:16 Sampling 14%|██████▏ | ETA: 0:00:16 Sampling 15%|██████▎ | ETA: 0:00:15 Sampling 15%|██████▌ | ETA: 0:00:15 Sampling 16%|██████▋ | ETA: 0:00:14 Sampling 16%|██████▉ | ETA: 0:00:14 Sampling 17%|███████ | ETA: 0:00:13 Sampling 17%|███████▎ | ETA: 0:00:13 Sampling 18%|███████▌ | ETA: 0:00:12 Sampling 18%|███████▋ | ETA: 0:00:12 Sampling 19%|███████▉ | ETA: 0:00:12 Sampling 19%|████████ | ETA: 0:00:11 Sampling 20%|████████▎ | ETA: 0:00:11 Sampling 20%|████████▍ | ETA: 0:00:11 Sampling 21%|████████▋ | ETA: 0:00:10 Sampling 21%|████████▉ | ETA: 0:00:10 Sampling 21%|█████████ | ETA: 0:00:10 Sampling 22%|█████████▎ | ETA: 0:00:10 Sampling 22%|█████████▍ | ETA: 0:00:09 Sampling 23%|█████████▋ | ETA: 0:00:09 Sampling 23%|█████████▊ | ETA: 0:00:09 Sampling 24%|██████████ | ETA: 0:00:09 Sampling 24%|██████████▎ | ETA: 0:00:08 Sampling 25%|██████████▍ | ETA: 0:00:08 Sampling 25%|██████████▋ | ETA: 0:00:08 Sampling 26%|██████████▊ | ETA: 0:00:08 Sampling 26%|███████████ | ETA: 0:00:08 Sampling 27%|███████████▏ | ETA: 0:00:07 Sampling 27%|███████████▍ | ETA: 0:00:07 Sampling 28%|███████████▋ | ETA: 0:00:07 Sampling 28%|███████████▊ | ETA: 0:00:07 Sampling 28%|████████████ | ETA: 0:00:07 Sampling 29%|████████████▏ | ETA: 0:00:07 Sampling 29%|████████████▍ | ETA: 0:00:06 Sampling 30%|████████████▌ | ETA: 0:00:06 Sampling 30%|████████████▊ | ETA: 0:00:06 Sampling 31%|████████████▉ | ETA: 0:00:06 Sampling 31%|█████████████▏ | ETA: 0:00:06 Sampling 32%|█████████████▍ | ETA: 0:00:06 Sampling 32%|█████████████▌ | ETA: 0:00:06 Sampling 33%|█████████████▊ | ETA: 0:00:06 Sampling 33%|█████████████▉ | ETA: 0:00:05 Sampling 34%|██████████████▏ | ETA: 0:00:05 Sampling 34%|██████████████▎ | ETA: 0:00:05 Sampling 35%|██████████████▌ | ETA: 0:00:05 Sampling 35%|██████████████▊ | ETA: 0:00:05 Sampling 35%|██████████████▉ | ETA: 0:00:05 Sampling 36%|███████████████▏ | ETA: 0:00:05 Sampling 36%|███████████████▎ | ETA: 0:00:05 Sampling 37%|███████████████▌ | ETA: 0:00:05 Sampling 37%|███████████████▋ | ETA: 0:00:05 Sampling 38%|███████████████▉ | ETA: 0:00:05 Sampling 38%|████████████████▏ | ETA: 0:00:04 Sampling 39%|████████████████▎ | ETA: 0:00:04 Sampling 39%|████████████████▌ | ETA: 0:00:04 Sampling 40%|████████████████▋ | ETA: 0:00:04 Sampling 40%|████████████████▉ | ETA: 0:00:04 Sampling 41%|█████████████████ | ETA: 0:00:04 Sampling 41%|█████████████████▎ | ETA: 0:00:04 Sampling 42%|█████████████████▌ | ETA: 0:00:04 Sampling 42%|█████████████████▋ | ETA: 0:00:04 Sampling 42%|█████████████████▉ | ETA: 0:00:04 Sampling 43%|██████████████████ | ETA: 0:00:04 Sampling 43%|██████████████████▎ | ETA: 0:00:04 Sampling 44%|██████████████████▍ | ETA: 0:00:04 Sampling 44%|██████████████████▋ | ETA: 0:00:03 Sampling 45%|██████████████████▉ | ETA: 0:00:03 Sampling 45%|███████████████████ | ETA: 0:00:03 Sampling 46%|███████████████████▎ | ETA: 0:00:03 Sampling 46%|███████████████████▍ | ETA: 0:00:03 Sampling 47%|███████████████████▋ | ETA: 0:00:03 Sampling 47%|███████████████████▊ | ETA: 0:00:03 Sampling 48%|████████████████████ | ETA: 0:00:03 Sampling 48%|████████████████████▎ | ETA: 0:00:03 Sampling 49%|████████████████████▍ | ETA: 0:00:03 Sampling 49%|████████████████████▋ | ETA: 0:00:03 Sampling 49%|████████████████████▊ | ETA: 0:00:03 Sampling 50%|█████████████████████ | ETA: 0:00:03 Sampling 50%|█████████████████████▏ | ETA: 0:00:03 Sampling 51%|█████████████████████▍ | ETA: 0:00:03 Sampling 51%|█████████████████████▌ | ETA: 0:00:03 Sampling 52%|█████████████████████▊ | ETA: 0:00:03 Sampling 52%|██████████████████████ | ETA: 0:00:03 Sampling 53%|██████████████████████▏ | ETA: 0:00:02 Sampling 53%|██████████████████████▍ | ETA: 0:00:02 Sampling 54%|██████████████████████▌ | ETA: 0:00:02 Sampling 54%|██████████████████████▊ | ETA: 0:00:02 Sampling 55%|██████████████████████▉ | ETA: 0:00:02 Sampling 55%|███████████████████████▏ | ETA: 0:00:02 Sampling 56%|███████████████████████▍ | ETA: 0:00:02 Sampling 56%|███████████████████████▌ | ETA: 0:00:02 Sampling 56%|███████████████████████▊ | ETA: 0:00:02 Sampling 57%|███████████████████████▉ | ETA: 0:00:02 Sampling 57%|████████████████████████▏ | ETA: 0:00:02 Sampling 58%|████████████████████████▎ | ETA: 0:00:02 Sampling 58%|████████████████████████▌ | ETA: 0:00:02 Sampling 59%|████████████████████████▊ | ETA: 0:00:02 Sampling 59%|████████████████████████▉ | ETA: 0:00:02 Sampling 60%|█████████████████████████▏ | ETA: 0:00:02 Sampling 60%|█████████████████████████▎ | ETA: 0:00:02 Sampling 61%|█████████████████████████▌ | ETA: 0:00:02 Sampling 61%|█████████████████████████▋ | ETA: 0:00:02 Sampling 62%|█████████████████████████▉ | ETA: 0:00:02 Sampling 62%|██████████████████████████▏ | ETA: 0:00:02 Sampling 63%|██████████████████████████▎ | ETA: 0:00:02 Sampling 63%|██████████████████████████▌ | ETA: 0:00:02 Sampling 63%|██████████████████████████▋ | ETA: 0:00:02 Sampling 64%|██████████████████████████▉ | ETA: 0:00:02 Sampling 64%|███████████████████████████ | ETA: 0:00:02 Sampling 65%|███████████████████████████▎ | ETA: 0:00:02 Sampling 65%|███████████████████████████▌ | ETA: 0:00:01 Sampling 66%|███████████████████████████▋ | ETA: 0:00:01 Sampling 66%|███████████████████████████▉ | ETA: 0:00:01 Sampling 67%|████████████████████████████ | ETA: 0:00:01 Sampling 67%|████████████████████████████▎ | ETA: 0:00:01 Sampling 68%|████████████████████████████▍ | ETA: 0:00:01 Sampling 68%|████████████████████████████▋ | ETA: 0:00:01 Sampling 69%|████████████████████████████▊ | ETA: 0:00:01 Sampling 69%|█████████████████████████████ | ETA: 0:00:01 Sampling 70%|█████████████████████████████▎ | ETA: 0:00:01 Sampling 70%|█████████████████████████████▍ | ETA: 0:00:01 Sampling 70%|█████████████████████████████▋ | ETA: 0:00:01 Sampling 71%|█████████████████████████████▊ | ETA: 0:00:01 Sampling 71%|██████████████████████████████ | ETA: 0:00:01 Sampling 72%|██████████████████████████████▏ | ETA: 0:00:01 Sampling 72%|██████████████████████████████▍ | ETA: 0:00:01 Sampling 73%|██████████████████████████████▋ | ETA: 0:00:01 Sampling 73%|██████████████████████████████▊ | ETA: 0:00:01 Sampling 74%|███████████████████████████████ | ETA: 0:00:01 Sampling 74%|███████████████████████████████▏ | ETA: 0:00:01 Sampling 75%|███████████████████████████████▍ | ETA: 0:00:01 Sampling 75%|███████████████████████████████▌ | ETA: 0:00:01 Sampling 76%|███████████████████████████████▊ | ETA: 0:00:01 Sampling 76%|████████████████████████████████ | ETA: 0:00:01 Sampling 77%|████████████████████████████████▏ | ETA: 0:00:01 Sampling 77%|████████████████████████████████▍ | ETA: 0:00:01 Sampling 77%|████████████████████████████████▌ | ETA: 0:00:01 Sampling 78%|████████████████████████████████▊ | ETA: 0:00:01 Sampling 78%|████████████████████████████████▉ | ETA: 0:00:01 Sampling 79%|█████████████████████████████████▏ | ETA: 0:00:01 Sampling 79%|█████████████████████████████████▍ | ETA: 0:00:01 Sampling 80%|█████████████████████████████████▌ | ETA: 0:00:01 Sampling 80%|█████████████████████████████████▊ | ETA: 0:00:01 Sampling 81%|█████████████████████████████████▉ | ETA: 0:00:01 Sampling 81%|██████████████████████████████████▏ | ETA: 0:00:01 Sampling 82%|██████████████████████████████████▎ | ETA: 0:00:01 Sampling 82%|██████████████████████████████████▌ | ETA: 0:00:01 Sampling 83%|██████████████████████████████████▊ | ETA: 0:00:01 Sampling 83%|██████████████████████████████████▉ | ETA: 0:00:01 Sampling 84%|███████████████████████████████████▏ | ETA: 0:00:01 Sampling 84%|███████████████████████████████████▎ | ETA: 0:00:01 Sampling 84%|███████████████████████████████████▌ | ETA: 0:00:01 Sampling 85%|███████████████████████████████████▋ | ETA: 0:00:01 Sampling 85%|███████████████████████████████████▉ | ETA: 0:00:00 Sampling 86%|████████████████████████████████████▏ | ETA: 0:00:00 Sampling 86%|████████████████████████████████████▎ | ETA: 0:00:00 Sampling 87%|████████████████████████████████████▌ | ETA: 0:00:00 Sampling 87%|████████████████████████████████████▋ | ETA: 0:00:00 Sampling 88%|████████████████████████████████████▉ | ETA: 0:00:00 Sampling 88%|█████████████████████████████████████ | ETA: 0:00:00 Sampling 89%|█████████████████████████████████████▎ | ETA: 0:00:00 Sampling 89%|█████████████████████████████████████▍ | ETA: 0:00:00 Sampling 90%|█████████████████████████████████████▋ | ETA: 0:00:00 Sampling 90%|█████████████████████████████████████▉ | ETA: 0:00:00 Sampling 91%|██████████████████████████████████████ | ETA: 0:00:00 Sampling 91%|██████████████████████████████████████▎ | ETA: 0:00:00 Sampling 91%|██████████████████████████████████████▍ | ETA: 0:00:00 Sampling 92%|██████████████████████████████████████▋ | ETA: 0:00:00 Sampling 92%|██████████████████████████████████████▊ | ETA: 0:00:00 Sampling 93%|███████████████████████████████████████ | ETA: 0:00:00 Sampling 93%|███████████████████████████████████████▎ | ETA: 0:00:00 Sampling 94%|███████████████████████████████████████▍ | ETA: 0:00:00 Sampling 94%|███████████████████████████████████████▋ | ETA: 0:00:00 Sampling 95%|███████████████████████████████████████▊ | ETA: 0:00:00 Sampling 95%|████████████████████████████████████████ | ETA: 0:00:00 Sampling 96%|████████████████████████████████████████▏ | ETA: 0:00:00 Sampling 96%|████████████████████████████████████████▍ | ETA: 0:00:00 Sampling 97%|████████████████████████████████████████▋ | ETA: 0:00:00 Sampling 97%|████████████████████████████████████████▊ | ETA: 0:00:00 Sampling 98%|█████████████████████████████████████████ | ETA: 0:00:00 Sampling 98%|█████████████████████████████████████████▏| ETA: 0:00:00 Sampling 98%|█████████████████████████████████████████▍| ETA: 0:00:00 Sampling 99%|█████████████████████████████████████████▌| ETA: 0:00:00 Sampling 99%|█████████████████████████████████████████▊| ETA: 0:00:00 Sampling 100%|██████████████████████████████████████████| ETA: 0:00:00 Sampling 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)
p = scatter(
data.att_age,
data.fraction,
markersize=data_weight,
alpha=0.5,
label="Experience data point (size indicates relative exposure quantity)",
xlabel="age",
ylim=(0.0, 0.25),
ylabel="mortality rate",
title="Parametric Bayseian Mortality"
)
# show n samples from the posterior plotted on the graph
n = 300
ages = sort!(unique(data.att_age))
for i in 1:n
p_posterior = sample(chain, 1)[:p][1]
hline!([p_posterior], label="", alpha=0.1)
end
p
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.029440898171244314, 0.02946157677851741)
2. Parametric model
In this example, we utilize a MakehamBeard parameterization because it’s already very similar in form to a logistic function. This is important because our desired output is a probability (ie the probablity of a death at a given age), so the value must be constrained to be in the interval between zero and one.
The prior values for a
,b
,c
, and k
are chosen to constrain the hazard (mortality) rate to be between zero and 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 = MortalityTables.hazard(m, age)
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 = sample(m2, NUTS(), 1000)
summarize(chain2)
Sampling 0%| | ETA: N/A ┌ Info: Found initial step size └ ϵ = 0.0125 Sampling 0%|▎ | ETA: 0:06:18 Sampling 1%|▍ | ETA: 0:03:08 Sampling 1%|▋ | ETA: 0:02:09 Sampling 2%|▊ | ETA: 0:01:40 Sampling 2%|█ | ETA: 0:01:21 Sampling 3%|█▏ | ETA: 0:01:09 Sampling 3%|█▍ | ETA: 0:01:00 Sampling 4%|█▋ | ETA: 0:00:55 Sampling 4%|█▊ | ETA: 0:00:50 Sampling 5%|██ | ETA: 0:00:45 Sampling 5%|██▏ | ETA: 0:00:42 Sampling 6%|██▍ | ETA: 0:00:39 Sampling 6%|██▌ | ETA: 0:00:36 Sampling 7%|██▊ | ETA: 0:00:34 Sampling 7%|███ | ETA: 0:00:32 Sampling 7%|███▏ | ETA: 0:00:30 Sampling 8%|███▍ | ETA: 0:00:29 Sampling 8%|███▌ | ETA: 0:00:27 Sampling 9%|███▊ | ETA: 0:00:26 Sampling 9%|███▉ | ETA: 0:00:25 Sampling 10%|████▏ | ETA: 0:00:23 Sampling 10%|████▎ | ETA: 0:00:23 Sampling 11%|████▌ | ETA: 0:00:22 Sampling 11%|████▊ | ETA: 0:00:21 Sampling 12%|████▉ | ETA: 0:00:20 Sampling 12%|█████▏ | ETA: 0:00:19 Sampling 13%|█████▎ | ETA: 0:00:19 Sampling 13%|█████▌ | ETA: 0:00:18 Sampling 14%|█████▋ | ETA: 0:00:17 Sampling 14%|█████▉ | ETA: 0:00:17 Sampling 14%|██████▏ | ETA: 0:00:16 Sampling 15%|██████▎ | ETA: 0:00:16 Sampling 15%|██████▌ | ETA: 0:00:15 Sampling 16%|██████▋ | ETA: 0:00:15 Sampling 16%|██████▉ | ETA: 0:00:15 Sampling 17%|███████ | ETA: 0:00:14 Sampling 17%|███████▎ | ETA: 0:00:14 Sampling 18%|███████▌ | ETA: 0:00:13 Sampling 18%|███████▋ | ETA: 0:00:13 Sampling 19%|███████▉ | ETA: 0:00:13 Sampling 19%|████████ | ETA: 0:00:13 Sampling 20%|████████▎ | ETA: 0:00:12 Sampling 20%|████████▍ | ETA: 0:00:12 Sampling 21%|████████▋ | ETA: 0:00:12 Sampling 21%|████████▉ | ETA: 0:00:12 Sampling 21%|█████████ | ETA: 0:00:11 Sampling 22%|█████████▎ | ETA: 0:00:11 Sampling 22%|█████████▍ | ETA: 0:00:11 Sampling 23%|█████████▋ | ETA: 0:00:11 Sampling 23%|█████████▊ | ETA: 0:00:10 Sampling 24%|██████████ | ETA: 0:00:10 Sampling 24%|██████████▎ | ETA: 0:00:10 Sampling 25%|██████████▍ | ETA: 0:00:10 Sampling 25%|██████████▋ | ETA: 0:00:09 Sampling 26%|██████████▊ | ETA: 0:00:09 Sampling 26%|███████████ | ETA: 0:00:09 Sampling 27%|███████████▏ | ETA: 0:00:09 Sampling 27%|███████████▍ | ETA: 0:00:09 Sampling 28%|███████████▋ | ETA: 0:00:09 Sampling 28%|███████████▊ | ETA: 0:00:08 Sampling 28%|████████████ | ETA: 0:00:08 Sampling 29%|████████████▏ | ETA: 0:00:08 Sampling 29%|████████████▍ | ETA: 0:00:08 Sampling 30%|████████████▌ | ETA: 0:00:08 Sampling 30%|████████████▊ | ETA: 0:00:08 Sampling 31%|████████████▉ | ETA: 0:00:08 Sampling 31%|█████████████▏ | ETA: 0:00:07 Sampling 32%|█████████████▍ | ETA: 0:00:07 Sampling 32%|█████████████▌ | ETA: 0:00:07 Sampling 33%|█████████████▊ | ETA: 0:00:07 Sampling 33%|█████████████▉ | ETA: 0:00:07 Sampling 34%|██████████████▏ | ETA: 0:00:07 Sampling 34%|██████████████▎ | ETA: 0:00:07 Sampling 35%|██████████████▌ | ETA: 0:00:07 Sampling 35%|██████████████▊ | ETA: 0:00:07 Sampling 35%|██████████████▉ | ETA: 0:00:07 Sampling 36%|███████████████▏ | ETA: 0:00:07 Sampling 36%|███████████████▎ | ETA: 0:00:06 Sampling 37%|███████████████▌ | ETA: 0:00:06 Sampling 37%|███████████████▋ | ETA: 0:00:06 Sampling 38%|███████████████▉ | ETA: 0:00:06 Sampling 38%|████████████████▏ | ETA: 0:00:06 Sampling 39%|████████████████▎ | ETA: 0:00:06 Sampling 39%|████████████████▌ | ETA: 0:00:06 Sampling 40%|████████████████▋ | ETA: 0:00:06 Sampling 40%|████████████████▉ | ETA: 0:00:06 Sampling 41%|█████████████████ | ETA: 0:00:06 Sampling 41%|█████████████████▎ | ETA: 0:00:06 Sampling 42%|█████████████████▌ | ETA: 0:00:06 Sampling 42%|█████████████████▋ | ETA: 0:00:05 Sampling 42%|█████████████████▉ | ETA: 0:00:05 Sampling 43%|██████████████████ | ETA: 0:00:05 Sampling 43%|██████████████████▎ | ETA: 0:00:05 Sampling 44%|██████████████████▍ | ETA: 0:00:05 Sampling 44%|██████████████████▋ | ETA: 0:00:05 Sampling 45%|██████████████████▉ | ETA: 0:00:05 Sampling 45%|███████████████████ | ETA: 0:00:05 Sampling 46%|███████████████████▎ | ETA: 0:00:05 Sampling 46%|███████████████████▍ | ETA: 0:00:05 Sampling 47%|███████████████████▋ | ETA: 0:00:05 Sampling 47%|███████████████████▊ | ETA: 0:00:05 Sampling 48%|████████████████████ | ETA: 0:00:05 Sampling 48%|████████████████████▎ | ETA: 0:00:05 Sampling 49%|████████████████████▍ | ETA: 0:00:05 Sampling 49%|████████████████████▋ | ETA: 0:00:04 Sampling 49%|████████████████████▊ | ETA: 0:00:04 Sampling 50%|█████████████████████ | ETA: 0:00:04 Sampling 50%|█████████████████████▏ | ETA: 0:00:04 Sampling 51%|█████████████████████▍ | ETA: 0:00:04 Sampling 51%|█████████████████████▌ | ETA: 0:00:04 Sampling 52%|█████████████████████▊ | ETA: 0:00:04 Sampling 52%|██████████████████████ | ETA: 0:00:04 Sampling 53%|██████████████████████▏ | ETA: 0:00:04 Sampling 53%|██████████████████████▍ | ETA: 0:00:04 Sampling 54%|██████████████████████▌ | ETA: 0:00:04 Sampling 54%|██████████████████████▊ | ETA: 0:00:04 Sampling 55%|██████████████████████▉ | ETA: 0:00:04 Sampling 55%|███████████████████████▏ | ETA: 0:00:04 Sampling 56%|███████████████████████▍ | ETA: 0:00:04 Sampling 56%|███████████████████████▌ | ETA: 0:00:04 Sampling 56%|███████████████████████▊ | ETA: 0:00:04 Sampling 57%|███████████████████████▉ | ETA: 0:00:04 Sampling 57%|████████████████████████▏ | ETA: 0:00:03 Sampling 58%|████████████████████████▎ | ETA: 0:00:03 Sampling 58%|████████████████████████▌ | ETA: 0:00:03 Sampling 59%|████████████████████████▊ | ETA: 0:00:03 Sampling 59%|████████████████████████▉ | ETA: 0:00:03 Sampling 60%|█████████████████████████▏ | ETA: 0:00:03 Sampling 60%|█████████████████████████▎ | ETA: 0:00:03 Sampling 61%|█████████████████████████▌ | ETA: 0:00:03 Sampling 61%|█████████████████████████▋ | ETA: 0:00:03 Sampling 62%|█████████████████████████▉ | ETA: 0:00:03 Sampling 62%|██████████████████████████▏ | ETA: 0:00:03 Sampling 63%|██████████████████████████▎ | ETA: 0:00:03 Sampling 63%|██████████████████████████▌ | ETA: 0:00:03 Sampling 63%|██████████████████████████▋ | ETA: 0:00:03 Sampling 64%|██████████████████████████▉ | ETA: 0:00:03 Sampling 64%|███████████████████████████ | ETA: 0:00:03 Sampling 65%|███████████████████████████▎ | ETA: 0:00:03 Sampling 65%|███████████████████████████▌ | ETA: 0:00:03 Sampling 66%|███████████████████████████▋ | ETA: 0:00:03 Sampling 66%|███████████████████████████▉ | ETA: 0:00:03 Sampling 67%|████████████████████████████ | ETA: 0:00:03 Sampling 67%|████████████████████████████▎ | ETA: 0:00:03 Sampling 68%|████████████████████████████▍ | ETA: 0:00:03 Sampling 68%|████████████████████████████▋ | ETA: 0:00:02 Sampling 69%|████████████████████████████▊ | ETA: 0:00:02 Sampling 69%|█████████████████████████████ | ETA: 0:00:02 Sampling 70%|█████████████████████████████▎ | ETA: 0:00:02 Sampling 70%|█████████████████████████████▍ | ETA: 0:00:02 Sampling 70%|█████████████████████████████▋ | ETA: 0:00:02 Sampling 71%|█████████████████████████████▊ | ETA: 0:00:02 Sampling 71%|██████████████████████████████ | ETA: 0:00:02 Sampling 72%|██████████████████████████████▏ | ETA: 0:00:02 Sampling 72%|██████████████████████████████▍ | ETA: 0:00:02 Sampling 73%|██████████████████████████████▋ | ETA: 0:00:02 Sampling 73%|██████████████████████████████▊ | ETA: 0:00:02 Sampling 74%|███████████████████████████████ | ETA: 0:00:02 Sampling 74%|███████████████████████████████▏ | ETA: 0:00:02 Sampling 75%|███████████████████████████████▍ | ETA: 0:00:02 Sampling 75%|███████████████████████████████▌ | ETA: 0:00:02 Sampling 76%|███████████████████████████████▊ | ETA: 0:00:02 Sampling 76%|████████████████████████████████ | ETA: 0:00:02 Sampling 77%|████████████████████████████████▏ | ETA: 0:00:02 Sampling 77%|████████████████████████████████▍ | ETA: 0:00:02 Sampling 77%|████████████████████████████████▌ | ETA: 0:00:02 Sampling 78%|████████████████████████████████▊ | ETA: 0:00:02 Sampling 78%|████████████████████████████████▉ | ETA: 0:00:02 Sampling 79%|█████████████████████████████████▏ | ETA: 0:00:02 Sampling 79%|█████████████████████████████████▍ | ETA: 0:00:02 Sampling 80%|█████████████████████████████████▌ | ETA: 0:00:01 Sampling 80%|█████████████████████████████████▊ | ETA: 0:00:01 Sampling 81%|█████████████████████████████████▉ | ETA: 0:00:01 Sampling 81%|██████████████████████████████████▏ | ETA: 0:00:01 Sampling 82%|██████████████████████████████████▎ | ETA: 0:00:01 Sampling 82%|██████████████████████████████████▌ | ETA: 0:00:01 Sampling 83%|██████████████████████████████████▊ | ETA: 0:00:01 Sampling 83%|██████████████████████████████████▉ | ETA: 0:00:01 Sampling 84%|███████████████████████████████████▏ | ETA: 0:00:01 Sampling 84%|███████████████████████████████████▎ | ETA: 0:00:01 Sampling 84%|███████████████████████████████████▌ | ETA: 0:00:01 Sampling 85%|███████████████████████████████████▋ | ETA: 0:00:01 Sampling 85%|███████████████████████████████████▉ | ETA: 0:00:01 Sampling 86%|████████████████████████████████████▏ | ETA: 0:00:01 Sampling 86%|████████████████████████████████████▎ | ETA: 0:00:01 Sampling 87%|████████████████████████████████████▌ | ETA: 0:00:01 Sampling 87%|████████████████████████████████████▋ | ETA: 0:00:01 Sampling 88%|████████████████████████████████████▉ | ETA: 0:00:01 Sampling 88%|█████████████████████████████████████ | ETA: 0:00:01 Sampling 89%|█████████████████████████████████████▎ | ETA: 0:00:01 Sampling 89%|█████████████████████████████████████▍ | ETA: 0:00:01 Sampling 90%|█████████████████████████████████████▋ | ETA: 0:00:01 Sampling 90%|█████████████████████████████████████▉ | ETA: 0:00:01 Sampling 91%|██████████████████████████████████████ | ETA: 0:00:01 Sampling 91%|██████████████████████████████████████▎ | ETA: 0:00:01 Sampling 91%|██████████████████████████████████████▍ | ETA: 0:00:01 Sampling 92%|██████████████████████████████████████▋ | ETA: 0:00:01 Sampling 92%|██████████████████████████████████████▊ | ETA: 0:00:01 Sampling 93%|███████████████████████████████████████ | ETA: 0:00:00 Sampling 93%|███████████████████████████████████████▎ | ETA: 0:00:00 Sampling 94%|███████████████████████████████████████▍ | ETA: 0:00:00 Sampling 94%|███████████████████████████████████████▋ | ETA: 0:00:00 Sampling 95%|███████████████████████████████████████▊ | ETA: 0:00:00 Sampling 95%|████████████████████████████████████████ | ETA: 0:00:00 Sampling 96%|████████████████████████████████████████▏ | ETA: 0:00:00 Sampling 96%|████████████████████████████████████████▍ | ETA: 0:00:00 Sampling 97%|████████████████████████████████████████▋ | ETA: 0:00:00 Sampling 97%|████████████████████████████████████████▊ | ETA: 0:00:00 Sampling 98%|█████████████████████████████████████████ | ETA: 0:00:00 Sampling 98%|█████████████████████████████████████████▏| ETA: 0:00:00 Sampling 98%|█████████████████████████████████████████▍| ETA: 0:00:00 Sampling 99%|█████████████████████████████████████████▌| ETA: 0:00:00 Sampling 99%|█████████████████████████████████████████▊| ETA: 0:00:00 Sampling 100%|██████████████████████████████████████████| ETA: 0:00:00 Sampling 100%|██████████████████████████████████████████| Time: 0:00:06
parameters mean std mcse ess_bulk ess_tail rhat e ⋯ Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯ a 0.0045 0.0008 0.0001 230.8300 355.1058 1.0061 ⋯ b 0.0326 0.0030 0.0002 236.6880 345.3864 1.0059 ⋯ k 2.0719 0.9777 0.0488 294.4018 261.2196 1.0044 ⋯ 1 column omitted
plot(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)
p = scatter(
data.att_age,
data.fraction,
markersize=data_weight,
alpha=0.5,
label="Experience data point (size indicates relative exposure quantity)",
xlabel="age",
ylim=(0.0, 0.25),
ylabel="mortality rate",
title="Parametric Bayseian Mortality"
)
# show n samples from the posterior plotted on the graph
n = 300
ages = sort!(unique(data.att_age))
for i 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)
plot!(ages, age -> MortalityTables.hazard(m, age), alpha=0.1, label="")
end
p
end
3. Parametric model
This model extends the prior to create a multi-level model. Each risk class (risklevel
) gets its own MakhamBeard
model. The prior for
@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 = MortalityTables.hazard(m, age)
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 = sample(m3, NUTS(), 1000)
summarize(chain3)
Sampling 0%| | ETA: N/A ┌ Info: Found initial step size └ ϵ = 0.0125 Sampling 0%|▎ | ETA: 0:06:41 Sampling 1%|▍ | ETA: 0:03:53 Sampling 1%|▋ | ETA: 0:02:59 Sampling 2%|▊ | ETA: 0:02:34 Sampling 2%|█ | ETA: 0:02:07 Sampling 3%|█▏ | ETA: 0:01:51 Sampling 3%|█▍ | ETA: 0:01:42 Sampling 4%|█▋ | ETA: 0:01:37 Sampling 4%|█▊ | ETA: 0:01:33 Sampling 5%|██ | ETA: 0:01:29 Sampling 5%|██▏ | ETA: 0:01:28 Sampling 6%|██▍ | ETA: 0:01:28 Sampling 6%|██▌ | ETA: 0:01:30 Sampling 7%|██▊ | ETA: 0:01:28 Sampling 7%|███ | ETA: 0:01:27 Sampling 7%|███▏ | ETA: 0:01:23 Sampling 8%|███▍ | ETA: 0:01:19 Sampling 8%|███▌ | ETA: 0:01:15 Sampling 9%|███▊ | ETA: 0:01:12 Sampling 9%|███▉ | ETA: 0:01:10 Sampling 10%|████▏ | ETA: 0:01:08 Sampling 10%|████▎ | ETA: 0:01:05 Sampling 11%|████▌ | ETA: 0:01:02 Sampling 11%|████▊ | ETA: 0:01:01 Sampling 12%|████▉ | ETA: 0:00:58 Sampling 12%|█████▏ | ETA: 0:00:57 Sampling 13%|█████▎ | ETA: 0:00:55 Sampling 13%|█████▌ | ETA: 0:00:54 Sampling 14%|█████▋ | ETA: 0:00:52 Sampling 14%|█████▉ | ETA: 0:00:51 Sampling 14%|██████▏ | ETA: 0:00:49 Sampling 15%|██████▎ | ETA: 0:00:48 Sampling 15%|██████▌ | ETA: 0:00:46 Sampling 16%|██████▋ | ETA: 0:00:45 Sampling 16%|██████▉ | ETA: 0:00:44 Sampling 17%|███████ | ETA: 0:00:43 Sampling 17%|███████▎ | ETA: 0:00:42 Sampling 18%|███████▌ | ETA: 0:00:41 Sampling 18%|███████▋ | ETA: 0:00:40 Sampling 19%|███████▉ | ETA: 0:00:39 Sampling 19%|████████ | ETA: 0:00:39 Sampling 20%|████████▎ | ETA: 0:00:38 Sampling 20%|████████▍ | ETA: 0:00:37 Sampling 21%|████████▋ | ETA: 0:00:36 Sampling 21%|████████▉ | ETA: 0:00:35 Sampling 21%|█████████ | ETA: 0:00:34 Sampling 22%|█████████▎ | ETA: 0:00:34 Sampling 22%|█████████▍ | ETA: 0:00:33 Sampling 23%|█████████▋ | ETA: 0:00:32 Sampling 23%|█████████▊ | ETA: 0:00:32 Sampling 24%|██████████ | ETA: 0:00:31 Sampling 24%|██████████▎ | ETA: 0:00:31 Sampling 25%|██████████▍ | ETA: 0:00:30 Sampling 25%|██████████▋ | ETA: 0:00:30 Sampling 26%|██████████▊ | ETA: 0:00:29 Sampling 26%|███████████ | ETA: 0:00:29 Sampling 27%|███████████▏ | ETA: 0:00:28 Sampling 27%|███████████▍ | ETA: 0:00:28 Sampling 28%|███████████▋ | ETA: 0:00:27 Sampling 28%|███████████▊ | ETA: 0:00:27 Sampling 28%|████████████ | ETA: 0:00:26 Sampling 29%|████████████▏ | ETA: 0:00:26 Sampling 29%|████████████▍ | ETA: 0:00:26 Sampling 30%|████████████▌ | ETA: 0:00:25 Sampling 30%|████████████▊ | ETA: 0:00:25 Sampling 31%|████████████▉ | ETA: 0:00:24 Sampling 31%|█████████████▏ | ETA: 0:00:24 Sampling 32%|█████████████▍ | ETA: 0:00:24 Sampling 32%|█████████████▌ | ETA: 0:00:23 Sampling 33%|█████████████▊ | ETA: 0:00:23 Sampling 33%|█████████████▉ | ETA: 0:00:23 Sampling 34%|██████████████▏ | ETA: 0:00:22 Sampling 34%|██████████████▎ | ETA: 0:00:22 Sampling 35%|██████████████▌ | ETA: 0:00:22 Sampling 35%|██████████████▊ | ETA: 0:00:22 Sampling 35%|██████████████▉ | ETA: 0:00:21 Sampling 36%|███████████████▏ | ETA: 0:00:21 Sampling 36%|███████████████▎ | ETA: 0:00:21 Sampling 37%|███████████████▌ | ETA: 0:00:20 Sampling 37%|███████████████▋ | ETA: 0:00:20 Sampling 38%|███████████████▉ | ETA: 0:00:20 Sampling 38%|████████████████▏ | ETA: 0:00:20 Sampling 39%|████████████████▎ | ETA: 0:00:19 Sampling 39%|████████████████▌ | ETA: 0:00:19 Sampling 40%|████████████████▋ | ETA: 0:00:19 Sampling 40%|████████████████▉ | ETA: 0:00:19 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 42%|█████████████████▉ | ETA: 0:00:17 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:16 Sampling 45%|███████████████████ | ETA: 0:00:16 Sampling 46%|███████████████████▎ | ETA: 0:00:16 Sampling 46%|███████████████████▍ | ETA: 0:00:16 Sampling 47%|███████████████████▋ | ETA: 0:00:16 Sampling 47%|███████████████████▊ | ETA: 0:00:15 Sampling 48%|████████████████████ | ETA: 0:00:15 Sampling 48%|████████████████████▎ | ETA: 0:00:15 Sampling 49%|████████████████████▍ | ETA: 0:00:15 Sampling 49%|████████████████████▋ | ETA: 0:00:15 Sampling 49%|████████████████████▊ | ETA: 0:00:15 Sampling 50%|█████████████████████ | ETA: 0:00:14 Sampling 50%|█████████████████████▏ | ETA: 0:00:14 Sampling 51%|█████████████████████▍ | ETA: 0:00:14 Sampling 51%|█████████████████████▌ | ETA: 0:00:14 Sampling 52%|█████████████████████▊ | ETA: 0:00:14 Sampling 52%|██████████████████████ | ETA: 0:00:13 Sampling 53%|██████████████████████▏ | ETA: 0:00:13 Sampling 53%|██████████████████████▍ | ETA: 0:00:13 Sampling 54%|██████████████████████▌ | ETA: 0:00:13 Sampling 54%|██████████████████████▊ | ETA: 0:00:13 Sampling 55%|██████████████████████▉ | ETA: 0:00:13 Sampling 55%|███████████████████████▏ | ETA: 0:00:12 Sampling 56%|███████████████████████▍ | ETA: 0:00:12 Sampling 56%|███████████████████████▌ | ETA: 0:00:12 Sampling 56%|███████████████████████▊ | ETA: 0:00:12 Sampling 57%|███████████████████████▉ | ETA: 0:00:12 Sampling 57%|████████████████████████▏ | ETA: 0:00:12 Sampling 58%|████████████████████████▎ | ETA: 0:00:11 Sampling 58%|████████████████████████▌ | ETA: 0:00:11 Sampling 59%|████████████████████████▊ | ETA: 0:00:11 Sampling 59%|████████████████████████▉ | ETA: 0:00:11 Sampling 60%|█████████████████████████▏ | ETA: 0:00:11 Sampling 60%|█████████████████████████▎ | ETA: 0:00:11 Sampling 61%|█████████████████████████▌ | ETA: 0:00:10 Sampling 61%|█████████████████████████▋ | ETA: 0:00:10 Sampling 62%|█████████████████████████▉ | ETA: 0:00:10 Sampling 62%|██████████████████████████▏ | ETA: 0:00:10 Sampling 63%|██████████████████████████▎ | ETA: 0:00:10 Sampling 63%|██████████████████████████▌ | ETA: 0:00:10 Sampling 63%|██████████████████████████▋ | ETA: 0:00:09 Sampling 64%|██████████████████████████▉ | ETA: 0:00:09 Sampling 64%|███████████████████████████ | ETA: 0:00:09 Sampling 65%|███████████████████████████▎ | ETA: 0:00:09 Sampling 65%|███████████████████████████▌ | ETA: 0:00:09 Sampling 66%|███████████████████████████▋ | ETA: 0:00:09 Sampling 66%|███████████████████████████▉ | ETA: 0:00:09 Sampling 67%|████████████████████████████ | ETA: 0:00:09 Sampling 67%|████████████████████████████▎ | ETA: 0:00:08 Sampling 68%|████████████████████████████▍ | ETA: 0:00:08 Sampling 68%|████████████████████████████▋ | ETA: 0:00:08 Sampling 69%|████████████████████████████▊ | ETA: 0:00:08 Sampling 69%|█████████████████████████████ | ETA: 0:00:08 Sampling 70%|█████████████████████████████▎ | ETA: 0:00:08 Sampling 70%|█████████████████████████████▍ | ETA: 0:00:08 Sampling 70%|█████████████████████████████▋ | ETA: 0:00:07 Sampling 71%|█████████████████████████████▊ | ETA: 0:00:07 Sampling 71%|██████████████████████████████ | ETA: 0:00:07 Sampling 72%|██████████████████████████████▏ | ETA: 0:00:07 Sampling 72%|██████████████████████████████▍ | ETA: 0:00:07 Sampling 73%|██████████████████████████████▋ | ETA: 0:00:07 Sampling 73%|██████████████████████████████▊ | ETA: 0:00:07 Sampling 74%|███████████████████████████████ | ETA: 0:00:06 Sampling 74%|███████████████████████████████▏ | ETA: 0:00:06 Sampling 75%|███████████████████████████████▍ | ETA: 0:00:06 Sampling 75%|███████████████████████████████▌ | ETA: 0:00:06 Sampling 76%|███████████████████████████████▊ | ETA: 0:00:06 Sampling 76%|████████████████████████████████ | ETA: 0:00:06 Sampling 77%|████████████████████████████████▏ | ETA: 0:00:06 Sampling 77%|████████████████████████████████▍ | ETA: 0:00:06 Sampling 77%|████████████████████████████████▌ | ETA: 0:00:05 Sampling 78%|████████████████████████████████▊ | ETA: 0:00:05 Sampling 78%|████████████████████████████████▉ | ETA: 0:00:05 Sampling 79%|█████████████████████████████████▏ | ETA: 0:00:05 Sampling 79%|█████████████████████████████████▍ | ETA: 0:00:05 Sampling 80%|█████████████████████████████████▌ | ETA: 0:00:05 Sampling 80%|█████████████████████████████████▊ | ETA: 0:00:05 Sampling 81%|█████████████████████████████████▉ | ETA: 0:00:05 Sampling 81%|██████████████████████████████████▏ | ETA: 0:00:05 Sampling 82%|██████████████████████████████████▎ | ETA: 0:00:04 Sampling 82%|██████████████████████████████████▌ | ETA: 0:00:04 Sampling 83%|██████████████████████████████████▊ | ETA: 0:00:04 Sampling 83%|██████████████████████████████████▉ | ETA: 0:00:04 Sampling 84%|███████████████████████████████████▏ | ETA: 0:00:04 Sampling 84%|███████████████████████████████████▎ | ETA: 0:00:04 Sampling 84%|███████████████████████████████████▌ | ETA: 0:00:04 Sampling 85%|███████████████████████████████████▋ | ETA: 0:00:04 Sampling 85%|███████████████████████████████████▉ | ETA: 0:00:03 Sampling 86%|████████████████████████████████████▏ | ETA: 0:00:03 Sampling 86%|████████████████████████████████████▎ | ETA: 0:00:03 Sampling 87%|████████████████████████████████████▌ | ETA: 0:00:03 Sampling 87%|████████████████████████████████████▋ | ETA: 0:00:03 Sampling 88%|████████████████████████████████████▉ | ETA: 0:00:03 Sampling 88%|█████████████████████████████████████ | ETA: 0:00:03 Sampling 89%|█████████████████████████████████████▎ | ETA: 0:00:03 Sampling 89%|█████████████████████████████████████▍ | ETA: 0:00:03 Sampling 90%|█████████████████████████████████████▋ | ETA: 0:00:02 Sampling 90%|█████████████████████████████████████▉ | ETA: 0:00:02 Sampling 91%|██████████████████████████████████████ | ETA: 0:00:02 Sampling 91%|██████████████████████████████████████▎ | ETA: 0:00:02 Sampling 91%|██████████████████████████████████████▍ | ETA: 0:00:02 Sampling 92%|██████████████████████████████████████▋ | ETA: 0:00:02 Sampling 92%|██████████████████████████████████████▊ | ETA: 0:00:02 Sampling 93%|███████████████████████████████████████ | ETA: 0:00:02 Sampling 93%|███████████████████████████████████████▎ | ETA: 0:00:02 Sampling 94%|███████████████████████████████████████▍ | ETA: 0:00:01 Sampling 94%|███████████████████████████████████████▋ | ETA: 0:00:01 Sampling 95%|███████████████████████████████████████▊ | ETA: 0:00:01 Sampling 95%|████████████████████████████████████████ | ETA: 0:00:01 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:00 Sampling 98%|█████████████████████████████████████████▍| ETA: 0:00:00 Sampling 99%|█████████████████████████████████████████▌| ETA: 0:00:00 Sampling 99%|█████████████████████████████████████████▊| ETA: 0:00:00 Sampling 100%|██████████████████████████████████████████| ETA: 0:00:00 Sampling 100%|██████████████████████████████████████████| Time: 0:00:22
parameters mean std mcse ess_bulk ess_tail rhat e ⋯ Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯ b 0.0378 0.0034 0.0003 148.7659 160.4862 0.9999 ⋯ ā 0.0087 0.0088 0.0004 446.5366 433.2182 1.0016 ⋯ a[1] 0.0023 0.0005 0.0000 158.2328 136.2125 1.0000 ⋯ a[2] 0.0031 0.0006 0.0000 156.3086 168.6002 1.0009 ⋯ a[3] 0.0048 0.0009 0.0001 163.6040 187.3317 1.0025 ⋯ k 2.0028 0.9763 0.0765 193.9719 176.6343 1.0001 ⋯ 1 column omitted
let data = data2
data_weight = data.exposures ./ sum(data.exposures)
data_weight = .√(data_weight ./ maximum(data_weight) .* 20)
color_i = data.risklevel
p = scatter(
data.att_age,
data.fraction,
markersize=data_weight,
alpha=0.5,
color=color_i,
label="Experience data point (size indicates relative exposure quantity)",
xlabel="age",
ylim=(0.0, 0.25),
ylabel="mortality rate",
title="Parametric Bayseian Mortality"
)
# show n samples from the posterior plotted on the graph
n = 100
ages = sort!(unique(data.att_age))
for r in 1:3
for 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)
if i == 1
plot!(ages, age -> MortalityTables.hazard(m, age), label="risk level $r", alpha=0.2, color=r)
else
plot!(ages, age -> MortalityTables.hazard(m, age), label="", alpha=0.2, color=r)
end
end
end
p
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.hazard(m, age)
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 = sample(m4, NUTS(), 1000)
summarize(chain4)
Sampling 0%| | ETA: N/A ┌ Info: Found initial step size └ ϵ = 0.0125 Sampling 0%|▎ | ETA: 0:08:35 Sampling 1%|▍ | ETA: 0:04:32 Sampling 1%|▋ | ETA: 0:03:31 Sampling 2%|▊ | ETA: 0:02:41 Sampling 2%|█ | ETA: 0:02:14 Sampling 3%|█▏ | ETA: 0:01:54 Sampling 3%|█▍ | ETA: 0:01:41 Sampling 4%|█▋ | ETA: 0:01:30 Sampling 4%|█▊ | ETA: 0:01:22 Sampling 5%|██ | ETA: 0:01:20 Sampling 5%|██▏ | ETA: 0:01:24 Sampling 6%|██▍ | ETA: 0:01:24 Sampling 6%|██▌ | ETA: 0:01:30 Sampling 7%|██▊ | ETA: 0:01:26 Sampling 7%|███ | ETA: 0:01:23 Sampling 7%|███▏ | ETA: 0:01:21 Sampling 8%|███▍ | ETA: 0:01:18 Sampling 8%|███▌ | ETA: 0:01:15 Sampling 9%|███▊ | ETA: 0:01:12 Sampling 9%|███▉ | ETA: 0:01:10 Sampling 10%|████▏ | ETA: 0:01:08 Sampling 10%|████▎ | ETA: 0:01:06 Sampling 11%|████▌ | ETA: 0:01:04 Sampling 11%|████▊ | ETA: 0:01:02 Sampling 12%|████▉ | ETA: 0:01:00 Sampling 12%|█████▏ | ETA: 0:00:59 Sampling 13%|█████▎ | ETA: 0:00:57 Sampling 13%|█████▌ | ETA: 0:00:56 Sampling 14%|█████▋ | ETA: 0:00:54 Sampling 14%|█████▉ | ETA: 0:00:53 Sampling 14%|██████▏ | ETA: 0:00:51 Sampling 15%|██████▎ | ETA: 0:00:50 Sampling 15%|██████▌ | ETA: 0:00:49 Sampling 16%|██████▋ | ETA: 0:00:48 Sampling 16%|██████▉ | ETA: 0:00:47 Sampling 17%|███████ | ETA: 0:00:46 Sampling 17%|███████▎ | ETA: 0:00:45 Sampling 18%|███████▌ | ETA: 0:00:44 Sampling 18%|███████▋ | ETA: 0:00:43 Sampling 19%|███████▉ | ETA: 0:00:42 Sampling 19%|████████ | ETA: 0:00:42 Sampling 20%|████████▎ | ETA: 0:00:41 Sampling 20%|████████▍ | ETA: 0:00:40 Sampling 21%|████████▋ | ETA: 0:00:40 Sampling 21%|████████▉ | ETA: 0:00:39 Sampling 21%|█████████ | ETA: 0:00:38 Sampling 22%|█████████▎ | ETA: 0:00:38 Sampling 22%|█████████▍ | ETA: 0:00:37 Sampling 23%|█████████▋ | ETA: 0:00:36 Sampling 23%|█████████▊ | ETA: 0:00:35 Sampling 24%|██████████ | ETA: 0:00:35 Sampling 24%|██████████▎ | ETA: 0:00:34 Sampling 25%|██████████▍ | ETA: 0:00:34 Sampling 25%|██████████▋ | ETA: 0:00:33 Sampling 26%|██████████▊ | ETA: 0:00:33 Sampling 26%|███████████ | ETA: 0:00:32 Sampling 27%|███████████▏ | ETA: 0:00:32 Sampling 27%|███████████▍ | ETA: 0:00:31 Sampling 28%|███████████▋ | ETA: 0:00:31 Sampling 28%|███████████▊ | ETA: 0:00:30 Sampling 28%|████████████ | ETA: 0:00:30 Sampling 29%|████████████▏ | ETA: 0:00:29 Sampling 29%|████████████▍ | ETA: 0:00:29 Sampling 30%|████████████▌ | ETA: 0:00:29 Sampling 30%|████████████▊ | ETA: 0:00:28 Sampling 31%|████████████▉ | ETA: 0:00:28 Sampling 31%|█████████████▏ | ETA: 0:00:28 Sampling 32%|█████████████▍ | ETA: 0:00:27 Sampling 32%|█████████████▌ | ETA: 0:00:27 Sampling 33%|█████████████▊ | ETA: 0:00:27 Sampling 33%|█████████████▉ | ETA: 0:00:26 Sampling 34%|██████████████▏ | ETA: 0:00:26 Sampling 34%|██████████████▎ | ETA: 0:00:26 Sampling 35%|██████████████▌ | ETA: 0:00:25 Sampling 35%|██████████████▊ | ETA: 0:00:25 Sampling 35%|██████████████▉ | ETA: 0:00:25 Sampling 36%|███████████████▏ | ETA: 0:00:25 Sampling 36%|███████████████▎ | ETA: 0:00:25 Sampling 37%|███████████████▌ | ETA: 0:00:24 Sampling 37%|███████████████▋ | ETA: 0:00:24 Sampling 38%|███████████████▉ | ETA: 0:00:24 Sampling 38%|████████████████▏ | ETA: 0:00:24 Sampling 39%|████████████████▎ | ETA: 0:00:23 Sampling 39%|████████████████▌ | ETA: 0:00:23 Sampling 40%|████████████████▋ | ETA: 0:00:23 Sampling 40%|████████████████▉ | ETA: 0:00:23 Sampling 41%|█████████████████ | ETA: 0:00:22 Sampling 41%|█████████████████▎ | ETA: 0:00:22 Sampling 42%|█████████████████▌ | ETA: 0:00:22 Sampling 42%|█████████████████▋ | ETA: 0:00:22 Sampling 42%|█████████████████▉ | ETA: 0:00:22 Sampling 43%|██████████████████ | ETA: 0:00:21 Sampling 43%|██████████████████▎ | ETA: 0:00:21 Sampling 44%|██████████████████▍ | ETA: 0:00:21 Sampling 44%|██████████████████▋ | ETA: 0:00:21 Sampling 45%|██████████████████▉ | ETA: 0:00:21 Sampling 45%|███████████████████ | ETA: 0:00:20 Sampling 46%|███████████████████▎ | ETA: 0:00:20 Sampling 46%|███████████████████▍ | ETA: 0:00:20 Sampling 47%|███████████████████▋ | ETA: 0:00:20 Sampling 47%|███████████████████▊ | ETA: 0:00:19 Sampling 48%|████████████████████ | ETA: 0:00:19 Sampling 48%|████████████████████▎ | ETA: 0:00:19 Sampling 49%|████████████████████▍ | ETA: 0:00:19 Sampling 49%|████████████████████▋ | ETA: 0:00:19 Sampling 49%|████████████████████▊ | ETA: 0:00:19 Sampling 50%|█████████████████████ | ETA: 0:00:18 Sampling 50%|█████████████████████▏ | ETA: 0:00:18 Sampling 51%|█████████████████████▍ | ETA: 0:00:18 Sampling 51%|█████████████████████▌ | ETA: 0:00:18 Sampling 52%|█████████████████████▊ | ETA: 0:00:17 Sampling 52%|██████████████████████ | ETA: 0:00:17 Sampling 53%|██████████████████████▏ | ETA: 0:00:17 Sampling 53%|██████████████████████▍ | ETA: 0:00:17 Sampling 54%|██████████████████████▌ | ETA: 0:00:17 Sampling 54%|██████████████████████▊ | ETA: 0:00:16 Sampling 55%|██████████████████████▉ | ETA: 0:00:16 Sampling 55%|███████████████████████▏ | ETA: 0:00:16 Sampling 56%|███████████████████████▍ | ETA: 0:00:16 Sampling 56%|███████████████████████▌ | ETA: 0:00:16 Sampling 56%|███████████████████████▊ | ETA: 0:00:15 Sampling 57%|███████████████████████▉ | ETA: 0:00:15 Sampling 57%|████████████████████████▏ | ETA: 0:00:15 Sampling 58%|████████████████████████▎ | ETA: 0:00:15 Sampling 58%|████████████████████████▌ | ETA: 0:00:15 Sampling 59%|████████████████████████▊ | ETA: 0:00:14 Sampling 59%|████████████████████████▉ | ETA: 0:00:14 Sampling 60%|█████████████████████████▏ | ETA: 0:00:14 Sampling 60%|█████████████████████████▎ | ETA: 0:00:14 Sampling 61%|█████████████████████████▌ | ETA: 0:00:14 Sampling 61%|█████████████████████████▋ | ETA: 0:00:14 Sampling 62%|█████████████████████████▉ | ETA: 0:00:13 Sampling 62%|██████████████████████████▏ | ETA: 0:00:13 Sampling 63%|██████████████████████████▎ | ETA: 0:00:13 Sampling 63%|██████████████████████████▌ | ETA: 0:00:13 Sampling 63%|██████████████████████████▋ | ETA: 0:00:13 Sampling 64%|██████████████████████████▉ | ETA: 0:00:12 Sampling 64%|███████████████████████████ | ETA: 0:00:12 Sampling 65%|███████████████████████████▎ | ETA: 0:00:12 Sampling 65%|███████████████████████████▌ | ETA: 0:00:12 Sampling 66%|███████████████████████████▋ | ETA: 0:00:12 Sampling 66%|███████████████████████████▉ | ETA: 0:00:12 Sampling 67%|████████████████████████████ | ETA: 0:00:11 Sampling 67%|████████████████████████████▎ | ETA: 0:00:11 Sampling 68%|████████████████████████████▍ | ETA: 0:00:11 Sampling 68%|████████████████████████████▋ | ETA: 0:00:11 Sampling 69%|████████████████████████████▊ | ETA: 0:00:11 Sampling 69%|█████████████████████████████ | ETA: 0:00:11 Sampling 70%|█████████████████████████████▎ | ETA: 0:00:10 Sampling 70%|█████████████████████████████▍ | ETA: 0:00:10 Sampling 70%|█████████████████████████████▋ | ETA: 0:00:10 Sampling 71%|█████████████████████████████▊ | ETA: 0:00:10 Sampling 71%|██████████████████████████████ | ETA: 0:00:10 Sampling 72%|██████████████████████████████▏ | ETA: 0:00:10 Sampling 72%|██████████████████████████████▍ | ETA: 0:00:09 Sampling 73%|██████████████████████████████▋ | ETA: 0:00:09 Sampling 73%|██████████████████████████████▊ | ETA: 0:00:09 Sampling 74%|███████████████████████████████ | ETA: 0:00:09 Sampling 74%|███████████████████████████████▏ | ETA: 0:00:09 Sampling 75%|███████████████████████████████▍ | ETA: 0:00:09 Sampling 75%|███████████████████████████████▌ | ETA: 0:00:08 Sampling 76%|███████████████████████████████▊ | ETA: 0:00:08 Sampling 76%|████████████████████████████████ | ETA: 0:00:08 Sampling 77%|████████████████████████████████▏ | ETA: 0:00:08 Sampling 77%|████████████████████████████████▍ | ETA: 0:00:08 Sampling 77%|████████████████████████████████▌ | ETA: 0:00:08 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:07 Sampling 80%|█████████████████████████████████▊ | ETA: 0:00:07 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:06 Sampling 83%|██████████████████████████████████▉ | ETA: 0:00:06 Sampling 84%|███████████████████████████████████▏ | 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:04 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 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 98%|█████████████████████████████████████████▍| ETA: 0:00:00 Sampling 99%|█████████████████████████████████████████▌| ETA: 0:00:00 Sampling 99%|█████████████████████████████████████████▊| ETA: 0:00:00 Sampling 100%|██████████████████████████████████████████| ETA: 0:00:00 Sampling 100%|██████████████████████████████████████████| Time: 0:00:32
parameters mean std mcse ess_bulk ess_tail rhat e ⋯ Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯ b 0.0498 0.0062 0.0005 172.0071 279.3448 1.0010 ⋯ ā 0.0030 0.0032 0.0002 303.4035 542.0983 1.0000 ⋯ a[1] 0.0007 0.0004 0.0000 160.0585 290.8747 1.0004 ⋯ a[2] 0.0011 0.0005 0.0000 162.1442 285.6050 1.0004 ⋯ a[3] 0.0019 0.0008 0.0001 168.3593 276.0461 1.0009 ⋯ c 0.0081 0.0025 0.0002 180.2294 274.7454 1.0005 ⋯ k 1.9699 0.9270 0.0440 354.5814 362.3854 0.9991 ⋯ 1 column omitted
PRECIS(DataFrame(chain4))
risk_factors4 = [mean(chain4[Symbol("a[$f]")]) for f in 1:3]
3-element Vector{Float64}:
0.0007411869483965575
0.0011369538569128979
0.0019479084625641859
let data = data2
data_weight = data.exposures ./ sum(data.exposures)
data_weight = .√(data_weight ./ maximum(data_weight) .* 20)
color_i = data.risklevel
p = scatter(
data.att_age,
data.fraction,
markersize=data_weight,
alpha=0.5,
color=color_i,
label="Experience data point (size indicates relative exposure quantity)",
xlabel="age",
ylim=(0.0, 0.25),
ylabel="mortality rate",
title="Parametric Bayseian Mortality"
)
# show n samples from the posterior plotted on the graph
n = 100
ages = sort!(unique(data.att_age))
for r in 1:3
for i in 1:n
s = sample(chain4, 1)
a = only(s[Symbol("a[$r]")])
b = only(s[:b])
k = only(s[:k])
c = 0
m = MortalityTables.MakehamBeard(; a, b, c, k)
if i == 1
plot!(ages, age -> MortalityTables.hazard(m, age), label="risk level $r", alpha=0.2, color=r)
else
plot!(ages, age -> MortalityTables.hazard(m, age), label="", alpha=0.2, color=r)
end
end
end
p
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×533×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]
internals =
Use `describe(chains)` for summary statistics and quantiles.