Replies: 8 comments 13 replies
-
Additional ruled out case
Additional clue
|
Beta Was this translation helpful? Give feedback.
-
|
Currently, the transitions are defined as follows in transitions:
# S --> I # New infections caused by contact with infecteds in age0to5
- source: [["S"], ["age0to5", "age5to15", "age15to65", "age65to120"]]
destination: [["I"], ["age0to5", "age5to15", "age15to65", "age65to120"]]
proportional_to: [
"source",
[
["I"],
[["age0to5"],["age0to5"],["age0to5"],["age0to5"]]
]
]
rate: [["beta"], ["c_age0to5_age0to5","c_age5to15_age0to5","c_age15to65_age0to5","c_age65to120_age0to5"]]
proportion_exponent: [
["1",["1","1","1","1"]],
["1",["1","1","1","1"]]]
# S --> I # New infections caused by contact with infecteds in age5to15
- source: [["S"], ["age0to5", "age5to15", "age15to65", "age65to120"]]
destination: [["I"], ["age0to5", "age5to15", "age15to65", "age65to120"]]
proportional_to: [
"source",
[
["I"],
[["age5to15"],["age5to15"],["age5to15"],["age5to15"]]
]
]
rate: [["beta"], ["c_age0to5_age5to15","c_age5to15_age5to15","c_age15to65_age5to15","c_age65to120_age5to15"]]
proportion_exponent: [
["1",["1","1","1","1"]],
["1",["1","1","1","1"]]] and similarly for new infections caused by contact with infecteds in age15to65 and age65to120. This means that we define every possible S --> I transition four times, i.e. four times S_age0to5 --> I_age0to5, etc. But are we sure this gets summed internally in flepi? Doesn't it get overwritten instead? |
Beta Was this translation helpful? Give feedback.
-
|
It is possible to confirm this by setting all elements of the contact matrix to 24.1, omitting any of the four S --> I transition blocks one at a time. In all cases the epidemic dies down, so I do think flepi does sum internally. |
Beta Was this translation helpful? Give feedback.
-
|
Summarizing my slack note here, I'm a bit confused about the FOI element in the reference model. But I also don't know how flepiMoP handles these sort of mobility / age-contact pattern matrices. I'll review my expected representation / calculation below (as random sample of modelers that might use this tool / what my feedback as a reviewer would have been on the original work). If I'm understanding correctly: FOI in a location (by age) should be the the fraction of that population in that location (visitors & non-visitors) that is infectious (* their generation of contacts & infection probability per contact). Then susceptible individuals should be exposed to that FOI weighted by their time in home vs away(s). It seems like the Then Then susceptibles of age a, location i would get f_v * FOI to age, at i + (1 - f_v)*sum(FOI age, at not i) exposure. I don't quite see if that's equivalent to what's in the reference model. |
Beta Was this translation helpful? Give feedback.
-
|
Hi Carl, some pointers that might help clarifying this. Happy to discuss if you see any mistakes in the reference, I'd be more-than-happy to have them pointed out.
|
Beta Was this translation helpful? Give feedback.
-
intriguing - this seems worth double checking for implementation equivalence. |
Beta Was this translation helpful? Give feedback.
-
|
@twallema have you tried back-interpretting what flepimop does into your python model? as in, instead of trying to reproduce what you've done in flepimop, reproduce what flepimop is doing by modifying your baseline code? |
Beta Was this translation helpful? Give feedback.
-
|
Also, re this - it seems like flepimop is setting the "self" mobility via complement? https://iddynamics.gitbook.io/flepimop/gempyor/model-description#population-structure |
Beta Was this translation helpful? Give feedback.
-
|
Thanks Carl, it seems that way indeed. I'll be sure to match mobility implementations today and get back to you. |
Beta Was this translation helpful? Give feedback.
-
|
I've matched the definition of the mobility matrices, in the reference mobility matrix rows now sum to one. The simulation result itself is largely unaltered (peaks a couple of days sooner). I'll see if I can reverse engineer the reference back to flepi, I like that suggestion.
|
Beta Was this translation helpful? Give feedback.
-
|
@pearsonca I've typed out an example computation of how the number of new infections in the reference model are computed. It's for two age groups and three spatial patches for the sake of illustration. Hope that helps to understand it better, I'm unsure how I would translate it into mathematics, but by typing that example out I know it does what I'd like it to do. If you would have suggestions on how this can be translated to maths I'd be open to them, because I'm a hack mathematician. |
Beta Was this translation helpful? Give feedback.
-
|
Sorry, I've been really slow and I really want to dig this, seems very important. Late and could not finish so here are some pointers on reproducing. Thank you @twallema for taking the time to inpect the model, very welcome after all this time. https://we.tl/t-C31BWkpNl5 this file is the "integration_dump.pkl" from the last config. I'll detail the steps on how to obtain it, but to use it do: import pickle
fn_dump = "integration_dump.pkl"
with open(fn_dump, 'rb') as pickle_file:
states, states_daily_incid, ncompartments, nspatial_nodes, ndays, parameters, \
dt, transitions, proportion_info, transition_sum_compartments, initial_conditions, \
seeding_data, seeding_amounts, mobility_data, mobility_row_indices, mobility_data_indices, \
population, stochastic_p, method = pickle.load(pickle_file)and you have everything that goes as argument to the rk4 integrator. if you want to retry integration, do import gempyor
output = gempyor.steps_rk4.rk4_integration(
ncompartments=ncompartments,
nspatial_nodes=nspatial_nodes,
ndays=ndays,
parameters=parameters,
dt=dt,
transitions=transitions,
proportion_info=proportion_info,
transition_sum_compartments=transition_sum_compartments,
initial_conditions=initial_conditions,
seeding_data=seeding_data,
seeding_amounts=seeding_amounts,
mobility_data=mobility_data,
mobility_row_indices=mobility_row_indices,
mobility_data_indices=mobility_data_indices,
population=population,
stochastic_p=stochastic_p)Will dig further as soon as I can |
Beta Was this translation helpful? Give feedback.
-
|
What is not entirely clear to me is why should the model with the mobility sum over the rows be equal to the one where each transition is modeled separately. right ? and the second model does: which yield vastly different results but that is expected (multiplication of sums != sum of multiplications), and the second model should have a much slower epidemic growth, as it does. Is that not the case PySODM ? As to why the correct model (the second one) yield different results compared to PySODM, it might because the mobility matrix is row count and scaled by total population to obtain rate (a debatable choice, like to census data, I would also prefer rate) ? I can dig that further when I come back from holidays. |
Beta Was this translation helpful? Give feedback.
-
|
Yes, you are right. The "rowsum" idea should have So for origin age group being destination-indiscriminate means (neglecting mobility), |
Beta Was this translation helpful? Give feedback.
-
|
In case, the documentation of the config syntax for transition is here which is not that good and we should add better examples. As mentioned, I find this syntax confusing and not that helpful, would appreciate thoughts on how to make it better. To further clarify how flepi translates this into equations (following @pearsonca suggestion) I made a little script to output equations from gempyor directly. Basically it runs the rhs of flepiMoP equation, but for each operation, it does the same with a string that represents the parameter. See below for implementation. It's still rough and does not output proper latex (yet). I have attached the transitions for your two configurations as understood by flepiMoP. If we look at the transition that have source and destination in the age0to5 age group, in the first config: ## transition 0 : Sage0to5 > Iage0to5
delta : (Sage0to5) * ((1*(Sage0to5)^(1*1=1.0)/(Sage0to5))*(((1-\sum_i (M_ij/pop_i)_i*p_away)*(Iage0to5+Iage5to15+Iage15to65+Iage65to120)_i^(1*1=1.0)_i/ H_i * (beta*24.8_i=0.7439999999999999))+((p_{away}*M_{ij}/H_i)*((Iage0to5+Iage5to15+Iage15to65+Iage65to120)_j^(1*1=1.0)_j/ H_j * (beta*24.8_j=0.7439999999999999)))))
(this was a multiple proportion transition, so mobility is activated)
## transition 4 : Iage0to5 > Rage0to5
delta : (Iage0to5) * (1*(gamma*1=0.2))
(this was a single proportion transition, no mobility)and in the second config: ## transition 0 : Sage0to5 > Iage0to5 proportional to Iage0to5
delta : (Sage0to5) * ((1*(Sage0to5)^(1*1=1.0)/(Sage0to5))*(((1-\sum_i (M_ij/pop_i)_i*p_away)*(Iage0to5)_i^(1*1=1.0)_i/ H_i * (beta*cage0to5age0to5_i=0.297))+((p_{away}*M_{ij}/H_i)*((Iage0to5)_j^(1*1=1.0)_j/ H_j * (beta*cage0to5age0to5_j=0.297)))))
(this was a multiple proportion transition, so mobility is activated)
## transition 4 : Sage0to5 > Iage0to5 proportional to Iage5to15
delta : (Sage0to5) * ((1*(Sage0to5)^(1*1=1.0)/(Sage0to5))*(((1-\sum_i (M_ij/pop_i)_i*p_away)*(Iage5to15)_i^(1*1=1.0)_i/ H_i * (beta*cage0to5age5to15_i=0.10800000000000007))+((p_{away}*M_{ij}/H_i)*((Iage5to15)_j^(1*1=1.0)_j/ H_j * (beta*cage0to5age5to15_j=0.10800000000000007)))))
(this was a multiple proportion transition, so mobility is activated)
## transition 8 : Sage0to5 > Iage0to5 proportional to Iage15to65
delta : (Sage0to5) * ((1*(Sage0to5)^(1*1=1.0)/(Sage0to5))*(((1-\sum_i (M_ij/pop_i)_i*p_away)*(Iage15to65)_i^(1*1=1.0)_i/ H_i * (beta*cage0to5age15to65_i=0.3090000000000001))+((p_{away}*M_{ij}/H_i)*((Iage15to65)_j^(1*1=1.0)_j/ H_j * (beta*cage0to5age15to65_j=0.3090000000000001)))))
(this was a multiple proportion transition, so mobility is activated)
## transition 12 : Sage0to5 > Iage0to5 proportional to Iage65to120
delta : (Sage0to5) * ((1*(Sage0to5)^(1*1=1.0)/(Sage0to5))*(((1-\sum_i (M_ij/pop_i)_i*p_away)*(Iage65to120)_i^(1*1=1.0)_i/ H_i * (beta*cage0to5age65to120_i=0.03000000000000002))+((p_{away}*M_{ij}/H_i)*((Iage65to120)_j^(1*1=1.0)_j/ H_j * (beta*cage0to5age65to120_j=0.03000000000000002)))))
(this was a multiple proportion transition, so mobility is activated)
## transition 16 : Iage0to5 > Rage0to5
delta : (Iage0to5) * (1*(gamma*1=0.2))
(this was a single proportion transition, no mobility)the mean numeric parameters value across subpop and times is shown as well. This confirms the above observation about the models being different, and I believe these equations match what would be expected from the configs as well. Implementation to print equation: very annoying I cannot attach a ipynb in github, see the notebook here: https://we.tl/t-PuRlCAq2Os(execute the first cells, then go to "JUMP THERE TO PRINT TRANSITIONS", it is dirty as it starts from my run_analysis notebook, I just deleted most of the cells). The goal would be to have it output proper latex, and then have a flag in gempyor that allows printing these transitions. |
Beta Was this translation helpful? Give feedback.
-
|
i someone wants to finish this and make it part of gempyor, that would be great otherwise I'll continue after vacations :) We could have a command |
Beta Was this translation helpful? Give feedback.
-
|
Thanks for digging @jcblemai! Having a script that outputs the equations is a helpful addition. I've typed up a step-by-step computation of the reference model along with the mathematical equations for the reference model: Will try to figure out today where the exact difference lies & what caused any confusion in the syntax. I was not able to deduce from the example in the documentation how the transitions present in both configs are currently coded up, so I actually had to rely heavily on @saraloo & @anjalika-nande to help with the code-up (indicating an overhaul with some more examples will likely be useful). |
Beta Was this translation helpful? Give feedback.
-
|
Quick question: What would |
Beta Was this translation helpful? Give feedback.
-
|
Closing the week off, the equations for the reference model are, for age group which is the sum of the infections happening on the "home" patch and the infections happening on the "visited" patch. The visiting populations are, with |
Beta Was this translation helpful? Give feedback.
-
|
@jcblemai What do the the H_i and H_j terms in Flepi's equations represent? |
Beta Was this translation helpful? Give feedback.
-
Most likely cause of the discrepancyA difference in normalisation in the FOIIn the reference model, for an age group with the visiting populations, Here, with,
The infectious population in every age group Verifying this is the likely cause: reverse engineeringThe model below implements the FlepiMoP FOI. class spatial_ODE_SIR(ODE):
"""
SIR model with FlepiMoP-like normalization
"""
states = ['S','I','R']
parameters = ['beta','gamma', 'f_v', 'N', 'M']
dimensions = ['age', 'location']
@staticmethod
def integrate(t, S, I, R, beta, gamma, f_v, N, M):
# Compute total population and visited populations
T = S + I + R
T_v = np.matmul(T, M)
I_v = np.matmul(I, M)
# sum total populations over the age groups
T = np.sum(T, axis=0)
T_v = np.sum(T_v, axis=0)
# divide every I_ij by T_j = sum_j T_ij
I_div_T = I / T[np.newaxis, :]
I_div_T_v = I_v / T_v[np.newaxis, :]
# compute force of infection
lambda_home = beta * (1-f_v) * np.einsum ('kj,ik -> ij' , I_div_T, N)
lambda_visit = beta * f_v * np . einsum ('jk, lk, il -> ij', M, I_div_T_v, N)
# calculate differentials
dS = - (lambda_home + lambda_visit) * S
dI = (lambda_home + lambda_visit) * S - 1/gamma*I
dR = 1/gamma*I
return dS, dI, dRAnd yields the following simulation output, which is very similar tot the output that triggered my further inquiries,
The implications of the FlepiMoP FOI: "heterogeneity-inconsistency in beta"To help study the implications it helps to set up a fictitious example with two age groups (young & old) and three spatial patches (A, B and C). Case 1: Homogeneous mixingAssuming Using the age groups (
Collapsing the age structure (
Conclusion: The reference model results in an identical number of new infections, irregardless if the age groups are present or collapsed. In the FlepiMoP implementation, the "apparent" Case 2: Heterogeneous mixingAssuming Using the age groups (
Collapsing the age structure (
Conclusion: The reference model results in 26 new infections if the age groups with heterogeneous contact rates are present, and 27 new infections when the age groups are collapsed. This is not identical, but it needn't be, a slight difference is even expected because collapsing the age groups removes heterogeneity and allows the disease to spread more easily. Opposed, when expanding the age structure from For lack of a better word, you could say the current FlepiMoP normalization makes beta "heterogeneity-inconsistent". Other implicationsImplications are not limited to the FOI, in Josh' flu model we can model the vaccination with an additional disease state for vaccination "V" and infection of a vaccinated individual "I_v". In this model, we transfer individuals from the S --> V compartment with a rate derived from a logistic model that results in 50% of the S being transferred to V.
When we try to model the vaccinations thru a strata, the rate of vaccination now gets normalized with the sum of (S,unvacc) and (S,vacc), which is the equivalent of normalizing the rate wit the sum of S and V in the original implementation. Because we are normalising the rate with more individuals than necessary, the apparent rate becomes lower,
Proposed changesChanges should be implemented in
|
Beta Was this translation helpful? Give feedback.
-
A test case: an SIR model with age- and space stratification
Over the past weeks I set up an SIR model for Belgium using a software I developed to simulate & calibrate dynamical systems in Python. My goal is to reproduce the model using flepiMoP in order to learn its syntax, gain insights in its code and test its speed. Additionally, having an age+space structured SIR with commuter mobility and sociall contact matrix could make a good template for flepi's examples directory.
This model has 581 spatial patches representing Belgian municipalities, connected with a commuter mobility matrix and four age groups whose contacts are governed by a contact matrix. One infected individual is seeded in the age group 65-120 yo in Arlon (bottom right-hand corner).
The "factual": pySODM variant
Model states
S,I, andRare 2D numpy matrices due to their dimensionsageandlocation.beta(float) is the infectivity parameter (0.03),gamma(float) is the rate of recovery (0.02),f_v(float) is the fraction of contactsN(2D np.ndarray) individuals have on the home patch (0.1).M(2D np.ndarray) is a square OD mobility matrix, diagonal elements have been set to zero to match flepi's definition of mobility.As you can see, some pretty barebones matrix multiplications, I'm fairly certain these are correct because of prior validation we did in Belgium (Please let me know if you see a flaw!). The model output looks like this,
The counterfactual: flepiMoP
I've been building the FlepiMoP equivalent in three steps ranging from least complex to identical to the factual.
1. An SIR model with spatial stratification
See
twallema/flepimop_influenza_BE/use_flepimop/wo_age_groups/SIR_BE_wo_age.yml,which produces the following output,
which is quite close to the factual case, especially considering the reduced heterogeneity (omitting age groups) should result in a more intense peak. So far so good.
2. An SIR model with age and space stratification, but contacts are averaged over the contact matrix's rows
See
twallema/flepimop_influenza_BE/use_flepimop/wo_age_groups/SIR_BE_w_age_rowsum.yml,which implements the transitions from
S_age0to5 --> I_age0to5proportional tobeta*\sum_j N_ij * S_age0to5 * (I_age0to5 + I_age5to15 + I_age15to65 + I_age65to120)(which is what I mean by row-sum of the contact matrix). This produces the following output,hell yeah! Getting there. Less people get infected, the peak is slightly lower, as you would expect when adding heterogeneity.
3. An SIR model with age and space stratification
Moment of truth, this one is equal to the pySODM variant. See
twallema/flepimop_influenza_BE/use_flepimop/wo_age_groups/SIR_BE_w_age.yml,which produces the following output,
Now, when I double beta from
beta=0.03tobeta = 0.06I get a consistent result,which is not likely a coincidence. My current hunch is "beta" is somehow sliced in two because it is broadcast over two dimensions. Any help figuring this one out would be appreciated.
Ruled-out causes
dt: does not impact simulation results (as it should).Reproducing the issue
All code needed to reproduce the model results shown in this discussion can be found in https://github.com/twallema/flepimop_influenza_BE.
The pySODM model can be found in the folder
use_pySODM, while its flepiMoP counterparts can be found inuse_flepiMoP. Instal pySODM usingpip install pySODM, which should install v0.2.4.A living document detailing the issue can be found in
twallema/flepimop_inlfuenza_BE/comparison.docx.Hoping for your most sincere un-TLDR insights, yours truly Tijs.
Beta Was this translation helpful? Give feedback.
All reactions