Bayesian Hierarchical Models for Sparse Data: A Practical Guide to Scale-Up Forecasting

Predicting scale-up to industrial production performance

Author — Paul Piho

Across process industries, drug development, and materials engineering, the same problem recurs: a long list of candidate configurations must be narrowed to a shortlist of ones worth committing to pilot and production. Full-scale evaluation of every candidate is prohibitive, so the shortlist has to be produced from bench-scale data alone, even when the candidates themselves (new cell lines, product formulations, compounds, ore samples, battery cell chemistries, alloys) lack data about their production scale feasibility.

Bench-scale experiments on a new candidate are fast and cheap. The goal is a model that learns the shared scale-up structure from this history and applies it to a new candidate's bench data, predicting feasibility before a dedicated pilot run.

Standard curve-fitting cannot do this. A least-squares fit on the new candidate's bench data returns a point estimate with a standard error, but offers no robust mechanism for predicting the pilot and production scale performance. Extrapolating to production yields either a wrong number or one with a misleading confidence interval. That gap between reported and real uncertainty is where wrong screening calls get made: greenlighting a candidate that will not scale, or discarding one that would. Committing the wrong candidate to pilot and production is expensive: each scale-up stage ties up equipment, materials, and months of runtime before the failure shows.

Bayesian hierarchical models close that gap. They return posterior distributions instead of point estimates, so the uncertainty in each candidate's forecast is visible. They let data-rich candidates inform data-poor ones without forcing them to agree, and the same mechanism lets dense small-scale data tighten a sparse large-scale forecast. In this post we introduce the methodology and apply these techniques to a scale-up problem inspired by cell culture in bioreactors, an application of Bayesian modeling to industrial scale-up.

What is a posterior distribution, and why does it beat a point estimate?

Bayes' rule combines a prior — what you knew before the data — with a likelihood that reweights it by how well each parameter pair explains the observations:

$$\underbrace{p(\theta \mid \text{data})}_{\text{posterior}} \;\propto\; \underbrace{p(\text{data} \mid \theta)}_{\text{likelihood}}\, \underbrace{p(\theta)}_{\text{prior}}$$

The output is a posterior distribution over the model parameters $\theta$ (for a textbook treatment of Bayesian modelling see Gelman et al. (2013)). The posterior in the equation above rarely has a closed form, so it is sampled rather than solved. Markov chain Monte Carlo (MCMC) draws a sequence of parameter values whose visiting frequency matches the posterior density: regions the data and prior favour get sampled often, unlikely regions rarely. The collection of draws summarises the posterior: histograms, medians, and highest-density intervals (HDIs) are read straight off it. An HDI is the narrowest parameter range holding a given share of the posterior mass; a 94% HDI is used throughout for parameter estimates.

In this post we use a first-order kinetic reaction model of saturating response:

$$y(t) = y_{\max}\,\bigl(1 - e^{-k\,t}\bigr) \cdot 100\%,\qquad y_i^{\text{obs}} = y(t_i) + \varepsilon_i,\; \varepsilon_i \sim \mathcal{N}(0,\sigma)$$

The model is parametrised by the rate of saturation $k$ and the maximum saturation value $y_{\max}$. The physical meaning is domain-specific. The parameter $k$ wraps everything that controls how fast the process saturates.

Note that the maximum saturation $y_{\max}$ and the rate $k$ are independent: any rate $k$ can reach any ceiling $y_{\max}$. However, data gathered while the curve is still climbing only informs the initial slope (the product $y_{\max}k$). A high ceiling approached slowly (large $y_{\max}$, small $k$) traces almost the same early points as a low ceiling approached quickly (small $y_{\max}$, large $k$).

To contrast the Bayesian fit with a least-squares maximum likelihood estimate (MLE), let us observe the response at three early times — while the curve is still climbing, before it reaches its plateau — and fit those three points two ways. Both methods produce a parametrisation of the curve that fits the data, but they disagree on uncertainty. Least squares returns a single best fit. The uncertainty can be gauged by a parametric bootstrap: estimate the noise from the fit residuals, draw fresh datasets from the fitted curve at that level, refit each, and read off the spread. A Bayesian fit on the same three points with deliberately uninformative prior (flat on the curve parameters over their physical ranges and a Jeffreys prior on the noise) returns instead a posterior distribution over curves. Least squares overestimates the maximum saturation with high confidence, while the Bayesian fit returns wide posterior distributions for the parameters. The Bayesian fit is consistent with our analysis: the early data cannot reliably pin the plateau.

Figure 1. The least-squares and Bayesian curve fit contrast; both methods are given the same physical constraint $y_{\max}\le 100\%$. Top row — a saturating response measured at three early times (dots), true curve dashed, extrapolation region shaded. Left: least squares estimates the ceiling $y_{\max}$ to be close to the physical boundary constraint (94%). Right: the Bayesian posterior (median + 94% band) is a fan of curves that widens where the early data cannot pin the plateau and contains the true curve. Bottom row — the fitted parameter distributions behind those bands. Left: the joint $(k, y_{\max})$ cloud. Least squares (bootstrap, red) collapses against the ceiling without seeing the curved rate–ceiling trade-off; the Bayesian posterior (blue) traces the ridge below it. Right: the ceiling marginal. The least-squares point estimate is pinned at the $100\%$ boundary and its bootstrap piles up against that wall (94% $[74\%,100\%]$) while the posterior spreads smoothly, mean $82\%$, covering the true $80\%$ (94% $[42\%,100\%]$).

For the same single-group setup we additionally ask how the posterior responds to data volume. A dense design places 12 points spanning $t = 3$ to $t = 150$; a sparse design keeps just three, at $t = 2.4, 30, 120$. Both span the full window, so the ceiling is observed either way. But the dense design tiles the rising limb and repeats near the plateau, averaging each noisy point down by its neighbours. The sparse design pins the rate on just two early points ($t = 2.4, 30$) and the whole ceiling on the single point at $t = 120$. Fewer points buy a wider posterior: the sparse $k$ posterior comes out roughly $2.4\times$ wider (94% HDI) than the dense one, sparse $y_{\max}$ roughly $2.1\times$ wider.

Figure 2. Posterior medians and 94% HDI for a dense (12-point) and sparse (3-point) fit of the same process. Bottom row overlays the prior (grey) on each posterior. The dense fit results in visibly tighter posterior distributions.

From a single fit to a hierarchy

Figure 2 describes a simplified problem: strip a fit down to three points and its posterior widens, and the only remedy shown so far is more data. However, experimental campaigns can be slow and costly to supply. There are two levers that remain when observations can't grow. The first is the prior itself: the fits so far used deliberately uninformative priors to keep the least-squares contrast fair, but domain expertise, such as a plausible rate range, a mechanistically bounded ceiling, a known noise floor, can be encoded directly. This tightens the posterior with no new data. The downside is that with sparse data the prior dominates, so an informative prior must be defensible; a wrong one buys false confidence.

The second lever, and the focus of the rest of this post, is structure shared across groups: a data-poor group borrows from data-rich ones that follow the same process. The forward model stays the first-order saturating response introduced previously with Gaussian observation noise, parametrised by the rate $k$ and the ceiling $y_{\max}$. What changes from here on is not the forward model but how information is shared: first hierarchically across many groups, then across measurement scales by treating scale as a covariate.

Hierarchical Partial pooling across groups

Take a toy example of group-level Bayesian analysis where eight clones run through the same cell-culture platform, each with its own specific productivity, growth rate, and metabolic burden. The simplest approach fits each independently (no pooling): a noisy outlier in a sparsely sampled experiment pulls the estimated titer rate for the clone far from where the rest of the population would put it. A hierarchy adds a layer saying that every clone's rate is drawn from a shared distribution, itself estimated from the data. Thus a data-poor clone is informed by the data-rich ones it resembles. This is the partial-pooling idea that goes back to the shrinkage estimators of Stein (1956) and James & Stein (1961), with the empirical-Bayes treatment of Efron & Morris (1975).

Figure 3. Pooling strategies as graphical models. (A) No-pool MLE fits each group by least squares — point estimates with no shared structure, so a sparse group overfits and its standard error understates the real uncertainty. (B) Complete pooling forces a single $\theta$ on every group, erasing genuine group-to-group variation. (C) Independent Bayes returns a per-group posterior but still learns nothing across groups. (D) The hierarchical model draws each $\theta_g = (k_g, y_{\max,g})$ from a shared hyperprior linked to group features, so sparse groups draw on dense ones and a new group can be scored from its features alone. Yellow: latent Bayesian parameters; blue: MLE point estimates; grey: observed data; pink: shared hyperprior.

The hierarchical model's parameters and priors:

Each group's rate is drawn from a feature regression on log-rate and its ceiling from a shared Beta hyperprior.

$$\begin{aligned} \log k_g &= \alpha + \mathbf{x}_g^\top \boldsymbol{\beta} + \varepsilon_g,\; \varepsilon_g \sim \mathcal{N}(0, \sigma_{\text{grp}}) \\ y_{\max,g} &\sim \text{Beta} \bigl(\mu_{y_{\max}}\kappa_{y_{\max}},\, (1-\mu_{y_{\max}})\kappa_{y_{\max}}\bigr) \\ y_i^{\text{obs}} &\sim \mathcal{N} \bigl(y(t_i;\, k_{g(i)},\, y_{\max,g(i)}),\,\sigma\bigr) \end{aligned}$$

$\sigma_{\text{grp}}$ controls how strong the pull toward the population is: small $\sigma_{\text{grp}}$ collapses toward complete pooling on the feature regression, large $\sigma_{\text{grp}}$ relaxes toward independent per-group fits.

The plate diagram below sets out the full structure of how the priors, structural equations, and likelihood connect.

The outer plate spans groups $g = 1,\ldots,G$; the inner plate spans observations $i = 1,\ldots,n_g$ within each group. Grey nodes are the observed responses and the green node holds the fixed group features $\mathbf{x}_g$; the remaining nodes are parameters drawn from the shared hyperpriors at the top.

Figure 4. Hierarchical pooling on eight clones, one of which has only three observations. No-pool least-squares MLE 94% CIs (purple) vs hierarchical 94% HDIs (blue); ground truth marked. Because the clone features vary near-independently, the covariate-effect panel recovers each regression coefficient $\beta_j$ within its 94% HDI (truth marked), so the fitted rate can be attributed to individual features rather than an unidentified blend of them.

Adding the shared population layer leaves the three-observation group's rate interval about $3.7\times$ tighter than a plain least-squares fit on its own data, and nudges its estimate from 0.096 to 0.074, toward the truth at 0.073. The tightening concentrates where it is needed: the dense groups barely move, their own data already pinning the rate, while the sparse group draws from the population to inform the fitted parameters.

Where Figure 4 compares interval widths, Figure 5 shows the same pull as a movement of the point estimate. Each clone's least-squares rate is pulled toward its feature-predicted value.

Figure 5. Shrinkage of the rate $k$. Because the hierarchy carries a feature regression, each clone's estimate is pulled toward its feature-predicted rate $\hat k_g = \exp(\alpha + \mathbf{x}_g^{\top} \boldsymbol{\beta})$ (the dashed diagonal), not a single grand mean. For every clone the no-pool least-squares MLE (purple) and the hierarchical mean (blue) share an $x$; the vertical drop between them is the shrinkage, green where pooling moved the estimate toward ground truth ($\diamond$). Dense clones already sit on the line and barely move; the sparse clone's MLE ($0.096$) sits well above it, and pooling pulls it back toward the prediction ($0.074$). Predictive strength is gained from its siblings, not from data it does not have. Clone H is sampled once near baseline and twice after it has saturated, so its three points barely constrain the rate: the MLE's parametric-bootstrap 94% CI (purple, computed with the known noise level $\sigma$, as the Bayesian fits are) runs to $[0.06, 0.23]$, its upper end clipped at the axis and labelled. The hierarchical 94% HDI (blue) is several times tighter and recentred near the truth, because the population supplies what H's own data cannot.

Once $\hat\alpha$ and $\hat{\boldsymbol{\beta}}$ are estimated, the same regression predicts $k$ for a new group directly from its feature vector, with no observations of its own:

$$\hat k_{\text{new}} \;=\; \exp \bigl(\hat\alpha + \mathbf{x}_{\text{new}}^{\top}\hat{\boldsymbol{\beta}}\bigr).$$

This is the mechanism that lets a candidate be scored before any dedicated campaign is run.

Forecasting production before production data

Continue with biologics cell culture as the running example. The same clone is grown at three scales, each layering new physics on the last: the bench-bioreactor isolates the clone's intrinsic growth and production kinetics under tightly controlled, well-mixed conditions; the pilot adds finite oxygen transfer and the first mixing and dissolved-gas gradients as the vessel grows; the production bioreactor (2000 L) adds imperfect mixing, dissolved-oxygen and pH gradients, and CO₂ accumulation. Each scale lowers the titer a clone can reach and slows how fast it gets there.

Generating data becomes more expensive and time consuming as the scale grows. For the candidate whose fate is being decided, production data is either absent or too sparse to answer the question directly. The historical record of past candidates carried through all three scales is what a defensible forecast has to lean on, and fitting each scale independently throws it away. A production-only fit cannot disentangle intrinsic kinetics from scale effects: a slower intrinsic rate and a harsher mass-transfer limitation explain the same observations equally well. A bench-bioreactor-only fit identifies the intrinsics but says nothing about what scale-up does to them.

Figure 6. The scale-covariate idea. One clone's saturating response at three scales: each scale climbs to its plateau more slowly and settles lower than the last. The move from one scale to the next is a shift on the rate and ceiling, applied cumulatively. Each shift splits into a population part $\gamma_s$ shared by all clones and a feature part $\mathbf{x}_g^{\top}\boldsymbol{\beta}_s$ that depends on the clone, so a held-out clone's production curve is predicted from its features without any production data of its own.

A hierarchical model with scale as a covariate

Write one model spanning all scales, partitioning parameters by where they belong:

$$p\bigl(\boldsymbol{\theta}_{\text{int}}, \{\boldsymbol{\gamma}_s\} \mid \{\text{data}_s\}\bigr) \;\propto\; p(\boldsymbol{\theta}_{\text{int}}) \prod_s p(\boldsymbol{\gamma}_s)\, p \bigl(\text{data}_s \mid \boldsymbol{\theta}_{\text{int}},\, \boldsymbol{\gamma}_{\le s}\bigr)$$

The intrinsic parameters $\boldsymbol{\theta}_{\text{int}} = (k_g^{\text{int}}, y_{\max,g}^{\text{int}})$ capture scale-invariant clone kinetics, shared across every scale's likelihood. The intrinsic rate is regressed on a vector of features $\mathbf{x}_g$ describing each clone; the intrinsic ceiling is partially pooled through a shared Beta hyperprior:

$$\log k_g^{\text{int}} = \alpha + \mathbf{x}_g^{\top}\boldsymbol{\beta} + \varepsilon_g,\qquad \varepsilon_g \sim \mathcal{N}(0, \sigma_{\text{grp}})$$

Scale itself enters as a covariate: each scale contributes a shift on the rate and on the ceiling, added cumulatively on their natural link scales. Each shift has a population part $\gamma_s$ and a feature-dependent part $\mathbf{x}_g^{\top}\boldsymbol{\beta}_s$, so the size of a scale's effect can vary with the same features that drive the intrinsic kinetics. The saturating curve shifts scale-to-scale: production climbs to its plateau more slowly and settles at a lower one.

$$\begin{aligned} \log k_g^{(s)} &= \log k_g^{\text{int}} + \sum_{s' \le s} \bigl(\gamma^{k}_{s'} + \mathbf{x}_g^{\top}\boldsymbol{\beta}^{k}_{s'}\bigr) \\[2pt] \operatorname{logit} y_{\max,g}^{(s)} &= \operatorname{logit} y_{\max,g}^{\text{int}} + \sum_{s' \le s} \bigl(\gamma^{y}_{s'} + \mathbf{x}_g^{\top}\boldsymbol{\beta}^{y}_{s'}\bigr) \end{aligned}$$

Bench-bioreactor is the reference level: $\gamma^{k}_{\text{bench}} \equiv \gamma^{y}_{\text{bench}} \equiv 0$. Each shift absorbs the lumped effect of whatever new physics that scale introduces (oxygen transfer, shear, mixing, dissolved-gas gradients) without committing to a detailed mechanistic model of these effects. The population part $\gamma_s$ (with features centred, the mean shift) is shared across clones; the feature part $\mathbf{x}_g^{\top}\boldsymbol{\beta}_s$ lets scale sensitivity depend on the clone. A clone predisposed to struggle at scale (high oxygen demand, say) is forecast to struggle rather than assigned the population average. Crucially there is no per-clone scale parameter to over-fit from thin production data. A held-out clone's shift is predicted from its features (metabolic burden, growth rate, specific productivity), not estimated from its own (absent) production data.

Each shift carries weakly-informative $\mathcal{N}(0, 0.5)$ priors on its link scale (log-rate for the $k$ shifts, logit-ceiling for the $y$ shifts), centred on "no shift" and wide enough to admit a substantial scale effect in either direction. This matters most for the production shifts, which are fit only from the clones that carry production data. The fewer those clones, the less negligible the prior, so defensible production-shift priors are part of the model, not an afterthought.

The scale-specific likelihood mean is then

$$\mu_{gi}^{(s)} = y_{\max,g}^{(s)}\, \bigl(1 - e^{-k_g^{(s)}\, t_{gi}}\bigr) \cdot 100 \%,$$

so no two scales share the same saturating curve. In particular, the bench-bioreactor sees the intrinsic profile, pilot sees it shifted in rate and ceiling by one scale shift, production by two.

Information flows asymmetrically by construction. Bench-bioreactor observations enter the likelihood for the intrinsic pair only, so dense bench-bioreactor data tightens both $k_g^{\text{int}}$ and $y_{\max,g}^{\text{int}}$. Pilot observations, available for every clone, identify the pilot shifts (population $\gamma_{\text{pilot}}$ and feature slopes $\boldsymbol{\beta}_{\text{pilot}}$) from all eight clones. Production observations, available only for the subset given production data, identify the production shifts $\gamma_{\text{prod}}, \boldsymbol{\beta}_{\text{prod}}$. A held-out clone contributes no production data: its production curve is predicted by applying those shifts (the population part plus its own feature-predicted part) on top of its own bench-bioreactor- and pilot-identified curve.

Figure 7. Scale progression for a seen clone (clone A, $\bullet$ production data used) and the unseen held-out clone (clone F, $\times$ validation). Production data is given to seven, and the eighth (clone F) is held out and predicted. Shaded bands are 94% posterior predictive and include measurement noise $\sigma$. Time is shown dimensionless as $t/T_{\max}$, with $T_{\max}$ the shared 14-day fed-batch horizon; all three scales run the same length and differ only in rate and ceiling, so their saturation speeds are directly comparable. The endpoint is $t/T_{\max} = 1$ (day 14).

Cross-validating the forecast. Figure 7 shows one leave-one-out fold. Clone F is held out and forecast from the seven clones that carry production data. Repeating this for every clone removes any dependence on which one is held out. Refit the joint model eight times. Each fold gives production data to seven clones and predicts the eighth clone's production endpoint from its bench-bioreactor and pilot data plus the shifts learned from the other seven. This asks the question the shortlist problem actually poses: can an unseen clone be forecast from the rest? Figure 8's right panel plots the answer: every clone's held-out endpoint against its truth, one point per fold.

Figure 8. Validation of the three-scale fit. Left: forest plot of the population per-scale shifts $\gamma_s$ on log-rate and logit-ceiling, with ground truth ($\times$) and the no-pooling MLE ($\circ$, curves fit independently per clone and scale, then differenced and averaged) marked. All four posteriors recover truth within tight HDIs (rate shifts identified from the densely-sampled rising limb, ceiling from the plateau), while the MLE, with no pooling and no feature regression, scatters wider. The feature slopes $\boldsymbol{\beta}_s$ are fit alongside but not shown. Right: every clone's held-out endpoint ($\blacktriangle$) pooled across all eight leave-one-out folds (each clone predicted from the fold in which it was held out), with the aggregate RMSE, bias, and coverage inset. Both panels report mean-curve credible intervals, which exclude observation noise $\sigma$.

Across the eight folds the predicted production endpoint tracks truth to an RMSE of $0.96\%$, spanning clones from $21\%$ to $51\%$ of target titer, with a slight $-0.49\%$ bias (predictions run marginally low). The mean-curve 94% intervals cover the true endpoint in seven of the eight folds; the one miss (clone F) lands about $0.5\%$ outside.

Conclusion

When many similar groups exist, a hierarchy lets them share a population-level prior so data-rich groups inform data-poor ones. The same idea extends along scale: a joint model that shares intrinsic parameters and treats scale as a covariate lets cheap dense small-scale data tighten an expensive sparse large-scale forecast — something per-scale fits cannot do, because they spend their data identifying intrinsic and scale-specific parameters at once.

Biologics cell culture was the running example (Xing et al., 2009), but the structure recurs wherever cheap-scale data isolates intrinsics and expensive-scale data adds scale-specific physics on top. For example, bench vs pilot column vs production heap in metal leaching, where intrinsic dissolution kinetics fix early and solution transport, channelling, and bed heterogeneity stack on with scale (Petersen, 2016); coin-cell vs pouch-cell vs pack in battery scale-up (Kwade et al., 2018); coupon vs component vs full-assembly fatigue in materials qualification (Tridello et al., 2021).

Frequently asked questions

What are the "features" and where do they come from? The features $\mathbf{x}_g$ are measured properties of each group known before scale-up. In the cell-culture example, these are a clone's metabolic burden, growth rate, and specific productivity. They must be available for a new candidate that has no large-scale data of its own, because predicting that candidate's shift from its features is the entire point. Good features are cheap to measure early and carry real signal about how the process behaves at scale. The modelling presented here reveals whether a feature carries that signal.

Does the forward model have to be this saturating curve? No. Everything that does the work — posteriors instead of point estimates, partial pooling across groups, a feature regression on the parameters, and scale entered as a cumulative covariate — applies to any parametric forward model. Swap the saturating curve for whatever your process follows (a sigmoid, a decay, a mechanistic rate law) and the same hierarchical machinery applies.

What if the groups are not actually similar? The hierarchy estimates the group spread $\sigma_{\text{grp}}$ from the data, so it does not force agreement: if the groups genuinely disagree, $\sigma_{\text{grp}}$ grows and pooling relaxes toward independent per-group fits.

Can pooling ever make things worse? Pooling groups can be wrongly treated as similar when the features do not capture how their performance differs. Then strength is borrowed from the wrong siblings and the estimate is biased, not merely wide. The fix is in the features and the model structure, not the pooling: give the model the covariates that explain the real differences.

References

• Stein, C. (1956). Inadmissibility of the usual estimator for the mean of a multivariate normal distribution. Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, Vol. 1, pp. 197–206. University of California Press. doi:10.1525/9780520313880-018
• James, W., & Stein, C. (1961). Estimation with quadratic loss. Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability, Vol. 1, pp. 361–379. University of California Press. doi:10.1007/978-1-4612-0919-5_30
• Efron, B., & Morris, C. (1975). Data analysis using Stein's estimator and its generalizations. Journal of the American Statistical Association, 70(350), 311–319. doi:10.1080/01621459.1975.10479864
• Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). Bayesian Data Analysis (3rd ed.). Chapman & Hall/CRC. doi:10.1201/b16018
• Kwade, A., Haselrieder, W., Leithoff, R., Modlinger, A., Dietrich, F., & Droeder, K. (2018). Current status and challenges for automotive battery production technologies. Nature Energy, 3, 290–300. doi:10.1038/s41560-018-0130-3
• Petersen, J. (2016). Heap leaching as a key technology for recovery of values from low-grade ores — a brief overview. Hydrometallurgy, 165, 206–212. doi:10.1016/j.hydromet.2015.09.001
• Tridello, A., Boursier Niutta, C., Berto, F., & Paolino, D. S. (2021). Size-effect in very high cycle fatigue: a review. International Journal of Fatigue, 153, 106462. doi:10.1016/j.ijfatigue.2021.106462
• Xing, Z., Kenty, B. M., Li, Z. J., & Lee, S. S. (2009). Scale-up analysis for a CHO cell culture process in large-scale bioreactors. Biotechnology and Bioengineering, 103(4), 733–746. doi:10.1002/bit.22287