crowdfightcovid19 / req-550-Syria

Repository to host the code of request 550 coming from the Pax Syriana Foundation
Other
3 stars 2 forks source link

Expressions for lambda and R0 #14

Closed Jennifer-Villers closed 4 years ago

Jennifer-Villers commented 4 years ago

@apascualgarcia @jordan-klein

Hi guys,

I had a look at the description of our model and derivation of R0 and compared it with the African modeling paper.

They define lambda a bit differently and I believe that they are right. In our case, I think that lambda should be defined as: image.

It would make sense because image represents the proportion of infectious people in the population group j, which influences the chances for a person i to become sick when meeting a person of group j.

If we assume that beta is the same for every infectious compartment (I don't remember whether we agreed on that or not), then we can simplify the expression: image

I also have one question regarding the derivation of R0: are you sure that 1/N is supposed to figure in it? Since it is associated with P, A, I and H to express a proportion, I would expect it not to be part of the expression that describes R0 but I am a bit agnostic about that.

Please let me know what you think.

Thanks!

jordan-klein commented 4 years ago

@Jennifer-Villers You're absolutely right! I updated my derivation to reflect this, and in the process, I made my derivation generalizable to a structured population. here

@apascualgarcia

apascualgarcia commented 4 years ago

@Jennifer-Villers and @jordan-klein Thanks both for pointing out the problem and finding the new expression. I reviewed Jordan's results and I think are correct, just please double check the definition of \Theta, I think it would be \Theta=C\diag(N)^(-1). I am going to review Eduard's proposal for the contacts matrix, note that the interpretation he does differs from your comment on the contact matrix "the rows of C describe the probabilities of contact with other classes and sum up to one", and it is related to issue #13, we could follow the discussion there or open a new issue.

jordan-klein commented 4 years ago

@apascualgarcia Done! Thanks for pointing out Theta.

apascualgarcia commented 4 years ago

@jordan-klein and @Jennifer-Villers

Hi guys,

I've been struggling with the computation of tau from lambda and R0 for a while and after reviewing the code, analytic expressions and literature I ended up with the conclusion that the problem comes from this expression for lambda:

imagen

I think it is incorrect to divide by N_j, and my argument comes as follows. The paper in which it is proposed this model the authors divide the population in different ages. But of course the thresholds to distinguish the population are often arbitrary, or at least could change under certain conditions. In our case, two classes may be distinguished by two areas in the camp.Therefore, I think that a basic requirement of a proper generalization of the basic SIR model into classes should be that if this distinction in classes is removed, we recover the model with a single class.

Consider a single infected compartment with two classes I1 and I2, with the same beta parameter. With the above formula lambda = beta* (I1/N1+I2/N2) and if I remove the distinction because the two populations mix, lambda = beta* I/N, how can be justified the formula in both cases if I1/N1+I2/N2 != I/N?

Moreover, working with that formula one ends up with a computation of the NGM that depends on N (and that was the source of all my problems). I think the correct way of modelling this is with the original formula that Judith proposed: lambda = beta/N* (I1+I2). This means that at steady state N=S (this requires assuming that there are mortality and birth rates, which we are neglecting but that if there is no infection should be there) and in the linearization N vanishes, so the matrix \Theta=C and the NGM becomes independent of N, as it must be (see for instance the first example here). Moreover, we end up with a tau<1 which would be impossible otherwise.

The most frustrating thing is that the authors of the two articles (who are the same) do not show the ODEs in the papers and they do not provide either the numerical values of the parameter u (equivalent to our tau) which I'm afraid may not have any sense. Let me know what you think.

jordan-klein commented 4 years ago

@apascualgarcia I disagree. The (P_j+A_j+I_j+H_j)/N_j term in lambda is the probability an individual is infected in each population class and the number of contacts between classes (Cij) needs to be multiplied by this probability as described in these papers. And if we were to collapse classes, wouldn't this be dealt with in the formula by the n in lambda's summation term \sum{j=1}^n changing (so for instance if we collapse all classes and assume a uniform contact rate, n=1)?

As far as tau derivation goes, I think this issue can be dealt with by fixing something that we might have overlooked in the linearization of the S_i x lambda_i term. At steady state, S_i=N_i, so we would linearize to tau x N_i x \sum C_ij/N_j. In matrix form we would express this as Theta = NCN^(-1) rather than CN^(-1). I tried calculating tau using some of the mean parameter values and this expression of Theta and got tau<1.

apascualgarcia commented 4 years ago

@jordan-klein thanks for the answer, I think I identified the source of discrepancy and we are all partially right. UPDATED, there was an error

The (P_j+A_j+I_j+H_j)/N_j term in lambda is the probability an individual is infected in each population class and the number of contacts between classes (C_ij) needs to be multiplied by this probability as described in these papers.

Please note that this is precisely the paper I was complaining about so yes, I noticed their interpretation. This is what they state: "Here, u is an individual’s susceptibility to infection upon contact with an infectious person, c ij,t is the number of age-j individuals contacted by an age-i individual per day at time t"

So u=tau and their C_ij and ours is also the same. In particular, since C_ij is the number of age-j individuals contacted by an age-i individual per day at time t then I assume that the way they build C_ij is the one we propose (but I don't know why they do not introduce it explicitly), namely: C_ij = c_i * N_j/N with c_i the average number of contacts and N_j/N the probability of finding an individual of class j by chance. In that case, it turns out that lambda_i = \sum_j tau * c_i * N_j/N * (P_j+A_j+I_j+H_j)/N_j = (1/N) * \sum_j tau * c_i* (P_j+A_j+I_j+H_j) which has the desired properties when compartments are merged/split. (This is, by the way, the additional simplification I was referring to in the Trello's card).

if we were to collapse classes, wouldn't this be dealt with in the formula by the n in lambda's summation term \sum_{j=1}^n changing (so for instance if we collapse all classes and assume a uniform contact rate, n=1)

There should be a natural way to present it mathematically (as above). If C_ij does not have this form I cannot see how. Proceeding as you suggest would make sense only if the classes were living in completely different areas and one susceptible individual walks in one area first and then in the other, or something like that.

As far as tau derivation goes, I think this issue can be dealt with by fixing something that we might have overlooked in the linearization of the S_i x lambda_i term. At steady state, S_i=N_i, so we would linearize to tau x N_i x \sum C_ij/N_j.

This was indeed my mistake, I had lost the index of S and then at steady state I assumed S=N.

In matrix form we would express this as Theta = NCN^(-1) rather than CN^(-1).

I guess you mean diag(N)*C*diag(N)^(-1)^1 but note that this is equal to C, which is what I used. However, following the above reasoning I see now that the result would be tau*diag(p)*c_i with p_i = N_i/N. Could you please double check that you agree?

I tried calculating tau using some of the mean parameter values and this expression of Theta and got tau<1

I obtain that the mean of tau is 0.025 with 95%CI [0.02, 0.03] after the correction is even lower: 0.006 [0.0045, 0.0075], do you obtain the same? I find this value somehow low though, it means that contacting with 100 infected will lead to an infection every 2-3 cases. This is not an issue because we are finding this value to obtain certain R0, but it suggests that looking for more realistic mechanistic models require to re-think assumptions.

Jennifer-Villers commented 4 years ago

Hi guys,

I am sorry for staying mainly silent the past week. The maths are getting very complicated for me and I don't have the right background to contribute at this stage.

I obtain that the mean of tau is 0.025 with 95%CI [0.02, 0.03] after the correction is even lower: 0.006 [0.0045, 0.0075], do you obtain the same? I find this value somehow low though, it means that contacting with 100 infected will lead to an infection every 2-3 cases. This is not an issue because we are finding this value to obtain certain R0, but it suggests that looking for more realistic mechanistic models require to re-think assumptions.

Regarding the low value for tau, it may have something to do with our extremely long infectious period for "hospitalized" people. Perhaps if we shorten that infectious period to a maximum of 10 days (one day above the maximum experimentally proven so far--the CDC suggested that it is very unlikely for people to be infectious longer than that), we could get a value for tau closer to our expectations.

jordan-klein commented 4 years ago

@apascualgarcia @Jennifer-Villers

Following up on my comment above, the same result is obtained whether Theta=diag(N)*C*diag(N)^(-1) or Theta=diag(p_i*c_i)U since they have the same spectral radius. @apascualgarcia your estimation of tau looks good!

I agree with @Jennifer-Villers that such a low tau could in part be due to our infectious period being so long. I'm not sure if this is a -need-to-change- because our parameters are sound and model is mathematically sound, but if we want to decrease the infectious period and see how much higher our tau is for external validity sake, that could make sense.

apascualgarcia commented 4 years ago

@jordan-klein

I don't think we can assume that C_ij=c_i * N_j/N universally and simplify lambda as you describe. This is only an expression for C_ij in a well-mixed population in which individuals from population class i are do not exhibit a differential preference for contact with any class j over another.

Assuming a well-mixed population is our null model, and the null model is what we must use to compute tau. You agreed with the parametrization of C_ij I proposed and, with the procedure to compute tau, and it is what I am using here.

I'm not sure I understand the logic here.

That justifying a departure from the null model requires a reasonable alternative model. If you have a better one please formulate it explicitly.

I don't agree with the result you derive.

In my tex file there is a command called \add{text here} that will allow you to highlight the error. Please go to the file, highlight the error, and provide the correct answer (also highlighted). Remember to add one more label after APG in the name of the file to track the version.

jordan-klein commented 4 years ago

@apascualgarcia Please disregard the previous comment, I got the same results as you, the derivation of tau looks good as is.