Discuss issues in specifying and inferring initial conditions, with particular reference to the polio example.
Suggest a possible improvement in the treatment of initial conditions here, code it up and make some preliminary assessment of its effectiveness. How will you decide if it is a substantial improvement?
Changing the initial conditions involves trying a new specification of rinit
.
After changing rinit
(or any other model component) we can re-run the entire analysis and see whether anything is substantially different.
If the new specification involves additional parameters, we must also add transformations for them, if appropriate, and random walk specifications. We may also want to add them to the scatterplots. Nothing else in the code needs to be changed.
A simplifying approximation in the original rinit
is to ignore initially infected babies, setting initial susceptibles to equal births and initial infected babies to zero.
Since the initial level of infection is low, and individuals do not stay long in the baby classes, this seems a reasonable approximation. Nevertheless, it is a testable one.
An alternative modeling approach is that babies might, to a first approximation, have the same prevalence as other individuals. IO_0/SO_0
is the current infection proportion for susceptible adults. If we suppose babies are exposed to this infection rate each month of their infancy, we can replace the previous rinit
accordingly. In practice, we use a more numerically stable alternative IO_0/(IO_0+SO_0)
which is practically equivalent in the anticipated situation with low prevalence, but is guaranteed not to exceed unity in pathological situations that might arise during maximization. Specifically, we replace
polio_rinit <- Csnippet("
SB1 = SB1_0;
SB2 = SB2_0;
SB3 = SB3_0;
SB4 = SB4_0;
SB5 = SB5_0;
SB6 = SB6_0;
IB = 0;
IO = IO_0 * P;
SO = SO_0 * P;
")
by
polio_rinit <- Csnippet("
double p=IO_0/(IO_0+SO_0);
SB1 = SB1_0 * pow(1-p,1);
SB2 = SB2_0 * pow(1-p,2);
SB3 = SB3_0 * pow(1-p,3);
SB4 = SB4_0 * pow(1-p,4);
SB5 = SB5_0 * pow(1-p,5);
SB6 = SB6_0 * pow(1-p,6);
IB = (SB1+SB2+SB3+SB4+SB5+SB6) * p/(1-p);
IO = IO_0 * P;
SO = SO_0 * P;
")
Here, there are no new parameters. We can just put the new rinit
into the R script and re-run it. The results are shown at initial-values-exercise/main.Rout.save. The key outcome is the maximized log likelihood. The search found a maximum of −794.7 for the modified model, whereas an identical search found −794.7 for the original model. The modification is certainly not a substantial improvement.
This version produced in R 4.3.2 on February 19, 2024.
Licensed under the Creative Commons Attribution-NonCommercial license. Please share and remix noncommercially, mentioning its origin.