GiulioRossetti / ndlib

Network Diffusion Library - (for NetworkX and iGraph)
http://ndlib.readthedocs.io/
BSD 2-Clause "Simplified" License
273 stars 79 forks source link

(non) probabilistic E > I transition (SEIR) #157

Closed ggrrll closed 4 years ago

ggrrll commented 4 years ago

hi,

if I understood correctly, here

https://github.com/GiulioRossetti/ndlib/blob/450216abd771cf945064a5f21659fa3eed36df1b/ndlib/models/epidemics/SEIRModel.py#L90

you have a 'sharp' transition (step function, w.r.t time from exposure)

I think a probabilistic version would be better (more realistic), like the one I am attempting here:

https://github.com/ggrrll/ndlib/commit/6fcec70b9948c62ae0935d061ab2451f31e113e2#diff-55a281f20b0ed52ec6ac2dc08e33ce64R91

that is basically the one mentioned in the paper you referenced

ggrrll commented 4 years ago

of course, if you want to keep it similar to the one of the paper, also this

https://github.com/ggrrll/ndlib/blob/6fcec70b9948c62ae0935d061ab2451f31e113e2/ndlib/models/epidemics/SEIRModel.py#L91

should be changed as ~ 1 - exp( - 𝛼 * t)

GiulioRossetti commented 4 years ago

Hi, I completely agree. The implementation was made by a group of students (as part of a small project) that opted for a step function to simplify their implementation.

If you like to address the discrepancy feel free to make a pull request, I'll be happy to merge it in the master branch!

Giulio

ggrrll commented 4 years ago

sorry, the implementation in this commit/line

https://github.com/ggrrll/ndlib/blob/6fcec70b9948c62ae0935d061ab2451f31e113e2/ndlib/models/epidemics/SEIRModel.py#L91

is wrong, as we should of course consider the time, since exposure

ggrrll commented 4 years ago

closed due to #164

mihaibanciu commented 4 years ago

@GiulioRossetti @ggrrll Sorry to revisit this, but according to the paper referenced at the top of the thread, both E->I and I->R transition probabilities follow exponential distributions. Here's what it says in section 2:

(iii) an exposed individual's probability of becoming infectious in a specified time interval is independent of time after initial contact, so that the probability of remaining in the exposed class at time t after initial contact is e^(-at), where 1/a is the mean latent period; (iv) after an individual becomes infective, the probability of that individual recovering after time t is e^(-gt), where 1/g is the mean infectious period.

Accordingly, I modified the code to reflect the same logic in transitioning from I to R. Didn't want to open a new issue, so I'm posting the suggestion below:

@ -62,6 +62,9 @@ class SEIRModel(DiffusionModel):
            else:
                return {"iteration": 0, "status": {},
                        "node_count": node_count.copy(), "status_delta": status_delta.copy()}
+            for u in self.graph.nodes:
+                if self.status[u] == 1:
+                    self.progress[u] = 0 # set timestamp t_0 = 0 for all infected nodes in fraction_infected

        for u in self.graph.nodes:

@ -87,14 +90,15 @@ class SEIRModel(DiffusionModel):

            elif u_status == 2:

-                # apply prob. of infection, after (t - t_i) 
-                if eventp < 1 - np.exp(- (self.actual_iteration - self.progress[u]) * self.params['model']['alpha']): 
+                # apply prob. of infection, after (t - t_i)
+                if eventp < 1 - np.exp(- (self.actual_iteration - self.progress[u]) * self.params['model']['alpha']):
                    actual_status[u] = 1  # Infected
-                    del self.progress[u]
+                    self.progress[u] =  self.actual_iteration # save time of infection t_i

            elif u_status == 1:
-                if eventp < self.params['model']['gamma']:
+                if eventp < 1 - np.exp(- (self.actual_iteration - self.progress[u]) * self.params['model']['gamma']):
                    actual_status[u] = 3  # Removed
+                    del self.progress[u]

        delta, node_count, status_delta = self.status_delta(actual_status)
        self.status = actual_status

Apologize if this is the wrong way to go about. I'm a noob at using Github and didn't want to break anything by accident!

ggrrll commented 4 years ago

@mihaibanciu sorry for the late reply

Actually, this is something I have implemented few moths ago in two v. https://github.com/ggrrll/ndlib/branches

and we actually quickly chatted with @GiulioRossetti about it, some time ago...

Indeed I was wondering whether we should have an 'exponential' implementation (like in branch: cont_time:SEIRModel ) or rather have all transitions in rate-like implementation (like in branch: SEIRModel:disc_time )

mihaibanciu commented 4 years ago

I think it should be done. The Exponential looks so different than the Uniform distribution and the models are so sensitive to initial parameters so that the final counts per each state could look very different if this is not addressed.

BTW, in your branch you mention an error on L96 in the seir model. Is that a KeyError? I'm seeing the same. I wonder if one needs to initialize self_progress with t_i = 0 for all nodes that are in fraction_infected at iteration 0. See my first three suggested additions, above.

GiulioRossetti commented 4 years ago

Sorry for the late reply.

Indeed, it make sense to have both implementations in NDlib since we got several requests (both here and via email) asking for them.

I think that, at this point, the optimal solution will be to add continuous versions for both SEIR and SEIS models.

If you agree I'll give a look to your implementations to understand the nature of the error and move forward with the inclusion.

GiulioRossetti commented 4 years ago

@mihaibanciu @ggrrll check out the implementations in the dedicated branch.

If they are fine, I'll push them into master and include the in the next release (that I would like to release within this week).

@ggrrll: The error you faced was due to a missing initialization of infection seeds progress variable, as correctly pointed out by @mihaibanciu. To make things more readable I splitted the progress tracking in two disjoint dict.

mihaibanciu commented 4 years ago

@GiulioRossetti As far as I can tell the cont branch looks good, I just looked at the *_ct.py files. I noticed one typo in a file name: DynamicCompostiteModel.py --> DynamicCompositeModel.py (there's an extra 't' in Composite that should disappear).

GiulioRossetti commented 4 years ago

Ok that's something that can be easily fixed :-)