Correct contact matrix processing?

[pratikunterwegs]

Hi @adamkucharski - something I noticed when looking at this repo as a reference, so thought I'd flag - the ODE models call .prepare_population() internally, which scales the contact matrix by its largest real eigenvalue (before scaling each row by the size of the corresponding demographic group).

This doesn't look right, and it's probably a holdover from a period when the models were parameterised in terms of the $R_0$ and infection duration, rather than the transmission rate $\beta$, so this step should probably be removed?

[adamkucharski]

There will need to be some normalisation at some point (either in definition of contact rates or transmission rate) to ensure the dominant eigenvalue of the next generation matrix is equal to R0 (i.e. beta matrix * duration of infection – as in the modelling course). So that would be a quick sense check to do to confirm that the current implemention is returning sensible R0 value – from the current plots, results don't look particularly strange.

[pratikunterwegs]

Yep, the plots look okay (e.g. Readme), but note the time-scale is 600 days, which seems like a lot for a single epidemic peak?

The package functions are not using $R0$ and infection duration, but $\beta$ directly. So I would expect force of infection to only be a function of $\beta$, $M{ij}$ and $Ij$, with $M{ij} = C_{ij} / N_j$, where $C$ is the raw contact matrix scaled by the demography.

Whereas in all ODE models, $M{ij} = (C{ij} / \text{max(EV}(C))) / N_{ij}$ - this seems not right? Just thought I'd flag in case this needs correcting.

[adamkucharski]

Although the README plot has interventions, so makes sense it takes longer to run the epidemic with R0=1.5. Note that beta is a function of R0 and contact rates (it's hard to define a meaningful beta a priori without either considering R0 and contacts, or fitting to data). So from what I can see, the beta = R0/duration definition in the readme is only valid because the matrix has already been normalised by dominant eigenvalue (and hence removed the effect of the overall magnitude of contact rates).

Effectively the model (as far as I can tell) is still fundamentally parametrised in terms of R0, but beta is now defined as R0/D and used as main input.

[adamkucharski]

Closing as don't think this needs further action (and alignment of {odin} implementation validates that it's doing what is expected).

[pratikunterwegs]

Great, thanks for taking a look.