Hierarchical Bayesian models
Parametric Bayesian framework
Currently, serosv only has models under parametric
Bayesian framework
Proposed approach
Prevalence has a parametric form π(ai, α) where α is a parameter vector
One can constraint the parameter space of the prior distribution P(α) in order to achieve the desired monotonicity of the posterior distribution P(π1, π2, ..., πm|y, n)
Where:
- n = (n1, n2, ..., nm) and ni is the sample size at age ai
- y = (y1, y2, ..., ym) and yi is the number of infected individual from the ni sampled subjects
Farrington
Refer to Chapter 10.3.1
Proposed model
The model for prevalence is as followed
$$ \pi (a) = 1 - exp\{ \frac{\alpha_1}{\alpha_2}ae^{-\alpha_2 a} + \frac{1}{\alpha_2}(\frac{\alpha_1}{\alpha_2} - \alpha_3)(e^{-\alpha_2 a} - 1) -\alpha_3 a \} $$
For likelihood model, independent binomial distribution are assumed for the number of infected individuals at age ai
yi ∼ Bin(ni, πi), for i = 1, 2, 3, ...m
The constraint on the parameter space can be incorporated by assuming truncated normal distribution for the components of α, α = (α1, α2, α3) in πi = π(ai, α)
αj ∼ truncated 𝒩(μj, τj), j = 1, 2, 3
The joint posterior distribution for α can be derived by combining the likelihood and prior as followed
$$ P(\alpha|y) \propto \prod^m_{i=1} \text{Bin}(y_i|n_i, \pi(a_i, \alpha)) \prod^3_{i=1}-\frac{1}{\tau_j}\text{exp}(\frac{1}{2\tau^2_j} (\alpha_j - \mu_j)^2) $$
Where the flat hyperprior distribution is defined as followed:
μj ∼ 𝒩(0, 10000)
τj−2 ∼ Γ(100, 100)
The full conditional distribution of αi is thus $$ P(\alpha_i|\alpha_j,\alpha_k, k, j \neq i) \propto -\frac{1}{\tau_i}\text{exp}(\frac{1}{2\tau^2_i} (\alpha_i - \mu_i)^2) \prod^m_{i=1} \text{Bin}(y_i|n_i, \pi(a_i, \alpha)) $$
Fitting data
To fit Farrington model, use
hierarchical_bayesian_model() and define
type = "far2" or type = "far3" where
type = "far2"refers to Farrington model with 2 parameters (α3 = 0)type = "far3"refers to Farrington model with 3 parameters (α3 > 0)
df <- mumps_uk_1986_1987
model <- hierarchical_bayesian_model(age = df$age, pos = df$pos, tot = df$tot, type="far3")
#>
#> SAMPLING FOR MODEL 'fra_3' NOW (CHAIN 1).
#> Chain 1: Rejecting initial value:
#> Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
#> Chain 1: Stan can't start sampling from this initial value.
#> Chain 1: Rejecting initial value:
#> Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
#> Chain 1: Stan can't start sampling from this initial value.
#> Chain 1: Rejecting initial value:
#> Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
#> Chain 1: Stan can't start sampling from this initial value.
#> Chain 1: Rejecting initial value:
#> Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
#> Chain 1: Stan can't start sampling from this initial value.
#> Chain 1: Rejecting initial value:
#> Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
#> Chain 1: Stan can't start sampling from this initial value.
#> Chain 1:
#> Chain 1: Gradient evaluation took 3.8e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.38 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 5000 [ 0%] (Warmup)
#> Chain 1: Iteration: 500 / 5000 [ 10%] (Warmup)
#> Chain 1: Iteration: 1000 / 5000 [ 20%] (Warmup)
#> Chain 1: Iteration: 1500 / 5000 [ 30%] (Warmup)
#> Chain 1: Iteration: 1501 / 5000 [ 30%] (Sampling)
#> Chain 1: Iteration: 2000 / 5000 [ 40%] (Sampling)
#> Chain 1: Iteration: 2500 / 5000 [ 50%] (Sampling)
#> Chain 1: Iteration: 3000 / 5000 [ 60%] (Sampling)
#> Chain 1: Iteration: 3500 / 5000 [ 70%] (Sampling)
#> Chain 1: Iteration: 4000 / 5000 [ 80%] (Sampling)
#> Chain 1: Iteration: 4500 / 5000 [ 90%] (Sampling)
#> Chain 1: Iteration: 5000 / 5000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 5.278 seconds (Warm-up)
#> Chain 1: 8.816 seconds (Sampling)
#> Chain 1: 14.094 seconds (Total)
#> Chain 1:
#> Warning: There were 1073 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: The largest R-hat is 1.26, indicating chains have not mixed.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#r-hat
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess
model$info
#> mean se_mean sd 2.5%
#> alpha1 1.394218e-01 7.184037e-04 5.370286e-03 1.302382e-01
#> alpha2 1.968698e-01 2.807780e-03 9.222811e-03 1.860160e-01
#> alpha3 6.810222e-03 1.940414e-03 7.288032e-03 3.211936e-04
#> tau_alpha1 2.660725e+00 1.707195e+00 7.244825e+00 6.015233e-06
#> tau_alpha2 8.709697e+00 6.243108e+00 1.429015e+01 5.594306e-06
#> tau_alpha3 5.438431e-01 2.934593e-01 1.298792e+00 9.791594e-06
#> mu_alpha1 8.465124e-01 2.759801e+00 3.146366e+01 -6.832796e+01
#> mu_alpha2 2.935836e+00 1.803454e+00 3.401768e+01 -5.828393e+01
#> mu_alpha3 1.987586e+00 8.635547e-01 2.807692e+01 -6.143534e+01
#> sigma_alpha1 1.161006e+02 4.388761e+01 1.389027e+03 1.887200e-01
#> sigma_alpha2 8.447457e+01 3.030156e+01 7.386139e+02 1.459998e-01
#> sigma_alpha3 4.514771e+01 7.045998e+00 2.357840e+02 4.568998e-01
#> lp__ -2.533498e+03 7.364915e-01 3.998086e+00 -2.542241e+03
#> 25% 50% 75% 97.5%
#> alpha1 1.361457e-01 1.375456e-01 1.435746e-01 1.510203e-01
#> alpha2 1.879245e-01 1.953378e-01 2.037840e-01 2.153396e-01
#> alpha3 4.607922e-04 4.078439e-03 1.123689e-02 2.404078e-02
#> tau_alpha1 4.648254e-03 1.770572e-02 2.365958e-01 2.807788e+01
#> tau_alpha2 2.614562e-03 1.079358e-01 1.570379e+01 4.691321e+01
#> tau_alpha3 2.674766e-03 8.677461e-03 1.110013e-01 4.790248e+00
#> mu_alpha1 -9.907055e-01 7.690655e-01 5.085120e+00 6.040674e+01
#> mu_alpha2 -1.827423e+00 1.519843e-01 1.248924e+00 9.193064e+01
#> mu_alpha3 -1.152008e+00 1.614805e+00 5.268640e+00 6.028124e+01
#> sigma_alpha1 2.055882e+00 7.515247e+00 1.466747e+01 4.077863e+02
#> sigma_alpha2 2.523468e-01 3.043809e+00 1.955693e+01 4.228959e+02
#> sigma_alpha3 3.001484e+00 1.073504e+01 1.933628e+01 3.195755e+02
#> lp__ -2.536046e+03 -2.532423e+03 -2.531240e+03 -2.526086e+03
#> n_eff Rhat
#> alpha1 55.880184 1.0618183
#> alpha2 10.789477 1.1607525
#> alpha3 14.106909 1.1302503
#> tau_alpha1 18.009003 1.1025413
#> tau_alpha2 5.239286 1.2809604
#> tau_alpha3 19.587688 1.0885712
#> mu_alpha1 129.975932 1.0018582
#> mu_alpha2 355.794542 1.0005630
#> mu_alpha3 1057.107638 1.0033144
#> sigma_alpha1 1001.700049 0.9997436
#> sigma_alpha2 594.162183 1.0016855
#> sigma_alpha3 1119.808606 1.0013388
#> lp__ 29.469229 0.9997289
plot(model)
#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.
Log-logistic
Proposed approach
The model for seroprevalence is as followed
$$ \pi(a) = \frac{\beta a^\alpha}{1 + \beta a^\alpha}, \text{ } \alpha, \beta > 0 $$
The likelihood is specified to be the same as Farrington model (yi ∼ Bin(ni, πi)) with
logit(π(a)) = α2 + α1log (a)
- Where α2 = log(β)
The prior model of α1 is specified as α1 ∼ truncated 𝒩(μ1, τ1) with flat hyperprior as in Farrington model
β is constrained to be positive by specifying α2 ∼ 𝒩(μ2, τ2)
The full conditional distribution of α1 is thus
$$ P(\alpha_1|\alpha_2) \propto -\frac{1}{\tau_1} \text{exp} (\frac{1}{2 \tau_1^2} (\alpha_1 - \mu_1)^2) \prod_{i=1}^m \text{Bin}(y_i|n_i,\pi(a_i, \alpha_1, \alpha_2) ) $$
And α2 can be derived in the same way
Fitting data
To fit Log-logistic model, use
hierarchical_bayesian_model() and define
type = "log_logistic"
df <- rubella_uk_1986_1987
model <- hierarchical_bayesian_model(age = df$age, pos = df$pos, tot = df$tot, type="log_logistic")
#>
#> SAMPLING FOR MODEL 'log_logistic' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 1.1e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.11 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 5000 [ 0%] (Warmup)
#> Chain 1: Iteration: 500 / 5000 [ 10%] (Warmup)
#> Chain 1: Iteration: 1000 / 5000 [ 20%] (Warmup)
#> Chain 1: Iteration: 1500 / 5000 [ 30%] (Warmup)
#> Chain 1: Iteration: 1501 / 5000 [ 30%] (Sampling)
#> Chain 1: Iteration: 2000 / 5000 [ 40%] (Sampling)
#> Chain 1: Iteration: 2500 / 5000 [ 50%] (Sampling)
#> Chain 1: Iteration: 3000 / 5000 [ 60%] (Sampling)
#> Chain 1: Iteration: 3500 / 5000 [ 70%] (Sampling)
#> Chain 1: Iteration: 4000 / 5000 [ 80%] (Sampling)
#> Chain 1: Iteration: 4500 / 5000 [ 90%] (Sampling)
#> Chain 1: Iteration: 5000 / 5000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 0.648 seconds (Warm-up)
#> Chain 1: 0.877 seconds (Sampling)
#> Chain 1: 1.525 seconds (Total)
#> Chain 1:
#> Warning: There were 1527 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: The largest R-hat is 2.02, indicating chains have not mixed.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#r-hat
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess
model$type
#> [1] "log_logistic"
plot(model)
#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.