A Patient-Anchored, Physics-Informed Gradient-Boosting Virtual Sensor for Arterial Oxygen Saturation Across a Nine-Hour Blood-Gas Gap

1. Introduction

Effective critical care depends on having timely and accurate knowledge of the patient’s physiological state, yet intensive care units operate under persistent structural constraints that systematically undermine this requirement [1]. High patient to nurse ratios, equipment saturation during surge periods, and the physical impossibility of continuous bedside observation for every admitted patient mean that clinically significant changes in physiology can go undetected for extended periods [2]. Laboratory based measurements remain the clinical standard for several critical physiological variables, including cardiac output [3] or arterial oxygen saturation ${\mathrm{SaO}}_2$ via co-oximetry [4,5], yet when the samples of intensive-care patients are analysed in a central laboratory the median turnaround time between blood collection and result availability is 56 minutes, so that in practice some clinical decisions are grounded in delayed information [6].

This growing divergence between the current state and the most recent measurement is not an edge case but a structural property of any monitoring system that relies on intermittent invasive sampling [7,8]. Retrospective analyses of large ICU databases have documented that missing or delayed values in key physiological variables are systematically associated with suboptimal alarm management and delayed recognition of clinical deterioration [7,9]. For patients receiving mechanical ventilation, however, a continuously available high-resolution data source already exists: the ventilator itself, which continuously records a rich stream of respiratory signals that reflect the patient’s respiratory mechanics and gas-exchange capacity in real time [2,10]. The predictive value of this stream has already been established for within-ventilator prediction tasks, Qiu et al. [11] demonstrated that a long short-term memory relation network trained on MIMIC-III can predict patient-specific ventilator parameters in real time. If ventilator signals are sufficiently informative to predict future ventilator parameters, they should also carry information relevant to predicting other physiological variables, such as ${\mathrm{SaO}}_2$, as demonstrated by prior work linking ventilator settings to patient oxygenation outcomes [12,13]. This follows the soft-sensor paradigm, where continuously available signals are used to estimate intermittently measured variables [14].

Recent advances in soft sensing have shown that capturing long-range temporal structure is essential, [14] propose a temporal convolutional network with both spatial and temporal attention mechanisms, showing that temporally aware models are preferred for such task.

In the medical domain, the soft-sensor paradigm has been successfully applied across a broad range of continuous physiological estimation problems, establishing a strong precedent for the clinical feasibility of the approach [15]. Beyond oxygenation-related variables, soft-sensor approaches have been successfully explored for glucose monitoring, cardiovascular assessment, respiratory monitoring, intracranial pressure estimation, and wearable health applications, consistently demonstrating the potential of machine-learning models to infer clinically relevant physiological states from routinely collected signals [4,16,17,18,19,20,21,22,23].

Recent work in ICU time-series modeling emphasizes the importance of explicitly representing temporal dynamics under irregular sampling. Lee et al. [8,24] show that incorporating inter-event time gaps into recurrent architectures improves predictive performance by allowing models to exploit temporal spacing as an informative signal rather than treating it as missingness. Similarly, adaptive recurrent formulations that condition hidden state updates on elapsed time further improve modeling of heterogeneous physiological trajectories in critical care settings [25]. Patharkar et al. [7] summarize these findings across biomedical time-series applications, highlighting that methods accounting for irregular sampling consistently outperform standard sequence models.

Among physiological targets in critical care, oxygenation estimation has received growing attention, yet existing approaches leave the core problem of this study unresolved. Woo et al. [26] trained Random Forest and XGBoost models on a large-scale multi-centre ICU dataset to predict ${\mathrm{SaO}}_2$ from ${\mathrm{SpO}}_2$ and clinical covariates, reducing estimation error substantially over the pulse oximetry baseline, but the approach corrects an available ${\mathrm{SpO}}_2$ signal rather than operating in its absence. Cabanas et al. [27] further highlight the regulatory and technical limitations of AI-enhanced pulse oximetry, reinforcing that ${\mathrm{SpO}}_2$-dependent correction methods have fundamental constraints in failure scenarios.

The closest line of work to ours estimates oxygenation directly from ventilator signals without relying on pulse oximetry. Where Ghazal et al. [28] trained a bootstrap-aggregated decision tree to predict ${\mathrm{SpO}}_2$ categories five minutes after a ventilator setting change, achieving an AUC below 0.75 a result that also reveals the core limitation: a five-minute classification horizon is far too narrow to address the multi-hour sensor-absence scenario we target. Dong et al. [5] take a closer approach, using a non-linear additive model to estimate blood gas pH and ${\mathrm{PaCO}}_2$ between invasive draws from ventilator-derived inputs, achieving clinically acceptable accuracy, but their target variable is not ${\mathrm{SaO}}_2$ and the gap length is never varied or evaluated as a modeling parameter. Hser et al. [29] incorporate oxygenation estimation within a broader ICU respiratory monitoring system using gradient boosted models, predicting respiratory failure events up to ten hours ahead, but their estimator still takes ${\mathrm{SpO}}_2$ as a continuous input, making it inapplicable when pulse oximetry is unavailable. Finally, Si et al. [30] estimate ${\mathrm{SpO}}_2$ from facial video using Random Forest and Support Vector Regression, demonstrating non-contact sensing potential but under controlled low-acuity conditions incompatible with the mechanically ventilated ICU patient.

Across these works, three gaps remain open: no approach models prediction accuracy as a function of elapsed time since the last measurement, none operates without assuming ${\mathrm{SpO}}_2$ availability, and none separates what a physical prior contributes from what the learned trees contribute under patient-level cross-validation. This paper addresses the three, building on the oxygen conservation law of our bioprocess soft sensors [31,32], on its transfer to human respiration [33,34], and on the maximum-relative-entropy weighting of our earlier estimators [32,35,36]. It bridges a nine-hour absence of the co-oximetry measurement from ventilator and gas-rig time series with gradient-boosted trees; it works without pulse oximetry; and it is reported as three explicit implementation scenarios of one code base, the unmodified library through our pipeline (X0), the library with our anchor blend and the cumulative oxygen uptake rate over the gap (XP), and the latter with a patient-level calibrated anchor (v4), so that every improvement is attributed to the code that produced it, every equation is matched to the line that implements it, and an exact decomposition of each prediction says how much of it is physics and how much is trees. The rest of the paper is organized as follows: Section 2 describes the data, the three scenarios, the shared boosted ensemble with a numeric trace of its leaves, the physics-informed terms with their derivations and code, and the evaluation protocol; Section 3 reports the three scenarios on identical folds, the anchor coefficients, the anatomy of the predictions and the ablations; Section 4 and Section 5 discuss and conclude.

2. Materials and Methods

The objective of this study is to predict $S_a{O}_{2}m$, the arterial oxygen saturation measured by co-oximetry, from hourly processed clinical variables, respiratory parameters and experimental gas-rig transducers, nine hours after the last available measurement. A principal characteristic of the study is the translation of soft-sensor technologies from bioprocess engineering to human physiological monitoring: as demonstrated by Kemesis et al. [34], industrial off-gas sensors can be adapted to the non-invasive continuous tracking of human respiratory gas exchange.

We define the dataset as the matrix $\mathbf{X}\in {\mathbb{R}}^{n_i\times n_m}$, where $n_i$ is the number of hourly samples and $n_m$ is the number of features that the boosted trees receive (its value depends on the scenario, Section 2.2), and the i-th row ${\mathbf{x}}_i\in {\mathbb{R}}^{n_m}$ is the feature vector of the i-th sample. The target vector is $\mathbf{y}\in {\mathbb{R}}^{n_i}$, the $S_a{O}_{2}m$ values of every patient-hour. The mapping is

$${\widehat{y}}_i=\varphi ({\mathbf{x}}_i,\theta )$$

(1)

where ${\widehat{y}}_i$ is the predicted $S_a{O}_{2}m$ of the i-th sample and $\varphi (·,\theta )$ is a regressor with learned parameters $\theta$. Three implementations of $\varphi$ are compared in this paper; they share one code base and differ in which of our additions are active.

2.1. Dataset and Preprocessing

The study cohort comprises $n_p=31$ mechanically ventilated patients, totaling $4,902$ recorded hourly entries. Of these, $n_i=4,017$ hours have a measured target; all of them are predicted out of fold under five-fold patient-grouped cross-validation, and the results are reported on the $3,730$ hours of 30 patients whose inputs contain no $S_a{O}_{2}m$ value younger than nine hours (Section 3.7), with the values on all $4,017$ hours listed alongside. In 322 patient-hours (644 rows, five patients) two blood gas measurements were documented within the same nominal hour. These rows are retained in the reported numbers; a sensitivity analysis keeping only the last measurement per hour is reported in Section 3.1.

The target variable is discrete and integer-valued, spanning 73 distinct levels from $20\%$ to $111\%$ (median $96\%$, 1st percentile $52\%$; 350 observations fall below $85\%$, and a single observation exceeding $100\%$ was retained as an authentic clinical recording). This discrete distribution is directly relevant to the mechanics of the persistence baseline (Section 3.6). All features entering the decision trees were standardized using Z-score normalization fitted strictly on the training folds; in contrast, all physiological quantities supplied directly to the analytical anchor (Section 2.4) were preserved in their raw, natural physical units.

Figure 1. Temporal problem formulation and prediction horizon. The model estimates the current oxygenation state $y_t$ by bridging a nine-hour absence of the co-oximetry measurement. The last available measurements $y_{t-9}$ and $y_{t-10}$ supply the physiological baseline; the continuous ventilator and gas-sensor channels supply the real-time information available at hour t.

Figure 1. Temporal problem formulation and prediction horizon. The model estimates the current oxygenation state $y_t$ by bridging a nine-hour absence of the co-oximetry measurement. The last available measurements $y_{t-9}$ and $y_{t-10}$ supply the physiological baseline; the continuous ventilator and gas-sensor channels supply the real-time information available at hour t.

We define a prediction horizon of nine hours. The $S_a{O}_{2}m$ history at $t-9$ and $t-10$ hours is included to establish the patient’s baseline prior to the simulated sensor gap, ensuring the model remains grounded in the patient’s physiological state without using future information. More precisely, the feature value provided at hour t represents the most recent observation recorded at or before$t-9\mathrm{h}$. Consequently, the actual elapsed time since that baseline reading was taken is at least nine hours and varies depending on clinical sampling irregularity; its empirical distribution is analyzed in Section 3.6.

2.1.1. Signals and Their Lags in Each Scenario

Six real-time signals are engineered from the ventilator and the gas rig and lagged; three target-derived and two static descriptors complete the input. Table 1 lists the six signals with the lag hours that each scenario builds and, of those, the ones the trees are allowed to see. Table 2 lists the remaining columns. The counts differ between the scenarios: the trees of X0 and XP receive $n_m=55$ columns (six signals at hour t, their 44 lagged copies, and the five columns of Table 2 that are tree inputs), whereas the trees of v4 receive $n_m=19$ (five signals at t with two lags each, plus four of those five columns). In v4 the pruning of the lag lists was decided by tolerance tests on the out-of-fold macro MAPE (Section 3); the signal O2diffxPower (a signal which is proportional to Oxygen Uptake Rate, OUR) is hidden from the trees in v4 because it feeds the physics terms instead.

Table 1. Dynamic signals and their lag hours. Each cell lists the hours t and $t-\ell$ at which the signal enters the data as a column (data_preparation.py); the count in brackets is the number of such columns. In X0 and XP every column is a tree input. For O2diffxPower the listed hours are also the sampling hours $\mathcal{L}$ of the cumulative oxygen uptake rate, Eq. (16): the available values at these hours are averaged and the average is multiplied by the gap D and divided by the body mass; in v4 these columns serve only that purpose and are not tree inputs. O2conc: inspired O₂ concentration at the mask inlet (%). GM102B_Volts: electrochemical gas-sensor output voltage (V). CO2xPower: expired CO₂ concentration times the ventilator power integral. O2xPress: O2conc times airway pressure. O2xPressplusPEEP: O2conc times pressure plus PEEP, a delivery composite. O2diffxPower: inlet-minus-expiratory O₂ concentration times the power integral, proportional to the oxygen uptake rate; it is the oxygen uptake rate signal $\mathrm{OUR}\left(t\right)$ of Section 2.4.

Signal	Hours in X0 and XP (all tree inputs)	Hours in v4	Of these, tree inputs in v4
O2conc	t, $t-1,3,4,5,6$ (6 cols)	t, $t-1,t-6$	t, $t-1,t-6$ (3 cols)
GM102B_Volts	t, $t-1,2,3,4,5,7,8$ (8 cols)	t, $t-1,t-8$	t, $t-1,t-8$ (3 cols)
CO2xPower	t, $t-1,2,3,4,5$ (6 cols)	t, $t-1,t-5$	t, $t-1,t-5$ (3 cols)
O2xPress	t, $t-1,2,3,4,5,6,8,11,12$ (10 cols)	t, $t-1,t-12$	t, $t-1,t-12$ (3 cols)
O2xPressplusPEEP	t, $t-1,2,4,5,6,7,8,9,10,11$ (11 cols)	t, $t-1,t-11$	t, $t-1,t-11$ (3 cols)
O2diffxPower	t, $t-1,2,4,5,7,8,10,11$ (9 cols; also the sampling hours of Eq. (16) in XP)	t, $t-1,2,4,7$ (the sampling hours of Eq. (16))	none
tree inputs from this table: 50			15

Table 1. Dynamic signals and their lag hours. Each cell lists the hours t and $t-\ell$ at which the signal enters the data as a column (data_preparation.py); the count in brackets is the number of such columns. In X0 and XP every column is a tree input. For O2diffxPower the listed hours are also the sampling hours $\mathcal{L}$ of the cumulative oxygen uptake rate, Eq. (16): the available values at these hours are averaged and the average is multiplied by the gap D and divided by the body mass; in v4 these columns serve only that purpose and are not tree inputs. O2conc: inspired O₂ concentration at the mask inlet (%). GM102B_Volts: electrochemical gas-sensor output voltage (V). CO2xPower: expired CO₂ concentration times the ventilator power integral. O2xPress: O2conc times airway pressure. O2xPressplusPEEP: O2conc times pressure plus PEEP, a delivery composite. O2diffxPower: inlet-minus-expiratory O₂ concentration times the power integral, proportional to the oxygen uptake rate; it is the oxygen uptake rate signal $\mathrm{OUR}\left(t\right)$ of Section 2.4.

Signal	Hours in X0 and XP (all tree inputs)	Hours in v4	Of these, tree inputs in v4
O2conc	t, $t-1,3,4,5,6$ (6 cols)	t, $t-1,t-6$	t, $t-1,t-6$ (3 cols)
GM102B_Volts	t, $t-1,2,3,4,5,7,8$ (8 cols)	t, $t-1,t-8$	t, $t-1,t-8$ (3 cols)
CO2xPower	t, $t-1,2,3,4,5$ (6 cols)	t, $t-1,t-5$	t, $t-1,t-5$ (3 cols)
O2xPress	t, $t-1,2,3,4,5,6,8,11,12$ (10 cols)	t, $t-1,t-12$	t, $t-1,t-12$ (3 cols)
O2xPressplusPEEP	t, $t-1,2,4,5,6,7,8,9,10,11$ (11 cols)	t, $t-1,t-11$	t, $t-1,t-11$ (3 cols)
O2diffxPower	t, $t-1,2,4,5,7,8,10,11$ (9 cols; also the sampling hours of Eq. (16) in XP)	t, $t-1,2,4,7$ (the sampling hours of Eq. (16))	none
tree inputs from this table: 50			15

Table 2. Target-derived and static columns. "Tree" says whether the boosted trees receive the column in that scenario; "anchor" says whether the column feeds the anchor of Section 2.4 instead. runmean_SaO2m exists only in v4.

Column	Description	X0, XP	v4
SaO2m_lag_9h	most recent $S_a{O}_{2}m$ measured at or before $t-9$ h (%); the persistence value $y\left({\tau }_i\right)$ of Eq. (17)	tree and anchor	tree and anchor
SaO2m_lag_10h	most recent $S_a{O}_{2}m$ measured at or before $t-10$ h (%)	tree	tree
SaO2m_first	the patient’s first $S_a{O}_{2}m$ measurement (%)	tree	tree
runmean_SaO2m	mean of all $S_a{O}_{2}m$ measurements at or before $t-9$ h, Eq. (20)	not built	anchor only
count_hours	elapsed treatment time (h)	tree	tree
Svoris	body mass (kg); also normalises the cumulative OUR term	tree	anchor scale only
tree inputs from this table		5	4
total tree inputs $n_m$		55	19

Table 2. Target-derived and static columns. "Tree" says whether the boosted trees receive the column in that scenario; "anchor" says whether the column feeds the anchor of Section 2.4 instead. runmean_SaO2m exists only in v4.

Column	Description	X0, XP	v4
SaO2m_lag_9h	most recent $S_a{O}_{2}m$ measured at or before $t-9$ h (%); the persistence value $y\left({\tau }_i\right)$ of Eq. (17)	tree and anchor	tree and anchor
SaO2m_lag_10h	most recent $S_a{O}_{2}m$ measured at or before $t-10$ h (%)	tree	tree
SaO2m_first	the patient’s first $S_a{O}_{2}m$ measurement (%)	tree	tree
runmean_SaO2m	mean of all $S_a{O}_{2}m$ measurements at or before $t-9$ h, Eq. (20)	not built	anchor only
count_hours	elapsed treatment time (h)	tree	tree
Svoris	body mass (kg); also normalises the cumulative OUR term	tree	anchor scale only
tree inputs from this table		5	4
total tree inputs $n_m$		55	19

2.2. Three Implementation Scenarios

Every number in this paper is attributed to one of three scenarios of one code base (Python pipeline and a patched XGBoost 2.1.4 C++ library). They are nested: each adds to the previous one and changes nothing else.

• X0—XGBoost. The XGBoost 2.1.4 library, trained by our Python pipeline with the tree inputs of Table 1Table 2 ($n_m=55$), the initial prediction equal to the persistence value $y\left({\tau }_i\right)$, the cost-sensitive weights of Section 2.3, and the hyper-parameters of Table 8. The two C++ terms added for XP are compiled into the library used here as well, but both coefficients are zero ($\beta ={\beta }_1=0$ in Eq. (37)), which makes them identities: X0 produces what unmodified XGBoost 2.1.4 produces from the same pipeline, up to float32 rounding.
• XP — XGBoost with our anchor blend and cumulative OUR term. X0 plus the two C++ terms of Section 2.6 with their configured coefficients $\beta =0.1418$ and ${\beta }_1=0.00139$: the final prediction is pulled toward the persistence value and receives the cumulative oxygen uptake rate over the gap per kilogram of body mass, Eq. (16). Same trees, same $n_m=55$.
• v4 — XP with the patient-level calibrated anchor. XP plus the anchor of Section 2.4.5, which replaces the raw persistence value by a least-squares combination of the persistence value and the patient’s own measurement history, fitted on each training fold; its trees receive $n_m=19$ inputs and the hyperparameters of Table 9. $\beta =0.3677$, ${\beta }_1=0.00124$.

Every code excerpt and equation in this paper follows a unified typographical convention.

In code listings: dark gray denotes unmodified upstream XGBoost (v2.1.4) source code; regular black font indicates our additions active in both XP and v4 (our Python pipeline and C++ extensions); and bold font denotes additions active exclusively in v4. Excerpts of original XGBoost source in Section 2.5 are therefore rendered entirely in gray.

In equations: our contributions carry the label [XP] if active in both XP and v4, and are set in bold with the label [v4] if unique to v4. Standard equations from baseline XGBoost or established physiological literature remain unmarked. Table 3 lists every mathematical symbol with its unit and corresponding code variable.

2.3. The Loss: Cost-Sensitive Sample Weights and the Relative-Entropy Criterion

Two fits in this paper minimise a weighted sum of squared errors over the training rows: the least-squares fit of the patient-level anchor of v4 (Section 2.4.5) and the boosted trees of all three scenarios (Section 2.5). Both use one dimensionless weight per row,

$$w_i={\displaystyle \frac{1-k_w}{|y_i{|}^{\nu }+\epsilon }}+k_w,$$

(2)

where $y_i$ is the target in %, $k_w\in [0,1]$ is a dimensionless floor (custom_kp), $\nu \ge 0$ a dimensionless exponent (custom_power) and $\epsilon ={10}^{-6}$ a stabiliser in ${\%}^{\nu }$ (custom_eps); the three are fields of config.py, class TrainConfig, and the weights are the array sw computed by compute_sample_weights (utils_prediction.py). The weights are the same in the three scenarios up to the values of $k_w$ and $\nu$ (Table 8 and Table 9).

The family (2) descends from our maximum-relative-entropy estimators [31,32,35], where it was written for the biomass concentration of a cultivation with a fixed square in place of the exponent $\nu$, and its single coefficient was named $K_{exp}$; in the notation of this paper that coefficient is $k_w$ and the observed quantity is $y_i$. For a Gaussian density $P_{post}$ of the prediction ${\widehat{y}}_i$ and a Gaussian likelihood $P_{lik}$ centred on the observation $y_i$ with variance ${\sigma }_i^2$ (in ${\%}^2$), the relative entropy of row i is [32]

$$S_i=-\int P_{post}\left(\widehat{y}\right)ln{\displaystyle \frac{P_{post}\left(\widehat{y}\right)}{P_{lik}\left(\widehat{y}\right)}}\mathrm{d}\widehat{y}=-{\displaystyle \frac{{\left({\widehat{y}}_i-y_i\right)}^2}{2{\sigma }_i^2}}+\mathrm{const},$$

(3)

dimensionless, with the variance model [32]

$${\sigma }_i^2={\sigma }_0^2{\displaystyle \frac{|y_i{|}^{\nu }}{1-k_w+k_w{\left|y_i\right|}^{\nu }}},$$

(4)

in which ${\sigma }_0^2$ is a constant scale in ${\%}^2$, $y_i$ enters with its numerical value in %, and $\nu =2$ gives the printed form of the cited papers. Inserting Eq. (4) into Eq. (3) and summing over the training rows gives

Table 3. Symbols, units and names in the code (dp: data_preparation.py; up: utils_prediction.py; lc: src/learner.cc). a.u.: arbitrary units of the gas signals.

Symbol	Meaning	Unit	Name in the code
i, $n_i$	row (one patient-hour), number of rows	—	row of the data frame
$p\left(i\right)$, $n_p$	patient of row i, number of patients	—	patient_id
$t_i$	hour of row i	h	LocalTime
$y_i$	measured $S_a{O}_{2}m$, the target	%	SaO2m
D	nominal gap	h	virtual_sensor_gap$=9$
${\tau }_i$, $D_i$	time and age of the last measurement available at $t_i$, Eq. (17)	h	inside _time_based_lag_last_known (dp)
$y\left({\tau }_i\right)$	last available measurement	%	SaO2m_lag_9h; raw_lag_tr, raw_lag_va (up)
${\mathcal{H}}_i$; ${\overline{y}}_{{\mathcal{H}}_i}$	measurement times available at $t_i-D$; mean of y over them, Eq. (20)	h; %	runmean_SaO2m
$m_i^{kg}$	body mass	kg	Svoris
$\mathrm{OUR}\left(t\right)$	oxygen uptake rate signal at the airway	a.u.	O2diffxPower
$\mathrm{CPR}\left(t\right)$	carbon dioxide production rate signal	a.u.	CO2xPower
$k_{OUR}$, $k_{CPR}$	coefficients of the balance, Eq. (13)	% kg h−1 a.u.−1	—
$\mathcal{L}$	lag set of the OUR signal	h	lags of O2diffxPower
${\mathrm{cOUR}}_i$	cumulative oxygen uptake rate over the gap per kg, Eq. (16)	a.u. h kg−1	integral_O2diffxPower; field custom_integral
${\widehat{y}}_{0,i}$	initial prediction, the anchor	%	raw_base_train, raw_base_val; base_margin; anchor (lc)
$c_0$; $c_1$, $c_2$	coefficients of the patient-level anchor, Eq. (20)	%; —	coef[0]; coef[1], coef[2] (up)
$w_i$	sample weight, Eq. (2)	—	sw
$k_w$; $\nu$; $\epsilon$	floor, exponent and stabiliser of $w_i$, Eq. (2)	—; —; ${\%}^{\nu }$	custom_kp, custom_power, custom_eps
b, $n_b$	boosting round, number of trees	—	rounds of xgb.train
$f_b$, $s_{b,j}$	tree b, score of its leaf j	%	XGBoost internal
$\eta$	learning rate	—	learning_rate
$\lambda$; $\alpha$; $\gamma$	L2, L1 and per-split penalties	—; %; ${\%}^2$	reg_lambda, reg_alpha, gamma
$g_{b,i}$, $h_{b,i}$	gradient and Hessian of row i in round b	%, —	XGBoost internal
${\widehat{y}}_{b,i}$; ${\widehat{y}}_{n_b,i}$	running and raw boosted prediction	%	prediction buffer; preds (lc)
${\tilde{y}}_i$	final prediction, Eq. (37)	%	pred_val, pred_tr (up)
$\beta$; ${\beta }_1$	anchor blend; weight of ${\mathrm{cOUR}}_i$	—; % kg h−1 a.u.−1	custom_b, custom_b1
$K_{L}a$; $p{\mathrm{O}}_2^{*}$, ${\mathrm{PaO}}_2^{*}$	precursor model: mass transfer capacity; maximal oxygen solubility	s−1; mmHg	—

Table 3. Symbols, units and names in the code (dp: data_preparation.py; up: utils_prediction.py; lc: src/learner.cc). a.u.: arbitrary units of the gas signals.

Symbol	Meaning	Unit	Name in the code
i, $n_i$	row (one patient-hour), number of rows	—	row of the data frame
$p\left(i\right)$, $n_p$	patient of row i, number of patients	—	patient_id
$t_i$	hour of row i	h	LocalTime
$y_i$	measured $S_a{O}_{2}m$, the target	%	SaO2m
D	nominal gap	h	virtual_sensor_gap$=9$
${\tau }_i$, $D_i$	time and age of the last measurement available at $t_i$, Eq. (17)	h	inside _time_based_lag_last_known (dp)
$y\left({\tau }_i\right)$	last available measurement	%	SaO2m_lag_9h; raw_lag_tr, raw_lag_va (up)
${\mathcal{H}}_i$; ${\overline{y}}_{{\mathcal{H}}_i}$	measurement times available at $t_i-D$; mean of y over them, Eq. (20)	h; %	runmean_SaO2m
$m_i^{kg}$	body mass	kg	Svoris
$\mathrm{OUR}\left(t\right)$	oxygen uptake rate signal at the airway	a.u.	O2diffxPower
$\mathrm{CPR}\left(t\right)$	carbon dioxide production rate signal	a.u.	CO2xPower
$k_{OUR}$, $k_{CPR}$	coefficients of the balance, Eq. (13)	% kg h−1 a.u.−1	—
$\mathcal{L}$	lag set of the OUR signal	h	lags of O2diffxPower
${\mathrm{cOUR}}_i$	cumulative oxygen uptake rate over the gap per kg, Eq. (16)	a.u. h kg−1	integral_O2diffxPower; field custom_integral
${\widehat{y}}_{0,i}$	initial prediction, the anchor	%	raw_base_train, raw_base_val; base_margin; anchor (lc)
$c_0$; $c_1$, $c_2$	coefficients of the patient-level anchor, Eq. (20)	%; —	coef[0]; coef[1], coef[2] (up)
$w_i$	sample weight, Eq. (2)	—	sw
$k_w$; $\nu$; $\epsilon$	floor, exponent and stabiliser of $w_i$, Eq. (2)	—; —; ${\%}^{\nu }$	custom_kp, custom_power, custom_eps
b, $n_b$	boosting round, number of trees	—	rounds of xgb.train
$f_b$, $s_{b,j}$	tree b, score of its leaf j	%	XGBoost internal
$\eta$	learning rate	—	learning_rate
$\lambda$; $\alpha$; $\gamma$	L2, L1 and per-split penalties	—; %; ${\%}^2$	reg_lambda, reg_alpha, gamma
$g_{b,i}$, $h_{b,i}$	gradient and Hessian of row i in round b	%, —	XGBoost internal
${\widehat{y}}_{b,i}$; ${\widehat{y}}_{n_b,i}$	running and raw boosted prediction	%	prediction buffer; preds (lc)
${\tilde{y}}_i$	final prediction, Eq. (37)	%	pred_val, pred_tr (up)
$\beta$; ${\beta }_1$	anchor blend; weight of ${\mathrm{cOUR}}_i$	—; % kg h−1 a.u.−1	custom_b, custom_b1
$K_{L}a$; $p{\mathrm{O}}_2^{*}$, ${\mathrm{PaO}}_2^{*}$	precursor model: mass transfer capacity; maximal oxygen solubility	s−1; mmHg	—

$$-2\sum _i{S}_i={\displaystyle \frac{1}{{\sigma }_0^2}}\sum _i\left[{\displaystyle \frac{1-k_w}{|y_i{|}^{\nu }}}+k_w\right]{\left(y_i-{\widehat{y}}_i\right)}^2+\mathrm{const},$$

(5)

whose bracket is Eq. (2) with $\epsilon =0$. The weighted squared loss used below is therefore the relative-entropy criterion of [32], generalised by the exponent $\nu$ and the stabiliser $\epsilon$ [36]; $k_w$ is the single coefficient that trades the error relative to $|y_i{|}^{\nu }$ (the MAPE end, $k_w\to 0$) against the squared error (the variance end, $k_w\to 1$). Its measured effect on this cohort is given in Section 3.5.

2.4. From the Oxygen Balance to the Physics-Informed Terms

This subsection derives, from the oxygen conservation law of our earlier work, the two terms that the scenarios add to the boosted trees of Section 2.5: the cumulative oxygen uptake rate of XP, Eq. (16), and the patient-level anchor of v4, Eq. (20). The derived columns integral_O2diffxPower and runmean_SaO2m are removed from the tree inputs (utils_prediction.py, _extract_base_margin) and reach the model only through these terms.

2.4.1. Precursors: the Oxygen Conservation Law of Our Earlier Work

In aerobic bioprocesses, at a stationary oxidative state the oxygen uptake rate of the culture, OUR, equals the oxygen transfer rate from the gas, OTR, both in amount of oxygen per unit volume and time [31,33]:

$$\mathrm{OUR}=\mathrm{OTR}.$$

(6)

In the time domain the same conservation law reads [33,37]: the dissolved oxygen, written as its partial pressure $p{\mathrm{O}}_2$ (mmHg), accumulates by what is transferred minus what is taken up,

$${\displaystyle \frac{\mathrm{d}p{\mathrm{O}}_2\left(t\right)}{\mathrm{d}t}}+\mathrm{OUR}\left(t\right)=\mathrm{OTR}\left(p{\mathrm{O}}_2\left(t\right),t\right),$$

(7)

and the transfer is proportional to the distance from the maximal oxygen solubility $p{\mathrm{O}}_2^{*}$ (mmHg), with the oxygen mass transfer capacity $K_{L}a$ (s⁻¹) [33]:

$$\mathrm{OTR}\left(p{\mathrm{O}}_2\left(t\right),t\right)=K_{L}a\left(p{\mathrm{O}}_2^{*}\left(t\right)-p{\mathrm{O}}_2\left(t\right)\right).$$

(8)

Substituting Eq. (8) into Eq. (7) and replacing the dissolved oxygen of the nutritional medium by the partial pressure of oxygen in arterial blood, ${\mathrm{PaO}}_2$ (mmHg), Survyla et al. [33] and Kemesis et al. [34] obtained the model of arterial oxygen they applied to a human volunteer and to ventilated children,

$${\displaystyle \frac{\mathrm{d}{\mathrm{PaO}}_2\left(t\right)}{\mathrm{d}t}}=K_{L}a\left({\mathrm{PaO}}_2^{*}\left(t\right)-{\mathrm{PaO}}_2\left(t\right)\right)-k_{OUR}\mathrm{OUR}\left(t\right),$$

(9)

in which $\mathrm{OUR}\left(t\right)$ is the airway signal of their gas-analysis device and $k_{OUR}$ its calibration coefficient, so that $k_{OUR}\mathrm{OUR}\left(t\right)$ is in mmHg s⁻¹ (in [33], Table I: $k_{OUR}=1.408$, $K_{L}a=0.4217$ s⁻¹, ${\mathrm{PaO}}_2^{*}=120.2$ mmHg for one volunteer). To use sampled data, Survyla et al. [33] integrated the balance between two neighbouring sample times $t_{i-1}$ and $t_i$ (s) and divided both sides by the interval $\Delta t_{i,i-1}=t_i-t_{i-1}$, a bar denoting the average over the interval (their Eqs. (8)–(9)); applied to Eq. (9) this step gives

$${\displaystyle \frac{{\mathrm{PaO}}_2\left(t_i\right)-{\mathrm{PaO}}_2\left(t_{i-1}\right)}{\Delta t_{i,i-1}}}+k_{OUR}\overline{\mathrm{OUR}}=K_{L}a\left(\overline{{\mathrm{PaO}}_2^{*}}-\overline{{\mathrm{PaO}}_2}\right).$$

(10)

Integrating the rate before using it is also the idea of the cumulative oxygen uptake rate cOUR [31]: the integral of OUR is less noisy than OUR itself. Two features of Eq. (9) do not fit the present task. It is written for partial pressure, and in the trial of [33] the haemoglobin saturation did not change noticeably, so the dissociation curve played no role, whereas here the target is the saturation itself. And it needs $K_{L}a$ and ${\mathrm{PaO}}_2^{*}$, which were identified per person from blood gases, whereas no per-patient calibration is available at prediction time here.

2.4.2. The Balance Used Here: Saturation, Whole-Body Amounts, Body Mass on the Left

Three physiological facts connect Eq. (7) to the saturation y (%). The oxygen content of arterial blood, ${\mathrm{C}}_{\mathrm{a}}{\mathrm{O}}_2$ (mL O₂ per dL of blood), is carried almost entirely by haemoglobin, of concentration Hb (g/dL) [38]:

$${\mathrm{C}}_{\mathrm{a}}{\mathrm{O}}_2=1.34\mathrm{Hb}{\displaystyle \frac{y}{100}}+0.003{\mathrm{PaO}}_2;$$

(11)

at Hb $=14$ g/dL, $y=95\%$ and ${\mathrm{PaO}}_2=100$ mmHg the dissolved term is $0.3$ against $17.8$ mL/dL bound to haemoglobin. The amount of oxygen in the circulating blood is the content times the blood volume $V_b$ (dL), and the blood volume is proportional to the body mass $m^{kg}$ (kg), $V_b=k_b{m}^{kg}$ with $k_b\approx 0.7$ dL/kg in adults [39]. Written for amounts instead of concentrations, Eq. (7) states that the oxygen held in the blood changes by the oxygen taken from the ventilator gas, ${\dot{V}}_{{\mathrm{O}}_2}^{air}$ (mL O₂ h⁻¹), which is what the airway signal measures (inspired minus expired oxygen times the gas flow [40]), minus the oxygen consumed by the tissues, which equals the carbon dioxide production ${\dot{V}}_{{\mathrm{CO}}_2}$ (mL CO₂ h⁻¹) divided by the dimensionless respiratory quotient RQ. With body mass on the left:

$$k_b{m}^{kg}{\displaystyle \frac{\mathrm{d}{\mathrm{C}}_{\mathrm{a}}{\mathrm{O}}_2\left(t\right)}{\mathrm{d}t}}={\dot{V}}_{{\mathrm{O}}_2}^{air}\left(t\right)-{\displaystyle \frac{{\dot{V}}_{{\mathrm{CO}}_2}\left(t\right)}{\mathrm{RQ}}}.$$

(12)

Eq. (12) is the conservation law of Eq. (7), as is Eq. (9) and its sampled form Eq. (10); Table 4 lists how each term of Eq. (10) corresponds to a term of Eq. (12) and of its integrated form, Eq. (15) below. The difference is the closure. Eq. (10) models the transfer from the gradient, $K_{L}a({\mathrm{PaO}}_2^{*}-{\mathrm{PaO}}_2)$, and uses the airway signal as the uptake, the two being equal in the stationary state of Eq. (6). Eq. (12) takes the transfer into the blood from the airway signal, which is what that signal physically measures, and takes the tissue consumption from the CO₂ production; it therefore needs neither $K_{L}a$ nor ${\mathrm{PaO}}_2^{*}$.

Table 4. Correspondence between the precursor balance, Eq. (10), and the balance of this paper, Eqs. (12) and (15).

Role in the balance	Precursor, Eq. (10)	This paper, Eqs. (12), (15)
state variable	${\mathrm{PaO}}_2$ (mmHg), per unit volume	y (%) through ${\mathrm{C}}_{\mathrm{a}}{\mathrm{O}}_2$, Eq. (11), whole body through $k_b{m}^{kg}$
accumulation	change of ${\mathrm{PaO}}_2$ over $\Delta t_{i,i-1}$	change of y over the gap D
transfer into the blood	$K_{L}a(\overline{{\mathrm{PaO}}_2^{*}}-\overline{{\mathrm{PaO}}_2})$, modelled	${\dot{V}}_{{\mathrm{O}}_2}^{air}$, measured at the airway as $\mathrm{OUR}\left(t\right)$
consumption	$k_{OUR}\overline{\mathrm{OUR}}$, from the airway signal	${\dot{V}}_{{\mathrm{CO}}_2}/\mathrm{RQ}$, from $\mathrm{CPR}\left(t\right)$; not used in the three scenarios ($k_{CPR}=0$)
time step	interval between two samples	gap $D=9$ h

Table 4. Correspondence between the precursor balance, Eq. (10), and the balance of this paper, Eqs. (12) and (15).

Role in the balance	Precursor, Eq. (10)	This paper, Eqs. (12), (15)
state variable	${\mathrm{PaO}}_2$ (mmHg), per unit volume	y (%) through ${\mathrm{C}}_{\mathrm{a}}{\mathrm{O}}_2$, Eq. (11), whole body through $k_b{m}^{kg}$
accumulation	change of ${\mathrm{PaO}}_2$ over $\Delta t_{i,i-1}$	change of y over the gap D
transfer into the blood	$K_{L}a(\overline{{\mathrm{PaO}}_2^{*}}-\overline{{\mathrm{PaO}}_2})$, modelled	${\dot{V}}_{{\mathrm{O}}_2}^{air}$, measured at the airway as $\mathrm{OUR}\left(t\right)$
consumption	$k_{OUR}\overline{\mathrm{OUR}}$, from the airway signal	${\dot{V}}_{{\mathrm{CO}}_2}/\mathrm{RQ}$, from $\mathrm{CPR}\left(t\right)$; not used in the three scenarios ($k_{CPR}=0$)
time step	interval between two samples	gap $D=9$ h

The assumptions are: (A1) the arterial blood is one well-mixed compartment whose volume $k_b{m}^{kg}$ does not change over the gap; (A2) Hb is constant over the gap and the dissolved term of Eq. (11) is neglected, so $\mathrm{d}{\mathrm{C}}_{\mathrm{a}}{\mathrm{O}}_2/\mathrm{d}t=(1.34\mathrm{Hb}/100)\mathrm{d}y/\mathrm{d}t$; (A3) the airway signal $\mathrm{OUR}\left(t\right)$, the column O2diffxPower (inlet-minus-expiratory O₂ concentration times the ventilator power integral, in arbitrary units a.u.), is proportional to the oxygen taken from the gas, ${\dot{V}}_{{\mathrm{O}}_2}^{air}=c_{OUR}\mathrm{OUR}\left(t\right)$ with $c_{OUR}$ in mL O₂ h⁻¹ a.u.⁻¹, the power integral standing for the gas flow; (A4) the signal $\mathrm{CPR}\left(t\right)$, the column CO2xPower (a.u.), is proportional to the CO₂ production, ${\dot{V}}_{{\mathrm{CO}}_2}=c_{CPR}\mathrm{CPR}\left(t\right)$ with $c_{CPR}$ in mL CO₂ h⁻¹ a.u.⁻¹, and RQ is constant over the gap; (A5) these constants are shared by the patients of a training fold; (A6) the hourly samples are missing at random. Substituting (A2)–(A4) into Eq. (12) and collecting the constants into $k_{OUR}=100c_{OUR}/\left(1.34\mathrm{Hb}{k}_b\right)$ and $k_{CPR}=100c_{CPR}/\left(1.34\mathrm{Hb}{k}_b\mathrm{RQ}\right)$, both in % kg h⁻¹ a.u.⁻¹, gives the balance with body mass on the left,

$$m^{kg}{\displaystyle \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}}=k_{OUR}\mathrm{OUR}\left(t\right)-k_{CPR}\mathrm{CPR}\left(t\right),$$

(13)

and dividing both sides by $m^{kg}$,

$${\displaystyle \frac{\mathrm{d}y\left(t\right)}{\mathrm{d}t}}=k_{OUR}{\displaystyle \frac{\mathrm{OUR}\left(t\right)}{m^{kg}}}-k_{CPR}{\displaystyle \frac{\mathrm{CPR}\left(t\right)}{m^{kg}}}.$$

(14)

This is why the oxygen uptake signal enters the code per kilogram of body mass. In the three scenarios of this paper the consumption term of Eq. (14) is not used: its coefficient $k_{CPR}$ is zero there, and only the uptake term is carried from here on.

2.4.3. Integration over the Gap: the Cumulative Oxygen Uptake Rate

Integrating Eq. (14) over the gap $(t_i-D,t_i]$, as Eq. (10) integrates between samples, with $t_i$ the hour of row i and $D=9$ h (config.py, virtual_sensor_gap):

$$y\left(t_i\right)-y(t_i-D)=k_{OUR}{\int }_{t_i-D}^{t_i}{\displaystyle \frac{\mathrm{OUR}\left(s\right)}{m_i^{kg}}}\mathrm{d}s,$$

(15)

where $m_i^{kg}$ is the last known body mass of the patient of row i (column Svoris; for a patient with no recorded mass the mean over the patients). The integral is the cumulative oxygen uptake rate of [31], taken over the gap and per kilogram. The code evaluates it with the rectangle rule from the OUR values recorded at the hours of the lag set $\mathcal{L}$ (h), the lags of Table 1 plus the current hour ($\mathcal{L}=\{0,1,2,4,5,7,8,10,11\}$ in X0 and XP, $\{0,1,2,4,7\}$ in v4): the available values are summed, the sum is divided by their number $|{\mathcal{L}}_i|$ (${\mathcal{L}}_i\subseteq \mathcal{L}$, the hours at which a value exists for row i), and the average is multiplied by the width D of the integration interval and divided by the body mass:

$${\mathrm{cOUR}}_i={\displaystyle \frac{D}{m_i^{kg}}}·{\displaystyle \frac{1}{|{\mathcal{L}}_i|}}\sum _{\ell \in {\mathcal{L}}_i}\mathrm{OUR}(t_i-\ell )\approx {\int }_{t_i-D}^{t_i}{\displaystyle \frac{\mathrm{OUR}\left(s\right)}{m_i^{kg}}}\mathrm{d}s\left[\mathsf{XP}\right],$$

(16)

in a.u. h kg⁻¹ (column integral_O2diffxPower, Listings 1 and 2; the mean over the available columns is mean(axis=1, skipna=True)). Under (A6) the sample mean is unbiased for the time average of $\mathrm{OUR}/m^{kg}$ over the gap, so Eq. (16) is unbiased for the integral [41]; the plain sum ${\sum }_{\ell }\mathrm{OUR}(t_i-\ell )·1\mathrm{h}$ would equal $|{\mathcal{L}}_i|/D$ times it and shrink with every missing hour. If $\mathrm{OUR}/m^{kg}$ is Lipschitz with constant $\Lambda$ (a.u. kg⁻¹ h⁻¹), the error of the rectangle rule on an hourly grid is at most $\frac{1}{2}\Lambda D·1\mathrm{h}$ [42]. The factor is the gap D and not the number of columns; the two coincide in X0 and XP when all nine values are present ($|{\mathcal{L}}_i|=D=9$) but not in v4, and the factor D keeps the units when lags are pruned. Only the row’s own signals up to $t_i$ are used.

2.4.4. The Last Available Measurement

The measurement of y that the sensor may use at hour $t_i$ is the most recent one taken at or before $t_i-D$:

$${\tau }_i=max\{s\le t_i-D:y\left(s\right)\mathrm{observed}\mathrm{for}\mathrm{patient}p\left(i\right)\},D_i=t_i-{\tau }_i\ge D,$$

(17)

where $D_i$ denotes the actual elapsed time (in hours). The term $y\left({\tau }_i\right)$ corresponds to the nine-hour lagged oxygen saturation (SaO2m_lag_9h), representing the most recent measurement recorded at or before $t_i-D$. During the initial hours of an ICU admission before a nine-hour history is established ($t_i<D$), this value defaults to the patient’s first recorded measurement. Substituting $y\left({\tau }_i\right)$ for $y(t_i-D)$ and substituting Eq. (16) into Eq. (15) yields the integrated oxygen mass balance,

$$y_i\approx y\left({\tau }_i\right)+k_{OUR}{\mathrm{cOUR}}_i,$$

(18)

whose two terms the scenarios use as listed in Table 5. In X0 and XP the initial prediction of the boosted trees is the last available measurement itself,

$${\widehat{y}}_{0,i}=y\left({\tau }_i\right)(\mathsf{X}\mathsf{0},\mathsf{XP}),$$

(19)

so the trees learn corrections to the nine-hour-old measurement rather than to a constant.

Table 5. How each scenario uses the terms of the integrated oxygen balance, Eq. (18).

Term	X0	XP	v4
$y\left({\tau }_i\right)$	initial prediction, Eq. (19)	initial prediction, Eq. (19)	recalibrated toward the patient’s history, Eq. (20)
$k_{OUR}{\mathrm{cOUR}}_i$	not used	${\beta }_1{\mathrm{cOUR}}_i$ in Eq. (37), ${\beta }_1=0.00139$	${\beta }_1{\mathrm{cOUR}}_i$ in Eq. (37), ${\beta }_1=0.00124$

Table 5. How each scenario uses the terms of the integrated oxygen balance, Eq. (18).

Term	X0	XP	v4
$y\left({\tau }_i\right)$	initial prediction, Eq. (19)	initial prediction, Eq. (19)	recalibrated toward the patient’s history, Eq. (20)
$k_{OUR}{\mathrm{cOUR}}_i$	not used	${\beta }_1{\mathrm{cOUR}}_i$ in Eq. (37), ${\beta }_1=0.00139$	${\beta }_1{\mathrm{cOUR}}_i$ in Eq. (37), ${\beta }_1=0.00124$

${\beta }_1$ is found by tuning, not by least squares, and in X0 and XP the trees also receive the OUR signal columns, so ${\beta }_1$ is not an estimate of $k_{OUR}$; both multiply the cumulative oxygen uptake per kilogram over the gap.

2.4.5. Scenario v4: the Patient-Level Anchor

Let ${\mathcal{H}}_i=\{s\le t_i-D:y\left(s\right)\mathrm{observed}\mathrm{for}p\left(i\right)\}$ be the times of all measurements of the patient available at $t_i-D$; the same information as $y\left({\tau }_i\right)$ uses, so no new leakage enters. In v4 the initial prediction recalibrates the last available measurement toward the patient’s own past mean:

$${\widehat{\mathbf{y}}}_{\mathbf{0},\mathbf{i}}={\mathbf{c}}_{\mathbf{0}}+{\mathbf{c}}_{\mathbf{1}}\mathbf{y}\left({\tau }_{\mathbf{i}}\right)+{\mathbf{c}}_{\mathbf{2}}{\overline{\mathbf{y}}}_{{\mathcal{H}}_{\mathbf{i}}},{\overline{y}}_{{\mathcal{H}}_i}={\displaystyle \frac{1}{|{\mathcal{H}}_i|}}\sum _{s\in {\mathcal{H}}_i}y\left(s\right)\left[\mathsf{v}\mathsf{4}\right]$$

(20)

with $c_0$ in % and $c_1$, $c_2$ dimensionless, and ${\overline{y}}_{{\mathcal{H}}_i}:=y\left({\tau }_i\right)$ when ${\mathcal{H}}_i$ is empty (the first D hours of a stay; utils_prediction.py). ${\overline{y}}_{{\mathcal{H}}_i}$ is the column runmean_SaO2m (data_preparation.py, compute_target_runmean). The three coefficients are the weighted least-squares solution on the training rows $\mathcal{T}$ of the fold, with the weights $w_i$ of Eq. (2):

$$({\mathbf{c}}_{\mathbf{0}},{\mathbf{c}}_{\mathbf{1}},{\mathbf{c}}_{\mathbf{2}})=argmin\sum _{\mathbf{i}\in \mathcal{T}}{\mathbf{w}}_{\mathbf{i}}{\left({\mathbf{y}}_{\mathbf{i}}-{\widehat{\mathbf{y}}}_{\mathbf{0},\mathbf{i}}\right)}^{\mathbf{2}}\left[\mathsf{v}\mathsf{4}\right]$$

(21)

with ${\widehat{y}}_{0,i}$ defined by Eq. (20). In the implementation, the three anchor coefficients ($c_0$, $c_1$, $c_2$, representing the intercept, persistence weight, and historical patient mean) are solved in Python and transferred via the base_margin field of the training and validation DMatrix objects. In the C++ core, this directly populates base_margin_ (${\widehat{y}}_{0,i}$), which acts as the initial reference score for each row in the tree recursion (Section 2.5) and the physiological anchor of the final ensemble prediction in Eq. (37).

Each cross-validation fold estimates these coefficients strictly using its training cohort and applies them out-of-fold to the held-out validation patients. The resulting fold-averaged coefficients are $c_0=3.368$, $c_1=0.2472$, and $c_2=0.7170$ (Section 3). Because $c_1+c_2=0.964$, Eq. (20) represents, to within four percent, a weighted combination between the last measured blood gas ($26\%$) and the patient’s long-term historical baseline ($74\%$). The theoretical rationale for assigning a weight strictly below unity to $y\left({\tau }_i\right)$ is regression attenuation [43,44]: framing the observed measurement as $y\left({\tau }_i\right)={\mu }_i+{\xi }_i$, where ${\mu }_i$ is the true underlying physiological state and ${\xi }_i$ represents independent measurement noise, the least-squares shrinkage slope on $y\left({\tau }_i\right)$ is

$$c_1^{*}={\displaystyle \frac{\mathrm{Var}\left(\mu \right)}{\mathrm{Var}\left(\mu \right)+\mathrm{Var}\left(\xi \right)}}<1,$$

(22)

and ${\overline{y}}_{{\mathcal{H}}_i}$ carries the same level with less noise, so Eq. (21) moves weight from $y\left({\tau }_i\right)$ to ${\overline{y}}_{{\mathcal{H}}_i}$, a shrinkage toward the patient’s own level [45].

2.5. The Boosted Tree Ensemble Shared by the Three Scenarios

The three scenarios share one boosted ensemble; what differs between them enters only through the initial prediction ${\widehat{y}}_{0,i}$ and the final prediction of Section 2.6. This subsection therefore describes the ensemble once. What is ours in this subsection. The ensemble is XGBoost 2.1.4 unchanged in all three scenarios; Section 2.5.1 states the four relations of its algorithm that our code touches and the three places where it does, without reproducing its derivations. Our contribution enters the ensemble at one place: the initial prediction ${\widehat{y}}_{0,i}$ from which the recursion of Eq. (31) starts, which is $y\left({\tau }_i\right)$ in X0 and XP, Eq. (19), and the patient-level anchor of Eq. (20) inv4. The final prediction of XP and v4, Eq. (37), is formed after the last tree and is not part of the recursion. Units: $y_i$, ${\widehat{y}}_{b,i}$, $f_b$, the leaf scores $s_{b,j}$ and the gradients $g_{b,i}$ are in %; $w_i$, the Hessians $h_{b,i}$ and their sums, $\lambda$ and $\eta$ are dimensionless; $\alpha$ is in %, the unit of the gradient sums; $\gamma$ and the split gain are in ${\%}^2$; the inputs $x_{i,m}$ are Z-scores. The library we patched and ran is XGBoost release 2.1.4; the excerpts of original XGBoost code in this subsection are copied, with their line numbers, from the release 2.1.4 files of the library we patched and ran, and are typeset in grey; the lines of our own additions are cited against the same sources (Section 2.4 and Section 2.6).

Boosting algorithms, such as the foundational AdaBoost [46], build ensemble models sequentially, where each new tree is specifically trained to correct the errors made by its predecessors [47]. XGBoost [47] improves upon traditional boosting by introducing a sparsity-aware split-finding algorithm for sparse data and a regularized learning objective to prevent overfitting, an approach since validated across many tabular clinical-prediction tasks [26,29]. Throughout this subsection the symbol ≡ denotes a definition. To keep the presentation uniform, every displayed equation is stated first and its newly introduced symbols are defined immediately afterwards in a where clause. The only exceptions are the few symbols that need a sentence of conceptual set-up (such as the routing $q_b$ and the leaf-score vector ${\mathbf{s}}_b$ below): these are introduced in the prose that motivates their equation and are then not repeated. In no case is a symbol defined twice. The fixed global counts are written as n with a subscript naming what is counted: $i\in \{1,\dots ,n_i\}$ indexes the $n_i$ training samples, $m\in \{1,\dots ,n_m\}$ indexes the $n_m$ features, and $b\in \{1,\dots ,n_b\}$ indexes the $n_b$ boosting rounds (equivalently, the $n_b$ trees). The one count that varies from tree to tree, the number of leaves, is written $T_b$, with leaf index $j\in \{1,\dots ,T_b\}$; it is a deliberately separate symbol from the tree count $n_b$ so the two can never be confused.

Input, target, and prediction.

Following the notation of Section 1, each sample i is one mechanically ventilated patient at one hour, described by the feature vector ${\mathbf{x}}_i\in {\mathbb{R}}^{n_m}$ of Table 1 and Table 2, where $n_m$ is the number of features. The value of $n_m$ is a property of the scenario (Section 2.2, Table 1Table 2): in X0 and XP the trees receive $n_m=55$ columns, the six real-time signals at hour t, their 44 lagged copies, and five target-derived or static descriptors; in v4 they receive $n_m=19$, five signals at t with two lags each and four descriptors. Grouping the columns by block, the input vector is

$${\mathbf{x}}_i={(\underset{\begin{array}{c}{n}_s\mathrm{signals}\\ \mathrm{at}\mathrm{hour}t\end{array}}{\underbrace{x_{i,1},\dots ,x_{i,n_s}}},\underset{\begin{array}{c}\mathrm{the}\mathrm{same}\mathrm{signals}\\ \mathrm{at}\mathrm{their}\mathrm{lag}\mathrm{hours}\end{array}}{\underbrace{x_{i,n_s+1},\dots ,x_{i,n_s+n_{\ell }}}},\underset{\begin{array}{c}\mathrm{target}\mathrm{history}\mathrm{and}\\ \mathrm{static}\mathrm{descriptors}\end{array}}{\underbrace{x_{i,n_s+n_{\ell }+1},\dots ,x_{i,n_m}}})}^{\top }\in {\mathbb{R}}^{n_m},$$

(23)

with $(n_s,n_{\ell },n_m)=(6,44,55)$ in X0 and XP and $(5,10,19)$ in v4, where each component $x_{i,m}\in \mathbb{R}$ (for $i=1,\dots ,n_i$ and $m=1,\dots ,n_m$) is a single scalar feature and the superscript ⊤ denotes transposition, so ${\mathbf{x}}_i$ is a column vector. The ground truth $y_i\in [0,100]$ is the co-oximetry oxygen saturation in percent at hour t, a single scalar; to make every symbol concrete enough to check by hand, a small running example throughout this subsection uses five hypothetical hourly samples with target vector $\mathbf{y}={(88,90,91,93,95)}^{\top }$. These five are a didactic device for hand-calculation and do not redefine the dataset size, which in our actual run is $n_i=4017$ in every scenario. Finally, the model’s prediction for sample i is a single scalar ${\widehat{y}}_i\in \mathbb{R}$: its estimate of the true target $y_i$, obtained by mapping the feature vector ${\mathbf{x}}_i$ through the ensemble. Input ${\mathbf{x}}_i$, target $y_i$, and prediction ${\widehat{y}}_i$ are the three quantities this paragraph names; the next paragraph constructs the map ${\mathbf{x}}_i\mapsto {\widehat{y}}_i$ explicitly.

From the input to a prediction.

Having fixed what a single input ${\mathbf{x}}_i$ from Eq. (23) is, we now build the function that turns it into the estimate ${\widehat{y}}_i$. Chen and Guestrin [47] construct this function as a sum of $n_b$ regression trees, each tree adding one correction to the running estimate:

$${\widehat{y}}_{n_b,i}=\varphi \left({\mathbf{x}}_i\right)={\widehat{y}}_{0,i}+\eta {\displaystyle \sum _{b=1}^{n_b}}{f}_b\left({\mathbf{x}}_i\right),f_b\in \mathcal{F},$$

(24)

where ${\widehat{y}}_{n_b,i}$ is the raw boosted estimate of $S_a{O}_{2}m$ for sample i after the last round, $\varphi$ denotes the full ensemble, ${\widehat{y}}_{0,i}$ is the per-row initial prediction of Section 2.4 (Eq. (19) in X0 and XP, Eq. (20) inv4), $\eta$ is the learning rate defined at Eq. (31), $n_b$ is the number of boosting rounds (capped at 2000 and chosen per fold by early stopping: on average $792.4$ trees in X0 and XP, $83.2$ in v4), and $f_b$ is the b-th regression tree, drawn from the tree space $\mathcal{F}$ defined next. In X0 the final prediction is ${\tilde{y}}_i={\widehat{y}}_{n_b,i}$; in XP and v4 it is the blend of Eq. (37).

What a single tree is.

A regression tree is a finite set of yes/no questions arranged as a binary tree, each question of the form "is feature m of the sample below the threshold $\vartheta$?", together with one number stored at every leaf. To evaluate the tree on a sample, one starts at the root, answers the question, moves to the left child on "yes" and to the right child on "no", and repeats until a leaf is reached; the tree’s output is the number stored at that leaf. Two objects therefore describe tree b completely. The first is the routing $q_b$, which is the branch structure written as a function: $q_b\left({\mathbf{x}}_i\right)\in \{1,\dots ,T_b\}$ is the index of the leaf that sample i reaches, and $q_b:{\mathbb{R}}^{n_m}\to \{1,\dots ,T_b\}$ is defined for any feature vector, including ventilator readings never seen during training. The second is the leaf-score vector ${\mathbf{s}}_b={(s_{b,1},\dots ,s_{b,T_b})}^{\top }\in {\mathbb{R}}^{T_b}$, one number per leaf. With $j=q_b\left({\mathbf{x}}_i\right)$ the leaf reached, the output of the tree is

$$f_b\left({\mathbf{x}}_i\right)=s_{b,j},j=q_b\left({\mathbf{x}}_i\right).$$

(25)

The same routing fixes which training samples share a leaf, and the objective derived below sums over those per-leaf groups; $q_b$ reappears in that role in Eq. (33).

$$\mathcal{F}=\left\{f_b:\mathbf{x}\mapsto s_{b,q_b\left(\mathbf{x}\right)}|q_b:{\mathbb{R}}^{n_m}\to \{1,\dots ,T_b\},{\mathbf{s}}_b\in {\mathbb{R}}^{T_b}\right\}.$$

(26)

The same routing fixes which training samples share a leaf, and the objective derived below sums over those per-leaf groups; $q_b$ reappears in that role in Eq. (33).

A concrete tree.

The purpose of this example is deliberately narrow: to exhibit the abstract pair $(q_b,{\mathbf{s}}_b)$ of the previous paragraph on a tree small enough to check by hand. It draws on two results that are only derived later, how a leaf score is obtained in closed form (Eq. (35), Section 2.5.1) and how a split is chosen (Section 2.5.1), and states them here without proof so that the objects can be seen in action before the machinery that produces them is built. This illustrate first, derive later ordering is common in expository treatments; a reader may take the two results on trust until those sections. Concretely, take $f_b$ to be a depth-1 tree whose root splits on the inspired oxygen concentration $O_2\mathrm{conc}$ at hour t, the first feature of the input, i.e. component $x_{i,1}$ of ${\mathbf{x}}_i$ in Eq. (23), producing $T_b=2$ leaves:

$$q_b\left({\mathbf{x}}_i\right)=\left\{\begin{array}{cc}1,\hfill & x_{i,1}\le 40\%,\hfill \\ 2,\hfill & x_{i,1}>40\%,\hfill \end{array}\right.{\mathbf{s}}_b={(s_{b,1},s_{b,2})}^{\top }.$$

(27)

The indices here are worth spelling out, because three different roles meet in this one equation. The value being tested, $x_{i,1}$, is a feature of samplei, so it carries the sample index i and the feature index 1; it does not carry the round index b, since the input does not depend on which tree reads it. What does belong to tree b is the split rule, the choice to branch on component 1 and the threshold $40\%$ are parameters of this particular $f_b$, which is exactly why the routing $q_b$ carries the subscript b. And the leaf index j does not appear on the input side at all: it is the output of the routing, $j=q_b\left({\mathbf{x}}_i\right)\in \{1,2\}$, fed into Eq. (25) to read the score. Concretely, sample i with $x_{i,1}=35\%$ gives $q_b\left({\mathbf{x}}_i\right)=1$, so by (25) the tree outputs $f_b\left({\mathbf{x}}_i\right)=s_{b,1}$; sample i with $x_{i,1}=55\%$ gives $q_b\left({\mathbf{x}}_i\right)=2$ and output $s_{b,2}$. Neither part of this tree is chosen by hand. XGBoost learns the leaf scores $s_{b,1},s_{b,2}$ from the data through the closed form () below, and it learns the split structure, which feature each node tests, and at what threshold, through the greedy search of Section 2.5.1, which scans the candidate thresholds of every feature and keeps the one that maximizes the split gain (). The number $40\%$ is a didactic placeholder chosen to keep the arithmetic legible; it is not a fitted value from our run. The quantity it stands for is, at training time, the C++ variable split_pt that EnumerateSplit selects for the node (the third argument of best.Update in the listing of Section 2.5.1); once training finishes, that value is stored as the node’s split condition inside the fitted booster and can be read back from the Python model with model.get_booster().trees_to_dataframe() (the Split column) or get_dump(). All split thresholds and leaf weights are determined automatically by the tree-growing algorithm via gradient-based gain maximization during training, rather than by manually specified decision rules.

In the source code, a tree’s prediction and the objective registered as reg:squarederror are realized by the LinearSquareLoss struct; the routing and leaf storage are handled by the tree builder discussed in Section 2.5.1.1

How the trees are chosen: the training objective.

The ensemble (24) is only as good as the trees placed in it, so the trees are selected to minimize a training objective that balances how well the model fits the data against how complex it becomes. The objective is assembled from two ingredients. The first is the data-fit term: inserting the ensemble prediction ${\widehat{y}}_i=\varphi \left({\mathbf{x}}_i\right)$ of Eq. (24) into the squared-error loss $\ell (y_i,{\widehat{y}}_i)={\displaystyle \frac{1}{2}}{(y_i-{\widehat{y}}_i)}^2$ (the factor ${\displaystyle \frac{1}{2}}$ is a convention that makes the gradient in Eq. (32) clean and does not change the minimizer) and summing over all $n_i$ samples gives the total prediction error ${\sum }_{i=1}^{n_i}\ell (y_i,{\widehat{y}}_i)$. The second is a complexity penalty charged to each individual tree $f_b$,

$$\Omega \left(f_b\right)=\gamma T_b+{\displaystyle \frac{1}{2}}\lambda {\displaystyle \sum _{j=1}^{T_b}}{s}_{b,j}^2,$$

(28)

which grows with that tree’s number of leaves $T_b$ (penalty rate $\gamma \ge 0$, the gamma parameter, left at its default 0 in our runs) and with the squared magnitude of its leaf scores $s_{b,j}\in \mathbb{R}$ (rate $\lambda \ge 0$, the $L_2$ strength reg_lambda $=5$); the index $j\in \{1,\dots ,T_b\}$ runs over the leaves of that one tree, and $s_{b,j}$ is the score at leaf j of tree b (the j-th component of ${\mathbf{s}}_b$). Now, because the ensemble of Eq. (24) is the sum of all $n_b$ trees, ${\widehat{y}}_i={\sum }_{b=1}^{n_b}{f}_b\left({\mathbf{x}}_i\right)$, the model’s total complexity is the sum of these per-tree penalties — one $\Omega \left(f_b\right)$ term for each of the $n_b$ trees — and this is precisely where the second summation below comes from. Adding the two ingredients yields the regularized objective of Chen and Guestrin [47]:

$$\mathcal{L}(f_1,\dots ,f_{n_b})=\underset{\mathrm{total}\mathrm{prediction}\mathrm{error}(\mathrm{sum}\mathrm{over}{n}_i\mathrm{samples})}{\underbrace{\sum _{i=1}^{n_i}\ell \left(y_i,{\widehat{y}}_i\right)}}+\underset{\mathrm{summed}\mathrm{complexity}(\mathrm{sum}\mathrm{over}{n}_b\mathrm{trees})}{\underbrace{\sum _{b=1}^{n_b}\Omega \left(f_b\right)}}.$$

(29)

The two summations run over different index sets and therefore cannot be merged: the first is indexed by the sample $i=1,\dots ,n_i$, the second by the tree $b=1,\dots ,n_b$. Nor can the complexity be pulled outside its own sum, because $\Omega \left(f_b\right)$ genuinely differs from tree to tree — the leaf count $T_b$ and the scores $s_{b,j}$ both depend on b, while only the rates $\gamma ,\lambda$ are shared — which is exactly why each carries the tree subscript b.

Making the loss cost-sensitive.

Finally, to make the predictions more sensitive to clinically critical low saturations, we insert the cost-sensitive sample weight $w_i$ of Eq. (2) into each loss term of the objective (29). This single substitution yields the weighted objective actually optimized in this study:

$$\mathcal{L}(f_1,\dots ,f_{n_b})={\displaystyle \sum _{i=1}^{n_i}}{w}_i·\ell \left(y_i,{\widehat{y}}_i\right)+{\displaystyle \sum _{b=1}^{n_b}}\Omega \left(f_b\right),$$

(30)

where $w_i=(1-k_w)/\left(\right|y_i{|}^{\nu }+\epsilon )+k_w$ from Eq. (2); the worked example uses the illustrative floor $k_w=0.15$, exponent $\nu =1.75$ and stabiliser $\epsilon ={10}^{-6}$ (the scenario values are in Table 8Table 9). For the five example targets these weights are all $w_i\approx 0.150$ (the targets are large and similar), whereas a hypothetical desaturation to $y_i=60$ would raise that sample’s weight, up-weighting its contribution to the loss. Equation (30) is written here in the compact summary form used in the literature; the remainder of this subsection derives, step by step, how the trees that minimize it are actually built.

2.5.1. Training the Trees: Original XGBoost, and the Three Places Where Our Code Enters

The trees are trained by the gradient-boosting algorithm of XGBoost 2.1.4 exactly as released. We reproduce none of its derivations: the additive recursion is that of Friedman [48], and the second-order objective, the closed-form leaf score and the split gain are those of Chen and Guestrin [47], to which we refer for their proofs. We state only the four relations of that algorithm that our code touches, in the notation of this paper, so that each of our interventions can be located (Figure 2). Starting from the initial prediction ${\widehat{y}}_{0,i}$, each round adds one tree scaled by the learning rate $\eta$; the tree of round b is fitted to the gradient and Hessian of the weighted squared loss at the previous prediction, summed over the training rows that its routing $q_b$ sends to leaf j; the leaf score and the gain of a split into children L and R follow in closed form with the penalties $\lambda$, $\alpha$ and $\gamma$:

$$\begin{array}{cc}\hfill {\widehat{y}}_{b,i}& ={\widehat{y}}_{b-1,i}+\eta f_b\left({\mathbf{x}}_i\right),b=1,\dots ,n_b,\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill g_{b,i}& =w_i({\widehat{y}}_{b-1,i}-y_i),h_{b,i}=w_i,\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\mathcal{I}}_{b,j}& =\{i:q_b\left({\mathbf{x}}_i\right)=j\},\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill G_{b,j}& =\sum _{i\in {\mathcal{I}}_{b,j}}{g}_{b,i},H_{b,j}=\sum _{i\in {\mathcal{I}}_{b,j}}{h}_{b,i},\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill s_{b,j}^{*}& =-{\displaystyle \frac{T_{\alpha }\left(G_{b,j}\right)}{H_{b,j}+\lambda }},T_{\alpha }\left(G\right)=\mathrm{sign}\left(G\right)max\left(\right|G|-\alpha ,0),\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill \mathrm{Gain}& ={\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{T_{\alpha }{\left(G_{b,L}\right)}^2}{H_{b,L}+\lambda }}+{\displaystyle \frac{T_{\alpha }{\left(G_{b,R}\right)}^2}{H_{b,R}+\lambda }}-{\displaystyle \frac{T_{\alpha }{(G_{b,L}+G_{b,R})}^2}{H_{b,L}+H_{b,R}+\lambda }}\right]-\gamma ,\hfill \end{array}$$

(36)

with ${\mathcal{I}}_{b,j}$ the set of training rows in leaf j, $s_{b,j}^{*}$ in %, G and $\alpha$ in %, H and $\lambda$ dimensionless, and the gain and $\gamma$ in ${\%}^2$; a node is split while the gain is positive, both children keep a summed Hessian of at least min_child_weight and the depth is below max_depth, and otherwise it becomes a leaf with the score of Eq. (35). In the library, Eq. (32) is src/objective/regression_obj.cu, Eq. (35) is CalcWeight in src/tree/param.h, Eq. (36) is CalcGain applied inside EnumerateSplit in src/tree/hist/evaluate_splits.h, and the scaling by $\eta$ of Eq. (31). Our code enters this algorithm at three places, all of them before the first tree is grown, and modifies nothing inside it. First, the starting value ${\widehat{y}}_{0,i}$ of Eq. (31) is our anchor, the last available measurement in X0 and XP, Eq. (19), and the patient-level anchor of Eq. (20) inv4, passed as the base_margin of both matrices (utils_prediction.py) and copied by the library into the prediction buffer of every row before the first tree (Predictor::InitOutPredictions, src/predictor/predictor.cc); consequently every gradient of Eq. (32) is evaluated relative to the anchor, and every leaf score of Eq. (35) is a correction to it rather than to a constant. Second, the weights $w_i$ of Eq. (2), our relative-entropy criterion of Section 2.3, are the weight of the training matrix and multiply every gradient and Hessian (regression_obj.cu). Third, the hyperparameters $\eta$, $\lambda$, $\alpha$, $\gamma$, max_depth, min_child_weight, the two subsampling ratios and the early-stopping patience are set per scenario (Table 8 and Table 9). The library’s evaluation of the validation rows after each round (src/learner.cc) sees the unblended ${\widehat{y}}_{b,i}$ and chooses the number of trees $n_b$ by early stopping (Section 2.8.1): on average $792.4$ trees per fold with $12.0$ leaves each in X0 and XP, and $83.2$ trees with $5.4$ leaves in v4 (Table 10). After the last tree, the final prediction of Section 2.6 is formed; it is ours and is not part of the training. The reason the anchor changes what the trees learn is visible in Eq. (32): with the raw persistence anchor, the first gradients are the errors of a noisy nine-hour-old value, and the trees spend their capacity on pulling it toward the patient’s level, which appears as a negative variance share of the trees in Section 3.3; with the patient-level anchor that work is done in closed form before boosting, and a far smaller ensemble suffices.

2.6. The Final Prediction

The trees give the raw boosted estimate ${\widehat{y}}_{n_b,i}$ of Eq. (24). After the last tree, our C++ function ApplyAnchorBlend forms the final prediction ${\tilde{y}}_i$ (%) of every row, training and validation rows alike:

$${\tilde{y}}_i={\widehat{y}}_{0,i}+(1-\beta )\left({\widehat{y}}_{n_b,i}-{\widehat{y}}_{0,i}\right)+{\beta }_1{\mathrm{cOUR}}_i\left[\mathsf{XP}\right]$$

(37)

where the dimensionless $\beta \in [0,1]$ (custom_b) pulls the prediction toward its initial value and ${\beta }_1$ (% kg h⁻¹ a.u.⁻¹; custom_b1) weighs the cumulative oxygen uptake rate of Eq. (16): $(\beta ,{\beta }_1)=(0,0)$ in X0, $(0.1418,0.00139)$ in XP (config.py) and $(0.3677,0.00124)$ in v4 (locked configuration file). Substituting Eq. (24), which gives ${\widehat{y}}_{n_b,i}-{\widehat{y}}_{0,i}=\eta {\sum }_b{f}_b\left({\mathbf{x}}_i\right)$, together with Eqs. (20) and (16), into Eq. (37) writes the prediction of v4 in its original quantities:

$${\tilde{y}}_i=\underset{{\widehat{y}}_{0,i}}{\underbrace{{\mathbf{c}}_{\mathbf{0}}+{\mathbf{c}}_{\mathbf{1}}\mathbf{y}\left({\tau }_{\mathbf{i}}\right)+{\mathbf{c}}_{\mathbf{2}}{\overline{\mathbf{y}}}_{{\mathcal{H}}_{\mathbf{i}}}}}+(1-\beta )\eta {\displaystyle \sum _{b=1}^{n_b}}{f}_b\left({\mathbf{x}}_i\right)+{\beta }_1{\displaystyle \frac{D}{m_i^{kg}}}{\displaystyle \frac{1}{|{\mathcal{L}}_i|}}\sum _{\ell \in {\mathcal{L}}_i}\mathrm{OUR}(t_i-\ell )\left[\mathsf{v}\mathsf{4}\right]$$

(38)

In X0 and XP the first term is $y\left({\tau }_i\right)$, Eq. (19). Table 6 gives, for every term of Eq. (38), the Python variable that holds it, how it reaches src/learner.cc, where it is read and when it is used.

Table 6. The terms of Eq. (38) in the code. File utils_prediction.py unless stated; C++ lines in src/learner.cc.

Term	Python	In C++	Used in
$c_0$, $c_1$, $c_2$	coef[0..2];	not present; they act through ${\widehat{y}}_{0,i}$	through ${\widehat{y}}_{0,i}$
$y\left({\tau }_i\right)$, ${\overline{y}}_{{\mathcal{H}}_i}$	columns SaO2m_lag_9h, runmean_SaO2m	not present	through ${\widehat{y}}_{0,i}$
${\widehat{y}}_{0,i}$	raw_base_train, raw_base_val; base_margin	base_margin_: start of the prediction buffer before the first tree; anchor	training rounds, early stopping, prediction
$\eta {\sum }_b{f}_b\left({\mathbf{x}}_i\right)$	the trees, trained by xgb.train	PredictRaw (training), (early stopping),(prediction)	training rounds, early stopping, prediction
${\mathrm{cOUR}}_i$	raw_int_train, raw_int_val;set_float_info("custom_integral",...),	our field custom_integral_ (declared in include/xgboost/data.h; filled by MetaInfo::SetInfo, src/data/data.cc; C API src/c_api/c_api.cc)	prediction
$\beta$, ${\beta }_1$	params["custom_b"], params["custom_b1"]; xgb.train(params=params, ...)	tparam_.custom_b, tparam_.custom_b1	prediction

Table 6. The terms of Eq. (38) in the code. File utils_prediction.py unless stated; C++ lines in src/learner.cc.

Term	Python	In C++	Used in
$c_0$, $c_1$, $c_2$	coef[0..2];	not present; they act through ${\widehat{y}}_{0,i}$	through ${\widehat{y}}_{0,i}$
$y\left({\tau }_i\right)$, ${\overline{y}}_{{\mathcal{H}}_i}$	columns SaO2m_lag_9h, runmean_SaO2m	not present	through ${\widehat{y}}_{0,i}$
${\widehat{y}}_{0,i}$	raw_base_train, raw_base_val; base_margin	base_margin_: start of the prediction buffer before the first tree; anchor	training rounds, early stopping, prediction
$\eta {\sum }_b{f}_b\left({\mathbf{x}}_i\right)$	the trees, trained by xgb.train	PredictRaw (training), (early stopping),(prediction)	training rounds, early stopping, prediction
${\mathrm{cOUR}}_i$	raw_int_train, raw_int_val;set_float_info("custom_integral",...),	our field custom_integral_ (declared in include/xgboost/data.h; filled by MetaInfo::SetInfo, src/data/data.cc; C API src/c_api/c_api.cc)	prediction
$\beta$, ${\beta }_1$	params["custom_b"], params["custom_b1"]; xgb.train(params=params, ...)	tparam_.custom_b, tparam_.custom_b1	prediction

The predicted saturation is formed in Predict: PredictRaw evaluates the trees from ${\widehat{y}}_{0,i}$ and gives ${\widehat{y}}_{n_b,i}$, and ApplyAnchorBlend overwrites it in place with ${\tilde{y}}_i$ . Python receives ${\tilde{y}}_i$ from bst.predict(dval) for the validation rows, which are the reported out-of-fold predictions, and from bst.predict(dtrain) for the training rows. $\beta$ and ${\beta }_1$ therefore act on the predictions of training and validation rows alike, but not on the training loss or the early-stopping criterion: UpdateOneIter and EvalOneIter never call ApplyAnchorBlend. For $\beta \in [0,1]$ and ${\beta }_1=0$, Eq. (37) is the convex combination $(1-\beta ){\widehat{y}}_{n_b,i}+\beta {\widehat{y}}_{0,i}$, hence by the triangle inequality

$$|{\tilde{y}}_i-y_i|\le (1-\beta )|{\widehat{y}}_{n_b,i}-y_i|+\beta |{\widehat{y}}_{0,i}-y_i|\le max\left(|{\widehat{y}}_{n_b,i}-y_i|,|{\widehat{y}}_{0,i}-y_i|\right),$$

(39)

which is why custom_b is declared on $[0,1]$ and why Eq. (37) is applied once, after the last tree: applied once per boosting round it would compound to ${(1-\beta )}^{n_b}$ and Eq. (39) would not hold.

2.7. Anatomy of a Prediction: How Much Is Physics and How Much Is Trees

Eq. (38) splits every final prediction exactly into three parts, written next to every out-of-fold row by the pipeline (anchorpart, treepart, integralpart; utilsprediction.py):

$${\tilde{y}}_i=A_i+B_i+C_i,A_i={\widehat{y}}_{0,i},B_i=(1-\beta )\eta {\displaystyle \sum _{b=1}^{n_b}}{f}_b\left({\mathbf{x}}_i\right),C_i={\beta }_1{\mathrm{cOUR}}_i,$$

(40)

all in %. For each part $X\in \{A,B,C\}$ the variance share and the work share are

$$\mathrm{variance}\mathrm{share}\left(X\right)={\displaystyle \frac{\mathrm{Cov}(X,\tilde{y})}{\mathrm{Var}\left(\tilde{y}\right)}},\mathrm{work}\mathrm{share}\left(X\right)={\displaystyle \frac{\overline{|X-X^0|}}{{\sum }_{X^{\prime }}\overline{|X^{\prime }-X^{\prime 0}|}}},$$

(41)

with $X^0=y\left({\tau }_i\right)$ for A and $X^0=0$ for B and C, and the error credit of the anchor is

$$\mathrm{error}\mathrm{credit}={\displaystyle \frac{\mathrm{MAPE}\left(y\left(\tau \right)\right)-\mathrm{MAPE}\left({\widehat{y}}_0\right)}{\mathrm{MAPE}\left(y\left(\tau \right)\right)-\mathrm{MAPE}\left(\tilde{y}\right)}},$$

(42)

with MAPE as in Eq. (43). By bilinearity of the covariance, substituting Eq. (40) into Eq. (41) gives $\mathrm{Cov}(A,\tilde{y})+\mathrm{Cov}(B,\tilde{y})+\mathrm{Cov}(C,\tilde{y})=\mathrm{Var}\left(\tilde{y}\right)$, so the three variance shares sum to 100%; a negative share means that the part moves against the prediction. The work share says how far each part moves the prediction away from the last available measurement, and the error credit which fraction of the improvement over it the anchor earns before any tree is grown.

2.8. Training and Evaluation Strategy

To avoid overly optimistic performance estimates caused by temporal correlation within the same patient trajectory, we implemented a 5-fold Grouped Cross-Validation strategy. In this configuration, data is partitioned such that all samples of a patient are kept together within a single fold. This ensures that no data from the same patient used during the training phase is ever present in the test set, effectively preventing data leakage. The patient-level fold assignment is shuffled with a fixed random seed (random_state $=42$), which also seeds the XGBoost estimator, so that the entire pipeline is reproducible; every scenario uses the same folds.

Target-history leakage and the strict gap.

Because two of the inputs are the target’s own past values ($S_a{O}_{2}m$ at $t-9$ and $t-10$), the validity of the reported error depends critically on which target lags are admitted as features. The blood-gas measurement that defines $S_a{O}_{2}m$ is itself delayed by roughly nine hours, so at deployment the freshest genuinely available value of the target is its 9-hour lag; any value of $S_a{O}_{2}m$ from the intervening eight hours ($t-1,\dots ,t-8$), and any feature derived from future samples, would simply not exist at prediction time. If such a near-past or future value of the target were allowed into the feature set, the model would in effect be shown a lightly time-shifted copy of the answer: the reported MAPE would collapse to an optimistically small figure, and the “virtual sensor” would be quietly reading the very signal it is meant to reconstruct rather than inferring it from ventilator and gas measurements. The pipeline constructs only target lags of nine hours or more. It does not, however, discard the rows that lack a nine-hour lag: in the first nine hours of a stay it fills the lag with the patient’s first measurement (fallback_to_first=True), and SaO2m_first and the patient mean carry the same value, so these 287 rows ($7.1\%$) see a value fresher than $t-9$, and 33 of them the current target itself. By the rule above these rows are invalid; all results of Section 3 are computed without them, so that the estimator at hour t sees the real-time gas and ventilator channels plus target history no fresher than $t-9$. It is important to stress that the patient-grouped cross-validation above does not on its own guard against this particular failure: grouping by patient eliminates between-patient leakage, but a near-past or future copy of a patient’s own target leaks within that patient’s trajectory, which grouped folds cannot detect. The headline accuracy must therefore be read as conditional on the strict-gap construction; relaxing it, admitting sub-gap target lags, or back-filling missing values from later rows, would lower the apparent error without any gain in genuine predictive power, and any such configuration should be treated as invalid for the nine-hour-ahead claim.

The primary evaluation metric utilized in this study is the Mean Absolute Percentage Error (MAPE), as it provides a relative measure of error that is clinically informative across varying oxygenation levels. The MAPE is defined as:

$$\mathrm{MAPE}={\displaystyle \frac{100\%}{n_i}}{\displaystyle \sum _{i=1}^{n_i}}\left|{\displaystyle \frac{y_i-{\widehat{y}}_i}{y_i}}\right|$$

(43)

Equation (43) is the pooled, row-weighted MAPE. Three patients (18, 3 and 21) hold about 85% of the squared error of every scenario, and five patients have hours with two readings, so the pooled MAPE mostly reports those patients. We therefore report, for every scenario, the per-patient MAPE

$${\mathrm{MAPE}}_p={\displaystyle \frac{100\%}{n_{i,p}}}\sum _{i:p\left(i\right)=p}\left|{\displaystyle \frac{y_i-{\tilde{y}}_i}{y_i}}\right|$$

(44)

and three summaries of it over the $n_p=31$ patients: the macro MAPE, its mean; the trimmed MAPE, its mean without the three worst patients; and its median. The macro MAPE is the objective that the pruning and the tuning of v4 used. Tracking within a patient is measured by

$$R_{\mathrm{within},p}^2=1-{\displaystyle \frac{\frac{1}{n_{i,p}}{\sum }_{i:p\left(i\right)=p}{(y_i-{\tilde{y}}_i)}^2}{{\mathrm{Var}}_p\left(y\right)}},$$

(45)

which is positive only when the model beats the patient’s own mean; its median over patients is reported together with the same quantity computed for persistence, the floor a virtual sensor must beat, and with the number of patients in which the model beats persistence. The pooled Pearson correlation r between $\tilde{y}$ and y is reported as in our earlier work, with the caveat that it is lifted by between-patient level differences.

2.8.1. Training Rows and Validation Rows

In each fold the rows of 24 or 25 patients are training rows and the rows of the other 6 or 7 patients are validation rows; the reported predictions are those of the validation rows. Table 7 lists every place where the two are treated differently.

Two consequences follow. The training loss and the early-stopping criterion never see $\beta$ and ${\beta }_1$: the blend of Eq. (37) exists only at prediction time. And the validation rows choose $n_b$ through early stopping, so they are not entirely unseen during training; the folds are identical in every scenario, so the comparison between scenarios is not affected.

Hyperparameters of each scenario.

X0 and XP use the configured operating point of our pipeline (config.py; Table 8); X0 differs from XP only in $\beta ={\beta }_1=0$. The hyperparameters of v4 (Table 9) were selected by a tree-structured Parzen estimator search (Optuna [49]) on the same patient-grouped folds for the 19-input tree with an earlier anchor, and were kept unchanged when the anchor of Eq. (20) was introduced; v4 itself was not re-tuned. Its lag lists (Table 1) were reduced from those of X0 and XP by lag-group ablation on the same folds; the column O2diffxPower is hidden from its trees because it feeds Eq. (16), and Svoris was hidden by a knockout test on the same folds. The configuration is fixed in the locked configuration file (Section 2.8.1.2).

Table 7. What differs between training and validation rows within a fold. File and line refer to the delivered sources.

Step	Training rows	Validation rows	Code
standardisation of the tree inputs	scaler fitted and applied	training scaler applied	data_preparation.py
$y\left({\tau }_i\right)$, ${\overline{y}}_{{\mathcal{H}}_i}$, ${\mathrm{cOUR}}_i$	from the row’s own past	the same formula	data_preparation.py
anchor coefficients, Eq. (21)	fitted	applied, not fitted	utils_prediction.py
sample weights $w_i$	in the loss	none	utils_prediction.py
trees $f_b$	fitted to the gradients	not used for fitting	utils_prediction.py
${\widehat{y}}_{b,i}$ after each round	enters the gradients	enters the early-stopping RMSE, unweighted and without the blend	learner.cc
number of trees $n_b$	—	decided by early stopping	utils_prediction.py
final prediction, Eq. (37)	computed	computed; reported out of fold	learner.cc

Table 7. What differs between training and validation rows within a fold. File and line refer to the delivered sources.

Step	Training rows	Validation rows	Code
standardisation of the tree inputs	scaler fitted and applied	training scaler applied	data_preparation.py
$y\left({\tau }_i\right)$, ${\overline{y}}_{{\mathcal{H}}_i}$, ${\mathrm{cOUR}}_i$	from the row’s own past	the same formula	data_preparation.py
anchor coefficients, Eq. (21)	fitted	applied, not fitted	utils_prediction.py
sample weights $w_i$	in the loss	none	utils_prediction.py
trees $f_b$	fitted to the gradients	not used for fitting	utils_prediction.py
${\widehat{y}}_{b,i}$ after each round	enters the gradients	enters the early-stopping RMSE, unweighted and without the blend	learner.cc
number of trees $n_b$	—	decided by early stopping	utils_prediction.py
final prediction, Eq. (37)	computed	computed; reported out of fold	learner.cc

Table 8. Hyperparameters of X0 and XP (config.py and the TrainConfig block). X0 sets $\beta ={\beta }_1=0$; everything else is identical.

Parameter (name in the code)	Value	Role
learning_rate ($\eta$)	$0.0150$	step applied to each tree’s correction, Eq. (31)
n_estimators ($n_b$ cap)	2000	ceiling on the tree count; early stopping selects $792.4$ on average
early_stopping_rounds	100	rounds without validation improvement before halting
max_depth	4	maximum tree depth
min_child_weight	$0.632$	minimum summed Hessian per leaf (the IsValid guard)
reg_lambda ($\lambda$)	$17.17$	$L_2$ penalty on leaf scores, Eq. (35)
reg_alpha ($\alpha$)	$0.328$	$L_1$ penalty on leaf scores
gamma ($\gamma$)	0	minimum gain to accept a split
subsample	$0.505$	fraction of rows per tree
colsample_bytree	$1.0$	fraction of columns per tree
tree_method	hist	histogram split finding
custom_b ($\beta$)	0 (X0) / $0.1418$ (XP)	anchor blend, Eq. (37)
custom_b1 (${\beta }_1$)	0 (X0) / $0.00139$ (XP)	integral weight, Eq. (37)
custom_power ($\nu$), custom_kp ($k_w$)	$3.656$, $0.148$	cost weights, Eq. (2); inert on this target range (Section 3.5)
random_state	42	fold shuffle and model seed

Table 8. Hyperparameters of X0 and XP (config.py and the TrainConfig block). X0 sets $\beta ={\beta }_1=0$; everything else is identical.

Parameter (name in the code)	Value	Role
learning_rate ($\eta$)	$0.0150$	step applied to each tree’s correction, Eq. (31)
n_estimators ($n_b$ cap)	2000	ceiling on the tree count; early stopping selects $792.4$ on average
early_stopping_rounds	100	rounds without validation improvement before halting
max_depth	4	maximum tree depth
min_child_weight	$0.632$	minimum summed Hessian per leaf (the IsValid guard)
reg_lambda ($\lambda$)	$17.17$	$L_2$ penalty on leaf scores, Eq. (35)
reg_alpha ($\alpha$)	$0.328$	$L_1$ penalty on leaf scores
gamma ($\gamma$)	0	minimum gain to accept a split
subsample	$0.505$	fraction of rows per tree
colsample_bytree	$1.0$	fraction of columns per tree
tree_method	hist	histogram split finding
custom_b ($\beta$)	0 (X0) / $0.1418$ (XP)	anchor blend, Eq. (37)
custom_b1 (${\beta }_1$)	0 (X0) / $0.00139$ (XP)	integral weight, Eq. (37)
custom_power ($\nu$), custom_kp ($k_w$)	$3.656$, $0.148$	cost weights, Eq. (2); inert on this target range (Section 3.5)
random_state	42	fold shuffle and model seed

Table 9. Hyperparameters of v4 (locked configuration file locked_model_v4.json, keys xgb_params, train and runtime), selected by the search described in the text.

Parameter (name in the code)	Value	Role
learning_rate ($\eta$)	$0.0986$	step applied to each tree’s correction
n_estimators ($n_b$ cap)	2000	ceiling; early stopping selects $83.2$ trees on average
early_stopping_rounds	50	patience
max_depth	3	maximum tree depth
min_child_weight	$97.96$	minimum summed Hessian per leaf
reg_lambda ($\lambda$)	$59.61$	$L_2$ penalty on leaf scores
reg_alpha ($\alpha$)	$11.19$	$L_1$ penalty on leaf scores
subsample, colsample_bytree	$0.975$, $0.927$	row and column fractions per tree
custom_b ($\beta$)	$0.3677$	anchor blend, Eq. (37)
custom_b1 (${\beta }_1$)	$0.00124$	integral weight, Eq. (37)
mass_balance=calib_pt	—	the anchor of Eq. (20), fitted per fold
tree_exclude	O2diffxPower, Svoris	columns hidden from the trees
custom_power ($\nu$), custom_kp ($k_w$)	$9.174$, $0.928$	cost weights; inert (Section 3.5)

Table 9. Hyperparameters of v4 (locked configuration file locked_model_v4.json, keys xgb_params, train and runtime), selected by the search described in the text.

Parameter (name in the code)	Value	Role
learning_rate ($\eta$)	$0.0986$	step applied to each tree’s correction
n_estimators ($n_b$ cap)	2000	ceiling; early stopping selects $83.2$ trees on average
early_stopping_rounds	50	patience
max_depth	3	maximum tree depth
min_child_weight	$97.96$	minimum summed Hessian per leaf
reg_lambda ($\lambda$)	$59.61$	$L_2$ penalty on leaf scores
reg_alpha ($\alpha$)	$11.19$	$L_1$ penalty on leaf scores
subsample, colsample_bytree	$0.975$, $0.927$	row and column fractions per tree
custom_b ($\beta$)	$0.3677$	anchor blend, Eq. (37)
custom_b1 (${\beta }_1$)	$0.00124$	integral weight, Eq. (37)
mass_balance=calib_pt	—	the anchor of Eq. (20), fitted per fold
tree_exclude	O2diffxPower, Svoris	columns hidden from the trees
custom_power ($\nu$), custom_kp ($k_w$)	$9.174$, $0.928$	cost weights; inert (Section 3.5)

Reproducibility.

Each scenario is fully specified by one configuration: X0 and XP by config.py with $\beta ={\beta }_1=0$ or with their configured values, and v4 by a locked JSON file that records its lag lists, the columns hidden from the trees, the anchor mode and every hyperparameter; XGB_LOCKED=<file> python train.py re-trains it and reproduces the metrics of Section 3 exactly (same folds, same seed). The three scenarios were evaluated in one run on identical folds, and the out-of-fold prediction of every row was stored together with its three parts of Eq. (40).

3. Results

All results are out-of-fold predictions under the patient-grouped five-fold cross-validation of Section 2.8, with identical folds for the three scenarios. Following the rule of Section 2.8, they are scored on the 3730 rows of 30 patients whose inputs contain no $S_a{O}_{2}m$ value younger than nine hours (Section 3.7); the values on all 4017 rows are listed alongside in Table 14. A patient counts as beating persistence when its $R_{\mathrm{within},p}^2$, Eq. (45), exceeds that of $y\left({\tau }_i\right)$.

3.1. Accuracy and Size of the Three Scenarios

Table 10 reports accuracy and size. The C++ terms of XP lower the macro MAPE by 0.032 percentage points: a paired patient bootstrap of XP − X0 (2000 resamples of the 30 patients) gives $-0.032$ with 95% interval $[-0.065,-0.001]$, while the intervals of the pooled ($-0.010$, $[-0.055,+0.027]$), trimmed ($-0.027$, $[-0.058,+0.001]$) and median per-patient MAPE ($-0.060$, $[-0.206,+0.053]$) include zero; 17 of the 30 patients improve. The effect of XP is therefore small and resolved only for the mean per-patient error. The patient-level anchor of v4 lowers, relative to XP, the pooled MAPE by 0.163 ($[-0.280,-0.021]$), the macro MAPE by 0.172 ($[-0.281,-0.056]$) and the trimmed MAPE by 0.142 ($[-0.241,-0.027]$), while the change of the median, $-0.073$ ($[-0.329,+0.125]$), is not resolved; 22 of the 30 patients improve. v4 raises the pooled correlation from 0.754 to 0.793 and the median within-patient $R^2$ from $-0.417$ to $-0.164$, beats persistence in 29 of the 30 patients, and does so with 19 inputs instead of 55, 83 trees instead of 792 and 450 leaves per fold instead of 9527. Keeping only the last of two readings in the same hour (3427 strict rows) changes the pooled MAPE to 3.919% (X0), 3.898% (XP) and 3.741% (v4) and leaves the ranking unchanged.

Table 10. Out-of-fold accuracy and size on the 3730 rows of 30 patients whose inputs contain no $S_a{O}_{2}m$ value younger than nine hours; identical patient-grouped folds. Trees and leaves are per-fold means. The persistence column scores $y\left({\tau }_i\right)$ itself, with no model.

Quantity	Persistence	X0	XP	v4
tree inputs $n_m$	—	55	55	19
trees per fold	—	792.4	792.4	83.2
leaves per fold	—	9527.4	9527.4	450.2
pooled MAPE (%)	4.308	3.900	3.890	3.727
macro MAPE (%)	3.499	3.321	3.289	3.117
trimmed MAPE (%)	2.038	1.878	1.850	1.709
median per-patient MAPE (%)	1.982	1.875	1.815	1.742
pooled Pearson r	0.728	0.754	0.754	0.793
median $R_{\mathrm{within}}^2$	$-0.800$	$-0.427$	$-0.417$	$-\mathbf{0}.\mathbf{164}$
patients beating persistence (of 30)	—	27	27	29

Table 10. Out-of-fold accuracy and size on the 3730 rows of 30 patients whose inputs contain no $S_a{O}_{2}m$ value younger than nine hours; identical patient-grouped folds. Trees and leaves are per-fold means. The persistence column scores $y\left({\tau }_i\right)$ itself, with no model.

Quantity	Persistence	X0	XP	v4
tree inputs $n_m$	—	55	55	19
trees per fold	—	792.4	792.4	83.2
leaves per fold	—	9527.4	9527.4	450.2
pooled MAPE (%)	4.308	3.900	3.890	3.727
macro MAPE (%)	3.499	3.321	3.289	3.117
trimmed MAPE (%)	2.038	1.878	1.850	1.709
median per-patient MAPE (%)	1.982	1.875	1.815	1.742
pooled Pearson r	0.728	0.754	0.754	0.793
median $R_{\mathrm{within}}^2$	$-0.800$	$-0.427$	$-0.417$	$-\mathbf{0}.\mathbf{164}$
patients beating persistence (of 30)	—	27	27	29

3.2. The Anchor Coefficients of v4

Table 11 lists the coefficients of Eq. (20) fitted on each training fold. $c_1$ and $c_2$ keep their sign in every fold, and their sum is $0.964$ on average, so the anchor is, to four per cent, a weighted average that gives the nine-hour-old measurement $26\%$ of the weight and the patient’s own history $74\%$; the intercept is small and changes sign across folds. The anchor alone, before any tree, has a pooled MAPE of $3.695\%$ against $4.130\%$ for the persistence value it recalibrates.

Table 11. Coefficients of the v4 anchor, Eq. (20), per training fold (locked configuration file, key anchor_coefficients_per_fold).

Fold	${\mathit{c}}_0$	${\mathit{c}}_1$ (last measurement)	${\mathit{c}}_2$ (patient history)	${\mathit{c}}_1+{\mathit{c}}_2$
1	$9.325$	$0.2429$	$0.6585$	$0.901$
2	$2.252$	$0.2531$	$0.7223$	$0.975$
3	$5.636$	$0.3056$	$0.6378$	$0.943$
4	$-2.883$	$0.1933$	$0.8350$	$1.028$
5	$2.508$	$0.2412$	$0.7313$	$0.972$
mean	$3.368$	$0.2472$	$0.7170$	$0.964$

Table 11. Coefficients of the v4 anchor, Eq. (20), per training fold (locked configuration file, key anchor_coefficients_per_fold).

Fold	${\mathit{c}}_0$	${\mathit{c}}_1$ (last measurement)	${\mathit{c}}_2$ (patient history)	${\mathit{c}}_1+{\mathit{c}}_2$
1	$9.325$	$0.2429$	$0.6585$	$0.901$
2	$2.252$	$0.2531$	$0.7223$	$0.975$
3	$5.636$	$0.3056$	$0.6378$	$0.943$
4	$-2.883$	$0.1933$	$0.8350$	$1.028$
5	$2.508$	$0.2412$	$0.7313$	$0.972$
mean	$3.368$	$0.2472$	$0.7170$	$0.964$

3.3. How Much of Each Prediction Is Physics

Table 12 applies the shares of Section 2.7 to the strict rows. In X0 and XP the anchor is the raw persistence value and the trees’ variance share is negative ($-13.7\%$, $-13.1\%$): the trees pull the noisy nine-hour-old value toward the patient’s level, which is the shrinkage of Eq. (22) done split by split with 792 trees. In v4 the anchor of Eq. (20) does this in closed form: alone it lowers the pooled MAPE from 4.308% (persistence) to 3.836%, which is 81% of the improvement of the full model (3.727%), it does 71% of the work, and the 83-tree ensemble corrects single rows while adding 1% of the variance. On the median patient the trees add nothing: the anchor alone has a median per-patient MAPE of 1.703%, the full v4 model 1.742%. The identity ${\tilde{y}}_i=A_i+B_i+C_i$ holds to $4\times {10}^{-6}$ on every row.

Table 12. Shares of Eqs. (41)–(42) in per cent on the strict rows. The variance shares of the three parts sum to 100.

Share (%)	X0	XP	v4
variance share of the anchor A	113.65	113.99	99.94
variance share of the trees B	$-13.65$	$-13.11$	1.02
variance share of the cumulative OUR term C	0.00	$-0.88$	$-0.97$
work share of the anchor A	0.00	0.00	71.33
work share of the trees B	100.00	84.76	20.44
work share of the cumulative OUR term C	0.00	15.24	8.23
error credit of the anchor	0.00	0.00	81.27
pooled MAPE of the anchor alone (%)	4.308	4.308	3.836

Table 12. Shares of Eqs. (41)–(42) in per cent on the strict rows. The variance shares of the three parts sum to 100.

Share (%)	X0	XP	v4
variance share of the anchor A	113.65	113.99	99.94
variance share of the trees B	$-13.65$	$-13.11$	1.02
variance share of the cumulative OUR term C	0.00	$-0.88$	$-0.97$
work share of the anchor A	0.00	0.00	71.33
work share of the trees B	100.00	84.76	20.44
work share of the cumulative OUR term C	0.00	15.24	8.23
error credit of the anchor	0.00	0.00	81.27
pooled MAPE of the anchor alone (%)	4.308	4.308	3.836

3.4. Ablation of the Two C++ Coefficients in v4

Table 13 evaluates the trained v4 model with $\beta$, ${\beta }_1$ or both set to zero in Eq. (37). The two coefficients act only on the predictions (Section 2.6), so the trees are the same in every row of the table, and each row measures what a coefficient contributes to the predictions of the validation rows. The predictions with a coefficient set to zero are not stored row by row, so the table is on all 4017 rows. Setting $\beta =0$ raises the pooled MAPE by 0.032 percentage points, the trimmed by 0.027 and the median by 0.061; setting ${\beta }_1=0$ raises the pooled, macro and trimmed MAPE by about 0.004 and lowers the median by 0.025. Both coefficients are kept by our rule that a physics coefficient is set to zero only when doing so changes the median per-patient MAPE by at most 0.001 and the macro by at most 0.005.

Table 13. Change of each metric when a coefficient of Eq. (37) is set to zero in the trained v4 model; all 4017 rows, percentage points.

Set to zero	$\mathit{\Delta }$ pooled	$\mathit{\Delta }$ macro	$\mathit{\Delta }$ trimmed	$\mathit{\Delta }$ median
${\beta }_1=0$	$+0.0043$	$+0.0037$	$+0.0041$	$-0.0254$
$\beta =0$	$+0.0318$	$+0.0071$	$+0.0270$	$+0.0614$
$\beta ={\beta }_1=0$	$+0.0249$	$+0.0027$	$+0.0245$	$+0.0038$

Table 13. Change of each metric when a coefficient of Eq. (37) is set to zero in the trained v4 model; all 4017 rows, percentage points.

Set to zero	$\mathit{\Delta }$ pooled	$\mathit{\Delta }$ macro	$\mathit{\Delta }$ trimmed	$\mathit{\Delta }$ median
${\beta }_1=0$	$+0.0043$	$+0.0037$	$+0.0041$	$-0.0254$
$\beta =0$	$+0.0318$	$+0.0071$	$+0.0270$	$+0.0614$
$\beta ={\beta }_1=0$	$+0.0249$	$+0.0027$	$+0.0245$	$+0.0038$

3.5. The Cost Weights Are Inert

The recorded target ranges from 20 to 111 % (Section 2.1). Over all 4017 rows, the largest relative deviation of $w_i$ from $k_w$ in Eq. (2) is $1.0\times {10}^{-4}$ in X0 and XP ($\nu =3.656$, $k_w=0.148$) and $9.0\times {10}^{-14}$ in v4 ($\nu =9.174$, $k_w=0.928$), so the weights are constant to that precision and the weighted least squares of Eq. (21) and the weighted objective of the trees are unweighted in effect: the tuned values sit at the variance end of the trade-off of Eq. (5). The squared-relative-error end would matter only if $(1-k_w)/{\left|y\right|}^{\nu }$ were comparable with $k_w$, which for y between 50 and 100 % requires $k_w$ of order ${10}^{-4}$ at $p=2$, below the range the tuning explored; exploring it is left for future work.

3.6. The Age of the Anchor and the Persistence Baseline

The persistence value is the most recent measurement at or before $t_i-9$ h, Eq. (17). Of the 4017 rows, $89.5\%$ have an anchor exactly nine hours old; $3.3\%$ have an older one (median age $22.0$ h, 90th percentile $74.8$ h, maximum 81 h) because no measurement fell in the preceding hour; and $7.1\%$ lie in the first nine hours of a stay and use the patient’s first measurement, whose age is below nine hours and is zero in 33 rows ($0.8\%$): the first row of each patient and, for patients 2 and 3, a second reading in that same hour. Because the target is integer-valued, the nine-hour-old value equals the target exactly in $17.1\%$ of the rows; in the rows that have a measurement at least nine hours old this is quantisation of the label, not leakage, and it sets the floor that Table 10 reports for persistence.

3.7. Strict Nine-Hour Check

The target-derived inputs use $S_a{O}_{2}m$ values younger than nine hours in one situation: when a patient has no measurement at or before $t_i-D$, that is, in the first D hours of a stay. There $y\left({\tau }_i\right)$ and SaO2m_lag_10h fall back to the patient’s first measurement (data_preparation.py, called with fallback_to_first=True), ${\overline{y}}_{{\mathcal{H}}_i}$ is filled with the same value (utils_prediction.py), and the tree input SaO2m_first is that first measurement too (data_preparation.py). This concerns 287 of the 4017 rows ($7.1\%$) in every scenario, and in 33 of them the value used is the current target itself; one patient has no row outside this window, so the strict rows cover 30 patients. No value from a later hour is used anywhere: the lags are taken with a backward merge, the patient mean stops at $t_i-D$, and SaO2m_first is the earliest reading. Because the validation patients of a fold never appear among its training patients, these inputs in training rows carry no information about the validation rows; removing the 287 rows at scoring therefore leaves out-of-fold predictions whose inputs contain no $S_a{O}_{2}m$ value younger than nine hours. Table 14 compares the two evaluations.

Table 14. Out-of-fold metrics on all 4017 rows (31 patients) and on the 3730 strict rows (30 patients).

	X0		XP		v4
Metric	all	strict	all	strict	all	strict
pooled MAPE (%)	3.791	3.900	3.766	3.890	3.603	3.727
macro MAPE (%)	3.272	3.321	3.209	3.289	3.029	3.117
trimmed MAPE (%)	1.962	1.878	1.906	1.850	1.758	1.709
median per-patient MAPE (%)	1.936	1.875	1.850	1.815	1.737	1.742
pooled Pearson r	0.763	0.754	0.764	0.754	0.801	0.793

Table 14. Out-of-fold metrics on all 4017 rows (31 patients) and on the 3730 strict rows (30 patients).

	X0		XP		v4
Metric	all	strict	all	strict	all	strict
pooled MAPE (%)	3.791	3.900	3.766	3.890	3.603	3.727
macro MAPE (%)	3.272	3.321	3.209	3.289	3.029	3.117
trimmed MAPE (%)	1.962	1.878	1.906	1.850	1.758	1.709
median per-patient MAPE (%)	1.936	1.875	1.850	1.815	1.737	1.742
pooled Pearson r	0.763	0.754	0.764	0.754	0.801	0.793

The ranking of the three scenarios is the same in both evaluations. A deployable sensor under the strict rule starts nine hours after a patient’s first measurement; rows before that point would need a population prior as the anchor, which is not part of this study.

4. Discussion

Three findings carry the paper, all measured on the rows whose inputs contain no $S_a{O}_{2}m$ value younger than nine hours. First, the physics pays through the anchor. The patient-level recalibration of v4 is a two-coefficient linear term, yet on its own it earns 81% of the improvement over persistence, and with it the pooled, macro and trimmed errors fall by 0.14 to 0.17 percentage points relative to XP, with bootstrap intervals that exclude zero, while the tree ensemble shrinks to 19 inputs and 450 leaves. Second, the anatomy of Table 12 explains why: with a raw persistence anchor the trees spend a negative variance share on shrinking a noisy measurement toward the patient’s level, a task that Eq. (22) performs in closed form once the patient’s own history is available; freed from it, a much smaller ensemble corrects single rows. Third, the direct contribution of the gas signals through the physics terms is small: the cumulative oxygen uptake rate and the blend of XP lower the mean per-patient error by 0.032 percentage points, an effect whose interval only just excludes zero and which the pooled, trimmed and median summaries do not resolve.

What the sensor does and does not do is now explicit. It reproduces each patient’s level and beats the last blood gas in 29 of 30 patients; it does not resolve hourly changes for typical patients, whose hourly standard deviation of about 1.9 percentage points is at the label’s integer quantisation, and the median within-patient $R^2$ stays negative in all scenarios. The median per-patient error is not resolved by either step, and on the median patient the anchor of v4 alone is as accurate as the full model.

The claims are bounded by the design. The cohort is small (30 patients on the strict rows), and three patients hold most of the squared error, which is why the macro, trimmed and median per-patient errors are reported next to the pooled one and every difference is tested with a paired patient bootstrap. The validation rows decide the number of trees through early stopping, which is common to all scenarios and therefore does not favour one of them. The training rows still include the first nine hours of the training patients’ stays; this carries no information about the validation patients, but a pipeline that drops these rows from training as well would be the natural next version.

Within the scope of this Special Issue, the model is a multivariate time-series application with an explicit account of what each component contributes, and its loss is the relative-entropy criterion of our earlier estimators; a natural next step is to add information-theoretic descriptors of the ventilator channels (permutation or sample entropy over the same lag windows) as tree inputs and to test whether they carry the within-patient dynamics that the present inputs do not.

5. Conclusions

A virtual sensor for arterial oxygen saturation across a nine-hour blood-gas gap was built as three nested scenarios of one code base and evaluated on identical folds on the rows whose inputs contain no saturation value younger than nine hours. Adding to XGBoost 2.1.4 a blend toward the anchor and the cumulative oxygen uptake rate over the gap, derived from an oxygen balance of the arterial blood (XP), gives a small improvement of the mean per-patient error. Replacing the raw nine-hour-old measurement by a patient-level calibrated anchor (v4) lowers the pooled, macro and trimmed errors with bootstrap intervals that exclude zero, beats the last blood gas in 29 of 30 patients, and reduces the ensemble to 19 inputs and 450 leaves; the anchor alone earns 81% of the gain over persistence. Every equation is matched to the line of C++ or Python that implements it, and each scenario is reproducible from one configuration file.