library(serosv)
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:
Refer to Chapter 10.3.1
Proposed model
The model for prevalence is as followed
π(a)=1−exp{α1α2ae−α2a+1α2(α1α2−α3)(e−α2a−1)−α3a}
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 N(μj,τj), j=1,2,3
The joint posterior distribution for α can be derived by combining the likelihood and prior as followed
P(α|y)∝m∏i=1Bin(yi|ni,π(ai,α))3∏i=1−1τjexp(12τ2j(αj−μj)2)
Where the flat hyperprior distribution is defined as followed:
μj∼N(0,10000)
τ−2j∼Γ(100,100)
The full conditional distribution of αi is thus P(αi|αj,αk,k,j≠i)∝−1τiexp(12τ2i(αi−μi)2)m∏i=1Bin(yi|ni,π(ai,α))
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)
<- 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: 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 6.9e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.69 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.602 seconds (Warm-up)
#> Chain 1: 0.839 seconds (Sampling)
#> Chain 1: 2.441 seconds (Total)
#> Chain 1:
#> Warning: There were 819 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: 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.388515e-01 3.143483e-04 5.686656e-03 1.283237e-01
#> alpha2 1.971229e-01 4.240475e-04 7.614690e-03 1.841599e-01
#> alpha3 7.639700e-03 4.500291e-04 6.817731e-03 3.534974e-04
#> tau_alpha1 1.663986e-01 2.691150e-02 3.429053e-01 1.791474e-06
#> tau_alpha2 7.574213e-02 2.141102e-02 1.909780e-01 3.895651e-06
#> tau_alpha3 2.082875e-01 7.744633e-02 6.028339e-01 2.965544e-06
#> mu_alpha1 -1.855848e-01 3.037000e+00 4.764718e+01 -1.097078e+02
#> mu_alpha2 2.225335e-01 2.954215e+00 4.037329e+01 -9.986680e+01
#> mu_alpha3 -4.196577e+00 2.551721e+00 4.143201e+01 -1.027110e+02
#> sigma_alpha1 1.434457e+02 6.087243e+01 1.397043e+03 8.863586e-01
#> sigma_alpha2 9.138533e+01 1.633538e+01 5.341490e+02 1.235928e+00
#> sigma_alpha3 9.076527e+01 1.643673e+01 4.931546e+02 6.803698e-01
#> lp__ -2.536086e+03 3.070347e-01 3.741862e+00 -2.544322e+03
#> 25% 50% 75% 97.5% n_eff
#> alpha1 1.351990e-01 1.384015e-01 1.424784e-01 1.506125e-01 327.25899
#> alpha2 1.921341e-01 1.963159e-01 2.016397e-01 2.133971e-01 322.45975
#> alpha3 2.424724e-03 5.710295e-03 1.043535e-02 2.551799e-02 229.50835
#> tau_alpha1 3.664167e-04 7.883421e-03 1.445835e-01 1.272861e+00 162.35759
#> tau_alpha2 3.012210e-04 3.573543e-03 4.029680e-02 6.546563e-01 79.55953
#> tau_alpha3 4.317701e-04 4.822130e-03 8.921679e-02 2.160280e+00 60.58899
#> mu_alpha1 -6.208494e+00 1.048231e-01 5.859708e+00 1.104209e+02 246.14159
#> mu_alpha2 -7.215519e+00 3.118702e-01 9.501639e+00 8.932361e+01 186.76865
#> mu_alpha3 -9.924458e+00 -1.992227e-01 5.134974e+00 8.396866e+01 263.63647
#> sigma_alpha1 2.629909e+00 1.126270e+01 5.224112e+01 7.474100e+02 526.71797
#> sigma_alpha2 4.981555e+00 1.672828e+01 5.761791e+01 5.068755e+02 1069.21787
#> sigma_alpha3 3.347934e+00 1.440060e+01 4.812534e+01 5.808039e+02 900.19228
#> lp__ -2.538429e+03 -2.535551e+03 -2.533332e+03 -2.529930e+03 148.52535
#> Rhat
#> alpha1 1.0103855
#> alpha2 1.0026672
#> alpha3 0.9998278
#> tau_alpha1 1.0071998
#> tau_alpha2 1.0304386
#> tau_alpha3 1.0395851
#> mu_alpha1 1.0195686
#> mu_alpha2 0.9999242
#> mu_alpha3 0.9997145
#> sigma_alpha1 1.0030298
#> sigma_alpha2 0.9998938
#> sigma_alpha3 0.9997435
#> lp__ 1.0038998
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
π(a)=βaα1+βaα, α,β>0
The likelihood is specified to be the same as Farrington model (yi∼Bin(ni,πi)) with
logit(π(a))=α2+α1log(a)
The prior model of α1 is specified as α1∼truncated N(μ1,τ1) with flat hyperprior as in Farrington model
β is constrained to be positive by specifying α2∼N(μ2,τ2)
The full conditional distribution of α1 is thus
P(α1|α2)∝−1τ1exp(12τ21(α1−μ1)2)m∏i=1Bin(yi|ni,π(ai,α1,α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"
<- 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 2.6e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.26 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.364 seconds (Warm-up)
#> Chain 1: 0.161 seconds (Sampling)
#> Chain 1: 0.525 seconds (Total)
#> Chain 1:
#> Warning: There were 843 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.