library(serosv)
Currently, serosv
only has models under parametric
Bayesian framework
Proposed approach
Prevalence has a parametric form \(\pi(a_i, \alpha)\) where \(\alpha\) is a parameter vector
One can constraint the parameter space of the prior distribution \(P(\alpha)\) in order to achieve the desired monotonicity of the posterior distribution \(P(\pi_1, \pi_2, ..., \pi_m|y,n)\)
Where:
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 \(a_i\)
\[ y_i \sim Bin(n_i, \pi_i), \text{ for } i = 1,2,3,...m \]
The constraint on the parameter space can be incorporated by assuming truncated normal distribution for the components of \(\alpha\), \(\alpha = (\alpha_1, \alpha_2, \alpha_3)\) in \(\pi_i = \pi(a_i,\alpha)\)
\[ \alpha_j \sim \text{truncated } \mathcal{N}(\mu_j, \tau_j), \text{ } j = 1,2,3 \]
The joint posterior distribution for \(\alpha\) 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:
\(\mu_j \sim \mathcal{N}(0, 10000)\)
\(\tau^{-2}_j \sim \Gamma(100,100)\)
The full conditional distribution of \(\alpha_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 (\(\alpha_3 = 0\))
type = "far3"
refers to Farrington model with 3
parameters (\(\alpha_3 >
0\))
<- mumps_uk_1986_1987
df <- hierarchical_bayesian_model(age = df$age, pos = df$pos, tot = df$tot, type="far3")
model #>
#> SAMPLING FOR MODEL 'fra_3' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 6.2e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.62 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: 1.526 seconds (Warm-up)
#> Chain 1: 0.798 seconds (Sampling)
#> Chain 1: 2.324 seconds (Total)
#> Chain 1:
#> Warning: There were 1725 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.1, 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
$info
model#> mean se_mean sd 2.5%
#> alpha1 1.401534e-01 2.109660e-04 4.274265e-03 1.305391e-01
#> alpha2 1.993244e-01 5.143642e-04 6.326270e-03 1.872606e-01
#> alpha3 7.323108e-03 9.902532e-04 6.397093e-03 5.341586e-04
#> tau_alpha1 3.977670e-02 9.654925e-03 1.324770e-01 1.078829e-05
#> tau_alpha2 5.456486e-02 2.034699e-02 1.589638e-01 2.356366e-05
#> tau_alpha3 4.083704e+00 2.159723e+00 4.344478e+00 1.883821e-05
#> mu_alpha1 3.136005e+00 2.064027e+00 3.553784e+01 -7.659224e+01
#> mu_alpha2 1.372515e+01 2.996959e+00 2.630097e+01 -4.643403e+01
#> mu_alpha3 -1.984782e-01 1.418436e+00 2.575022e+01 -6.249099e+01
#> sigma_alpha1 5.201257e+01 1.071658e+01 5.132115e+02 1.525223e+00
#> sigma_alpha2 4.007395e+01 5.558736e+00 1.561176e+02 1.259397e+00
#> sigma_alpha3 3.939620e+01 1.220862e+01 4.353595e+02 2.598392e-01
#> lp__ -2.533385e+03 1.184709e+00 3.443622e+00 -2.542190e+03
#> 25% 50% 75% 97.5%
#> alpha1 1.389235e-01 1.398818e-01 1.414889e-01 1.503498e-01
#> alpha2 1.973405e-01 1.984498e-01 1.997835e-01 2.192120e-01
#> alpha3 4.159475e-03 5.546748e-03 6.725417e-03 2.906667e-02
#> tau_alpha1 3.030719e-03 4.559242e-03 9.759346e-03 4.298828e-01
#> tau_alpha2 1.121173e-03 2.577967e-03 5.997993e-03 6.304918e-01
#> tau_alpha3 1.137083e-02 3.010172e+00 7.191923e+00 1.481121e+01
#> mu_alpha1 -6.094159e-01 3.579482e+00 1.099871e+01 8.154943e+01
#> mu_alpha2 1.133719e+00 1.804835e+01 2.442407e+01 5.371687e+01
#> mu_alpha3 -2.301439e-01 1.983247e-02 3.253283e-01 6.389640e+01
#> sigma_alpha1 1.012255e+01 1.480995e+01 1.816466e+01 3.044661e+02
#> sigma_alpha2 1.291214e+01 1.969524e+01 2.986526e+01 2.060331e+02
#> sigma_alpha3 3.728872e-01 5.763740e-01 9.377868e+00 2.304189e+02
#> lp__ -2.535333e+03 -2.531686e+03 -2.531029e+03 -2.530101e+03
#> n_eff Rhat
#> alpha1 410.485703 1.0148006
#> alpha2 151.270414 1.0376582
#> alpha3 41.732339 1.0326831
#> tau_alpha1 188.270864 1.0003399
#> tau_alpha2 61.037429 1.0049624
#> tau_alpha3 4.046495 1.0010162
#> mu_alpha1 296.450085 1.0022898
#> mu_alpha2 77.016167 0.9997180
#> mu_alpha3 329.565805 1.0092678
#> sigma_alpha1 2293.404989 0.9997162
#> sigma_alpha2 788.772413 0.9998269
#> sigma_alpha3 1271.635495 0.9998229
#> lp__ 8.449041 1.0013864
plot(model)
#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.
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 (\(y_i \sim Bin(n_i, \pi_i)\)) with
\[ \text{logit}(\pi(a)) = \alpha_2 + \alpha_1\log(a) \]
The prior model of \(\alpha_1\) is specified as \(\alpha_1 \sim \text{truncated } \mathcal{N}(\mu_1, \tau_1)\) with flat hyperprior as in Farrington model
\(\beta\) is constrained to be positive by specifying \(\alpha_2 \sim \mathcal{N}(\mu_2, \tau_2)\)
The full conditional distribution of \(\alpha_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 \(\alpha_2\) can be derived in the same way
Fitting data
To fit Log-logistic model, use
hierarchical_bayesian_model()
and define
type = "log_logistic"
<- rubella_uk_1986_1987
df <- hierarchical_bayesian_model(age = df$age, pos = df$pos, tot = df$tot, type="log_logistic")
model #>
#> SAMPLING FOR MODEL 'log_logistic' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 1.8e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.18 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.379 seconds (Warm-up)
#> Chain 1: 0.387 seconds (Sampling)
#> Chain 1: 0.766 seconds (Total)
#> Chain 1:
#> Warning: There were 589 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: 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
$type
model#> [1] "log_logistic"
plot(model)
#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.