Post-Inference Model Chaining: A Problem Statement and a Model-based Bootstrap for Training End-to-End Digital Twins

Post-Inference Model Chaining: A Problem Statement and a Model-based Bootstrap for Training End-to-End Digital Twins Lukas Madenach*^, Timo Bailer*, Tobias Eilert*°^ * CMC Statistics BioPharma, Boehringer Ingelheim Pharma GmbH & Co. KG, Ingelheim, Germany ° Corresponding author ^These authors contributed equally.

Abstract End-to-end digital twins can be constructed by fitting separate stage models and connecting them after inference. We investigate the Frequentist prediction properties of a commonly used Monte Carlo simulation scheme that forwards complete predictive draws through such a chain when measured stage outputs contain independent normal errors, but the latent stage means are deterministic. Analytical two-stage examples show persistent asymptotic interval and coverage inflation and, with quadratic propagation, mean bias dependent on the parameters to be estimated. Simulations of four- and eight-stage linear–quadratic univariate chains with 7, 15 or 30 training runs demonstrate pronounced input-dependent distortions. A diagnostic using error-free latent training predictors shows that the effect is mainly not attributable to errors in variables, but to the MCSim approach even showing underestimation of coverage. We propose a model-based bootstrap: an initial chain forwards fitted means, bootstrap chains are refitted on fixed-center bootstrap data while forwarding their mean prediction, and a target-stage future error enters a studentized prediction root. Across the investigated scenarios, relative mean bias is substantially reduced, and coverage is closer to the nominal 90%. At sample size 30, the eightstage bootstrap has mean width ratios of 1.16–1.42 relative to the known-parameter normal reference, compared with 4.11–4.97 for observed-to-observed MCSim. These results support improved endpoint prediction in the tested setting. Further investigations of the bootstrap scheme in the supplement for a bivariate 6-stage process show similar satisfactory results.

1 Introduction As industry moves towards greater automation and digitalization, digital twins — virtual representations of physical processes — are becoming essential tools for monitoring, control, steering and optimization [1, 2, 3]. While digital twins are increasingly implemented at individual process stages, their true potential lies in predicting how changes at one stage affect later stages and ultimately the final product. However, constructing a single, comprehensive model that captures the entire process, i.e., an end-to-end digital twin (E2E-DT), is often infeasible due to the system’s complexity, data heterogeneity, data scarcity, computational constraints, and the need to model different stages with different model types. A pragmatic alternative is to model each stage separately using available input and output data and then reconnect these models sequentially also termed the modular approach [4, 5]. The statistical assessment distinguishes

bias of the mean predictor, pointwise coverage for one new measured observation and prediction-interval width relative to a known-parameter reference [6]. Biopharmaceutical production motivates this work, but the results in this manuscript apply to all fields in which separately trained models of whichever type need to be chained post-inference to yield a complete model with approximate proper statistical properties for example in metrology [7, 8]. The main task in the modular approach is how to reconnect the proper single models so that the entire process is described by one model that satisfies the classical frequentist statistical properties. We do not discuss the Bayesian or the GUM approach via the law of propagation of uncertainty [9, 10]. We call this task holistic process modelling. In the literature there is one popular approach known for measurement error propagation, which uses Monte Carlo simulations abbreviated here as MCSim [11, 12, 13, 14, 15]. The idea is simple: given a future input for the first stage, the first stage model is utilized to sample from its predictive distribution, i.e., we obtain a collection of samples representing the uncertainty of the output. These samples are forwarded as inputs to the second model and for each input sample one output sample is drawn from the predictive distribution of the second model. This is now repeated until the last stage model. The lastly produced samples (or also for some stage in between) are then used to compute a mean prediction and quantiles that serve as prediction intervals. We examine the bias, coverage and width of this particular predictive-draw propagation scheme. The conclusions concern how fitted predictive distributions are connected after inference, not the numerical validity of Monte Carlo integration. The following illustration distinguishes variability of a physical quantity from measurement error around a fixed underlying quantity. Both are represented by random variables, but their prediction targets differ. Composing the known physical variables gives the distribution of a transformed physical quantity; applying the same operation to a distribution around an estimate need not yield an unbiased estimate of the underlying transformed value. For example, let us assume we have a solution of spherical nanoparticles whose distribution of radii 𝑅 is known to be normal, i.e. 𝑅 ~ 𝒩(𝜇, 𝜎) with known mean radius 𝜇 and standard deviation 𝜎 (ignoring, for simplicity, the fact that radii cannot be negative). We emphasize again that there is no uncertainty here, only a stochastic process of random radii which is completely known since its type of probability distribution and its parameters are known. Let us assume that we are now interested in the distribution of the surface area of the nanoparticles giving rise to the random variable 𝑆 = 4𝜋𝑅 2 = 4𝜋(𝜇 + 𝜎𝑧)2 = 4𝜋(𝜇2 + 2𝜇𝜎𝑧 + 𝜎 2 𝑧 2 ), where we used that 𝑅 = 𝜇 + 𝜎𝑧 with a standard normal random variable 𝑧. The expected value is 𝐸(𝑆) = 4𝜋𝜇2 + 4𝜋𝜎 2 , because of its linearity as well as 𝐸(𝑧) = 0 and 𝐸(𝑧 2 ) = 1. We see that we have an offset of 4𝜋𝜎 2 to the result of 4𝜋𝜇2 . This is the impact of the randomness in our stochastic process of random radii on the average surface area, but not a “mistake”! Since the composition of random variables is effectively an integration, Monte Carlo simulations are the perfect tool here since they are a class of numerical methods to solve integrals. We now change the target from a varying radius to a fixed but unknown radius measured with error. The issue is whether directly transforming a distribution around its estimate preserves the desired target, rather than whether Monte Carlo can evaluate the transformation numerically. As an introductory example, let us now assume a different situation: we still have nanoparticles in a solution, but they all have exactly the same radius 𝜇 , which is unknown. However, the

analytical method to measure the radius has a normally distributed error with (for simplicity here) known standard deviation 𝜎. If we now measure once, i.e. the measured radius is given by the random variable 𝑅 = 𝜇 + 𝜎𝑧, the uncertainty distribution estimator for the mean radius 𝜇 is 𝜇̂ = 𝑅 − 𝜎𝑧 ′ , where 𝑧 and 𝑧′ are independently standard normally distributed [16, 17]. Following the MCSim approach analytically, we plug the predictive uncertainty distribution estimator 𝜇̂ into the surface area formula to get 2

𝑆̂ = 4𝜋𝜇̂ 2 = 4𝜋(𝑅 − 𝜎𝑧′)2 = 4𝜋(𝜇 + 𝜎𝑧 − 𝜎𝑧′)2 = 4𝜋(𝜇 + √2𝜎𝑧̅) = 4𝜋(𝜇2 + 2√2𝜇𝜎𝑧̅ + 2𝜎 2 𝑧̅ 2 ), with 𝑧̅ standard normally distributed, and the expectation value becomes 𝐸(𝑆̂) = 4𝜋𝜇2 + 8𝜋𝜎 2 .

Here, the offset of 8𝜋𝜎 2 is not correct anymore, because it is obvious that the true surface area of any of the nanoparticles is 𝑆 = 4𝜋𝜇2 . The ordinary squared plug-in estimator already has a nonlinearity-induced offset; adding and transforming an independent uncertainty draw doubles that offset. Thus, the additional distortion comes from the statistical construction, not from Monte Carlo integration itself. We will elaborate the effect of the MCSim approach on other statistical metrics in the context of holistic process modelling with more realistic examples in much more detail in Section 3. For post-inference model chaining, the objective is therefore to quantify uncertainty for a clearly specified endpoint target rather than merely to compose fitted predictive distributions. Stage estimators may be dependent because they share measured training variables. Fitting the stages separately does not make these random estimators independent, but propagating independent future predictive errors through fixed fitted models does not reproduce their joint sampling distribution. A sequential bootstrap can instead repeat the fitted-chain estimator on surrogate datasets, making its sources of variation explicit. Whether that surrogate approximation is adequate remains a model-specific question. However, even before the task of chaining, we are already faced with hurdles in the estimation of the single models. First, transformed predictors can become nearly collinear or yield numerically rank-deficient stage designs. Second, it is inherent to the problem that the input data for at least the second and the following stages are subject to random measurement errors. Near collinearity can make coefficient estimates unstable, but usually does not affect prediction performance drastically. The second point is known as the error-in-variables problem; it is more severe and needs to be dealt with [18]. Measurement error in the predictors can bias ordinary regression fits even when the latent stage model is correctly specified; subsequent nonlinear transformations can compound the resulting discrepancy. The latent-input diagnostic in this study isolates this issue from the separate consequences of the predictive propagation rule for the investigation of the MCSim approach. The problem statement first specifies the deterministic latent chain and measured stage outputs. The investigation then develops analytical intercept-to-linear and intercept-toquadratic examples, followed by simulations of four- and eight-stage chains alternating between linear and quadratic functions. The same simulation design includes observed-to-observed and latent-to-observed training for MCSim, so that predictor measurement error is not conflated with post-inference propagation.

The proposed alternative is a fixed-center studentized model-based bootstrap. Its initial estimator uses fitted preceding-stage means as predictors, and its pseudo-responses retain fixed initial stage centers while the bootstrap designs vary. A separate target-stage future error supplies the predictive component of the studentized root. We evaluate the resulting endpoint mean bias, prediction coverage and interval width on the same underlying four- and eight-stage processes as for the MCSim scheme. Additionally, we investigate its performance on a bivariate 3-stage chain in the supplement.

2 Problem statement and Notation Industrial production processes are inherently iterative, i.e., consisting of consecutive steps, referred to hereafter as unit operations (UOs) or stages. Since the number of unit operations, say 𝑛, usually range up to 8-10, the complexity of each UO is high and the data is usually scarce, it is often not feasible or advisable to estimate one model describing the whole production line at once. A common workaround is that each unit operation is inferred separately and then post-hoc re-connected to yield a model describing the complete production process. The first step to understand and then overcome the problems that arise in this approach is to define a statistical model for the production process further denoted by 𝐹. We start with the description of a mean production line originating from mean functional models for each unit operation 𝑖 given by the functions 𝑓𝑖̅ : ℝ𝑛𝑖−1 ⟶ ℝ𝑛𝑖 each having 𝑛𝑖−1 real inputs and 𝑛𝑖 real outputs. The complete mean production line is then given by the composition of the single unit operations 𝐹̅ ≔ 𝐹̅1,..,𝑛 ≔ Ο𝑛𝑖=1 𝑓𝑖̅ : ℝ𝑛0 ⟶ ℝ𝑛𝑛 , where we used the convention ̅ ̅ (⋯ 𝑓2̅ (𝑓1̅ (𝑥0 )))), ∘ ⋯ ∘ 𝑓2̅ ∘ 𝑓1̅ )(𝑥0 ) = 𝑓𝑛̅ (𝑓𝑛−1 (Ο𝑛𝑖=1 𝑓𝑖̅ )(𝑥0 ) = (𝑓𝑛̅ ∘ 𝑓𝑛−1 for a starting configuration 𝑥0 ∈ ℝ𝑛0 . After each UO, the output material is measured in many different attributes so that we add a homoscedastic centered error vector 𝜖𝑖 to the mean functions to account for this. Thus, the statistical model for the unit operation 𝑖 can be expressed by 𝑓𝑖 (𝑥𝑖−1 ) = 𝑓𝑖̅ (𝑥𝑖−1 ) + 𝜖𝑖 for all 𝑥𝑖−1 ∈ ℝ𝑛𝑖−1 . The measurement errors are not propagated to the next unit operation; they are added to the measured output only. This results in a statistical model for the complete production line constructed as 𝐹(𝑥0 ) = 𝐹̅1,..,𝑛 (𝑥0 ) + 𝜀𝑛 = (Ο𝑛𝑖=1 𝑓𝑖̅ )(𝑥0 ) + 𝜀𝑛 . Note that this statistical model represents the truth, which is not known to the engineers in reality. We assume that the true process is constructed from deterministic functions, i.e., we restrict ourselves to measurement errors and do not consider process variability at the unit operations. Measurement errors are always present, and separating them from process variability keeps the propagation of pure measurement uncertainty clear. However, correct propagation of statistical uncertainty including process variability is part of forthcoming work. Given that the 𝑘 -th execution of the process, also called a run or batch, starts at a starting configuration 𝑥0𝑘 ∈ ℝ𝑛0 for 𝑘 ∈ {1, … , 𝑚}, we in total collect 𝑚 data sets, where 𝑚 is the number of executions of the production line. These consist of independent random vectors defined by 𝑘 𝑖 𝑋𝑖𝑘 ≔ 𝐹1,..,𝑖 𝑓𝑗̅ )(𝑥0𝑘 ) + 𝜀𝑖𝑘 , i.e., the random outputs of each unit operation given that (𝑥0𝑘 ) ≔ (Ο𝑗=1 𝑇

𝑇

𝑇

𝑇 𝑇

we started at 𝑥0𝑘 , building up the random data vector 𝑋 𝑘 ≔ (𝑥0𝑘 , 𝑋1𝑘 , 𝑋2𝑘 , ⋯ , 𝑋𝑛𝑘 ) for execution 𝑘 ∈ {1, … , 𝑚} and yielding the random data matrix 𝑋 = (𝑋1 , … , 𝑋 𝑚 )𝑇 . The error vectors

𝜀𝑖𝑘 are assumed to be iid between executions, i.e., independent copies of 𝜀𝑖 , as well as independent between unit operations. From each realization of the random data submatrix 𝑋𝑖 ≔ 𝑇 𝑇

𝑇

1 2 𝑚 ((𝑋𝑖−1 , 𝑋𝑖−1 , … , 𝑋𝑖−1 ) , (𝑋𝑖1 , 𝑋𝑖2 , … , 𝑋𝑖𝑚 ) ) , we have information to build a predictive estimator 𝑓̃𝑖

for the statistical model 𝑓𝑖 of the 𝑖-th unit operation [19]. An important observation is that the 𝑇

estimators of consecutive unit operations 𝑖 − 1 and 𝑖 share the random data (𝑋𝑖1 , 𝑋𝑖2 , … , 𝑋𝑖𝑚 ) . This means that the estimator of a unit operation depends at least on the directly preceding unit operation, but usually also on earlier ones. Now the ultimate goal in holistic process modelling is to estimate 𝐹, i.e., the entire production line, by one estimator, say 𝐹̃ . Since this is not feasible, the post-inference approach needs to somehow use the collection of trained models {𝑓̃𝑖 }𝑖∈{1,…,𝑛} to obtain a suitable approximation of 𝐹̃ . The MCSim approach takes the intuitive way of chaining these distribution estimators to obtain an approximation by 𝐹̃𝑀𝐶𝑆𝑖𝑚 (𝑥0 ) ≔ (Ο𝑛𝑖=1 𝑓̃𝑖 )(𝑥0 ). The next section investigates this approach in a frequentist setting.

3 Investigation of MCSim Approach 3.1 Analytical Univariate Case As an example – that can be solved purely analytically – we look at a production line with only two simplified unit operations which have a univariate input as well as univariate output, i.e., 𝑛0 = 𝑛1 = 𝑛2 = 1. The first mean function is set to be an intercept model, that is 𝑓1̅ (𝑥0 ) = 𝑎1 , which is independent of the starting value 𝑥0 . The random output data points are then given by 𝑋1𝑘 = 𝑓1𝑘 (𝑥0 ) = 𝑎1 + 𝜀1𝑘 with independent and normally distributed measurement errors 𝜀1𝑘 ~ 𝑁(0, 𝜎1 ) and positive and finite standard deviation 𝜎1 . The random data is denoted by 𝑋1 = {𝑋11 , … , 𝑋1𝑚 }. It follows that a predictive distribution estimator for 𝑓1 is given in the form 𝑓̃1 (𝑥0 ) = 𝑋̅ + 𝑡𝑚−1 𝑆√1 +

1 𝑚

with 𝑡𝑚−1 ~ 𝑇𝑚−1 , the 𝑡-distribution with 𝑚 − 1 degrees of freedom, 𝑋̅ =

𝑚 𝑚 1 1 ∑ 𝑋1𝑘 = ∑ 𝑓1𝑘 (𝑥0 ) ~ 𝑁(𝑎1 , 𝜎1 /√𝑚) 𝑚 𝑚 𝑘=1 𝑘=1

and 𝑆2 =

𝑚 𝑚 1 1 2 2 ∑ (𝑋1𝑘 − 𝑋̅) = ∑ (𝑓1𝑘 (𝑥0 ) − 𝑋̅) = 𝜎12 𝜃/(𝑚 − 1). 𝑚 − 1 𝑘=1 𝑚 − 1 𝑘=1

The latter two are the estimators for the sample mean and the unbiased sample variance of the random data. Here, the random variable 𝜃 follows a 𝜒 2 -distribution with 𝑚 − 1 degrees of freedom. It is readily seen that the distribution estimator 𝑓̃1 (𝑥0 ) is consistent, unbiased and shows correct coverage behavior for 𝑓1 (𝑥0 ).

As a first holistic example, we choose an affine function for the second unit operation that is 𝑓2 (𝑥1 ) = 𝑎2 + 𝑏2 𝑥1 + 𝜀2 with 𝜀2 ~ 𝑁(0, 𝜎2 ) , 𝜎2 > 0 and 𝑏2 ≠ 0 , where we assume that the coefficients are known, i.e., 𝑓̃2 = 𝑓2 . This is not realistic, but it is useful to derive an analytic solution, and it is a simplification of the real situation. Given the true intercept model of the first unit operation and the true affine function of the second UO, the true holistic mean function is 𝐹̅ (𝑥0 ) = 𝑎2 + 𝑏2 𝑎1. Applying the MCSim approach analytically we obtain 𝐹̃𝑀𝐶𝑆𝑖𝑚 (𝑥0 ) = 𝑓̃2 (𝑓̃1 (𝑥0 )) = 𝑓2 (𝑓̃1 (𝑥0 )) + 𝜀2 = 𝑎2 + 𝑏2 (𝑋̅ + 𝑡𝑚−1 𝑆√1 +

1 ) + 𝜀2 𝑚

= 𝑎2 + 𝑏2 𝑋̅ + 𝑏2 𝑡𝑚−1 𝑆√1 +

1 + 𝜀2 . 𝑚

The target of the prediction is a new observation 𝐹(𝑥0 ) = 𝐹̅ (𝑥0 ) + 𝜀2 with target mean 𝐹̅ (𝑥0 ) = 𝐸(𝐹(𝑥0 )) = 𝑎2 + 𝑏2 𝑎1. Since 𝐸(𝑡𝑚−1 ) = 0 (the draw of 𝑡𝑚−1 is independent of the data), 𝐸(𝑋̅) = 𝑎1 and 𝐸(𝜀2 ) = 0, the expected value of the MCSim estimator is 𝐸 (𝐹̃𝑀𝐶𝑆𝑖𝑚 (𝑥0 )) = 𝑎2 + 𝑏2 𝑎1 = 𝐹̅ (𝑥0 ) i.e., the MCSim estimator is unbiased for every 𝑚, and its expectation trivially converges to the target mean for 𝑚 → ∞ (consistency). For the coverage, the central 1 − 𝛼 interval of the MCSim distribution has to contain a new observation with probability 1 − 𝛼, i.e., 𝑃 (𝑄𝐹̃𝑀𝐶𝑆𝑖𝑚 (𝑥0 )|𝑋1 (𝛼/2) ≤ 𝐹(𝑥0 ) ≤ 𝑄𝐹̃𝑀𝐶𝑆𝑖𝑚(𝑥0 )|𝑋1 (1 − 𝛼/2)) = 1 − 𝛼 where 𝑄 denotes the quantile function of the MCSim distribution given the data 𝑋1 , which is, thus, a random quantity that depends on the data only. Unfortunately, no closed-form solution exists for finite 𝑚, because the MCSim distribution is the sum of a t-distributed and a normally distributed term. We therefore investigate the limit 𝑚 → ∞. Since 𝑋̅ → 𝑎1 and 𝑆 → 𝜎1 almost surely, and 𝑡𝑚−1 converges in distribution to a standard normally distributed variable 𝑧𝑡 , we obtain 𝐹̃𝑀𝐶𝑆𝑖𝑚 (𝑥0 ) = 𝑎2 + 𝑏2 𝑋̅ + 𝑏2 𝑡𝑚−1 𝑆√1 +

1 + 𝜀2 → 𝑎2 + 𝑏2 𝑎1 + 𝑏2 𝜎1 𝑧𝑡 + 𝜎2 𝑧2 𝑚

= 𝑎2 + 𝑏2 𝑎1 + √𝑏22 𝜎12 + 𝜎22 𝑧 with a standard normally distributed variable 𝑧 (the variables 𝑧𝑡 and 𝑧2 are independent). The quantile function of this limit distribution, which no longer depends on the data, is 2 2 2 𝑄𝐹∞ ̃𝑀𝐶𝑆𝑖𝑚 (𝑥0 ) (𝛼) = 𝑎2 + 𝑏2 𝑎1 + √𝑏2 𝜎1 + 𝜎2 𝑄𝑧 (𝛼)

The asymptotic coverage is therefore

∞ 𝑃 (𝑄𝐹∞ ̃𝑀𝐶𝑆𝑖𝑚 (𝑥0 ) (𝛼/2) ≤ 𝐹(𝑥0 ) ≤ 𝑄𝐹̃𝑀𝐶𝑆𝑖𝑚 (𝑥0 ) (1 − 𝛼/2))

𝜎22 = 𝑃 (𝑄𝑧 (𝛼/2) ≤ √ 2 2 𝑧 ≤ 𝑄𝑧 (1 − 𝛼/2)) 𝑏2 𝜎1 + 𝜎22 2 Hence, the central quantity is the standard normally distributed variable 𝑧2 scaled by the factor 𝜎2 √ 2 2 2 2. 𝜎2 + 𝑏2 𝜎1 The asymptotic coverage equals 1 − 𝛼 if 𝑏22 𝜎12 = 0, and it is approximately 1 − 𝛼 if 𝜎22 ≫ 𝑏22 𝜎12, because this factor is then (approximately) one. For 𝑏22 𝜎12 > 0 the factor is smaller than one and the coverage is larger than the nominal level, i.e., it is overestimated; for 𝜎22 = 0 the coverage would be 100%. In the limit, fortunately, the coverage of the MCSim interval is therefore never below the target. However, how strongly it is overestimated depends on the unknown quantities 𝜎1 , 𝜎2 and 𝑏2 , i.e. the extent is not known to the engineers in a modelling scenario. Since this is an asymptotic result, we see that even with an infinite amount of data the MCSim approach does not perform as it is expected from a well calibrated statistical model. The same holds for the width of the interval. The central 1 − 𝛼 interval of the distribution of the target random variable 𝐹(𝑥0 ) is [𝑄𝐹(𝑥0 ) (𝛼/2), 𝑄𝐹(𝑥0 ) (1 − 𝛼/2)] and its width is 𝑄𝐹(𝑥0 ) (1 − 𝛼/2) − 𝑄𝐹(𝑥0 ) (𝛼/2) = 𝜎2 (𝑄𝑧 (1 − 𝛼/2) − 𝑄𝑧 (𝛼/2)) In the limit, the average width of the MCSim interval is ∞ 𝐸 (𝑄𝐹∞ ̃𝑀𝐶𝑆𝑖𝑚 (𝑥0 ) (1 − 𝛼/2) − 𝑄𝐹̃𝑀𝐶𝑆𝑖𝑚 (𝑥0 ) (𝛼/2))

= √𝑏22 𝜎12 + 𝜎22 (𝑄𝑧 (1 − 𝛼/2) − 𝑄𝑧 (𝛼/2)) Hence, the ratio of the average width of the MCSim interval to the true average width is 𝑏2𝜎2 + 𝜎 2 √ 2 12 2 𝜎2 This ratio is one if 𝑏22 𝜎12 = 0 and approximately one if 𝜎22 ≫ 𝑏22 𝜎12 , but it is larger than one for 𝑏22 𝜎12 > 0: the average width of the MCSim interval is then overestimated by this factor. As for the coverage, the amount of the overestimation depends on the unknown quantities 𝜎1 , 𝜎2 and 𝑏2 , and again this result implies that collecting more data does not bring convergence to proper behavior. As an even stronger counter example, the second unit operation is chosen to be a quadratic function with intercept, i.e. 𝑓2 (𝑥1 ) = 𝑎2 + 𝑐2 𝑥12 + 𝜀2 with 𝜀2 ~ 𝑁(0, 𝜎2 ) with 𝜎2 > 0 and 𝑐2 ≠ 0, which coefficients are known as well, i.e. 𝑓̃2 = 𝑓2 . Given the true intercept model of the first UO

and the true quadratic polynomial of the second UO, the true holistic mean function is 𝐹̅ (𝑥0 ) = 𝑓2̅ (𝑓1̅ (𝑥0 )) = 𝑎2 + 𝑐2 𝑎12 . Again, analytically applying the MCSim approach, we obtain 𝐹̃𝑀𝐶𝑆𝑖𝑚 (𝑥0 ) = 𝑓̃2 (𝑓̃1 (𝑥0 )) = 𝑎2 + 𝑐2 𝑋̅ 2 + 2𝑐2 𝑋̅𝑡𝑚−1 𝑆√1 +

1 1 2 + 𝑐2 𝑡𝑚−1 𝑆 2 (1 + ) + 𝜀2 . 𝑚 𝑚

2 An interesting observation is that 𝑡𝑚−1 is 𝐹-distributed with 1 and 𝑚 − 1 degrees of freedom, i.e., 𝐹̃𝑀𝐶𝑆𝑖𝑚 (𝑥0 ) has no finite mean for 𝑚 ≤ 3 and no finite variance for 𝑚 ≤ 5. That means that it is

possible that even though the intercept estimator has finite moments, the MCSim approach might not. Taking expectations gives 𝐸 (𝐹̃𝑀𝐶𝑆𝑖𝑚 (𝑥0 )) = 𝑎2 + 𝑐2 𝑎12 + 𝑐2

𝜎12 𝑚−1 1 + 𝑐2 𝜎12 ( ) (1 + ) 𝑚 𝑚−3 𝑚

= 𝐹̅ (𝑥0 ) + 𝑐2 𝜎12 𝑘𝑚 , 𝑚2 +𝑚−4

with 𝜅𝑚 ≔ 𝑚(𝑚−3) > 1 for all 𝑚 ≥ 4. Above, we used that 2 )= 𝐸(𝑡𝑚−1

𝑚−1 , 𝑚−3

and 𝐸(𝑋̅ 2 ) = 𝑎12 +

𝜎12 , 𝑚

where 𝑧 is standard normally distributed. Therefore, we have a bias of the MCSim mean prediction relative to the true holistic mean 𝐸 (𝐹̃𝑀𝐶𝑆𝑖𝑚 (𝑥0 )) − 𝐹̅ (𝑥0 ) = 𝑐2 𝜎12 𝑘𝑚 which depends on the variance of measurement error and quadratic coefficient. Its sign follows the sign of 𝑐2 , i.e. the direction of the bias is unknown to the engineers in reality. The finite sample size impact via 𝜅𝑚 is great when 𝑚 is small, for example 𝜅4 = 4 or even still 𝜅15 ≈ 1.3, which is oftentimes the case in holistic process modelling. But also for 𝑚 → ∞, the bias 𝑐2 𝜎12 𝜅𝑚 ⟶ 𝑐2 𝜎12 does not vanish, i.e. the problem cannot be fixed with more data. The coverage and width behavior cannot be investigated analytically nor asymptotically anymore, because the quantiles of 𝐹𝑀𝐶𝑆𝑖𝑚 (𝑥0 ) are not accessible. These properties are evaluated in the simulation study in Section 3.2. The linear example shows that MCSim can preserve the mean while still overestimating asymptotic coverage and interval width. The quadratic example shows an additional failure: nonlinear propagation introduces a mean bias that persists as the sample size increases. Table 1 summarizes the results established analytically.

TABLE 1: S UMMARY OF EFFECTS DUE TO CHAINING WITH MCSIM APPROACH Known secondstage model

Chaining method

Unbiased point predictor

Consistent point predictor

Coverage as Width relative sample size to the true increases interval

Affine

MCSim

Yes

Quadratic

MCSim

No

Yes (almost At or above surely) the nominal level in the limit No Not derived here; evaluated by simulation

Ratio at least 1 in the limit Not derived here; evaluated by simulation

We note that both examples are highly idealized, because the second unit operation is assumed to be known. If we would need to infer these coefficients as well or for even more consecutive UOs than two, the situation is expected to synergistically get worse. These scenarios are addressed by Monte Carlo simulation of the entire experiment, where this is the appropriate tool.

3.2 Simulation Studies To investigate post-inference chaining beyond the analytical two-stage examples, we simulated a deterministic eight-stage chain and evaluated its outputs after stages 4 and 8. These endpoint models are denoted LQLQ and LQLQLQLQ, respectively, where L and Q indicate linear and quadratic stage means. The four-stage model is the prefix of the same eight-stage chain, not a separately specified data-generating process. The true stage means are: 𝑓1̅ (𝑧) = 𝑧 𝑓2̅ (𝑧) = 0.5𝑧 + 0.1𝑧 2 𝑓3̅ (𝑧) = 5 − 𝑧 𝑓4̅ (𝑧) =

20 20 𝑧 − 𝑧2 7 49

𝑓5̅ (𝑧) = 𝑧 𝑓6̅ (𝑧) = 5 − 0.6𝑧 − 0.08𝑧 2 𝑓7̅ (𝑧) = 5 − 𝑧 𝑓8̅ (𝑧) =

200 60 20 + 𝑧 − 𝑧2 49 49 49

Each stage maps [0, 5] onto [0, 5]. The linear stages alternate between increasing and decreasing functions. The quadratic stages alternate between monotone functions (stages 2 and 6) and functions with an interior maximum (stage 4 at z = 3.5; stage 8 at z = 1.5). The latent output at stage 𝑖 is 𝐹̅𝑖 (𝑥) = 𝑓𝑖̅ (𝐹̅𝑖−1 (𝑥)), with 𝐹̅0 (𝑥) = 𝑥. Only latent stage means enter the next true stage; measurement errors are added to the measured outputs, not propagated through the true chain.

For each sample size 𝑚 ∈ {7,15,30}, we performed 𝑅 = 1000 independent outer repetitions. The 𝑚 initial inputs were equally spaced over [0, 5], including both endpoints. Each run 𝑘 provided one measured output at each of the eight stages: 𝑋𝑖𝑘 = 𝐹̅𝑖 (𝑥0𝑘 ) + 𝜀𝑖𝑘 , 𝜀𝑖𝑘 ∼ 𝑁(0, 𝜎), 𝜎 = 0.5, where the measurement errors were independent across runs and stages, with the same standard deviation at every stage. Thus, the measurement-error standard deviation was 10% of the stage signal range of 5. Predictions were evaluated on the grid 𝐺 = {0,0.25, … ,5}. The true endpoint mean on this grid is 𝜇𝑖 (𝑔) = 𝐹̅𝑖 (𝑔), for 𝑖 = 4 or 8 and 𝑔 ∈ 𝐺. Each stage was fitted separately by ordinary least squares, using an intercept and a linear term for L stages, and an intercept, linear term and squared input for Q stages. Stage 1 used the known initial inputs. In the observed-to-observed variant shown in the main figures, each later stage regressed its measured output on the measured output of the preceding stage. The latent-toobserved variant instead used the true preceding latent outputs as training inputs. This latter variant is an oracle diagnostic, not a feasible procedure with unknown latent values: it removes measurement error from the training predictors while retaining the same MCSim propagation rule. It therefore distinguishes an errors-in-variables effect from the effect of forwarding predictive draws after inference. Both variants used the same generated training data within each outer repetition. The endpoint prediction interval was the empirical 5th to 95th percentile interval of the 𝐵 endpoint draws (nominal level 90%, 𝛼 = 0.1 ). No predictive draw was truncated to [0, 5] ; consequently, propagated inputs could leave the training range, where the fitted polynomials and leverage terms were still evaluated. As we already saw in Section 3.1 the mean in the MCSim approach might not exist. Thus, the code tracks a predictive tail index 𝜏𝑖 . The predictive mean exists only when the tail index exceeds 1; a finite variance requires an index exceeding 2. Thus, a finite 𝐵 always yields a numerical sample mean, but that value need not estimate an existing population mean. Even when the mean exists, small tail indices can make it unstable. For each endpoint and outer repetition 𝑟, the code selected the mean of the 𝐵 draws if the tail index exceeded 1; otherwise it selected their median. Denote this selected statistic by 𝜇̂ 𝑖,𝑟 (𝑔). Its relative deviation from the true endpoint mean was 𝛿𝑖,𝑟 (𝑔) =

𝜇̂ 𝑖,𝑟 (𝑔) − 𝜇𝑖 (𝑔) , 𝛥𝑖 = 𝑚𝑎𝑥𝜇𝑖 (𝑔) − 𝑚𝑖𝑛𝜇𝑖 (𝑔), 𝑔 ∈ 𝐺. 𝛥𝑖

The denominator is the range of the true endpoint means over the evaluation grid, not automatically 5. Coverage concerns one new measured endpoint observation at a fixed grid point. The reported coverage is the mean over the 1000 outer repetitions. Interval width was likewise averaged over repetitions after division by the true central 90% normal-interval width. A width ratio of 1 corresponds to the interval width with known endpoint mean and measurement-error standard deviation. It is an oracle reference, not a finite-sample optimum for intervals based on fitted models. For the summary curves, pointwise 90% percentile bootstrap confidence intervals were calculated from 2.000 resamples of the outer repetitions.

Figure 1 shows that MCSim point predictions deviate substantially from the true endpoint means. For LQLQ, relative mean bias ranges from −22.1% to −16.9% at 𝑚 = 7, from −13.0% to +18.8% at 𝑚 = 15 and from −11.5% to +29.3% at 𝑚 = 30 over the grid. A larger training set therefore does not uniformly remove the discrepancy. For LQLQLQLQ, the median-based relative deviations range from −44.8% to +41.4% at 𝑚 = 7 and from −38.2% to +39.5% at 𝑚 = 15. At 𝑚 = 30, where the predictive mean exists, relative mean bias ranges from −45.0% to +34.1%. We must note that wherever the mean prediction does not exist anymore, the MCSim approach is not applicable anymore. Consequently, their ordering must not be interpreted as a monotone sample-size effect on mean bias. The pronounced input-dependent deviations demonstrate that the distortion is not merely a common constant offset.

FIGURE 1: RELATIVE ENDPOINT DEVIATIONS FROM THE TRUE MEAN UNDER OBSERVED -TO-OBSERVED MCS IM TRAINING . PANELS A AND B SHOW LQLQ ( FOUR STAGES ) AND LQLQLQLQ ( EIGHT STAGES ). B LUE, ORANGE AND GREEN CORRESPOND TO M = 7, 15 AND 30. VALUES ARE NORMALIZED BY THE TRUE ENDPOINT -MEAN RANGE OVER THE 21- POINT INPUT GRID . M EAN- BASED CELLS REPORT THE MEAN ACROSS 5,000 REPETITIONS OF THE PREDICTIVE -MEAN DEVIATION . F OR STAGE 8 AT M = 7 AND 15, THE PREDICTIVE MEAN DOES NOT EXIST ; THE STATISTIC IS INSTEAD THE MEDIAN ACROSS REPETITIONS OF THE PREDICTIVE - MEDIAN DEVIATION . RIBBONS ARE POINTWISE 90% OUTER REPETITION BOOTSTRAP CONFIDENCE INTERVALS (2,000 RESAMPLES ); B = 5,000 PREDICTIVE CHAINS PER GRID POINT . Reliable prediction also requires appropriate coverage for a new measured endpoint observation. Figure 2 shows that the observed-to-observed MCSim intervals exceed the nominal 90% target at every evaluated grid point by far. Coverage ranges from 95.4% to 100% for LQLQ and from 94.1% to 100% for LQLQLQLQ across the three sample sizes. This is conservative behavior in these particular simulations, not evidence of well-calibrated 90% intervals. Increasing 𝑚 does not bring the curves uniformly towards the nominal target.

FIGURE 2: P OINTWISE COVERAGE OF NOMINAL 90% MCS IM PREDICTION INTERVALS FOR A NEW MEASURED ENDPOINT OBSERVATION , USING OBSERVED -TO -OBSERVED TRAINING . P ANELS A AND B SHOW LQLQ AND LQLQLQLQ; COLORS IDENTIFY M = 7, 15 AND 30 AS IN F IGURE 1. C OVERAGE IS THE TRUE NORMAL PROBABILITY INTEGRATED OVER EACH FITTED INTERVAL AND THEN AVERAGED OVER 5,000 OUTER REPETITIONS . THE BLACK HORIZONTAL LINE IS THE 90% TARGET . RIBBONS ARE POINTWISE 90% OUTER -REPETITION BOOTSTRAP CONFIDENCE INTERVALS . B = 5,000 PREDICTIVE CHAINS PER GRID POINT .

Figure 3 quantifies the loss of precision behind this overcoverage. For LQLQ, average interval widths are 3.66–11.16 times the true normal-interval width at 𝑚 = 7, 2.19–6.30 times at 𝑚 = 15 and 1.87–4.97 times at 𝑚 = 30 . For LQLQLQLQ, the corresponding ranges are 10.26– 656.78, 4.83–12.93 and 4.11–4.97. Thus, even m = 30 leaves endpoint intervals approximately four to five times the oracle width for the eight-stage chain. The especially large values for 𝑚 = 7 occur in a regime of severe tail amplification. Widths substantially above the oracle reference are not by themselves proof of an invalid finite-sample interval, but the combination of very large widths, systematic point-prediction deviations and coverage far above the requested level makes this implementation poorly calibrated and imprecise in the investigated designs.

FIGURE 3: MEAN MCS IM INTERVAL -WIDTH RATIO UNDER OBSERVED -TO-OBSERVED TRAINING . PANELS A AND B SHOW LQLQ AND LQLQLQLQ; COLORS IDENTIFY SAMPLE SIZES AS IN FIGURE 1. EACH 90% PREDICTION -INTERVAL WIDTH IS DIVIDED BY THE TRUE CENTRAL 90% NORMAL -INTERVAL WIDTH AND AVERAGED OVER 1000 OUTER REPETITIONS . T HE PERCENTAGE SCALE EXPRESSES A RATIO OF 1 AS 100%, WHICH IS THE DASHED ORACLE REFERENCE . RIBBONS ARE POINTWISE 90% OUTER REPETITION BOOTSTRAP CONFIDENCE INTERVALS ; B = 5,000 PREDICTIVE CHAINS PER GRID POINT . The latent-to-observed diagnostic (Figure S1, Figure S2 and Figure S3) separates the inputmeasurement-error issue from post-inference propagation. In that variant, every stage is fitted using the true preceding latent output, so the training predictors are measurement-error free. Nevertheless, for LQLQ at m = 30, relative mean bias still ranges from −14.5% to −5.9%, and mean interval widths range from 1.42 to 7.03 times the oracle width. For LQLQLQLQ at m = 30, widths range from 5.56 to 27.83 times the reference and coverage falls to about 83.1% at some inputs despite being near 100% elsewhere. These results show that correcting the training-input measurement error alone does not repair the forwarding of full predictive draws. They also caution against extending the analytical linear two-stage overcoverage result to arbitrary nonlinear multi-stage chains. The simulation supports the need for a different post-inference construction, evaluated next using measurement-error setting and models.

4 Holistic Retraining The main hurdles in holistic process modelling are the error-in-variable problem and the missed interdependencies due to separate training as discussed in more detail in the introduction. In this section we propose an algorithm that aims to holistically retrain the single models to reconnect them to a holistic model that has drastically reduced bias, fulfills the target coverage approximately and has a smaller interval width in comparison with the MCSim approach. To test this algorithm, we apply it to the two chains introduced in the previous section as well as additionally to a bivariate 6-stage model shown in the supplement.

4.1 Model-based Bootstrap The algorithm consists of two phases: the first constructs the initial holistic or multi-stage model that is used for the holistic mean prediction and the second phase generates a set of bootstrap holistic models which are used to compute the uncertainty of the holistic prediction. The core of the algorithm is a model-based bootstrap, i.e., we train models on the data and reuse them as data-generating processes. This already highlights the core requirement for the algorithm: the model type needs to know its own estimated error distribution. This is usually given for (mixed effect) linear parameter models or whenever the single models can be represented by a bootstrap model (parametric, residual or other bootstraps) the algorithm is applicable. The two phases of the algorithm to train the holistic model can be found in Pseudo Code 1 and 2 and the algorithm to compute its predictions with uncertainty in Pseudo Code 3. Its variables are explained in Table 2. TABLE 2: VARIABLES USED IN THE PSEUDO CODES Variable 𝑛 𝑖 𝐷𝑖 𝑀𝑖 𝑀𝑖 (𝑥, "𝑚") 𝑀𝑖 (𝑥, "𝑠") ̅𝑖 𝐷 𝐵 𝑏 ∈ {1, … , 𝐵} 𝐷𝑖𝑏 𝑀𝑖𝑏 𝑏 𝑀𝑖 (𝑥, "𝑚") 𝑀𝑖𝑏 (𝑥, "𝑠") ̅𝑖𝑏 𝐷

Explanation number of stage-models in the multi-stage model index to (initial or bootstrap) stage-model and mean predictions (𝑖 ∈ {1, … , 𝑛}) as well as measured data (𝑖 ∈ {0, 1, … , 𝑛}) output data from stage 𝑖 and input data to stage 𝑖 + 1 𝑖-th initial stage-model generating mean predictions from 𝑖-th initial stage-model for input data 𝑥 drawing samples from the 𝑖-th initial stage-model error distribution for input data 𝑥 mean predictions from 𝑖-th initial stage-model number of bootstrap loops and bootstrap (multi-)stage models index to bootstrap iteration and bootstrap (multi-)stage model bootstrap output data from stage 𝑖 of the 𝑏-th bootstrap loop 𝑖-th bootstrap stage-model of the 𝑏-th bootstrap loop generating mean predictions from 𝑖 -th bootstrap stage-model of the 𝑏-th bootstrap loop for input data 𝑥 drawing samples from the 𝑖-th bootstrap stage-model uncertainty distribution of the 𝑏-th bootstrap loop for input data 𝑥 mean predictions from 𝑖 -th bootstrap stage-model of the 𝑏 -th bootstrap loop

PSEUDO C ODE 1: For 𝑖 in 1 till 𝑁: If 𝑖 = 1 Else

Train 𝑀1 on (𝐷0 , 𝐷1) ̅1 to 𝑀1 (𝐷0 , "𝑚") Set 𝐷 ̅𝑖−1 , 𝐷𝑖 ) Train 𝑀𝑖 on (𝐷 ̅𝑖 to 𝑀𝑖 (𝐷 ̅𝑖−1 , "𝑚") Set 𝐷

PSEUDO C ODE 2: For 𝑏 in 1 till 𝐵: For 𝑖 in 1 till 𝑁: If 𝑖 = 1

Else:

Set 𝐷1𝑏 to 𝑀1 (𝐷0 , "𝑠") Train 𝑀1𝑏 on (𝐷0 , 𝐷1𝑏 ) ̅1𝑏 to 𝑀1𝑏 (𝐷0 , "𝑚") Set 𝐷 ̅𝑖−1 , "𝑠") Set 𝐷𝑖𝑏 to 𝑀𝑖 (𝐷 𝑏 𝑏 ̅𝑖−1 Train 𝑀𝑖 on (𝐷 , 𝐷𝑖𝑏 ) 𝑏 𝑏 𝑏 ̅𝑖 to 𝑀𝑖 (𝐷 ̅𝑖−1 , "𝑚") Set 𝐷

In prose, the algorithm reads: in the first phase, we build an initial multi-stage model by training the first stage model 𝑀1 on the input data 𝐷0 (that is assumed to be deterministic) and output ̅1, i.e. the data 𝐷1. In the next step, the second model is not trained on 𝐷1 as input data, but on 𝐷 mean predictions from the first stage model 𝑀1 applied to 𝐷0 . With this, we approximately mitigate the error-in-variable problem, because the mean prediction excludes the measurement error as best as possible. This gets better, the more data is available, i.e. the more precise the estimation is. This is now repeated until the last stage. This initial holistic model is now used to establish multi-stage bootstrap uncertainty in the second training phase. For each bootstrap loop iteration 𝑏, we apply the first initial stage model 𝑀1 on 𝐷0 as inputs to sample a bootstrap data set 𝐷1𝑏 that has the same dimensions as the original data set 𝐷1 and train a bootstrap model 𝑀1𝑏 on 𝐷0 and 𝐷1𝑏 . This model is used to provide a ̅1𝑏 and, in the second stage, model 𝑀2 is used to generate a bootstrap data set mean prediction 𝐷 𝑏 ̅1 as inputs. On this generated pair (𝐷 ̅1𝑏 , 𝐷2𝑏 ) we then train a second stage bootstrap 𝐷2 with 𝐷 model 𝑀2𝑏 . This is now repeated until the last stage and for all bootstrap iterations. Here, we used the same strategy as in the training of the initial multi-stage model to approximately eliminate the error-in-variable problem. Another reason for the strategy of forwarding the bootstrap mean predictions is that the bootstrap data dependency is carried through the whole chain reintroducing the dependencies between model coefficients, i.e., the second problem of holistic process modelling is addressed. The next step is the prediction of a future state by the holistic process model. This is described in Pseudo Code 3 which also uses the notations of Table 2. The algorithm is again divided into two phases. Given a future first stage input 𝐷0 , the first phase computes the mean prediction by simply executing the initial model chain. The second phase evaluates all bootstrap multi-stage models at 𝐷0 and approximates any desired statistical property from the obtained bootstrap samples in whichever way seems suitable. In our case we used a studentized prediction root [20]. PSEUDO C ODE 3: For 𝑖 in 1 till 𝑁: If 𝑖 = 1 ̅1 to 𝑀1 (𝐷0 , "𝑚") Set 𝐷 Else ̅𝑖 to 𝑀𝑖 (𝐷 ̅𝑖−1 , ”𝑚”) Set 𝐷 For 𝑏 in 1 till 𝐵: For 𝑖 in 1 till 𝑁: If 𝑖 = 1

Else:

̅1𝑏 to 𝑀1𝑏 (𝐷0 , ”𝑚”) Set 𝐷 Set 𝐷1𝑏 to 𝑀1𝑏 (𝐷0 , "𝑠") 𝑏 ̅𝑖𝑏 to 𝑀𝑖𝑏 (𝐷 ̅𝑖−1 Set 𝐷 , "𝑚") 𝑏 𝑏 ̅𝑏 Set 𝐷𝑖 to 𝑀𝑖 (𝐷𝑖−1 , "𝑠")

It is noteworthy that the mean predictions with statistical uncertainty intervals cannot only be computed at the last stage 𝑛 of the holistic process, but also at any other intermediate stage 𝑖 ∈ {2, … , 𝑛 − 1}. The next subsection investigates the just presented algorithm on the LQLQ and LQLQLQLQ models.

4.2 Simulation Studies applying the model-based Bootstrap We evaluated the holistic bootstrap on the LQLQ and LQLQLQLQ endpoints from Section 3.2, using the same true stage functions, measurement-error standard deviation 0.5, sample sizes 𝑚 = 7, 15 and 30. There were 1000 outer repetitions and B = 1000 refitted bootstrap chains per repetition. The point predictor was the noise-free initial fitted-to-observed composition, not an average of bootstrap predictions. Relative mean bias and interval-width ratios use the same true endpoint ranges and normal-interval reference as the MCSim study. Figure 4 shows substantially smaller relative mean biases than the MCSim results. For LQLQ, the ranges over the grid are −1.34% to +1.11% at 𝑚 = 7, −0.73% to +1.01% at 𝑚 = 15 and −0.38% to +0.57% at 𝑚 = 30 . For LQLQLQLQ, they are −6.98% to +7.87%, −4.97% to +4.84% and −3.16% to +2.65%, respectively. The longer chain therefore retains visible, input-dependent bias, although its magnitude decreases over the investigated sample sizes. These finite-sample results are compatible with improved estimation as 𝑚 increases, but do not establish general consistency. Figure 5 shows coverage much closer to the 90% target than the main MCSim study. For LQLQ, pointwise coverage ranges are 88.26–90.02%, 88.76–90.76% and 89.30–90.72% for 𝑚 = 7, 15 and 30. For LQLQLQLQ, the corresponding ranges are 88.18–92.16%, 86.04–90.64% and 88.32– 90.68%. Figure 6 shows considerably less width inflation than MCSim. For LQLQ, mean width ratios range from 1.43 to 1.73 at 𝑚 = 7, 1.14 to 1.38 at 𝑚 = 15 and 1.06 to 1.23 at 𝑚 = 30. For LQLQLQLQ, they range from 1.97 to 2.78, 1.36 to 1.65 and 1.16 to 1.42. Widths decrease towards the knownparameter normal reference as 𝑚 increases, and the longer chain generally requires wider intervals. A ratio above 1 is expected to some extent because the oracle reference excludes parameter-estimation uncertainty; it does not alone imply conservatism. Indeed, Figure 5 contains local undercoverage despite ratios above 1. At 𝑚 = 30, the eight-stage bootstrap ratios of 1.16–1.42 contrast with 4.11–4.97 for observed-to-observed MCSim, while achieving coverage far closer to the requested level. This is evidence of a useful improvement within the tested design.

FIGURE 4: RELATIVE MEAN BIAS OF THE INITIAL FITTED -TO-OBSERVED CHAIN USED BY THE FIXED CENTER STUDENTIZED BOOTSTRAP . P ANELS A AND B SHOW LQLQ AND LQLQLQLQ; BLUE , ORANGE AND GREEN DENOTE M = 7, 15 AND 30. BIAS IS NORMALIZED BY THE TRUE ENDPOINT -MEAN RANGE OVER THE EVALUATION GRID AND AVERAGED OVER 1000 OUTER REPETITIONS . R IBBONS ARE POINTWISE 90% PERCENTILE BOOTSTRAP CONFIDENCE INTERVALS FROM 2,000 RESAMPLES OF THE OUTER REPETITIONS .

FIGURE 5: P OINTWISE COVERAGE OF NOMINAL 90% FIXED -CENTER STUDENTIZED BOOTSTRAP PREDICTION INTERVALS . P ANELS A AND B SHOW LQLQ AND LQLQLQLQ; COLORS IDENTIFY M = 7, 15 AND 30 AS IN FIGURE 4. COVERAGE IS THE PROPORTION OF 1000 REPETITIONS IN WHICH ONE FRESH NORMAL ENDPOINT OBSERVATION LIES INSIDE THE FITTED INTERVAL . T HE BLACK LINE IS THE 90% TARGET. RIBBONS ARE POINTWISE 90% OUTER -REPETITION BOOTSTRAP CONFIDENCE INTERVALS ; B = 5,000 REFITTED CHAINS PER REPETITION .

FIGURE 6: MEAN WIDTH RATIOS OF THE FIXED -CENTER STUDENTIZED BOOTSTRAP PREDICTION INTERVALS . P ANELS A AND B SHOW LQLQ AND LQLQLQLQ; COLORS IDENTIFY SAMPLE SIZES AS IN

FIGURE 4. W IDTHS ARE DIVIDED BY THE TRUE CENTRAL 90% NORMAL-INTERVAL WIDTH AND AVERAGED OVER 1000 OUTER REPETITIONS . T HE DASHED LINE IS THIS KNOWN - PARAMETER REFERENCE . R IBBONS ARE POINTWISE 90% CONFIDENCE INTERVALS FROM 2,000 RESAMPLES OF THE OUTER REPETITIONS ; B = 5,000 BOOTSTRAP REFITS . A bivariate 6-stage illustration is retained in the supplement: Figures S4–S9 describe its bivariate stage functions, while Figures S9–S12 contain results from the holistic bootstrap. The example with 𝑛 = 60 batches shows the same tendencies for the final output, but also for the intermediate output after UO2: minimal relative bias (below 15% even after UO6), coverage is approximately met and the interval width is close to 100% except at the edges of the design space (Figures S6-8). The latter is a wanted behavior since models tend to be less precise at the edges to extrapolation, i.e., it is also seen for UO1, where no chaining has happened yet.

5 Conclusion Chaining separately fitted stage models by forwarding complete predictive draws can confuse uncertainty in the estimated latent stage means with measurement error in a future observed output. The analytical intercept-to-linear example shows finite but also asymptotic overcoverage and width inflation, even though its mean remains unbiased. Introducing a quadratic second stage adds a finite but also asymptotic mean bias. The four- and eight-stage simulations show that the resulting distortions can be substantial, input-dependent and not eliminated by larger training samples over the investigated range. In both investigations, a crucial observation was that the first as well as second (and at some point also higher) moments can get undefined. The latent-input diagnostic demonstrates that training-input measurement error alone does not explain the failure. It also shows local undercoverage, so the special-case linear overcoverage result must not be generalized to arbitrary nonlinear chains. The holistic retraining via a model-based bootstrap provides a substantially better-performing alternative for these simulated chains. It first estimates a fitted-to-observed mean chain, then refits that estimator on stage-specific pseudodata centered at fixed initial fitted values. Relative mean bias in the four-stage chain is approximately within −1.34% to +1.11% at 𝑚 = 7 and −0.38% to +0.57% at 𝑚 = 30. For the eight-stage chain, the corresponding ranges are −6.98% to +7.87% and −3.16% to +2.65%. At 𝑚 = 30, mean width ratios of 1.16–1.42 for the eight-stage bootstrap compare with 4.11–4.97 for observed-to-observed MCSim. These gains concern the combined change in point estimation and interval construction, not studentization in isolation. The bootstrap is nevertheless not uniformly exact. Pointwise coverage across the evaluated scenarios ranges from 86.04% to 92.16% at a nominal 90%, with visible local undercoverage in the eight-stage chain. Widths remain above the known-parameter reference, especially for small samples and longer chains. Model-form misspecification, correlated or heteroscedastic errors, other designs and other model families have not been evaluated. Fitted-mean regressors remain estimated and can therefore mitigate rather than generally eliminate errors-in-variables effects. Within the stated scope, the results support the holistic bootstrap as a promising basis for more informative endpoint prediction.

6 Literature [1] K. Iranshahi, J. Brun, T. Arnold, T. Sergi and U. C. Müller, "Digital twins: Recent advances and future directions in engineering fields," Intelligent Systems with Applications, 2025. [2] M. Shariatifar, M. S. Rizi, R. Sotudeh-Gharebagh, R. Zarghami and N. Mostoufi, "On digital twins in bioprocessing: Opportunities and limitations," Process Biochemistry, 2025. [3] J. Smiatek, A. Jung and E. Bluhmki, "Towards a Digital Bioprocess Replica: Computational Approaches in Biopharmaceutical Development and Manufacturing," Trends in Biotechnology, 2020. [4] C. Taylor, L. Marschall, M. Kunzelmann, M. Richter, F. Rudolph, J. Vajda, B. Presser, T. Zahel, J. Studts and C. Herwig, "Integrated Process Model Applications Linking Bioprocess Development to Quality by Design Milestones," bioengineering, 2017. [5] A. S. Rathore, N. Kateja and D. Kumar, "Process integration and control in continuous bioprocessing," Current Opinion in Chemical Engineering, 2018. [6] L. Rüschendorf, Mathematische Statistik, Springer, 2015. [7] H. Xia, Y. Ding, W. Chen, J. Zhang, H. Li, H. Zhao and R. Song, "Measurement uncertainty evaluation with small samples: A review and prospect," Measurement, 2026. [8] B. D. Hall, "Modelling Metrological Traceability," Metrology , 2025. [9] L. M. Herrera, T. Eilert, I.-T. Ho, M. Matysik, M. Laussegger, R. Guderlei, B. Schrantz, A. Jung, E. Bluhmki and J. Smiatek, "Holistic Process Models: A Bayesian Predictive Ensemble Method for Single and Coupled Unit Operation Models," processes, 2022. [10] J. C. f. G. i. Metrology, Guide to the expression of uncertainty in measurment, 2023. [11] N. Y. S. Garg and D. Aswal, "Monte Carlo Simulation in Uncertainty Evaluation: Strategy, Implications and Future Prospects," MAPAN, 2019. [12] A. Razumić, B. Runje, V. Alar, B. Štrbac and Z. Trzun, "A Review of Methods for Assessing the Quality of Measurement Systems and Results," Applied Sciences, 2025. [13] C. Taylor, B. Pretzner, T. Zahel and C. Herwig, "Architectural and Technological Improvements to Integrated Bioprocess Models towards Real-Time Applications," Bioengineering, 2022. [14] T. Oberleitner, T. Zahel, B. Pretzner and C. Herwig, "Holistic Design of Experiments Using an Integrated Process Model," Bioengineering, 2022. [15] J. Kager and C. Herwig, "Monte Carlo-Based Error Propagation for a More Reliable Regression Analysis across Specific Rates in Bioprocesses," Bioengineering, 2021.

[16] M.-g. Xie and K. Singh, " Confidence Distribution, the Frequentist Distribution Estimator of a Parameter: A Review," International Statistical Review, 2013. [17] T. Schweder and N. Hjort, Confidence, Likelihood, Probability, Cambridge University Press, 2016. [18] R. J. Carroll, D. Ruppert, L. A. Stefanski and C. Crainiceanu, Measurement Error in Nonlinear Models: A Modern Perspective, CRC Press, 2006. [19] J. Shen, R. Y. Liu and M.-g. Xie, "Prediction with confidence—A general framework for predictive inference," Journal of Statistical Planning and Inference, 2017. [20] L. Pan and D. N. Politis, "Bootstrap prediction intervals for linear, nonlinear and nonparametric autoregressions.," Journal of Statistical Planning and Inference, no. 177, pp. 1-27, 2016.

Supplements S1.

MCSim diagnostic with latent training inputs

This diagnostic uses the same simulated process, sample sizes, measured responses, prediction grid and MCSim propagation scheme as the main study. Only the training predictors change: later-stage models use the true preceding latent output rather than its measured value. The diagnostic therefore removes the errors-in-variables component from the training inputs without removing the predictive noise that MCSim subsequently propagates. Its purpose is not to propose an implementable oracle method, but to test whether training-input measurement error is the sole source of the observed distortions. Figure S1 shows persistent discrepancies even with latent training inputs. For LQLQ, mean bias remains negative over the grid: −38.9% to −13.8% at 𝑚 = 7, −19.9% to −7.5% at 𝑚 = 15 and −14.5% to −5.9% at 𝑚 = 30. For LQLQLQLQ at 𝑚 = 30, mean bias reaches −378.9% near the right boundary. Figure S2 shows LQLQ coverage of approximately 90.7–100%. The eight-stage diagnostic instead ranges from approximately 82.8% to 100% across sample sizes, with a localized trough near the upper end of the input range. Figure S3 shows that undercoverage can coexist with very large widths: at 𝑚 = 7, eight-stage mean widths range from 11.11 to 1733.35 times the reference; even at 𝑚 = 30 they range from 5.56 to 27.83 times. These are distinct grid-dependent features, not a claim that the maximum width occurs where coverage is minimal. Nonlinearity changes the location, shape and tails of the predictive distribution; wider intervals alone need not restore correct coverage. Taken together, these results rule out training-input measurement error as the sole explanation of the MCSim failure in these scenarios.

FIGURE S1: RELATIVE ENDPOINT DEVIATIONS UNDER LATENT -TO-OBSERVED TRAINING . PANELS A AND B SHOW LQLQ AND LQLQLQLQ; COLORS DENOTE 𝒎 = 𝟕, 𝟏𝟓 AND 𝟑𝟎. N ORMALIZATION , AGGREGATION AND POINTWISE 90% BOOTSTRAP BANDS ARE AS IN FIGURE 1.

FIGURE S2: P OINTWISE COVERAGE OF NOMINAL 90% MCS IM PREDICTION INTERVALS WITH LATENT TO - OBSERVED TRAINING . P ANELS A AND B SHOW LQLQ AND LQLQLQLQ; COLORS IDENTIFY 𝒎 = 𝟕, 𝟏𝟓 AND 𝟑𝟎. THE BLACK LINE IS THE 90% TARGET. RIBBONS ARE POINTWISE 90% OUTER REPETITION BOOTSTRAP CONFIDENCE INTERVALS . T HE LOCALIZED UNDERCOVERAGE IN PANEL B SHOWS THAT THE OVERCOVERAGE OF THE MAIN OBSERVED -TO- OBSERVED STUDY IS NOT UNIVERSAL .

FIGURE S3: MEAN MCS IM INTERVAL -WIDTH RATIOS WITH LATENT -TO-OBSERVED TRAINING . PANELS A AND B SHOW LQLQ AND LQLQLQLQ; COLORS IDENTIFY 𝒎 = 𝟕, 𝟏𝟓 AND 𝟑𝟎. W IDTHS ARE NORMALIZED BY THE TRUE 90% NORMAL -INTERVAL WIDTH . R IBBONS ARE POINTWISE 90% OUTER REPETITION BOOTSTRAP CONFIDENCE INTERVALS . O THER SIMULATION SETTINGS ARE THE SAME AS IN FIGURE 3.

S2 Bivariate 3-UO Example and Simulation Results

FIGURE S4: FUNCTION OF U NIT OPERATION 1 MAPPING THE INPUT VARIABLES TO THE OUTPUT VARIABLES .

FIGURE S5: FUNCTION OF U NIT OPERATION 2 MAPPING THE INPUT VARIABLES TO THE OUTPUT VARIABLES .

FIGURE S6: FUNCTION OF U NIT OPERATION 3 MAPPING THE INPUT VARIABLES TO THE OUTPUT VARIABLES .

FIGURE S7: FUNCTION OF U NIT OPERATION 4 MAPPING THE INPUT VARIABLES TO THE OUTPUT VARIABLES .

FIGURE S8: FUNCTION OF U NIT OPERATION 4 MAPPING THE INPUT VARIABLES TO THE OUTPUT VARIABLES .

FIGURE S9: FUNCTION OF U NIT OPERATION 4 MAPPING THE INPUT VARIABLES TO THE OUTPUT VARIABLES .

FIGURE S10: RELATIVE BIAS FOR ALL UO S AND VARIABLES FOR THE APPLICATION OF THE BOOTSTRAP ALGORITHM AND 𝒏 = 𝟔𝟎 INPUT DATA POINTS ON AN EQUALLY SPACED GRID OVER THE DESIGN SPACE .

FIGURE S11: C OVERAGE FOR ALL UO S AND VARIABLES FOR THE APPLICATION OF THE BOOTSTRAP ALGORITHM AND 𝒏 = 𝟔𝟎 INPUT DATA POINTS .

FIGURE S12: PI WIDTH RATIO FOR ALL UO S AND VARIABLES FOR THE APPLICATION OF THE BOOTSTRAP ALGORITHM AND 𝒏 = 𝟔𝟎 INPUT DATA POINTS .

Related documents

• Conformal ordering: Distribution-free predictive inference for orders and rank#1491948
• GBTSVR: Twin Support Vector Regression via Granular Ball Computing#1354732
• StepCOPS: Closed-Testing Lower-Tail Certificates for Language-Model Policy Selection#1077477
• 50 Years of Successful Predictive Modeling Should Be Enough: Lessons for Philosophy of Science#701684
• Enhancing Language Model Reasoning with Structured Multi-Level Modeling#1536446