Estimating Physiological Dose Optima in Plant Radiohormesis Through Censoring-Aware Hierarchical Bayesian Modeling

1.1 Current Predictive Frameworks and Their Limitations

A meta-regression model was recently introduced to predict the empirical optimum $x^{\mathrm{\text{*}}}$ (the observed $D^{\mathrm{\text{*}}}$) from high-level covariates such as variety, physiological endpoint, and measurement time 31. Although this model showed that radiosensitivity patterns are structured and predictable, the approach fits those optima by ordinary least squares (OLS). OLS assumes that the response entering the regression is the quantity of interest, observed exactly, with errors that are independent and homoscedastic and, for the usual interval estimates, normal. Censoring is a property of how the optimum is observed. When a fitted optimum is reported at the edge of the tested range, the recorded number is a bound on the quantity of interest rather than the quantity itself, and the exact-observation assumption is the one that fails.

In practice, seed-priming experiments are commonly constrained to a narrow dose range; the experiment re-analyzed here used 0–20 Gy 32. When the true optimum lies below or above the tested range, curve-fitting procedures return a boundary value (0 Gy or 20 Gy), yielding a censored estimate of $D^{\mathrm{\text{*}}}$. The OLS model treats these censored values as precise observations, leading to regression coefficients that are biased downward or upward depending on the censoring direction and frequency. These biases propagate into dose recommendations.

1.2 Censoring as a Structural Feature of Irradiation Experiments

Boundary optima are an intrinsic consequence of experimental design. For physiologically vigorous genotypes, the true optimum may exceed 20 Gy, while sensitive genotypes may peak at doses below the minimum tested. Writing $D^{\mathrm{\text{*latent}}}$ for the dose at which the underlying response surface peaks, $[D_L,D_U]$ for the tested range, and $D^{\mathrm{\text{*}}}$ for the optimum a curve fit reports, these outcomes are

An optimum recorded at $D_U$ is right-censored and one recorded at $D_L$ is left-censored; Section 3 gives the full observation model and the likelihood it implies. The distribution of censoring is itself informative: varieties whose optima reach the upper boundary more often may have higher intrinsic oxidative resilience or slower lesion accumulation. Censoring therefore has to be incorporated into inference.

1.3 A Statistical Framework for Latent Dose Estimation

To address these limitations, this work introduces: (i) a Tobit model that accounts for censoring in a frequentist framework, adjusting estimates of latent optima and regression parameters for boundary effects; and (ii) a hierarchical Bayesian model (HBM) that incorporates censoring directly into the likelihood, partially pools information across varieties and endpoints, and delivers posterior distributions of latent radiosensitivity. The Bayesian formulation represents variability across varieties and endpoints explicitly, quantifies posterior uncertainty, and remains consistent with redox kinetic theory.

1.4 Biological and Agronomic Relevance

Improved estimation of $D^{\mathrm{\text{*}}}$ has direct practical relevance because seed-priming protocols must balance stimulation against the risk of entering the inhibitory region of the dose–response curve. A statistical model that accounts for boundary optima can help distinguish between an experimentally observed maximum and a maximum that may have been truncated by the dose window. This distinction is especially important when different physiological endpoints respond over different time scales. Endpoint identity already accounts for part of the variation in the empirical optima of this experiment 31, and the analysis below places the early germination and seedling traits lower in the tested range than the biomass traits. A censoring-aware model can therefore support more cautious, endpoint-specific experimental designs.

The remainder of the manuscript sets out the steps from the censored observations to latent-dose inference and endpoint-specific dose summaries.

2.1 Dataset Structure and Preprocessing

The dataset comprises fitted optima ($D^{\mathrm{\text{*}}}$) from multiple dose–response curves across tobacco cultivars: Burley (BY), Dark (DK), and Virginia (VA). The underlying irradiation experiment and its dose–response measurements are reported in full by Alves et al. 32; the present work re-analyzes those measurements and introduces no new observations. The same measurements were previously fitted per variety and endpoint with quadratic and Brain–Cousens models, and the resulting optima summarized by the meta-regression described below 31; the optima used here are recomputed by quadratic fits alone, as set out in Section 3.

The physiological endpoints are the developmental traits recorded in that experiment, summarized into the nine endpoints analyzed here. Three of the traits are time-course measurements: germination, hypocotyl emergence, and cotyledon expansion, each recorded on days 3, 5, 7, and 9 after sowing as a count out of the 50 seeds sown per replicate. Each is summarized in two ways, giving six endpoints: the day-9 count, which is the terminal value of that trait, and the sum of the four counts recorded for it. The remaining three are plantlet biomass at two weeks, plantlet biomass at four weeks, and the increment between them. Table 1 lists all nine.

All optima were extracted by least-squares quadratic fits to the dose-wise means over a common dose range of 0–20 Gy, by the procedure set out in Section 3.

After constructing dose–response summaries as described below, the meta-dataset used for regression comprises $N\mathrm{\text{=}}27$ optimal-dose estimates $D^{\mathrm{\text{*}}}$, each corresponding to a specific combination of variety, endpoint, and time scale within the dose window $[0,20]$ Gy. As detailed in Section 3, two of these optima are formally right-censored at 20 Gy, while the remainder lie strictly inside the experimental range.

2.2 Baseline OLS Model

The existing predictive model 31 is:

where $\epsilon \sim N(0,{\sigma }^2)$ and the predictors code varieties, physiological endpoints, and measurement time.

OLS assumes complete observations and ignores censoring, which structurally biases the estimates when boundary optima are present.

2.3 Mechanistic Context: Redox–Lesion Kinetics

Generalized kinetic models describing reactive oxygen species (ROS) production ($R$), antioxidant pools ($A$), and lesion load ($X$) yield a scaling relation of the form

in the regime $D^{\mathrm{\text{*}}}\ll D_f$, where $D_f$ is the fidelity-collapse threshold and $(a,q,p)$ summarize adaptive and damaging processes. The exponent is negative because raising adaptive capacity $a$ relative to damage $q$ moves the optimum to a higher dose, so the ratio $D_f\mathrm{\text{/}}{D}^{\mathrm{\text{*}}}$ decreases. The lesion load $X$ is the effective burden of persistent oxidative lesions, that is, the macromolecular damage left after repair; that model keeps it deliberately unresolved, aggregating the damage classes rather than tracking DNA strand breaks, oxidized bases, and lipid peroxidation products separately 12. A predictive model for $D^{\mathrm{\text{*}}}$ should remain compatible with this mechanistic structure so that its parameters can, at least qualitatively, be related to redox–lesion kinetics 12, 31, 33.

3.1 Construction of Endpoint-Specific Optima ${\mathbf{D}}^{\mathrm{\text{*}}}$

3.2 Quadratic Approximation and Preliminary ${\mathbf{D}}^{\mathrm{\text{*}}}$

For each combination $(v,e)$ the quadratic model

is fitted by ordinary least squares over the available doses $D_j$. The fitted quadratic has vertex

The sign of $c_{v,e}$ is tested first. If $c_{v,e}\ge 0$ the quadratic has no interior maximum and the optimum falls on a boundary; the two boundary predictions are compared, $D_{v,e}^{\mathrm{\text{*}}}$ is set to the boundary with the higher predicted response, and the observation is marked as censored.

If $c_{v,e}\mathrm{\text{<}}0$ the quadratic is concave and $D_{v,e}^{\mathrm{\text{quad}}}$ is its maximum. When $D_{v,e}^{\mathrm{\text{quad}}}\in [D_L,D_U]$ the vertex is taken as the preliminary optimum $D_{v,e}^{\mathrm{\text{*}}}$ and labeled uncensored. When it falls outside the interval the optimum is placed at the nearer boundary and marked as censored:

This yields a meta-dataset in which each row corresponds to a variety $v$, endpoint $e$, and time scale, with an associated preliminary optimum $D_{v,e}^{\mathrm{\text{*}}}$ and an indicator flag specifying whether the optimum is censored at $D_L$ or $D_U$. This table serves as input for the subsequent censored regression analysis.

3.3 Censored Regression (Tobit Model)

3.4 Observation Mechanism and Censored Likelihood

The observed optimum $D_i^{\mathrm{\text{*}}}$ is related to the latent value via

With $f$ and $F$ denoting the normal density and cumulative distribution function with mean $x_i^{\top }\beta$ and variance ${\sigma }^2$, the likelihood contributions are

for uncensored, left-censored, and right-censored sets $U,L,R$. The full likelihood is

3.5 OLS as a Naïve Benchmark

A naïve OLS regression that treats $D_i^{\mathrm{\text{*}}}$ as fully observed ignores censoring. Boundary optima are interpreted as exact, which induces bias: right-censored optima at $D_U$ are treated as if the true optimum were exactly $D_U$, rather than somewhere beyond that point. The Tobit model explicitly accounts for this by treating boundary values as inequalities in the likelihood.

3.6 Hierarchical Bayesian Model (HBM)

The classical Tobit model does not explicitly represent hierarchical structure across varieties and endpoints, and it is not naturally framed for prediction with partial pooling. A hierarchical Bayesian model that generalizes the Tobit formulation is therefore adopted.

3.7 Hierarchical Structure and Linear Predictor

Let $v[i]$ denote the variety associated with observation $i$ and $g[i]$ denote the physiological endpoint. Random effects for varieties and endpoints and a time effect are introduced:

where ${\alpha }_{v[i]}$ is a variety-specific intercept, ${\zeta }_{g[i]}$ is an endpoint-specific intercept, ${\mathrm{\text{Time}}}_i^{\mathrm{\text{(scaled)}}}$ is the time covariate after centering and scaling, and $\gamma$ is a fixed slope.

Because each endpoint is observed on a single time scale, the time covariate is a deterministic function of endpoint identity. The coefficient $\gamma$ is therefore not identified as a dose–time relationship and is not interpreted as one; it acts as a grouping of endpoints by measurement scale, which the hierarchical prior shrinks less strongly than the individual endpoint effects.

3.8 Logit–Normal Mapping to the Dose Interval

The physiological optimum dose is constrained to lie in $[D_L,D_U]$. ${\eta }_i\mathbb{\in }ℝ$ is therefore mapped to the dose scale by

with

This yields a logit–normal prior structure for the latent optimum, so that $D_i^{\mathrm{\text{*latent}}}\in (D_L,D_U)$ by construction.

3.9 Observation Model with Additive Noise

Conditional on $D_i^{\mathrm{\text{*latent}}}$, the preliminary optimum $D_i^{\mathrm{\text{*}}}$ is modeled as

subject to the same censoring mechanism at $D_L$ and $D_U$:

where $\phi$ and $\Phi$ are the standard normal density and distribution functions.

3.10 Hierarchical Priors

Hierarchical priors are assigned to the random effects:

with hyperpriors

The time effect $\gamma$ is assigned a weakly informative prior

after centering and scaling the time covariate. The residual standard deviation is modeled as

3.11 Full Joint Distribution

Collecting the model parameters in $\Theta \mathrm{\text{=}}\{{\alpha }_v,{\zeta }_g,{\mu }_{\alpha },{\mu }_{\zeta },{\tau }_{\alpha },{\tau }_{\zeta },\gamma ,\sigma \}$ and denoting $D^{\mathrm{\text{*}}}\mathrm{\text{=}}(D_1^{\mathrm{\text{*}}},\dots ,D_N^{\mathrm{\text{*}}})$, the joint posterior is

with ${\eta }_i\left(\Theta \right)\mathrm{\text{=}}{\alpha }_{v[i]}\mathrm{\text{+}}{\zeta }_{g[i]}\mathrm{\text{+}}\gamma \cdot {\mathrm{\text{Time}}}_i^{\mathrm{\text{(scaled)}}}$.

3.12 Posterior Inference via Hamiltonian Monte Carlo

Posterior inference is performed using the No-U-Turn Sampler (NUTS) implemented in PyMC. Four independent Markov chains were run, each with 1500 tuning iterations at a target acceptance probability of 0.95 followed by 2000 post-tuning draws, yielding 8000 posterior samples. Sampling was well behaved: no divergent transitions were recorded, the split-$\hat{R}$ statistic equaled 1.00 for every parameter, and the smallest bulk effective sample size was 2901 of the 8000 draws. The sampler produced an Inference Data object containing posterior draws, sampling statistics, and the observed data.

Posterior summaries for variety and endpoint effects are obtained by transforming draws of ${\eta }_i$ through the logistic function to obtain draws of $D_i^{\mathrm{\text{*latent}}}$ for each observation and aggregating across the relevant groups.

3.13 Model Comparison and Information Criteria

Three models are compared:

• 1. a naïve OLS regression that ignores censoring;

• 2. a classical Tobit regression with right-censoring at $D_U\mathrm{\text{=}}20$ Gy; and
• 3. the hierarchical Bayesian model described above.

Maximized log-likelihoods were computed for the OLS and Tobit models. The OLS regression yielded a maximized log-likelihood of $\mathrm{\text{-}}59.1402$, whereas the right-censored Tobit model yielded a maximized log-likelihood of $\mathrm{\text{-}}59.1506$, so

The two fits carry the same regression parameters and the same residual scale and differ only in how the two boundary observations enter the likelihood; they are therefore not nested, and no likelihood-ratio test applies. The comparison establishes only that the difference is of the order of one hundredth of a log-likelihood unit: with this sample size and this censoring pattern, treating the two boundary values as inequalities rather than as exact observations leaves the maximized likelihood essentially unchanged. The case for the censored formulation therefore rests on its structure rather than on a measured gain in fit.

For the hierarchical Bayesian model the Watanabe–Akaike Information Criterion (WAIC) and approximate leave-one-out cross-validation using Pareto-smoothed importance sampling (PSIS–LOO) were computed from the pointwise censored-normal log-likelihood. Writing $logp(D_i^{\mathrm{\text{*}}}|{\Theta }^s)$ for the log-likelihood at observation $i$ and posterior draw $s\mathrm{\text{=}}1,\dots ,S$, the log pointwise predictive density is

from which WAIC estimates the expected log predictive density by subtracting an effective number of parameters $p_{\mathrm{\text{WAIC}}}$, obtained from the posterior variance of the pointwise log-likelihoods:

On the log (elpd) scale reported by ArviZ, the hierarchical model gave

so $lppd\mathrm{\text{=-}}61.47$ and $\mathrm{\text{WAIC=}}134.48$ on the deviance scale. ArviZ issued a diagnostic warning that for four observations the posterior variance of the log predictive density exceeded 0.4, the largest being 1.16, so WAIC should be read with caution for those points. They are the two censored biomass optima, the two-week biomass optimum of DK, and the low hypocotyl-AUC optimum of VA.

PSIS–LOO approximates exact leave-one-out cross-validation using importance weights stabilized by a generalized Pareto distribution. For the hierarchical model, on the same elpd scale,

consistent with the WAIC estimate ($\mathrm{\text{-}}67.24$ versus $\mathrm{\text{-}}67.68$, well within either standard error). The Pareto shape parameter $k$ exceeded 0.70 for one of the 27 observations, the two-week biomass optimum of VA at $k\mathrm{\text{=}}0.84$, with no value above 1. The importance-sampling approximation is therefore reliable apart from that single point, whose leave-one-out contribution should be treated as approximate.

The effective number of parameters is small relative to the sample size ($p_{\mathrm{\text{WAIC}}}\mathrm{\text{=}}5.8$ and $p_{\mathrm{\text{LOO}}}\mathrm{\text{=}}6.2$ against $N\mathrm{\text{=}}27$), which is what partial pooling of three variety effects and nine endpoint effects should produce, and the two criteria agree to within their standard errors. They indicate that the model is not overfitting the 27 optima, subject to the single influential observation noted above. Both flags fall on the censored and the most extreme optima, which is where the sample is thinnest. WAIC and PSIS–LOO were not computed for the OLS or Tobit models, and no numerical ranking of the three models is made. Table 2 collects the values.

3.14 Reproducibility and Software

All analyses were performed in Python using the PyMC ecosystem for Bayesian modeling and ArviZ for posterior diagnostics and information criteria. The code used to construct the endpoint tables, fit the censored hierarchical models, and generate the figures is available from the corresponding author upon reasonable request and can be readily adapted to other dose–response datasets. The posterior reported here was drawn with a fixed random seed, so the figures and every quoted value regenerate exactly.

4.1 Effects of Censoring on Regression Structure

A central objective of this study is to understand how censoring affects the regression structure linking varieties, physiological endpoints, and time to the dose of maximal stimulation. In the meta-dataset, 27 optima are available and two are formally right-censored at 20 Gy. Several additional optima lie close to the upper boundary and behave, in practice, like censored points.

In a naïve OLS analysis, all boundary values are treated as exact optima. This forces the regression surface to accommodate artificial constraints and tends to compress apparent differences between groups, particularly for endpoints whose fitted optima are at or near 20 Gy. The Tobit likelihood treats a boundary observation as an inequality on the dose scale,

In the hierarchical logit–normal model, by contrast, $D_i^{\mathrm{\text{*latent}}}$ is constrained to the experimental interval $(D_L,D_U)$, and right-boundary observations contribute through the upper-tail probability of the censored observation model. Thus, in both cases, a boundary value is not treated as an exact interior optimum.

In the hierarchical model, variety and endpoint effects are shrunk toward a global mean, but in a way that respects censoring: groups associated with boundary optima are not forced to share the same latent optimum as groups whose optima lie well inside the experimental range. This improves interpretability of differences across endpoints relative to the OLS benchmark.

4.2 Latent vs. Observed Optima

Observed optima $D_i^{\mathrm{\text{*}}}$ arise from quadratic fits to noisy dose–response curves over a finite set of doses. Curves whose maxima are not well resolved within the tested range may therefore produce apparent optima near the upper boundary. By contrast, the hierarchical Bayesian model returns a posterior distribution for the underlying optimum within the experimental interval, integrating information across varieties and endpoints.

Figure 1 shows the effect. The latent estimates separate by endpoint class in a way the observed optima do not: the 18 early-vigor observations occupy 9.7 to 11.3 Gy and the nine biomass observations 14.7 to 16.7 Gy, a gap of 3.4 Gy with no overlap, whereas the observed optima of the two groups overlap heavily (3.6–13.0 Gy against 11.2–20.0 Gy). Within each band the latent value varies little with the observed one, so a least-squares line through the 27 pairs has slope 0.55 rather than 1: partial pooling replaces each individual optimum with the value implied by its variety and endpoint. The extreme values move most. The 3.6 Gy hypocotyl-AUC optimum of VA is placed at 9.7 Gy, and the two optima censored at 20 Gy at 15.6 and 16.0 Gy. The censored pair falls in the upper quarter of the tested window, which is the closest the logit–normal mapping can place an optimum that the data bound only from below.

4.3 Posterior Predictive Checks

Posterior predictive checks (PPCs) assess absolute model adequacy. For each posterior draw a replicated dataset $D_i^{\mathrm{\text{*rep}}}$ was generated by sampling latent optima from the logit–normal structure and applying the censored observation model. Four statistics of each replicate were then compared with the observed data: the mean of the 27 optima, their standard deviation, the number right-censored at 20 Gy, and the height of the tallest 2 Gy histogram bin. The first two test location and spread, the third tests whether the model reproduces boundary censoring at the observed rate, and the fourth tests how tightly the interior optima are concentrated.

The model reproduces both the location and the overall spread of the optima. The observed $D^{\mathrm{\text{*}}}$ had mean 11.84 Gy and standard deviation 3.38 Gy; across 1000 replicated datasets the mean of the replicate means was 11.95 Gy, with standard deviation 0.75 Gy across replicates, against the $3.38\mathrm{\text{/}}\sqrt{27}\mathrm{\text{=}}0.65$ Gy that the observed sample implies for a mean of 27 optima. The within-replicate standard deviation averaged 3.85 Gy, with 95% interval 2.67 to 5.17 Gy, so the observed 3.38 Gy falls inside the predictive range.

The frequency of boundary optima is reproduced as well. Two optima were right-censored at 20 Gy in the data; across replicates the censored count had mean 0.68 (standard deviation 0.93; central 95% interval $[0,3]$), two or more arose in 15.1% of replicated datasets and none in 54.2% (Figure 3), so the observed count of two lies inside the predictive distribution. This follows from treating those two observations as lower bounds rather than as exact values. The censored likelihood places their latent optima at 15.6 and 16.0 Gy. At the posterior mean of the observation noise, $\sigma \mathrm{\text{=}}2.86$ Gy, a latent optimum that close to the boundary carries a replicated draw past 20 Gy with appreciable probability.

The model does not reproduce the concentration of the interior optima. Eleven of the 27 observed optima fall in the single 2 Gy interval from 10 to 12 Gy (Figure 2), whereas the tallest bin of a replicated dataset holds 7.1 optima on average, with 95% interval 5 to 10, and reaches eleven in 1.9% of replicates. The observation noise required to accommodate the extremes of this sample is therefore larger than the dispersion of its center, and the interior optima are spread more evenly than in the data. The interval estimates in Tables 3 and 4 are consequently more likely to be too wide than too narrow.

One limitation is structural, the logit–normal mapping constrains every latent optimum to lie strictly inside $(0,20)$ Gy, so an endpoint whose true optimum lies beyond the tested range can be placed near the upper boundary but never past it. For the biomass endpoints, with posterior means from 15.0 to 16.2 Gy, the model should therefore not be read as evidence that the optimum lies within 20 Gy; the cautious reading is that the tested window may be incomplete and that a wider dose range is needed to resolve the upper tail.

4.4 Empirical Structure of BY Optima

Before the hierarchical summaries, the empirical optima for variety BY are examined on their own. Figure 4 displays the distribution of observed $D^{\mathrm{\text{*}}}$ values obtained from quadratic fits, showing a concentration around mid-range doses and no left-censored values, with one right-censored observation at 20 Gy.

Figure 5 gives one dose–response curve in full, cotyledon expansion at day 9 in BY, and shows the basis of a single extracted optimum. The dose-wise means occupy a range of four seedlings out of 50, they fall and rise again above the fitted maximum, and the quadratic accounts for a quarter of their variation ($R^2\mathrm{\text{=}}0.24$), so the vertex at 10.3 Gy is a weakly determined summary of that curve. This is the uncertainty the hierarchical model carries, and the reason the individual optima are pooled rather than interpreted separately.

4.5 Qualitative Patterns Across Endpoints and Varieties

Figure 6 presents posterior estimates of the latent optimal dose $D_{\mathrm{\text{latent}}}^{\mathrm{\text{*}}}$ for each physiological endpoint, with posterior means and 94% highest-density intervals, and Table 3 gives the same summaries numerically. Early-stage traits (germination, hypocotyl, cotyledons) have posterior means from 10.08 to 10.55 Gy, close to the midpoint of the 0–20 Gy window, whereas biomass-related endpoints ($W_2$, $W_4$, $\Delta W$) sit higher, from 14.97 to 16.23 Gy, and include right-censored contributions. Averaged across varieties, the biomass optima exceed the early-vigor optima by 5.17 Gy, with 94% highest-density interval 2.89 to 7.47 Gy, and the posterior probability that the biomass group carries the higher optimum is 0.9998.

At the variety level, Figure 7 and Table 4 show that the three Nicotiana cultivars have overlapping posterior intervals, with latent optima between 11.63 and 12.64 Gy. Genotypic differences in the location of the hormetic peak are therefore modest compared with endpoint-driven variation.

Finally, Figure 8 displays the observed optima for each endpoint in BY. Vigor-related endpoints tend to have interior maxima at lower doses, whereas the biomass endpoints shift upward and include the right-censored $\Delta W$ optimum at 20 Gy. The two sets interleave in this variety: the hypocotyl-AUC optimum at 13.0 Gy lies above the two-week biomass optimum at 12.9 Gy. The raw pattern therefore agrees with the hierarchical inference in direction. Table 5 sets out what these posterior patterns imply for the design of subsequent experiments.

4.6 Summary of Main Patterns

The censored hierarchical Bayesian model characterizes the optimal dose estimates of multiple traits under boundary censoring. It corrects the main conceptual deficiency of OLS (treating boundary optima as exact), reveals endpoint-specific differences that agree with the raw data, and yields posterior distributions that can be used in downstream decision-support tools. At the same time, the small sample size and the limited number of censored points indicate that model-based conclusions should be interpreted cautiously, especially for extrapolation beyond the observed dose range. The results should therefore be viewed as a proof-of-concept application of censored hierarchical modeling to plant radiohormesis.

5.1 Limitations

The present analysis should be interpreted in light of several limitations. First, the meta-dataset contains only 27 variety-by-endpoint optimal-dose estimates, which restricts the precision with which variety, endpoint, and time effects can be separated. Second, only two observations are formally right-censored at the upper experimental boundary, although several additional optima lie close to 20 Gy. Those two observations carry most of the information about the upper tail: they raise the biomass optima from about 13 Gy, where an uncensored treatment leaves them, to 15.0–16.2 Gy. An estimate resting on two observations is correspondingly fragile, and the biomass optima in particular should be re-examined on a design that brackets the peak.

Third, the logit–normal mapping used in the hierarchical model constrains latent optima to the experimental interval $(0,20)$ Gy. This is useful for stable inference within the tested range, but it also means that the model cannot directly estimate true optima beyond 20 Gy. Boundary-supported posterior estimates should therefore be interpreted as evidence that the tested window may be incomplete for some endpoints.

Fourth, although the information criteria are well behaved for this fit, two diagnostics still fall on the sparsest part of the sample: one observation, the two-week biomass optimum of VA, has a Pareto shape parameter of 0.84, and four observations have a pointwise log-predictive variance above the WAIC threshold. Both flags sit on the censored and the most extreme optima, so the fit is least reliable exactly where the censoring information enters; larger datasets with broader dose ranges will be needed to evaluate the stability and generality of the inferred endpoint-specific patterns.

Fifth, the separation between the early-vigor and biomass optima depends on the grouping of endpoints by measurement scale described in Section 3. Omitting that term reduces the separation from 5.17 Gy to 2.69 Gy, with a 94% interval that then includes zero, although its direction is unchanged (posterior probability 0.97). The overlap of the variety optima is unaffected.

5.2 Future Directions

On the experimental side, future studies should use wider and denser dose ranges, especially for biomass-related endpoints whose optima approached the upper boundary in the present dataset. A design including additional treatments above 20 Gy, for example 25 and 30 Gy, would help determine whether the biomass response truly peaks near the upper part of the current interval or whether the full decline occurs beyond it. Additional dose levels below 10 Gy would also be useful for early germination and seedling endpoints.

On the biological side, future work should measure intermediate physiological markers, especially oxidative stress indicators, antioxidant enzyme activity, and early growth dynamics. These measurements would allow the statistical optimum $D^{\mathrm{\text{*}}}$ to be linked more directly to redox adaptation and lesion accumulation. On the modeling side, the present framework could be extended to full dose–response curves rather than extracted optima alone, allowing uncertainty from curve fitting and censoring to be propagated in a single model. Broader genotype panels and validation across species would also allow the hierarchical structure to exploit partial pooling across a larger biological space.