chaining confidence intervals from samples (with covid19 application)

Question

I've been dealing a lot with covid19 related data and it seems a lot of the work seems to calculate confidence intervals and pValue just for the step of the work they did while taking the previous results they rely on as exact. I'm hoping we can do better. And get better bounds by propagating uncertainty from one step to the next.

Specific motivating example (which of obvious interest in it's own right):

Serological test by BioMedics sensitivity - 88.66% specificity - 90.63% results came from: 397 PCR confirmed COVID-19 patients and 128 negative patients.

Chelsea study: 200 participants (generally appeared healthy off the street): 64 out of 200 tested positive for antibodies using the Biomedics test

Population of the town 40160 deaths at relevant time point 41

Can we chain everything to get X% confidence intervals for lethality rate? deaths/contracted.

(obviously we will still have to assume populations are similar and representative but that is another issue). Can we get good error bounds? when is it reasonable to ignore errors from previous step?

Answer 1

Firstly, I don't feel qualified to comment on the details of the COVID-19 models, so my comments are more about the general principles of using outcomes from one model for subsequent models.

Your intuition is certainly right that using point estimates ignores uncertainty and can thus lead to very wrong results, particularly, if unlikely but possible parameter values have very important implications E.g. if some parameter values conceivable based on the data, but away from the point estimate, would imply extreme outcomes that are extremely unlikely with the point estimates, then this would be important to know.

In principle, when using outputs of previous experiments / analyses / trials etc. into subsequent analyses, one of the most logical approach is using a Bayesian approach. I.e. you take the posterior distribution of one analysis as the prior for the next analysis.

In particular, if

1. a estimate +- standard error description is appropriate (or is appropriate after some distribution e.g. for the logit of a proportion or for log viral load or ...) for the first analysis,
2. if the knowledge about the other parameters in a model and the other parameters int he other models is independent (for example, it is likely a problem if the same data are used in both models), and
3. if there is no differences/generalization issues, then a Normal(mean=estimate, SD=standard error) distribution would be appropriate for the parameter of the next analysis.

Obviously, any of the above can be relaxed. E.g. allowing non-normal distributions can be quite easy in some cases (e.g. exploiting conjugate distributions, for a proportion one might assume a Beta(0.5, 0.5) prior and get a Beta(0.5+number of events, 0.5 + number of non-events) posterior). Using a joint model that can deal with the same data going into multiple of your chained models. When there's a question of whether data from one setting generalizes (e.g. laboratory testing versus what happens out in the real-world), there methods that attempt to deal with prior-data conflict or one can try to explicitly elicit the uncertainty from experts (and adjust the distribution accordingly).

However, the technical/statistical/implementation difficulty/potential for errors increases the more complicated we have to make things. However, there are a lot of tools (like the R and Python interfaces to Stan - rstan and pystan) that allow hand-crafted specification of almost arbitrarily complex models.

Answer 2

So I took this opportunity to learn Stan and some Bayesian modelling. This allows filling everything in into a model and let it do it's thing Here is the Stan model I used:

data {
  int<lower=0> Npos; 
  int<lower=0> Nneg;
  int<lower=0,upper=1> testPos[Npos];
  int<lower=0,upper=1> testNeg[Nneg];

  int<lower=0> Nc;
  int<lower=0,upper=1> chelsea[Nc];

  int<lower=0> Npop;
  int<lower=0,upper=1> deaths[Npop];


}
parameters {
    real<lower=0,upper=1> tpr;
    real<lower=0,upper=1> fpr;
    real<lower=0,upper=1> infectionRate;
    real<lower=0,upper=1> lethality;

}
model {   
  testPos ~ bernoulli(tpr);
  testNeg ~ bernoulli(fpr);
  chelsea ~ bernoulli(infectionRate*tpr + (1-infectionRate)*fpr);
  deaths ~ bernoulli(infectionRate*lethality);
}

This models the true positive rate (tpr) and False Positive rate (fpr) based on the Biomedics paper. We use these together with the chelsea study date to estiamte the infectionRate. And lastly we use this to estimate the mortality rate. We get confidence intervals for all paramaters estimated.

If we would have just used the mean estimates from previous step we would get similar estimation for all the means but the 97.5% quantile estimation for lethality rate would have been 0.48% without propagating uncertainty vs 0.63% with. So not as huge of a difference as I thought, but a nice effect and I learned a new tool. Thanks Stan.(and Björn for the reference)