Bayesian sequential approach to monitor COVID-19 variants through test positivity rate from wastewater

WW settled solids for the Davis wastewater treatment plant (WWTP) were collected daily from the primary clarifier, transported on the same day of collection to the analytical laboratory, and processed within 24 h as previously described (36). WW settled solids were obtained from daily composite influent samples from UC Davis WWTP. Composite influent samples were collected using a refrigerated autosampler (Hach Sigma 900 MAX) located at the WWTP headworks and programmed to collect flow-weighted influent sample volumes every 20 min for a total volume of 19 L in 24 h. Composite influent samples were then transferred to one or two 4 L low-density polyethylene containers (LDPE Cubitainers, Thermo Scientific I-Chem) and stored at 4°C prior to settling (up to 6 d of storage). Each 4 L sample was pasteurized in a 60°C water bath for 45 min immediately prior to settling. Pasteurized influent samples were inverted to mix, poured into a 3-gallon high-density polyethylene conical vessel equipped with a sampling port (FF3G, FastFerment), and left to settle for 2 h. Settled solids were obtained from either 4 L or 8 L of influent from a single day (two 4 L samples were combined into one settling vessel when 8 L was used). Settled solids were collected by dispensing from the bottom of the settling vessel into one or two 50 mL polypropylene centrifuge tubes (for 4 L or 8 L initial volume, respectively). If two tubes of settled solids were obtained, the supernatants were carefully decanted from each, and the remaining settled solids were combined. Between sampling episodes, the settling tanks were emptied, and the tank and sampling port valves were bleached (10% commercial bleach for 1 h), rinsed with deionized water, and left to air dry. Samples of settled solids were stored at 4°C and subsequently transported on ice in a cooler by a courier to the laboratory. Sample processing was completed within 8 d of initial sample collection. Sample RNA extraction, purification, and droplet digital reverse transcriptase PCR (ddRT-PCR) followed the same protocol as for the Davis samples. These protocols are described in detail elsewhere (37, 38).

We normalize the SARS-CoV-2 RNA concentration determined (N gene copies per gram dry weight solids) by the concentration of mild pepper mottle virus (PMMoV gene copies per gram dry weight solids) to yield the dimensionless metric, N/PMMoV. The N gene is present in all variants of the virus. PMMoV is a highly abundant RNA virus detected broadly in WW (39, 40). PMMoV serves as a process control such that normalization of N by PMMoV for each sample helps correct SARS-CoV-2 concentrations for virus extraction efficiency. PMMoV is also often used to account for variations in population size, rainfall, and water usage between different WW collection areas (41). We expect these latter factors to have less of an effect on N gene concentrations determined herein because water is removed from the WW settled solids samples prior to sample analysis and concentrations are reported in terms of the dry weight of the dewatered solids. Fig. 5 shows the normalized N gene concentration.

Given that WW signals are often noise corrupted, we applied a 7-d trimmed mean for daily WW data (smoothed N/PMMoV) to reduce uncertainty and minimize daily fluctuations. The smoothed data are later correlated to raw TPR (Fig. 6).

Table 2 illustrates Pearson’s correlation coefficient between the WW data (smoothed N gene and smoothed N/PMMoV) and clinical data (cases and TPR) in different waves. An improvement in correlation is observed when using TPR, compared with other cases. In general, an increase in correlation is also observed when using the normalized signal (N/PMMoV). Additionally, Table 2 highlights a notable change (from Omicron to BA.2–5 wave) in the correlation between WW and case data, while the correlation between WW data and TPR showed a relatively small change. These observations suggest that one of the primary factors responsible for the shift in the WW-to-cases relationship could be the reduction in the number of tests conducted during the second wave (see Fig. 5).

We model TPR as a beta distribution using WW viral loads as a covariate, assuming a Bayesian approach. We assume $Y_i$ as the TPR at day $i$ , defined as the ratio of the number of new positive cases among the number of tests performed at day $i$ . Beta regression is a good choice of model for continuous data with response variables expressed as proportions.

The beta distribution is reparametrized as $Y_i\sim ℬ({\mu }_i,\varphi )$ , with mean ${\mu }_i$ and variance ${\sigma }_i^2={\mu }_i(1-{\mu }_i)/(1+\varphi )$ , $\varphi$ is known as the precision parameter since, for fixed ${\mu }_i$ , the larger the $\varphi$ value, the smaller the variance of ${\displaystyle Y_i;{\varphi }^{-1}}$ as a dispersion parameter (42). The mean ${\mu }_i$ can be expressed as a function of the linear predictor ${\eta }_i={\beta }^T{\mathbf{x}}_i$ , where $\beta ={({\beta }_0,{\beta }_1,\mathrm{\dots },{\beta }_p)}^T$ is a $(p+1)$ -dimensional vector of unknown regression coefficients (including the intercept), and ${\mathbf{x}}_i={(1,x_1,\mathrm{\dots },x_p)}^T$ is the vector of covariates plus a one for the intercept. In this study, only the 7 d trimmed mean of the WW data (smoothed N/PMMoV), denoted as $C_i$ , is included in the linear predictor ( $p=1$ ). Thus, the linear predictor is given by ${\eta }_i={\beta }_0+{\beta }_1{C}_i$ , and the logit link, the inverse of the logistic function, is used in the beta regression (i.e., ${\displaystyle logit({\mu }_i)=\log \left(\frac{{\mu }_i}{1-{\mu }_i}\right)={\eta }_i}$ ).

The aim of this inference problem is to estimate $\theta =({\beta }_0,{\beta }_1)$ from measurements of WW data, ${\displaystyle \mathbf{C}=(C_1,C_2,\dots ,C_n)}$ , and COVID-19 TPR, ${\displaystyle \mathbf{Y}=(y_1,y_2,\dots ,y_n)}$ . Thus, the likelihood function for the previous model is given by:

where $a_i={\mu }_i\varphi$ , $b_i=\varphi -a_i$ , ${\mu }_i=1/(1+\exp (-{\eta }_i)$ is the mean (logistic function), ${\eta }_i={\beta }_0+{\beta }_1{C}_i$ is the linear predictor, and $\varphi$ is the dispersion parameter.

We adopt a Bayesian statistical approach, which is well suited to model multiple sources of uncertainty and allows for incorporating background knowledge of the model’s parameters. In this framework, a prior distribution, ${\pi }_{\mathrm{\Theta }}(\theta )$ , is required to account for the unknown parameter $\theta$ in order to obtain the posterior distribution. Having specified the likelihood and the prior, we use Bayes’ rule to calculate the posterior distribution,

where ${\displaystyle Z(\mathbf{Y})=\int {\pi }_{\mathrm{\Theta }}(\theta )L(\theta |\mathbf{C},\mathbf{Y})d\theta }$ is the normalization constant. The posterior distribution is simulated using an existing Markov chain Monte Carlo (MCMC) method, the t-walk algorithm (43).

We adapt the sequential approach proposed in reference 3 to our model to update forecasts over time. The aim is to train the model using only a subset of the most recent data. The forecast is then updated sequentially in a sliding window of data.

We let $L$ be the length in days of the period used to train the model. The data window is then moved forward every $n$ day as new data become available. We set t ₀ as the first initial time to start the analysis and the subsequent initial times as $t_{k+1}=t_k+n$ . The training period is taken as $[t_k,t_k+L]$ and the forecasting period as $[t_k+L,t_k+L+F]$ , see Fig. 7.

We denote ${\theta }_k=({\beta }_0^{(k)},{\beta }_1^{(k)})$ the model parameters to be inferred and the vectors of data as ${\displaystyle {\mathbf{C}}_{k,n}=(C_k,\dots ,C_{k+n}),{\mathbf{Y}}_{k,n}=(y_k,\dots ,y_{k+n})}$ at period $k$ . Note that, from the beginning, ${\theta }_k$ is assumed to change in time within each forecast window. If $k=0$ , we postulate a prior distribution ${\pi }_{{\mathrm{\Theta }}_k}({\theta }_k)$ and a likelihood $L({\theta }_k|{\mathbf{x}}_{k,n},Y_{k,n})$ previously described. The probabilistic prediction of y _t , in the forecasting period $t\in [t_k+L,t_k+L+F]$ , is obtained by using the estimated parameters through the MCMC method and the WW data concentration ( ${\displaystyle {\mathbf{C}}_{k+L,F}}$ ).

Afterward, the forecasting window is updated by setting $t_{k+1}=t_k+n$ , with $n$ as the number of days until the next forecast. In the new training window $[t_{k+1},t_{k+1}+L]$ , we propose a new prior distribution ${\pi }_{{\mathrm{\Theta }}_{k+1}}({\theta }_{k+1})$ for the model parameters ${\theta }_{k+1}$ using samples from the posterior distribution obtained in the previous forecast. Finally, we set $k=k+1$ and repeat the process described above to create a new forecast, see reference (3) for implementation details.

Regarding the elicitation of the parameters’ prior distribution for the first forecast, at $k=0$ , we assume a normal distribution for the parameters ${\beta }_0^{(k)}$ , ${\beta }_1^{(k)}$ , with mean and standard deviation $(0,0)$ and (1,000, 1,000), respectively (i.e., noninformative priors). We set $L$ to twice the length from symptoms onset to mild disease clinical outcome, namely 30 d, and $F$ is chosen to be 10 d. The forecasting was updated every $n=10$ d.

We let $Y$ be the TPR, and $C$ corresponds to the SARS-CoV-2 RNA concentrations measured in WW. Then, the cumulative distribution function of $Y$ given $C$ is defined as $F_{Y|C}(y|c):=P(Y\le y|C=c),$ which represents the probability that TPR is less than or equal to $y$ , given that the SARS-CoV-2 RNA concentration in WW is $c$ . Henceforth for simplicity, we will use $F(y|c)$ instead of $F_{Y|C}(y|c)$ without loss of generality. The quantile function $F^{-1}$ of $Y$ given $C=c$ is defined by

The $p$ -quantile of a data set is defined as the value where a $p$ fraction of the data is below that value and $(1-p)$ fraction of the data is above that value (e.g., the $0.5$ quantile is the median).

Using the CDC thresholds for TPR values corresponding to low ( $Y\le 0.05$ ), moderate ( $Y\in (0.05,0.08)$ ), substantial ( ${\displaystyle Y\in (0.08,0.1)}$ ), and high ( $Y\ge 0.1$ ) transmission levels of SARS-CoV-2 and the parameter estimates from the assumed beta regression model, we propose a methodology for estimating WW concentrations associated with TPR thresholds at a given point in time. We find the value of WW concentrations, $c$ , such that with probability $1-\alpha$ , the TPR is less than the CDC threshold $y\in \{0.5,0.8,0.1\}$ (i.e., to find $c$ such that $F_c^{-1}(\alpha )=y$ ). Note that $\alpha$ is the precision we set to estimate the threshold. With $\alpha =0.05$ , we are being conservative and choose lower bounds of WW concentration values associated with TPR thresholds.

The proposed method for finding these thresholds, assuming that we have simulations of the posterior distribution for the parameter $\theta =({\beta }_0,{\beta }_1)$ , is given in Algorithm 1. First, we suggest a search grid for the WW viral load concentration. Then, for each concentration, we simulate the predicted posterior distribution of the TPR. Lastly, we calculate the $\alpha$ -quantile for each concentration level and find the concentration level that yields the quantile closest to the value of the desired TPR threshold.

Input: A sample, ${\theta }_1,{\theta }_2,\mathrm{\dots },{\theta }_N$ , from the posterior distribution of $\mathrm{\Theta }$ , where       ${\displaystyle {\theta }_i=({\beta }_{0,i},{\beta }_{1,i}),i=1,\dots ,N}$ . The precision parameter $\varphi$ , the threshold for the       TPR, $T_{pr}$ , and the probability $\alpha$ ;

Output: Threshold for WW concentration $T_{ww}$ ;

Step 1. Generate a grid of WW concentration, $c_1,c_2,\mathrm{\dots },c_L$ ;

for ${\displaystyle l}$  in 1 to $L$ do

  for $i$  in 1 to $N$ do

   Step 2. Simulate $Y_{l,i}\sim Beta({\mu }_{l,i},\varphi )$ , where ${\mu }_{l,i}=\frac{1}{1+\exp \left(-{\eta }_{l,i}\right)}$ and      ${\eta }_{l,i}={\beta }_{0,i}+{\beta }_{1,i}{c}_l$ ;

  Step 3. Compute the $\alpha$ -quantile of ${\displaystyle {\mathbf{Y}}_l=(Y_{l,1},\dots ,Y_{l,N})}$ , namely $Q_l$ ;

Step 4. Find $l^{*}$ such that $l^{*}=\underset{x}{\mathrm{argmin}}\mathit{\hspace{1em}}|Q_l-T_{pr}|$ ;

Step 5. Set $T_{ww}=c_{l^{*}}$

The number of people in a population who are susceptible to infection by an infected individual at any particular time is denoted by $R_e$ , the effective reproductive number. This dimensionless quantity is sensitive to time-dependent variation due to reductions in susceptible individuals, changes in population immunity, and other factors. $R_e$ can be estimated by the ratio of the number of new infections ( $I_t$ ) generated at time $t$ , to the total infectious individuals at time $t$ , given by ${\sum }_{s=1}^t{I}_{t-s}{w}_s$ , the sum of new infections up to time step $t-1$ , weighted by the infectivity function w _s . We implement the method proposed in reference 44 to calculate the $R_e$ from the TPR estimated with the WW data.

Note that the TPR is an estimation of the proportion of infected persons. Therefore, if we multiply the TPR by a T value, representing the total number of tests carried out in the study period, we will have an estimate of the incidence, which we can use to compute the $R_e$ . Since $R_e$ is a scale-free metric, we should get similar results for different $T$ values. We use $T$ as the average of the tests carried out in Davis or UC Davis.