4.3 Application to German rotavirus incidence data
4.3.4 Statistical insights
ever, note that these averaged cumulative distribution functions are primarily displayed to illustrate the methodological possibilities as the main purpose of our analysis was to inves-tigate the effects of fixing parameters, not to obtain an averaged posterior among models with partially fixed parameters. One could have obtained the same posterior considering one single model with a mixed prior for the parameters (µ, ω) consisting of degenerated distributions and the respective skew-normal distributions defined above – each weighted equally.
1.0 1.5 2.0 2.5 3.0 3.5
0.00.20.40.60.81.0
µ
cdf
Θ(µ, ω)−model Θ(µ)−model Θ()−model averaged cdf prior cdf
0.5 1.0 1.5 2.0
0.00.20.40.60.81.0
ω
cdf
Θ(µ, ω)−model Θ(µ)−model Θ()−model averaged cdf prior cdf
Figure 4.7: Averaged posterior cumulative distribution function for the parameters µand ω according to the models using contact pattern C6.
relation is of particular interest considering the consequences these estimates might have on health policy making or subsequent epidemiological studies.
Parameter spaces and collinearities
Looking at the median estimates for the relative infectiousness of symptomatically infected individualspas shown in Figure 4.5, it becomes clear that certain parameter estimates react very sensitively to changes in model structure and parameter space. Of particular interest is the impact of fixing parameter values in advance as often done in the literature (Pitzer et al., 2009; Atchison et al., 2010; de Blasio et al., 2010). In our case, the models using configuration C6,Θ(µ,ω) and C6,Θ(µ), where one or both of the parameters µ = 7/3 and ω = 7/8 were fixed, the posterior median of variable p was 10.0 and 9.0, respectively.
In contrast the model using the full space C6, Z() yielded a median estimate for p at 12.6.
Adding only ω to the parameter vector ϑ did yield a slight decrease in the estimated value forp, while the posterior median forωitself at 1.93 was higher than its corresponding fixed value. Adding also µ, we obtained notably different estimates. In the corresponding model C6,Θ(), the posterior medians for ω and µ were given at 1.01 and 1.02, signif-icantly lower than the estimates from model C6,Θ(µ), while the median for p at 12.6 increased. With a closer look at the corresponding posterior correlations, it seems that the model prefers parameter vectors yielding the same µ/ω ratio and compensates variations in these mean durations of infection by adjusting the contact ratesc1, c2, c3 andp. In order to prevent large variations of some parameter’s estimates when adding other correlated parameters to the vector which is subject of estimation, one approach might be to define stricter priors instead of fixing a parameter entirely. In that regard, we have already seen that the prior distribution for the contact parameters appears to be too vague, especially in the light of the large amounts of data which dominate the posterior through the likelihood.
The above insights are also useful when thinking about potential extensions of the model structure. One modification could be the introduction of an infectious period preceding the symptomatic phase, which was left out in the original model because epidemiological studies (Anderson and Weber, 2004; Pickering et al., 1988) have shown this period to be rather short. However, while the estimated durations of symptomatic and asymptomatic phase might change under such a model variation, the relative infectiousnesspis likely to be robust. This can be seen from the duration parametersµandωvarying significantly among the compared models, while the estimate for parameter p remains relatively stable, since it is primarily affected by the mean infectiousness of a childhood case (many symptomatic infections) versus an adult case (almost no symptomatic infections), but only partially by the duration of symptomatic and asymptomatic period.
Two step optimization
Also of interest is the impact of the second step optimization on the posterior, in which we adjusted for the estimated cumulative autocorrelation in the data as given in equation (4.4). Using the modified likelihood LLCA had significant effects on the final estimation results, as we can see by, e.g., analysing the posterior estimation for parameter h(e).
As shown in Figure 4.4 in model C6,Θ(µ,ω) the second step posterior median for h(e) was 19.0% with 95% CI (17.6-20.5%) whereas after the first optimization step using only the posterior density as defined in equation (4.3) the median was estimated at 19.6%
(19.3-19.9%). By reducing the impact of the likelihood by a factor determined by the cumulative autocorrelation the posterior density flattened around the posterior mode. As a consequence the posterior median estimates moved closer to the prior medians, but more importantly we obtained a larger posterior variance and hence also wider credibility regions.
Furthermore, by downscaling the likelihood the marginal likelihoods of each of the selected models, as computed via equation (4.9), moved closer together. This can also be observed in the case of vertical averaging with respect to contact structure C5, where the largest difference of the marginal likelihoods was 7.2 on the log scale after the second step procedure while the minimal difference after only one estimation step was 25.4, which would lead to nearly degenerated model weights.
Normal approximation of the posterior
The accuracy of the normal approximation of the posterior distribution was already in-vestigated in Section 2.3.4, although for a much simpler epidemic model. Here, we aim to check whether the normal approximation holds also for the higher dimensional setting of the rotavirus transmission model.
As it is not feasible to compare the whole 17-dimensional posterior against its approx-imation, we instead checked the match of the respective conditional log posteriors with respect to one or two parameter components, conditioned on all other parameter compo-nents being equal to their posterior mean. The conditional log posterior density of the parameter component ϑi is given by
logπCAϑiϑ−i =ϑ∗−i,D= logπCAϑi,ϑ∗−i|D−logπCA(ϑ∗|D),
where ϑ−i denotes the vector of components in ϑ which are not ϑi, and ϑ∗ denotes the second step posterior mode ϑ∗CA. Hereby, the term −logπCA(ϑ∗|D) functions as an additive constant such that the conditional log posterior’s maximum is at zero (this does not imply that logπCA(ϑ∗|D) is the normalizing constant of the log posterior). Analogously the conditional log posterior according to the normal approximation as defined in (4.5) is given by
logπN
ϑiϑ−i =ϑ∗−i,D= logφϑ∗,Σϑi,ϑ∗−i−logφϑ∗,Σ(ϑ∗)
=−1 2
ϑi,ϑ∗−i−ϑ∗T Σ−1ϑi,ϑ∗−i−ϑ∗,
whereΣ is the approximated posterior covariance matrix based on the Fisher information (see Equation (4.5)) and (ϑi,ϑ∗−i) denotes ϑ∗ with itsi-th component being replaced byϑi. Thus, the conditional log densities each represent a slice of the corresponding unconditioned log posterior density. To check whether the two conditional densities are approximately equal, we evaluated them over a range determined by the 95% credibility region (CrR) of the unconditioned posterior. Based on the approximated posterior the random variable
(ϑ−ϑ∗)TΣ−1(ϑ−ϑ∗)
is χ2d distributed with 17 degrees of freedom. Thus, the 95% CrR is given by CrR95% =
(
ϑ ∈Θ
logφϑ∗,Σ(ϑ)−logφϑ∗,Σ(ϑ∗)>−q0.95(χ217)
2 ≈ −13.79
)
.
We calculated the conditional posteriors for the parametersµandωaccording to model
C6,Θ(). Note that the densities were computed on the respective log transformed space.
The conditional log densities according to the true posterior and its normal approximation together with the 95% CrR thresholds are given in Figure 4.8 (top row). Additionally, we computed the two-dimensional conditional log posterior for the parameter vector (µ, ω) and plotted the contour lines where the true posterior and the normal approximation cross the thresholds corresponding to the 5%-, 25%-, 50%-, 75%-, and 95% CrRs (bottom figure).
One can see that the normal approximation matches the actual log posterior density quite well such that all qualitative statements about the parameter estimates, e.g. those from Section 4.3.2, remain true. However, the match is certainly not perfect as the pa-rameter values at which the densities cross the critical 95% CrR threshold slightly deviate.
Moreover, the approximation is well capable to capture the parameter correlations as dis-played by the 2-dimensional distribution in Figure 4.8. It should be noted that the relative deviation of posterior approximation increases with growing distance to the posterior mode.
This can be easily explained since the normal approximation is based on a second order Taylor approximation (Abramowitz and Stegun, 1964) of the log posterior around its mode and therefore, becomes increasingly inaccurate within distant regions.
Convergence of the optimization procedure
One of the main reasons for computing the normal approximation of the posterior was to save many costly model evaluation which would be otherwise required if, e.g., the posterior had been approximated by an MCMC-sample (see Section 3.2). However, computing the posterior mode using the Nelder-Mead algorithm as described in Section 4.2 already required up to 100,000 posterior evaluations for some considered models, which yielded a computation time of a few days for one model. Moreover, it is not possible to shorten the procedure at the cost of approximation accuracy as the approximation approach requires the Hessian at the computed mode to be negative definite, which is not guaranteed if the optimization is stopped too early.
0.00 0.05 0.10 0.15 0.20 0.25
−15−10−50
unnormalized conditional log posterior
log µ log posterior normal approximation 95% CI threshold
0.085 0.090 0.095 0.100
−20−15−10−50
unnormalized conditional log posterior
log ω log posterior normal approximation 95% CI threshold
Contour lines of conditional log posterior
log µ
log ω
5%
25%
50%
75%
95%
95%
0.00 0.05 0.10 0.15 0.20 0.25
0.0850.0900.0950.100
5%
25%
50%
75%
95% log posterior
normal approximation
Figure 4.8: Comparison of the conditional log posterior density and its corresponding normal approximation for the single parametersµ,ω and the joint parameter vector. Also displayed are the approximate 95% credibility region borders corresponding to the full posterior distribution.