Lack of epidemiological metrics on netdiffuseR

[aoliveram]

netdiffuseR currently lacks metrics for epidemiological analysis of contact networks. The diffnet class treats adoption as permanent (like a SI model), hazard_rate() is network-blind, there is no way to store a transmission tree, and stochastic transmission (already in progress on the stochastic-transmission branch) supports only a sigmoid kernel (not well-suited for epidemiological simulations).

What we should do

1. Time of Disadoption. Add an optional $tod slot to diffnet parallel in shape to $toa. Both should be matrices N × K to handle k reinfection cycles, and be a vector when one episode per node.

2. Reproductive number. We can construct a native function repr_number() that mirrors epiworldR::get_reproductive_number() exactly: offspring count aggregated by (source, source_exposure_date) from the transmission tree. No other R_0 flavour in this PR.

3. Supporting epi metrics. secondary_attack_rate(), generation_time(), survival_curve(), peak_prevalence(), peak_time(). All native, and relatively simple once with the adoption tree.

4. Transmission tree. We should add an optional $transmission slot on diffnet (schema identical to epiworldR::get_transmissions()). Tested first with the Epigames histories.csv and checked if it fits well with the one provided by Andres's scripts.

5. Stochastic transmission. 1) Generalise the existing stochastic-transmission work: add a pluggable link_fun parameter to allocate a preferred function (linear, sigmoid, Wells-Riley, or a function given by the user); and 2) add a adoption_model = "logit" path in rdiffnet() with beta0 + beta_expo parameters.

Key aspects

• Zero new dependencies.
• Backward compatible: When $tod and $transmission are NULL, every existing function returns identical results.

Origin of prior work

This PR integrates two existing branches:

• stochastic-transmission (link-level Bernoulli, sigmoid kernel).
• issue-75-epigames-dynamic-attrs (Epigames dataset, dynamic attrs, collapse_timeframes).

Both will be merged into a new issue-XX-epidemiological-metrics branch.

[aoliveram]

hey @gvegayon, I've created a plan for implementing all the stuff. I can share a detailed summary if you want to double-check each point.

I'll open a draft PR and start implementing it. Please let me know if you want to discuss any point more deeply ;)

[gvegayon]

@aoliveram, my responses below

netdiffuseR currently lacks metrics for epidemiological analysis of contact networks. The diffnet class treats adoption as permanent (like a SI model), hazard_rate() is network-blind, there is no way to store a transmission tree, and stochastic transmission (already in progress on the stochastic-transmission branch) supports only a sigmoid kernel (not well-suited for epidemiological simulations).

What we should do

1. Time of Disadoption. Add an optional $tod slot to diffnet parallel in shape to $toa. Both should be matrices N × K to handle k reinfection cycles, and be a vector when one episode per node.

I think this feature (or something similar) already exists, and you implement it! Perhaps this could be simplified by having a matrix indicator of the current status of the agent or something? Instead of having multiple matrices, you could have just one that we can use to compute all the others.

3. Reproductive number. We can construct a native function repr_number() that mirrors epiworldR::get_reproductive_number() exactly: offspring count aggregated by (source, source_exposure_date) from the transmission tree. No other R_0 flavour in this PR.

For that, as you suggest, you would need to have a transmission tree. This is irrelevant in the current version of the model, since transmission is based on exposure, not individual contacts. You would need to have a full model that indicates who transmitted to who (like epiworldR), but that requires a different theoretical approach.

5. Supporting epi metrics. secondary_attack_rate(), generation_time(), survival_curve(), peak_prevalence(), peak_time(). All native, and relatively simple once with the adoption tree.

I think these generally should be easy to compute, based on the adopt counts.

7. Transmission tree. We should add an optional $transmission slot on diffnet (schema identical to epiworldR::get_transmissions()). Tested first with the Epigames histories.csv and checked if it fits well with the one provided by Andres's scripts.

Perhaps you could apply something like inheritance here. This could be a new class of diffnet object called, for example, diffnet_epi, that allows adding transmission trees, for example, and would have other associated methods. In R, you can have objects with multiple classes, so this would be something like class(x) <- c("diffnet_epi", "diffnet") so dispatch would be first checking for diffnet_epi class and then for diffnet.

9. Stochastic transmission. 1) Generalise the existing stochastic-transmission work: add a pluggable link_fun parameter to allocate a preferred function (linear, sigmoid, Wells-Riley, or a function given by the user); and 2) add a adoption_model = "logit" path in rdiffnet() with beta0 + beta_expo parameters.> > ## Key aspects> * Zero new dependencies.> * Backward compatible: When $tod and $transmission are NULL, every existing function returns identical results.

I think it could be even more general. The deterministic model (current) is just a degenerate random variable, so you could have instead a adoption_mechanism which would be a function that receives (a) the current network object, (b) the set of adopters, (c) etc. We could have some predefined like adoptmech_threshold (for Tom's model), adoptmech_logit for the logit model, etc.

Origin of prior work

This PR integrates two existing branches:

• stochastic-transmission (link-level Bernoulli, sigmoid kernel).
• issue-75-epigames-dynamic-attrs (Epigames dataset, dynamic attrs, collapse_timeframes).

Both will be merged into a new issue-XX-epidemiological-metrics branch.

[aoliveram]

Thank you @gvegayon! Let me comment on some of your points:

@aoliveram, my responses below

netdiffuseR currently lacks metrics for epidemiological analysis of contact networks. The diffnet class treats adoption as permanent (like a SI model), hazard_rate() is network-blind, there is no way to store a transmission tree, and stochastic transmission (already in progress on the stochastic-transmission branch) supports only a sigmoid kernel (not well-suited for epidemiological simulations).

What we should do

1. Time of Disadoption. Add an optional $tod slot to diffnet parallel in shape to $toa. Both should be matrices N × K to handle k reinfection cycles, and be a vector when one episode per node.

I think this feature (or something similar) already exists, and you implement it! Perhaps this could be simplified by having a matrix indicator of the current status of the agent or something? Instead of having multiple matrices, you could have just one that we can use to compute all the others.

Yes, there is a disadopt entry on rdiffnet, but it only helps during simulations. So, to generalize it to the as_diffnet constructor, I think the tod slot is the simplest way to do it. Under the hood, the information could be fused into $status, from which everything is derived. It should be possible for the user to give the status matrix as well:

# Giving tod
toa_vec <- c(1, 2, 3, NA, 4)
tod_vec <- c(5, 4, NA, NA, 6) 

dn <- as_diffnet(graph = g, toa = toa_vec, tod = tod_vec)
# status[i, t, 1] = 1 ssi toa[i] <= t <= tod[i] - 1

# Giving status matrix (n = 5, T = 10, Q = 1 behaviors)
status <- array(0L, dim = c(5, 10, 1))
status[1, 1:4,   1] <- 1L
status[2, 2:3,   1] <- 1L
status[4, c(4:5, 7:9), 1] <- 1L  # multi-cycle

dn <- as_diffnet(graph = g, status = status)
# $toa y $tod are derived from status

This should be possible both for epidemiological and behavioral sims.

2. Reproductive number. We can construct a native function repr_number() that mirrors epiworldR::get_reproductive_number() exactly: offspring count aggregated by (source, source_exposure_date) from the transmission tree. No other R_0 flavour in this PR.

For that, as you suggest, you would need to have a transmission tree. This is irrelevant in the current version of the model, since transmission is based on exposure, not individual contacts. You would need to have a full model that indicates who transmitted to who (like epiworldR), but that requires a different theoretical approach.

4. Transmission tree. We should add an optional $transmission slot on diffnet (schema identical to epiworldR::get_transmissions()). Tested first with the Epigames histories.csv and checked if it fits well with the one provided by Andres's scripts.

Perhaps you could apply something like inheritance here. This could be a new class of diffnet object called, for example, diffnet_epi, that allows adding transmission trees, for example, and would have other associated methods. In R, you can have objects with multiple classes, so this would be something like class(x) <- c("diffnet_epi", "diffnet") so dispatch would be first checking for diffnet_epi class and then for diffnet.

Interesting. This could avoid adapting the working rdiffnet to the new epi necessities, without having to modify several lines. I'll update the milestones table in the PR.

5. Stochastic transmission. 1) Generalise the existing stochastic-transmission work: add a pluggable link_fun parameter to allocate a preferred function (linear, sigmoid, Wells-Riley, or a function given by the user); and 2) add a adoption_model = "logit" path in rdiffnet() with beta0 + beta_expo parameters.> > ## Key aspects> * Zero new dependencies.> * Backward compatible: When $tod and $transmission are NULL, every existing function returns identical results.

I think it could be even more general. The deterministic model (current) is just a degenerate random variable, so you could have instead a adoption_mechanism which would be a function that receives (a) the current network object, (b) the set of adopters, (c) etc. We could have some predefined like adoptmech_threshold (for Tom's model), adoptmech_logit for the logit model, etc.

I get you. So, inside rdiffnet, the entry adoption_mechanism would be a mirror of disadopt for adoption. This is what I have now:

rdiffnet(..., adoption_model = c("deterministic", "stochastic"),
adoption_pars = NULL)

# Inside rdiffnet ...
if (adoption_model == "stochastic") {
p <- plogis(adoption_pars$beta0 + adoption_pars$beta_expo * expo[, , q])
whoadopts <- which(runif(length(p)) < p & is.na(toa[, q]))
} else {
whoadopts <- which((expo[, , q] >= thr[, q]) & is.na(toa[, q]))
}

Following your comment, it could become something like this:

# pre-defined mechanism

adoptmech_threshold <- function(expo, thresholds, not_adopted, time, pars) { 
which(expo >= thresholds & not_adopted)
}
adoptmech_logit <- function(expo, thresholds, not_adopted, time, pars) { 
p <- plogis(pars$beta0 + pars$beta_expo * expo) 
which(runif(length(p)) < p & not_adopted)
}
adoptmech_probit <- function(expo, thresholds, not_adopted, time, pars) { 
p <- pnorm(pars$beta0 + pars$beta_expo * expo) 
which(runif(length(p)) < p & not_adopted)
}

# user written:

my_mechanism <- function(expo, thresholds, not_adopted, time, pars) {
... a function...
}

# rdiffnet now receives a function:

rdiffnet(..., adoption_mechanism = adoptmech_logit,
adoption_pars = list(beta0 = -2, beta_expo = 5))

# Inside rdiffnet ...
whoadopts <- adoption_mechanism(
expo = expo[, , q], 
thresholds = thr[, q], 
not_adopted = is.na(toa[, q]), 
time = t, 
pars = adoption_pars
)

3. Supporting epi metrics. secondary_attack_rate(), generation_time(), survival_curve(), peak_prevalence(), peak_time(). All native, and relatively simple once with the adoption tree.

I think these generally should be easy to compute, based on the adopt counts.

I agree, these shouldn't be problematic.

[aoliveram]

@gvegayon Let me know your thoughts!

[gvegayon]

Agree with all your comments, @aoliveram, except the first one. Throughout the years, I have been naive in thinking that people would be happy with many ways of doing the same thing. But I have realized that it is way better to have a single way of doing things like this. Particularly, I think the idea of having a status argument would be the best, let's avoid adding the other td vector, since it can create confusion and more room for bugs.