Support multiple sampling sites

[seabbs]

Do you have a clear view of how hard it would be to support multiple sites of data at once here?

Looking at the model code it looks like doing something fairly simple where you assume the difference is just some kind of observation error (potentially AR structured or similar) would be relatively easy?

[adrian-lison]

Especially when assuming a shared Rt and infection trajectory, this would be mostly straightforward.

A small challenge when not being tightly integrated with case count data is to find the right weights for the different sites. One can try to choose them such that the load per case in the total population (all catchments summed up) for each site is proportional to the relative population size of the catchment. This however makes the assumption that the recovery of genetic material is identical across sites. If catchment-specific case count data is provided one could calibrate to that.

Another slightly complication is that when using the dPCR-specific observation model, one either has to assume identical lab parameters (i.e. if all sites are processed by the same lab) or also stratify across these.

Now that the single-site model is rather mature, I am motivated to support this not too far in the future.

[seabbs]

Especially when assuming a shared Rt and infection trajectory, this would be mostly straightforward.

This is my preferred approach to thinking about this in the first instance. @kaitlynjohnson and co had a variant of their model that looked like this but it was dropped due to lack of resource.

. If catchment-specific case count data is provided one could calibrate to that.

Agree that this would be useful but the linkage this implies is beyond the simple feature I was suggesting here.

Another slightly complication is that when using the dPCR-specific observation model, one either has to assume identical lab parameters (i.e. if all sites are processed by the same lab) or also stratify across these.

Either of which would make sense as options right? Both of these approaches would be possible in the current framework without new features?

Now that the single-site model is rather mature, I am motivated to support this not too far in the future.

Nice. My personal preference is to fit the weights as a intercept (all the sites are smaller than the total pop) and a random effect (sites have varying catchments) and strongly suggest these are given as informative priors. Then as you say you need support for different labs and to additionally allow sites to randomly vary to unmodelled factors.

[kaitejohnson]

A small challenge when not being tightly integrated with case count data is to find the right weights for the different sites. One can try to choose them such that the load per case in the total population (all catchments summed up) for each site is proportional to the relative population size of the catchment. This however makes the assumption that the recovery of genetic material is identical across sites. If catchment-specific case count data is provided one could calibrate to that.

I think I might be missing something here, but regarding this point, what is the need for the relative weights of each catchment if the assumption is that all sites have the same Rt and infection trajectory? Couldn't you just estimate a site-specific recovery of genetic material which would account for differences in the magnitude of observed concentrations across the sites?

I really like the idea of stratifying across the lab parameters. We've thought about this, as we currently treat each lab-site pair as an independent group. Would be nice to learn lab parameter specific properties in labs processing samples from different sites.

[Orthomyxoviridae]

Hey from Chicago - we're also dipping our toes into sewer water collection.

We might be large enough where just performing this type of analysis for each individual sewershed makes sense, but I think we're also interested in constructing a regional picture and weighting appropriately across sites.

[adrian-lison]

I think I might be missing something here, but regarding this point, what is the need for the relative weights of each catchment if the assumption is that all sites have the same Rt and infection trajectory? Couldn't you just estimate a site-specific recovery of genetic material which would account for differences in the magnitude of observed concentrations across the sites?

@kaitejohnson The reason why I was thinking about weights is that I would like to distinguish between recovery and population sizes. If one of the sites makes up a large share of the population but for some reason has pretty bad recovery and therefore low concentrations with lots of noise, I would still like the total infection trajectory to be dominated by that site. In this example, that would probably mean that Rt estimates remain rather uncertain even if the other much smaller sites have better recovery with stronger signal. This is because I think such Rt estimates would be interpreted as representative of the overall population, so we should account for the representativeness of the sites... Does that make sense?

Nice. My personal preference is to fit the weights as a intercept (all the sites are smaller than the total pop) and a random effect (sites have varying catchments) and strongly suggest these are given as informative priors.

@seabbs When doing this, do you think it would be important to use a simplex to enforce that weights sum up to 100% of infections or rather not?

[kaitejohnson]

Ah yes this does make sense, I had assumed that the relative weights of each catchment would be by population size of the catchment if you were to sum them to get a total infection trajectory.

From each catchment's infections per capita, I think you could still estimate a catchment specific recovery and have that be part of the observation process. I think it just depends on at what level you are calculating the total infections, but I don't think it should be from the total genomes, right?

So in your model, you are suggesting that if there are $K$ subpops of population sizes $n_k = {n_1, n_2,... n_k}$ then:

$$I_{total} = \sum{I_k}$$

where

$$I_k = n_k \frac{I_{total}}{\sum{n_k}}$$

Since we you would impose the same infections per capita on each catchment area. Then I think your recovery is independent so you would have:

$$g_k(t)= r_k\sum_{\tau = 0}{I_k(t-\tau)s(\tau)}$$

Where $r_k$ is independently estimated "recovery" (genomes per infection that get captured in wastewater), $g_k$ is the total genomes collected in that catchment on that day, $s(\tau)$ is normalized shedding kinetics. This is for the simplified version where you assume number of genomes shed per infection is constant. Then you could also add a catchment level observation error term, similar to what we do https://github.com/CDCgov/ww-inference-model/blob/main/model_definition.md#viral-genome-concentration-observation-model (where $ij$ is lab-catchment-area combination, which is our observation unit, though I think the better way to do this would be to estimate catchment area and lab specific recovery and noise parameters...)

[seabbs]

When doing this, do you think it would be important to use a simplex to enforce that weights sum up to 100% of infections or rather not?

No I think its not in most settings unless you have some certainty its likely to be the case.

[seabbs]

Especially when assuming a shared Rt and infection trajectory, this would be mostly straightforward.

@adrian-lison and others any interest in picking this up as a potential feature? In principle I am happy to take a first pass at implementing.