Example: serial dilution assay

1. Experimental setup and problem

A serial dilution assay measures compound concentration using an optical color readout.
The color response is informative only in a limited concentration range:

• At low concentration: essentially no color change
• At high concentration: color saturates

To get usable measurements, each sample is measured at several dilutions.

A typical plate has 96 wells:

• First columns: standard sample with known concentration and its dilutions
• Remaining columns: unknown samples and their dilutions

For the example plate:

• The standard occupies the first two columns (many closely spaced dilutions)
• Ten unknown samples occupy the remaining columns (fewer, more widely spaced dilutions)

The goal is to estimate concentrations of the 10 unknown samples, using the standard data for calibration and properly using all information, including measurements labeled “below detection limit” by standard software.

---

2. Problems with the standard analysis

Standard analysis typically:

1. Fits a calibration curve using the standard wells
2. Inverts this curve to turn each optical reading from an unknown into a concentration estimate
3. For each unknown, averages the back-calculated concentrations (after adjusting for dilution)
4. Treats low readings as “below detection limit” and often discards them

In the example:

• For Unknown 8, all readings are flagged “below detection limit,” so no estimate is produced.
  Yet the measurements decrease systematically across dilutions (1, 1/3, 1/9, …), which indicates a real signal.
• For Unknown 9, only four of eight measurements are above detection limit. The others are discarded, even though they also show a decreasing pattern with dilution.

Discarding these measurements wastes information, especially for low concentrations that are scientifically important (e.g., allergens).

A Bayesian model can use all measurements, including those below the arbitrary detection limit, and properly quantify uncertainty.

---

3. Model specification

3.1 Parameters and notation

• Unknown concentrations: $\theta_1, \dots, \theta_{10}$​ (one for each unknown sample)
• Known standard concentration: $\theta_0$​
• Plate wells: indexed by $i = 1, \dots, 96$
  • $x_i$​: true concentration in well $i$
  • $y_i$​: observed optical intensity in well $i$
  • $j(i) \in \{0,1,\dots,10\}$: which sample (standard or unknown) well $i$ corresponds to
  • $d_i$: dilution factor for well $i$ relative to that sample’s “initial” concentration

3.2 Calibration curve $E(y \mid x, \beta)$

Use the standard four-parameter logistic-type curve:$$E(y \mid x, \beta) = g(x, \beta) = \beta_1 + \frac{\beta_2}{1 + (x / \beta_3)^{\beta_4}} \quad (19.1)$$

Interpretation:

• $\beta_1$​: baseline intensity as concentration $x \to 0$
• $\beta_2$​: increase from baseline to saturation
• $\beta_3$​: concentration where the curve bends (midpoint / turning point)
• $\beta_4$​: steepness of transition (how quickly saturation occurs)

All $\beta_k \ge 0$. This is equivalent to a shifted/scaled logistic function of $\log x$.
The fitted curve matches the standard data well (as seen in the original figure).

3.3 Measurement error model

Intensity measurements are noisy, with larger variance at higher intensity. Model:$$y_i \sim N\big(g(x_i, \beta),\ \sigma_y^2 \, (g(x_i, \beta)/A)^{2\alpha}\big) \quad (19.2)$$

• $\sigma_y$​: typical error SD
• $\alpha \in [0, 1]$: controls how variance grows with intensity
• $A$: fixed scale constant (set to 30, a mid-range value of $y$)

Special cases:

• $\alpha = 0$: constant variance (homoscedastic normal errors)
• $\alpha = 1$: approximately constant variance on the log scale of $y$

Getting this variance pattern approximately correct matters because many observations are at low intensities; the model must use them without overstating precision.

3.4 Dilution and initial dilution error

There are two sources of potential error:

1. Error in the initial preparation of the standard solution
2. Error in the serial dilutions

For this dataset:

• Serial dilution error is small and ignored
• Initial dilution error for the standard is explicitly modeled

Let:

• $\theta_0$​: true concentration of the stock standard
• $d_0^{\text{init}}$​: known initial dilution factor used to prepare the standard’s first working solution
• $x_0^{\text{init}}$: actual (unknown) concentration after this initial dilution

Model the initial standard concentration on the log scale:$$\log x_0^{\text{init}} \sim N\big(\log(d_0^{\text{init}} \, \theta_0),\ \sigma_{\text{init}}^2\big) \quad (19.3)$$

For unknown samples:

• No initial dilution is applied in this model
• Their initial concentration equals their true concentration:

$$x_j^{\text{init}} = \theta_j, \quad j = 1,\dots,10.$$

Then every well concentration is:$$x_i = d_i \, x_{j(i)}^{\text{init}} \quad (19.4)$$

where $d_i$​ is the known dilution factor for well $i$, and $j(i)$ identifies which sample (standard or which unknown) that well belongs to.

$\sigma_{\text{init}}$​ is not identifiable from a single plate, so it is fixed at 0.02 (2% SD) based on prior analyses.

3.5 Prior distributions

Use weak, noninformative priors for the calibration parameters:

• $\log \beta_k \sim \text{Uniform}(-\infty, \infty)$ for $k = 1,\dots,4$
• $\sigma_y \sim \text{Uniform}(0, \infty)$
• $\alpha \sim \text{Uniform}(0, 1)$

For unknown concentrations:

• No pooling: $p(\log \theta_j) \propto 1$, for each $j = 1,\dots,10$

A hierarchical model such as $\log \theta_j \sim N(\mu_\theta, \sigma_\theta^2)$, or a mixture model that allows some true concentrations to be exactly zero, would be more realistic. For simplicity, this example uses the no-pooling model.

---

4. Inference strategy and computation

The model is implemented in BUGS (predecessor of Stan):

• Two parallel Gibbs chains
• 50,000 iterations each
• Convergence monitoring via potential scale reduction $\hat{R} < 1.1$ for all parameters
• To reduce storage, every 20th draw is kept

To improve convergence:

• Work with log-parameters: $\log \beta_k$​, $\log \theta_j$
• Reparameterize in terms of $\gamma_j = \theta_j / \beta_3$ and sample $\log \gamma_j$​ instead of $\log \theta_j$
  This reduces posterior correlation between unknown concentrations and $\beta_3$​ (curve’s x-location).

---

5. Results for calibration and variance

Posterior medians (50% intervals):

• $\beta_1 = 14.7 \, [14.5, 14.9]$
• $\beta_2 = 99.7 \, [96.8, 102.9]$
• $\beta_3 = 0.054 \, [0.051, 0.058]$
• $\beta_4 = 1.34 \, [1.30, 1.38]$

Variance parameters:

• $\sigma_y = 2.2 \, [2.1, 2.3]$
• $\alpha = 0.97 \, [0.94, 0.99]$

The calibration curve with these parameters fits the standard data well.
The high precision of the assay (clear replicates) allows precise estimation of curve and variance parameters from one plate.

---

6. Results for unknown concentrations

Posterior medians and central intervals for the 10 unknown concentrations are obtained directly from the posterior samples of $\theta_j$​. All 10 unknowns have estimates, including:

• Unknown 8: all measurements nominally “below detection limit” under the standard method, but the Bayesian model still infers a nonzero concentration using the decreasing pattern across dilutions.

Residual diagnostics (standardized residuals vs expected values) show no major lack of fit.

---

7. Comparison to the standard method

Standard (non-Bayesian) pipeline:

1. Fit calibration curve using the standard data
2. Convert each unknown’s measurement into a concentration by inverting the curve
3. For each unknown, divide each concentration estimate by its dilution and average
4. Discard or flag “below detection limit” values and do not use them

Example for Unknown 9 (standard method):

• Four dilutions produce measurable values above detection limit
• Four higher dilutions are “below detection limit” and ignored
• A typical calculation: Estimated concentration
  = ¼ (0.040 + 0.031 + 3·0.005 + 3·0.005) = 0.025

Issues with that approach:

• Wastes information from low but informative readings
• Produces no estimate for unknowns with all readings below detection limit (e.g., Unknown 8)

Bayesian approach:

• Uses all measurements, including very low readings
• Propagates uncertainty through the full joint model
• Produces concentration estimates for all unknowns
• Shows better reliability (agreement between estimates from top vs bottom half of the plate) than the standard method

A broader comparison (across multiple plates and conditions) confirms that the Bayesian method is more robust and informative, especially when concentrations are low and detection limits are a serious issue.