A Bayesian framework to model variance of grain yield response to plant density

Abstract

Background

The expected grain yield response to plant density in winter wheat (Triticum aestivum L.) follows a diminishing returns function. To our knowledge, all previous studies dealing with plant density have assumed constant variance. The gap relies on quantifying the optimum plant density that optimizes grain yield at the lowest risk. Here, we propose a Bayesian hierarchical framework to model the variance of grain yield response to plant density. We demonstrate our framework by identifying the plant density in each seed size, seed treatment and environment combination that maximizes the expected yield and minimizes yield variance.

Results

To fit the model, we used data from field experiments conducted in the Canadian Prairies to identify informative priors and Kansas experiments to demonstrate and validate our framework. Kansas experiments were conducted in 25 environments and consisted of a complete factorial combination of three seed cleaning methods leading to three different seed sizes (light, moderate, heavy), two or three seeding rates, and two seed chemical treatments (insecticide + fungicide vs. none). We described both expected yield and variance of yield in response to plant density. The proposed model allowed us to quantify the minimum risk plant density (minRPD), which represents the minimum plant density at which grain yield variance becomes constant. Plant density at the minRPD was always greater than the agronomic optimum plant density (AOPD, i.e.: the plant density that maximizes expected yield); thus, minRPD could be used to estimate the minimum plant density that maximizes expected yield and minimizes yield variance. When compared at the AOPD, four seed cleaning × chemical treatments combinations resulted in similar yield advantages over the control under high and low yielding environments. However, in low-yielding environments, only two cleaning × chemical treatments combinations resulted in smaller variance when compared at the minRPD against the control. All seed cleaning × chemical treatments combinations resulted in similar AOPD. However, two cleaning × chemical treatments combinations had greater minRPD in low-yield environments compared to the control.

Conclusion

Modeling grain yield response to plant density with the proposed framework is suitable for heteroscedastic data scenarios. Future research may focus on exploring how genotypes, environments and their interaction modulate the difference between AOPD and minRPD and, extend the framework to a variety of processes involving crop management decisions.

Background

The response of winter wheat to plant density has been widely studied [5, 11, 16, 17]. Previous studies focused on the average effects (i.e., expected value) of plant density on yield [4, 44, 46] in interaction with cultivars [2, 7, 30, 43, 44, 47], seeding management strategies (e.g., seed size and seed treatment; [6]), sowing date [39] and environment [26]. These studies have always used similar models that assume constant variance (i.e., ‘spread of the data’) for all plant densities, which is oftentimes overlooked [12]. Despite this fact, no attempts have been made to model yield variance in response to plant density thus far.

Variance models were developed for heteroscedastic data scenarios and have been used for decades [13, 14, 18, 40] in fields as diverse as macroeconomics [13], epidemiology [23] and ecology [37]. More recently, Bayesian modeling has been successfully implemented for modeling heteroscedastic data [34], with multiple advantages compared to frequentist approaches allowing to integrate multiple sources of information (e.g.: selecting informative priors; [3]) and providing the means for making inference based on probabilities [20].

Adopting a fully Bayesian analysis requires specifying objective priors reflecting the current state of knowledge. In a Bayesian setting, informative priors are based on data that weights model parameters. Similarly, in hierarchical modelling structures, which are widely accepted and applied in the scope of agronomy [25, 29], model parameters are informed based on a given process model. Hence, using informative priors supported on previous data not only reflects the appropriate use of Bayesian models but also can increase parameter precision (i.e., narrowing credible intervals) without reducing accuracy (e.g., targeted median minRPD), when specified correctly [33].

The expected wheat yield response to plant density follows a diminishing returns function [17, 48] which is supported by crop eco-physiological processes [46]. The mathematical equation of expected yield response to plant density can vary from quadratic [26, 30], quadratic-plateau [48], linear-plateau [4, 28] or exponential rise to maximum [27, 39, 48]. Across proposed mathematical models, light interception limits yields under low plant densities [17, 46]. At higher plant densities yield reaches a plateau [4, 17, 41]. The nature of the plateau relies on compensation mechanisms to maintain grain yield while increasing plant density, including (i) reduction in the number of shoots per plant, (ii) reduction in the green area per shoot, (iii) lower radiation use efficiency through less evenly light distribution across the vertical canopy profile and (iv) changes in the resource partitioning within the plant related to sink size [46].

Environmental factors interact with crop management practices to influence the optimum plant density required to maximize grain yield [6, 21, 22]. The inherent uncertainty of environmental conditions within a cropping season requires us to think about alternative approaches that account for the uncertainty on optimum plant density quantification. Thus, we developed a Bayesian hierarchical framework to model the variance of grain yield response to plant density. We demonstrate the model by combining information in two datasets of similar characteristics (see Sects. "Developing informative priors: Canadian Prairies data", "Data analysis: Central Great Plains, Kansas US") aiming to identify the plant density in each seed size, seed treatment and environment combination that maximizes expected yield and minimizes yield variance.

Methods

Modeling grain yield response to plant density

To model the density-yield relationship we implement a hierarchical Bayesian framework as follows:

$${y}_{ij}\sim N\left({\mu }_{ij},{\sigma }_{ij}^{2}\right),$$

(1)

$${\mu }_{ij}=\left\{\begin{array}{cc}{\beta }_{1j}{s}_{ij}+{\beta }_{2j}{s}_{ij}^{2},& {s}_{ij}<{AOPD}_{j}\\ {\beta }_{1j}{AOPD}_{j}+{\beta }_{2j}{AOPD}_{j}^{2},& \text{ otherwise}\end{array}\right.,$$

(2)

$${\sigma }_{ij}^{2}=\left\{\begin{array}{cc}0,& {s}_{ij}=0\\ {e}^{{\alpha }_{1j}{s}_{ij}+{\alpha }_{2j}{s}_{ij}^{2}},& {s}_{ij}<{minRPD}_{j}\\ {e}^{{\alpha }_{1j}{minRPD}_{j}+{\alpha }_{2j}{minRPD}_{j}^{2}},& {s}_{ij}>{minRPD}_{j}\end{array}\right..$$

(3)

In Eq. 1${y}_{ij}$ represents an observation of the response variable grain yield at the $ith$ plot ($i=\text{1,2},\dots n$) in the $jth$ combination of seed size, seed treatment and yielding environment ($j=\text{1,2},\dots 12$) and follows a normal distribution with an expected value ${\mu }_{ij}$ and variance ${\sigma }_{ij}^{2}$. The physiological mechanisms described in Sect. "Background" support that ${\mu }_{ij}$ is affected by plant density (${s}_{ij}$) following a quadratic plateau function described in Eq. 2. The quadratic-plateau expression in Eq. 2 is biologically meaningful and allows to interpret model parameters from an agronomic perspective. Here, ${\beta }_{1j}$ and ${\beta }_{2j}$ are the linear terms specific to the jth seed size × seed treatment × environment and ${AOPD}_{j}$ is the agronomic optimum planting density, which is the smallest value of ${s}_{ij}$ that maximizes ${\mu }_{ij}$. In other words, ${AOPD}_{j}$ can be interpreted as the lowest plant density required to reach yield plateau. We note that the intercept in this process model (Eq. 2) is set to zero, as we assume ${\mu }_{ij}$ is equal to zero when ${s}_{ij}$ is zero.

The knowledge gap lies on how yield variance (${\sigma }_{ij}^{2}$) is affected by ${s}_{ij}$. Thus, we propose Eq. 3 to describe the non-linear process of ${\sigma }_{ij}^{2}$ as a function of ${s}_{ij}$. First, ${\sigma }_{ij}^{2}$ is equal to zero when ${s}_{ij}$ is zero. Secondly, ${\sigma }_{ij}^{2}$ increases with ${s}_{ij}$, following an exponential-quadratic function until a maximum. Last, ${\sigma }_{ij}^{2}$ decreases until a point where further increases in plant density do not modify ${\sigma }_{ij}^{2}$ (i.e., ${\sigma }_{ij}^{2}$ is constant). The value of ${s}_{ij}$ where ${\sigma }_{ij}^{2}$ becomes constant is defined by the parameter ${minRPD}_{j}$ (minimum risk plant density). In other words, ${minRPD}_{j}$ represents the minimum value of plant density that minimizes yield variance on the right-hand side of the maximum variance. Likewise, ${minRPD}_{j}$ is a separate value from the ${AOPD}_{j}$, that identifies the minimum plant density required to minimize risk while still within the optimum yield range.

The model described in Eqs. 1–3 can be easily adapted to commonly used experimental designs of agricultural experiments, such as experiments testing the interaction of plant density and genotype (i.e. [44]:) or row spacing (i.e. [16]:). To facilitate the implementation and extensions of our model, we provide a detailed description of the model fitting process in Supplementary material and its implementation within the R programming language (https://github.com/giordanon/yield_variance_plant_density,[35]).

To fully specify our Bayesian model, we need to define the priors for the parameters in Eqs. 1 to 3. Instead of explicitly specifying our priors for ${\beta }_{2j}$ and ${\alpha }_{2j}$, we chose to consider ${AOPD}_{j}$ and ${maxVPD}_{j}$ respectively because they are more interpretable quantities (i.e., their units are in plants m⁻²), thus resulting in a more conscious prior selection. Assuming ${\mu }_{ij}$ and ${s}_{ij}$ are constants, then the ${AOPD}_{j}$ can be expressed as a function of ${\beta }_{1j}$ and ${\beta }_{2j}$ (Eq. 4):

$$\frac{\partial {\mu }_{ij}}{\partial {s}_{ij}}=0={\beta }_{1j}+2{\beta }_{2j}{AOPD}_{j}, {AOPD}_{j}=\frac{-{\beta }_{1j}}{2{\beta }_{2j}}.$$

(4)

Likewise, ${maxVPD}_{j}$, maximum variance plant density, is the value of plant density that maximizes ${\sigma }_{ij}^{2}$ (Eq. 5) and can be expressed as a function of ${\alpha }_{1j}$ and ${\alpha }_{2j}$.

$$\frac{\partial {\sigma }_{ij}^{2}}{\partial {s}_{ij}}=0 \left({\alpha }_{1j}+2{\alpha }_{2j}{maxVPD}_{j}\right) {e}^{{\alpha }_{1j}{maxVPD}_{j}+{\alpha }_{2j}{maxVPD}_{j}^{2}}, {maxVPD}_{j}=\frac{-{\alpha }_{1j}}{2{\alpha }_{2j}}.=$$

(5)

Developing informative priors: Canadian Prairies data

Prior hyperparameters were defined using data of winter wheat response to density in the Canadian Prairies (see details below, Eqs. 6-10) and previous literature [17]. We acknowledge that the Canadian Prairies and Kansas are widely differing environments; however, priors used for the purpose of this study have a support wide enough for model parameters accepetable to both regions. The procedure consisted of: (1) defining flat priors based on previous literature [17], see Supplementary material for details on flat priors); (2) fitting the proposed model (Eq. 1–3) to the Canadian Prairies dataset adopting flat priors; (3) derive quantities (e.g., the expected value, standard deviation, and quantiles) from the fitted model in step 2 and; (4) parametrize priors in the desired model.

Detailed procedures for this trial comprised of field experiments conducted in 26 environments across Canada during three cropping seasons between 2010 and 2013 are cited in Beres et al. [6]. Briefly, this study consisted of a factorial combination of two seeding rates (200 and 400 seeds m⁻²), three seed size levels (small, moderate, or heavy), and two seed treatment levels (check-no seed treatment and fungicide). Seeding rate treatments were imposed to generate a range of plant density (i.e., plants emerged) wide enough to study yield response to plant density. Treatments were applied to the experimental units (ith plot) which were laid out in a randomized complete block design. Seed sizes were determined by passing the bulk seed lot over a no. 9 round top screen and a no. 8 round bottom screen. Scalped seeds from the top of the no. 9 round screen were categorized as “heavy,” while those that scalped off of the no. 8 round screen were deemed “light” and while remaining seeds between both screens resulted in “moderate” seed size. Cultivar CDC Buteo, a hard-red winter wheat cultivar widely adapted to the Canadian Prairies, was sown in all trials. Seed germination was approx. 98% (data not shown). Stand counts were performed during late October-early November in each season by counting 2 linear row meters. Once the crop reached physiological maturity and grain moisture was lower than 15%, plots were harvested with a plot combine. Grain yields were corrected to 13% moisture content.

In order to build the priors indicated in Eqs. 6–10, we fitted the model described in Eqs. 1–3 across all of the environments in the Canadian dataset. For this instance, we defined flat priors based on previous literature [17] allowing a reasonable range of values for our parameter of interest. These flat priors specifications are described in Supplementary Material Sect. "Background".

Equations 6–10 show the resulting priors combining the Canada dataset and previous literature [17] that were used to analyze data from Kansas (see Sect. "Data analysis: Central Great Plains, Kansas US" below):

$${\beta }_{1j}\sim\uppsi \left(\text{0.08,0.01,0},0.20\right),$$

(6)

$${AOPD}_{j}\sim\uppsi \left(10\text{5.29,8.72,10,200}\right),$$

(7)

$${\alpha }_{1j}\sim\uppsi \left(\text{0.031,0.004,0},0.04\right),$$

(8)

$${maxVPD}_{j}\sim\uppsi \left(\text{91.88,7.76,0},111.87\right),$$

(9)

$${minRPD}_{j}\sim\uppsi \left(\text{160.19,15.04,0},198.94\right),$$

(10)

We selected the truncated normal (ψ) distribution for ease of analysis and its capacity to clearly define the support of acceptable values for each of our priors. The ψ distribution is defined by four parameters. In the order listed hyperparameters of the ψ distribution are: the expected value (µ), the precision (σ) and the lower (a) and upper (b) truncation points. Evidence for the lower trunctation point for the ${AOPD}_{j}$ of 10 plants m⁻² was obtained from Fischer et al. [17]. This study was conducted under irrigated conditions with plants sown in raised beds in a honeycomb design in the Yaqui Valley of Mexico (27°22′13", −109°55′39"). The combination of low-latitude environments where freeze damage is unlikely, and agronomic management adopted allows wheat crops to express its tillering potential and maximize yields with a plant stand as low as 10 plants m⁻² (see Fig. 3 in [17]). While environmental conditions and agronomic management explored in this study differ from that of Fischer et al. [17], we utilize their evidence to inform the model about the lowest plausible plant density that can maximize yield in Kansas, where environments are much more restrictive for tillering and stand establishment compared to the Yaqui Valley.

Data analysis: central great plains, Kansas US

Our data to demonstrate the model comprised of rainfed field experiments conducted in Kansas, US during three cropping seasons (2018–2021) in a total of 25 environments (site-year combinations). Each experimental unit (ith plot) was 9 m long by seven rows (0.19 m spaced). A split-split plot design with three to four replications was used to establish 18 treatments, which were created by combining three seeding rates (148, 222 and 297 seeds per m⁻²), three levels of seed cleaning [none (light), air screen (moderate), and gravity table + color sorting (heavy)], and two seed treatments (none and insecticide + fungicide). In order to generate a wide range of plant density (i.e., plants emerged), we imposed seeding rate treatments. Seed treatment consisted of 5 oz ac⁻¹ of Cruiser Maxx (thiamethoxam and mefenoxam–fludioxonil) and 0.75 oz ac⁻¹ Cruiser 5FS (thiamethoxam). The wheat variety evaluated at all locations was 'SY Monument,' which was the most widely sown winter wheat variety in Kansas during the study period [42]. Weeds, insects, and diseases were controlled as needed using standard practices so that they were non-limiting to grain yield. Other practices (e.g., sowing date, tillage, etc.) were specific to each environment (Supplementary Table 1 and 2), not allowing us to explore their interactions with the treatments. Stand counts were done 3–4 weeks after sowing on 1 linear row-meter. Once the crop reached physiological maturity and grain moisture was lower than 15%, experimental units were harvested using a Massey Ferguson XP8 small-plot combine. To avoid border effect, the plot ends were trimmed at harvest time. Grain yields were corrected to 13% moisture content.

Using a k-means clustering algorithm, we classified all 25 environments based on attained yield as high or low yielding environments. The number of clusters was set to two (k = 2). The individual environment classification comprised the following rules: in a given environment if more than 50% of the plot level yield observations fell in the high-yielding cluster the environment was classified as high, or vice versa for low-yielding environments [4, 31, 36].

In order to rapidly quantify our parameters of interest (AOPD and minRPD) in each seed size × seed treatment × environment combination, we adopted a divide and conquer approach [45]. This approach consisted of fitting the preceding model (Eq. 1–3) to 12 subsets of our data corresponding to each permutation of seed size, seed treatment, and environment. Note that the divide and conquer approach does not modify the statistical model but rather serves as a computational advantage. It is easily applicable in this case, since each response curve (for both ${\mu }_{ij}$ and ${\sigma }_{ij}^{2}$) is specific to each seed size × seed treatment × environment combination.

In order to analyze how different treatments, seed sizes and environments impacted yield response to plant density, we compared all treatments against the control treatment (untreated, light seed) in each environment cluster (high or low-yielding environments). For the ease of interpretation of the variable units, the variance of yield values was converted to standard deviation. From the posterior distritbution of model parameters and predictions we calculated the median (q_.50), the upper (q_.025) and lower (q_.975) limits of the 95% credible interval. Comparisons between treatments imposed were calculated as:

$$\Delta {x}_{j}={x}_{j}-{x}_{Control},$$

(11)

where $x$ is one of the following ${\mu }_{j}$, ${\sigma }_{j}$, ${AOPD}_{j}$, ${minRPD}_{j}$.

We evaluated the economic impacts of the following plant density recommendations by the ${minRPD}_{j}$ against ${AOPD}_{j}$:

$${Net Income}_{j} = {y}_{ij} * 220 - {PD}_{j}* 6.16 * {10}^{-5},$$

(12)

$$\Delta {Net Income}_{j}={Net Income}_{{minRPD}_{j}}-{Net Income}_{{AOPD}_{j}},$$

(13)

where the ${NetIncome}_{j}$ is calculated as the difference between income and cost of seed. Income is calculated as the product between the posterior distribution of grain yield ${y}_{ij}$ and the price of wheat grain assumed to be 220 $ Mg ⁻¹. Seed costs are expressed by the product between, a given plant density reccomentation ${PD}_{j}$, where ${PD}_{j}$ can be either ${minRPD}_{j}$ or ${AOPD}_{j}$, and $6.16 * {10}^{-5}$ $ seed ⁻¹. Price per unit seed was calculated by assuming 40 $ seed bag ⁻¹, a bag weight of 22.7 kg and grain weight of 35 mg grain ⁻¹. $\Delta {NetIncome}_{j}$ represents the net income difference from reccomendations following ${minRPD}_{j}$ against ${AOPD}_{j}$.

Results

Informative priors: Canadian prairies

Yields in the Canadian Prairies averaged 4.42 Mg ha⁻¹ (q_0.025 = 1.35 Mg ha⁻¹; q_0.975 = 8.18 Mg ha⁻¹). Plant density ranged from 4 to 501 plants m⁻². Across all treatments and environments, the AOPD ranged from 79 to 128 plants m⁻². Meanwhile, minRPD was consistently higher than the AOPD (posterior probability = 0.0001) and ranged from 138 to 199 plants m⁻². Capturing the whole range of probable values for AOPD and minRPD parameters allowed to design informative priors reflecting the current knowledge and its associated uncertainty (Fig. 1).

Fig. 1

Posterior distributions of yield response to plant density model parameters (Eqs. 2 and 3) pooling all 12 combinations of seed treatment, seed size and environment. Derived quantities from these distributions (e.g., the expected value, the standard deviation, and quantiles) were used to build informative priors described in Eqs. 6–10. Purple bars refer to model parameters related to the expected value ${(\mu }_{ij})$ and yellow bars to model parameters related to the variance ${(\sigma }_{ij}^{2})$ of grain yield. Abbreviations: sample size (n), linear term of expected yield response to plant density (${\beta }_{1}$), agronomic optimum plant density (AOPD), maximum variance plant density (maxVPD), minimum risk plant density (minRPD), linear term of variance of yield response to plant density (${\alpha }_{1}$)

Yield response to plant density under varying seed size, seed treatment and environment in Kansas

Across all treatments and Kansas environments, grain yield averaged 4.10 Mg ha⁻¹ (q_.025 = 1.87 Mg ha⁻¹; q_.975 = 7.25 Mg ha⁻¹) and plant density ranged from 8 to 275 plants m⁻². In low- yielding environments grain yield was on average 3.37 Mg ha⁻¹ (q_.025 = 1.62 Mg ha⁻¹; q_.975 = 4.64 Mg ha⁻¹). In high yield environments grain yield averaged 5.59 Mg ha⁻¹ (q_.025 = 3.79 Mg ha⁻¹; q_.975 = 7.42 Mg ha⁻¹).

Plant density effects on expected yield and yield variance can be visually inspected (as a first check) for winter wheat seeded in Kansas US under untreated, light seeds treatments in low-yielding environments cluster (Fig. 2). All the other seed treatments × seed size × environments combinations can be found in Supplementary Fig. 2. Moreover, to compare how seed treatment, seed size and environment affected expected yield, yield variance and model parameters, we derived quantities from models presented in Fig. 2 and Supplementary Fig. 2 using Eq. 11 and explored in Sects. "Seed treatment and size impacts expected yield and yield variance" and "Seed treatment and size effects on AOPD and minRPD.".

Fig. 2

Panel A shows expected grain yield (${\mu }_{ij}$) in response to plant density (${s}_{ij}$, purple); panel B shows standard deviation of yield (derived from yield variance,${\sigma }_{ij}^{2}$) in response to plant density (yellow); panel C shows the posterior predictive distribution of grain yield (${y}_{ij}$) response to plant density (brown) for the jth treatment × environment combination (in this example untreated, light seeds in low yield environment cluster). Solid black lines represent the median and ribbons the 95% credible interval. Purple refers to parameters or quantities related to modelling the expected value ${(\mu }_{ij})$ and yellow to parameters or quantities related to modelling the variance ${(\sigma }_{ij}^{2})$ of grain yield

Seed treatment and size impacts expected yield and yield variance

Heavy and treated seeds conferred greater yields at the AOPD compared to the control seed management practice untreated light (Fig. 3). Across environments, seed treatment and heavy seed size resulted in similar yield advantage over the control (Fig. 3), averaging 0.27 Mg ha⁻¹ (q_.025 = 0.06 Mg ha⁻¹; q_.975 = 0.47 Mg ha⁻¹). Results from the rest of seed treatment and size combinations are shown in Supplementary Fig. 3. Three treatments (i.e., None_Moderate, None_Heavy, Treated_Moderate in Supplementary Fig. 3) resulted in similar yield advantage over the control under high (q_.50 = 0.238 Mg ha⁻¹; q_.025 = −0.026 Mg ha⁻¹; q_.975 = 0.498 Mg ha⁻¹) and low (q_.50 = 0.167 Mg ha⁻¹; q_.025 = 0.001 Mg ha⁻¹; q_.975 = 0.331 Mg ha⁻¹) yielding environments (Supplementary Fig. 3).

Fig. 3

A Posterior distributions of expected grain yield response to seed management evaluated at the AOPD (Eq. 11); B standard deviation of yield response to seed management evaluated at the minRPD (Eq. 11) in high (circles, opaque) and low (diamond, transparent) yielding environments. Contrasts represent heavy, treated seeds against light seed, untreated control treatment (Eq. 11). Points represent the median and whiskers are the 95% credible interval; p represents the posterior probability of $\Delta {\mu }_{j}$, $\Delta {\sigma }_{j}$ being greater than zero. Purple refers to derived quantities related to modelling the expected value ${(\mu }_{ij})$ and yellow to derived quantities related to modelling the variance ${(\sigma }_{ij}^{2})$ of grain yield

In low-yielding environments, seeding heavy or moderate size seeds and treated seeds resulted in lower yield variance (more stable yields) compared to the control (Fig. 3B, Supplementary Fig. 3), reducing yield variance on average by −0.13 Mg ha⁻¹ (q_.025 = 0.02 Mg ha⁻¹; q_.975 = −0.31 Mg ha⁻¹).

Seed treatment and size effects on AOPD and minRPD

Regardless of the environment, treated, heavy seeds resulted in similar AOPD as compared to the untreated, light seeds control (Figs. 4, Fig. 5). All other seed treatments and sizes resulted in similar AOPD compared to the control (Supplementary Fig. 5). Seed treatment influenced the minRPD, particularly in low-yielding environments where two treatments [Treated_Moderate (posterior probability = 0.880) and Treated_Heavy (posterior probability = 0.863)] had a higher minRPD compared to the control untreated, light seeds (posterior probability = 0.863; Fig. 5B, Supplementary Fig. 5).

Fig. 4

Posterior distributions of AOPD (A) and minRPD (B) for contrasting seed management treatments (no treatment and light seeds, “None_Light”; treated and heavy seeds, “Treated_Heavy”) under high (circles, opaque) and low (diamond, transparent) yielding environments. Points represent the median and whiskers are the 95% credible interval. Purple refers to model parameters related to the expected value ${(\mu }_{ij})$ and yellow to model parameters related to the variance ${(\sigma }_{ij}^{2})$ of grain yield

Fig. 5

Posterior distributions of AOPD (A) and minRPD (B) as affected by seed management treatments under high (circles, opaque) and low (diamond, transparent) yielding environments. Contrasts represent heavy, treated seeds against light, untreated (control) seeds. Points represent the median and whiskers are the 95% credible interval. Equation 11 describes the calculations used for contrasts between treatments. Purple refers to derived quantities related to modelling the expected value ${(\mu }_{ij})$ and yellow to derived quantities related to modelling the variance ${(\sigma }_{ij}^{2})$ of grain yield. Abbreviations: agronomic optimum plant density (AOPD), minimum risk plant density (minRPD)

Economic assessment of plant densiy recommendations

When comparing net income resulting from plant density recommendations following minRPD against AOPD, across seed chemical treatments and cleaning methods net income was on average 37 $ ha⁻¹ higher (q_.025 = 6 $ ha⁻¹; q_.975 = 77 $ ha⁻¹) in low yielding environment (Fig. 6A). In high yielding environments net income was 19 $ ha⁻¹ higher (q_.025 = −8 $ ha⁻¹; q_.975 = 54 $ ha⁻¹) when following plant density recommendations based on the minRPD compared to the AOPD (Fig. 6B).

Fig. 6

Posterior distribution of the difference in net income from plant density reccomendations following minRPD compared to AOPD, under (A) low and (B) high yielding environments and varying seed chemical treatment and cleaning methods. Positive ΔNet Income values represent greater net income by following minRPD compared to AOPD (Eq. 12, 13). Points represent the mean and whiskers the 95% credible interval

Discussion

We introduced a model describing both expected yield and yield variance in response to plant density using comprehensive datasets from two hard red winter wheat growing areas. In addition, assessed how seed management affected the optimal planting density. We determined the optimum planting density that maximizes yield and minimizes variance (minRPD) and showed how our framework could be used to compare yield-plant density relations between a combination of treatments (seed size and seed treatment) and environments. The novelty on the proposed framework lies on the assumptions behind the variance of yield response to plant density, which allows quantifying the minimum plant density that minimizes grain yield variance.

Limitations

One limitation of this study is the use of k-means clustering algorithm to perform environmental classification. First, k-means does not account for uncertainty associated with yield environment classifications [1, 19]. Second, k-means classification approach requires post hoc yield data to classify environments, restricting its applicability in real world scenarios where growers have no knowledge of the conditions for the upcoming season. A feasible alternative to overcome the latter limitation is to perform long-term assessment of drought patterns and yield environments in the region (e.g., [8, 9]), subsequently adopting seeding rate management strategies that maximize yield and minimize variance for the most common environment in a region.

Another limitation of this study is that we evaluated a single cultivar, thus ignoring cultivar × plant density interactions that can occur in wheat [2, 15, 41, 44]. Previous studies showed the optimum plant density depended on the cultivar’s tillering potential [4, 44], yield of a cultivar at very low planting density (i.e.: plant yield efficiency; [41]), and to lesser extent differences in leaf insertion angle (i.e.: planophile or erectophile; [16]). Future research can focus on estimating the minRPD for cultivars varying tillering potential or other physiological characteristics.

Modeling assumptions

In our model, yield response to plant density is defined by models for two processes (one for expected yield and one for yield variance, Eq. 2. and Eq. 3.) and are based on mechanisms of wheat physiology [46]. From the expected yield perspective, we acknowledge multiple processes were proposed in the literature to model the yield-plant density relation [4, 26, 28, 30, 39, 48]. Critically thinking about the underlying biological basis behind the tested hypothesis (i.e., the model selected) is crucial to provide meaningful inference [24]). Thus, our designed process model relies on parameters and derived quantities that are agronomically meaningful and interpretable to guide producers in their plant density decisions (e.g., minRPD and AOPD). Contrastingly, model comparison and selection approaches according to the empirical support of a model (e.g., selecting a model based on the lowest Akaike or Deviance Information Criterion, etc.) should only be utilized when the data shows support for more than one hypothesis [24].

Agronomic and economic implications behind minRPD

Our research shows minRPD was always greater than the AOPD (Fig. 4). This phenomenon could be partially explained by plant density × environment interactions within each environmental cluster (i.e., low- or high- yielding environment), resulting in each environment having its own AOPD, as AOPD depends on the target environment [4]. Environmental factors influencing the optimal plant density (i.e.: AOPD), remain to be tested. Improved predictive performance of plant density recommendation models can seek to incorporate static (i.e., those with lower variability within the same cropping season such as soil organic matter, and soil water holding capacity) environmental covariates. While dynamic environmental factors (those that have greater variability within the same cropping season) might be prone to better explain changes in yield response to plant population (e.g., fall precipitation, average temperature, and cumulative radiation during the Fall and Winter [4]), their inherent uncertainty limits their use on predictive models.

There are multiple issues associated with modeling grain yield response to plant density for each individual environment. First, the range of plant density is constrained by densities explored in a single environment. Second, subgrouping sample size is limited to observations within a single environment resulting in less precise estimators (i.e., greater uncertainty on model parameters). Finally, it has no predictive power as environmental conditions resulting in a given AOPD are largely unpredictable.

Heavy seed sizes and treated seeds resulted in more stable yields compared with light untreated seeds only under low yielding environments (Fig. 3B). Evidence in the Canadian Praries suggest similar behavior, where the benefits of seed size and chemical treatment were more apparent under conditions conducive to thinner wheat stands [6]. Mian and Nafziger [32] found that the use of smaller seed sizes can translate into thinner wheat stands and even provide a yield advantage over heavier seed sizes under G x E combinations conducive to lodging. However, genetic progess in lodging resistance makes lodging conditions unlikely in more recently released wheat germplasm [49]. A probabilistic assessment of genotypes with varying tillering potential, and planting dates can contribute to classify G x E combinations where intensified seed management strategies can provide yield advantages.

Guiding plant density reccomendations following minRPD comprise higher seed costs when compared against AOPD. However, our data shows that more stable yields resulting from adopting minRPD translate in higher net income compared to those following AOPD. Economic benefits from minRPD are higher in low yielding environments, as plant density modulates spikes per m⁻², a coarse regulator of grain yield [38]. An independent dataset shows empirical evidence from field trials conducted in the Central Great Plains region showing that winter wheat yield response to plant density is higher as environmental mean yield decreses [4], aligning with our findings. Future models can consider incorporating uncertainty in seed to grain price ratio as a component within the model [10], which would allow to quantify the impacts of plant density following minRPD on long term net returns. Therefore, delivering plant density recommendations following minRPD can mitigate, at the expense of higher seeding costs, risks associated with plant density decisions in unpredictable environments.

Conclusion

We described both expected yield and yield variance in response to plant density. Our model consisted of process models generated with our expert knowledge on eco-physiological processes and the prior distributions were specified using previous data and past studies. We used this model to answer agronomically relevant questions regarding seed management (e.g., yield and variance of yield response to seed size and treatment; AOPD and minRPD response to seed size and treatment). To facilitate the usability of our model, we developed a user interface (https://ksuwheat.shinyapps.io/minrpd) which allows informing seed management decisions based on probabilistic inference of the minimum risk plant density (minRPD). Under heteroscedastic data scenarios the proposed framework can be used to model grain yield response to plant density. Future research can (i) explore how genotypes, environments and their interaction modulate the difference between AOPD and minRPD and, (ii) extend the framework to a variety of processes involving crop management decisions.