A Bayesian approach to assessing uncertainty in the effect of fertilization strategies on paddy rice yield via multiple on-farm experiments

ABSTRACT

The surge in global agricultural input prices, particularly for nitrogen (N) fertilizers, has presented challenges for farmers facing stagnant crop prices and rising production costs. This study explored a system designed to sustain or enhance farmers’ income while reducing production costs, with a specific focus on implementing practices that involve cultivating crops with lower rates of fertilizer than conventional rates, particularly in smaller land settings, through on-farm experimentation (OFE). OFE was conducted over three years across 11 fields in Gifu, Japan, utilizing a strip trial design to investigate the effect of fertilizer input rates on rice yield. To collect yield data, unmanned aerial vehicle (UAV)-mounted multispectral cameras were deployed. A Gaussian process model was used to address the impact of fertilizer input rates on yield while accounting for spatial yield variation. The study further assessed the number of field trials achieving positive profits at certain probability levels by integrating outcomes from multiple OFEs. The results highlighted little rice yield response to fertilizer input in most of the fields, while several fields showed significant positive yield responses. Notably, the use of lower fertilizer input rates in certain fields proved to be more profitable whereas higher fertilizer input rates were less profitable in most fields. The results further emphasized the economic viability of lower rates, especially when fertilizer prices are high, thus contributing valuable insights for sustainable and profitable rice production in Japanese cultivation.

GRAPHICAL ABSTRACT

KEYWORDS:

• Reduced fertilizer rate
• small field
• strip trial
• UAV
• revenue analysis

Introduction

The escalation of global agricultural input prices, particularly for nitrogen (N) fertilizers (Hellegers, 2022; Rabbi et al., 2023; Randive et al., 2021), has posed challenges for farmers facing stagnant crop prices and rising production costs (Weerahewa & Dayananda, 2023). Furthermore, there is an urgent need to reduce N fertilizer inputs to mitigate ground water pollution (Bijay-Singh & Craswell, 2021; Gao et al., 2012; Ju et al., 2006; Zhang et al., 1996) and global warming (Koneswaran & Nierenberg, 2008; Tian et al., 2012). Overcoming these challenges involves promoting practices where farmers can grow crops with less N fertilizer while adopting sustainable cropping practices. This approach seeks to balance economic viability with ecological sustainability, offering benefits for farmers, the broader agricultural sector, and the natural environment.

To facilitate an understanding of crop physiological responses to crop management practices, i.e. N fertilizer input rates, agricultural experimentation is traditionally conducted at research centers through randomized and replicated small-plot designs (Piepho et al., 2011). However, the inadequacy of small-plot experiments arises from the limitations in capturing the spatial variability inherent in large and diverse fields (Sadeghi et al., 2006; Stefanova et al., 2023). Small-plot experiments might oversimplify this complexity, leading to results that are not representative of broader field conditions (Eltaib et al., 2002). To obtain a more accurate and nuanced understanding of optimal treatments across fields, the on-farm experimentation (OFE) method has the potential to expedite the implementation of this system (Evans et al., 2020; P. B. Hegedus, Maxwell, et al., 2023; Roques et al., 2022).

OFE is a collaborative and innovative approach in agricultural research involving a dynamic partnership between farmers and professional researchers (Lacoste et al., 2022; Marchant et al., 2019). In OFEs, farmers actively participate in agronomic trials conducted directly in their fields, providing valuable real-world insights. OFE is designed to produce scientifically significant outcomes while remaining highly relevant to local farming contexts, acknowledging the diversity of environmental and management conditions (Lacoste et al., 2022). Collaboration ensures that experiments reflect the details of actual farming practices, enhancing the translational impact of research findings. This active engagement promotes a two-way exchange, with farmers contributing practical knowledge and researchers providing technical expertise. The holistic and inclusive nature of OFE promotes a deeper understanding of the local context and empowers farmers in decision-making, resulting in more effective and contextually relevant agricultural strategies.

Strip trials are among the popular experimental designs for OFE. This design involves dividing a field into strips, with each strip representing a distinct experimental unit where a specific treatment is applied based on the working width of agricultural machinery (e.g. seeder, spreader, and sprayer). Field-length strips serve as experimental units, and multiple replications of treatment and control strips are arranged within a field. Although designs such as chequerboard and eggboxes can achieve a finer treatment scale by using variable rate technology (Cook et al., 2018), strip trials remain more cost effective (Li et al., 2023) and easily implementable for farmers. Furthermore, the emergence and widespread adoption of yield monitoring systems and remote sensing technologies allows for data collection within strips and reveals spatial variability caused by diverse factors (Istiak Md et al., 2023; Javaid et al., 2023; Marchant et al., 2019; Roques et al., 2022). However, the traditional statistical method such as ordinary least squares regression or analysis of variance is not suitable for data analysis in OFEs because it may not effectively address the inherent spatial variability in specific response variables (Tanaka, 2021). Several studies have explored appropriate data analytics methods for strip trials via various statistical models, considering the large spatial dataset of yields per strip rather than a single value for the entire strip (Alesso et al., 2021; Marchant et al., 2019; Tanaka, 2021). Those studies demonstrated that a linear mixed model accounting for spatial variations in residuals could be used to effectively assess treatment effects, even in straightforward OFE designs, namely, strip trials.

OFEs are commonly found on expansive field sizes (Marchant et al., 2019) and typically exceed 10 hectares. The advantage of larger land sizes lies in their capacity to facilitate a more comprehensive evaluation of diverse agricultural factors, providing a broader and more representative understanding of real-world conditions. With the global adoption of precision agricultural technologies, even small-scale or smallholder farmers can reap the benefits of OFEs (Tanaka, 2021). However, designing experiments for small-scale fields requires consideration of unique challenges distinct from those encountered in large-scale fields (Tanaka, 2021). The principal hurdle faced by small-scale farming systems is the potential constraint on the economic benefits of OFEs due to their limited scale, which hinders their expansion via farmer networks. Economies of scale become pivotal; larger field sizes yield more substantial benefits from a single OFE. However, smaller field sizes may compromise the generalizability of OFE outcomes for the broader farming community. Aggregating and sharing data from multiple OFEs can overcome this challenge and foster a collaborative and integrated approach. This collaborative effort ensures that recommendations for small-scale farmers are robust, inclusive, and reflective of the diverse conditions encountered in small-scale farming systems, thus mitigating the challenges posed by individual OFE scales. The incorporation of multiple experimental sites offers a significant advantage in research and analysis, particularly in the context of OFEs, thus facilitating a heightened level of precision in cross-site analyses and allowing for the discernment of intricate patterns, trends, and relationships across diverse environments and conditions (Roques et al., 2022).

One approach for analyzing the data from multiple experimental sites is the Bayesian approach. The Bayesian approach is preferred in various fields, including medical (Achim et al., 2001; Natanegara et al., 2014; Spiegelhalter et al., 1999; Sutton & Abrams, 2001), ecological (Carl & Kühn, 2007; Choy et al., 2009; Ellison, 1996; Hooten & Hobbs, 2015), and even agricultural (Donald et al., 2011; P. B. Hegedus, B. D. Maxwell, et al., 2023; Shirley et al., 2020; Wang et al., 2012) fields, due to its robust and flexible analysis system that effectively addresses uncertainty. Unlike traditional frequentist statistics, Bayesian analysis incorporates prior knowledge, expressed as a probability distribution, and updates it with observed data to obtain a posterior distribution. In the context of spatial modeling, Bayesian analysis is particularly advantageous because it can also account for spatial autocorrelation (Besag & Higdon, 1999; Dendoncker et al., 2007; Donald et al., 2011; Morrison et al., 2012; Paccioretti et al., 2021; Wang et al., 2012). Bayesian data analysis framework can explicitly use the posterior probability for quantifying uncertainty in inference (Gelman et al., 2013). Therefore, Bayesian analysis enables decision-makers to quantify the treatment effect and its uncertainty for each OFE. Furthermore, pooling statistical information from multiple OFEs was proposed to facilitate an understanding of the overall treatment effect and its uncertainty synthesized from multiple OFEs (Kyveryga et al., 2018). While Bayesian analysis and pooling approach has great potential, to the best of our knowledge, few studies have focused on analyzing multiple OFEs via Bayesian analysis while accounting for spatial autocorrelation of each OFE.

The primary objective of this study was to develop a Bayesian data analytical framework to integrate outcomes from multiple strip trials to provide recommendations for fertilizer input management in small-scale paddy rice production systems in Japan. To date, such systems have not been adequately studied, nor have their unique circumstances been addressed. The specific goals of this study were to 1) assess the effect of fertilizer application rates on rice yield for each field trial via a Gaussian process model, 2) comprehensively explore the effect of fertilizer application rates and the associated probability of profit gain from multiple field trials, and 3) provide insights into fertilization strategies that can balance crop productivity and financial sustainability. Therefore, on-farm strip trials with different fertilizer input rates were conducted over three years across 11 paddy fields in Gifu, Japan. The yield data were collected using a UAV-based yield prediction model. The question tested in this study was whether reducing or increasing the fertilizer input rate from the conventional rate would result in enhancing farmers’ economic returns.

Materials and methods

Experimental design of OFE

OFEs were conducted over a three-year period across 11 farmers’ fields in Gifu Prefecture, Japan, from 2020 to 2022 (as detailed in Table 1). The investigation focused on the high-yielding rice cultivar Hoshijirushi. This cultivar is tolerant of lodging even at a high N input rate; therefore, farmers tend to apply overfertilization to achieve higher yields. However, several farmers had the question of whether it was still profitable to increase or decrease fertilization under local circumstances in which the fertilizer price continued to rise while the rice price rarely changed. This socioeconomic situation facilitated the recruitment of farmers who were willing to become involved in OFEs with the aid of a local crop advisory service. The experiment was carried out across various locations involving multiple farmers who utilized different transplanting machineries and types of fertilizer. Furthermore, each farmer has different experiences and strategies for crop management. These factors presented challenges in maintaining consistent experimental conditions across the fields. Through the discussion with the farmers and local crop advisory service, three different basal fertilizer input rates were primarily applied to each farmers’ field by increasing or decreasing the basal fertilizer input from farmers’ conventional rates. The fertilizer input rates ranged from 300 to 500 kg ha⁻¹, resulting in an increase or reduction of 10–14% compared with conventional rates. However, some farmers also experimented with much lower or higher application rates to test various hypotheses. For instance, in fields 5 to 7 and 8 to 11, other fertilizer treatments of 0 kg ha⁻¹ and 250 kg ha⁻¹, respectively, were additionally applied, as depicted in Figure 1. For this study, variations in the experimental treatments of fertilizer rates were only considered within the range of 10–14% of the conventional rates due to the limitation of complete dataset; therefore, rates either lower (<10%) or higher (>14%) than the specified ranges were excluded from further analysis. To evaluate the effect of different basal fertilizer input rates on rice yield, a series of strip trials were systematically performed (Figure 1). The trials involved the arrangement of long strips alongside a field. Each of these strips received different rates of fertilizer application, a difference determined based on the working width of the rice transplanter, which was approximately 1.9 m. Both types of slow-releasing basal fertilizer (NPK = 25:6:6 or 24:9:9) were used by the local farmers. There was no top-dressing application for any of the fields, and weed and pest management were performed according to the local crop advisory recommendation.

Figure 1. Yield and treatment maps of strip trials. Within the depicted strips, the numerical values denote the corresponding fertilizer application rates. The gradient colors represent the distribution of the predicted yield.

Table 1. Basic information of the OFE fields.

Field ID	Latitude	Longitude	Transplanting date	UAV imagery acquisition date	Fertilizer application rate, kg ha^−1			The number of yield observations
					Conventional rate	Lower rate	Higher rate
1	35°13′N	136°40′E	11 May 2020	12 August 2020	500	450	550	2,360
2	35°13′N	136°40′E	11 May 2020	6 August 2020	500	450	550	1,984
3	35°16′N	136°35′E	11 May 2020	6 August 2020	400	350	450	2,307
4	35°16′N	136°35′E	23 May 2020	6 August 2020	400	350	450	2,832
5	35°16′N	136°35′E	25 April 2021	13 July 2021	350	300	400	532
6	35°16′N	136°35′E	19 April 2021	13 July 2021	350	300	400	356
7	35°15′N	136°36′E	14 May 2021	11 August 2021	350	300	400	877
8	35°13′N	136°40′E	11 May 2021	11 August 2021	350	300	400	834
9	35°15′N	136°36′E	16 May 2022	9 August 2022	350	300	400	515
10	35°16′N	136°35′E	23 May 2022	9 August 2022	350	300	400	334
11	35°16′N	136°35′E	2 May 2022	9 August 2022	350	300	400	263

Data collection

A yield map was generated according to the method of Mia et al. (2023). Briefly, a total of 894 samples was surveyed across 22 farmers’ fields in Japan over six years (2017–2022) to develop globally applicable rice yield prediction models. The fields of interest were surveyed using a commercial multispectral camera (Rededge Altum, MicaSense, Seattle, USA) mounted on an unmanned aerial vehicle (UAV) (Inspire 2, DJI, Shenzhen, China) during the rice flowering stage. This aerial imagery, coupled with ground truth yield data and weekly summations of weather data, facilitated the development of a yield prediction model, specifically a convolutional neural network (CNN) model. Mia et al. (2023) compared rice yield prediction accuracy among 18 CNN models, and the best model (Model 8 of their work) was utilized in the present study. The model consists of AlexNet feature extractor with weekly weather data, and one additional hidden layer with 64 neurons between concatenation and output layers. The performance of the developed CNN model exhibited a R² value of 0.68 and a root mean square error of 0.859 t ha⁻¹ (Mia et al., 2023). The prediction model operated on a 1.0 m × 1.0 m unit, aligning with the input data resolution. To visualize the predicted yield values, maps were generated using QGIS 3.22.4, wherein gradient-colored points represented centroids of 1.0 m × 1.0 m prediction units (Figure 1). Subsequently, the predicted yield data were stratified based on the polygons of the treatment strips.

Gaussian process model

To assess the effect of fertilizer application rates on rice yield, a Gaussian process model was established for each field trial as follows:

(1) $\mathit{y}\sim \mathrm{multvariate}\mathrm{normal}\left(\mathit{\beta }\left(\mathit{X}\right),\mathit{K}\right)$(1)

where y is a vector of length n of the response variable rice yield, n is the number of observations, β(X) is an n vector of the mean function of input design matrix X, and $\mathit{K}$ is an n × n covariance matrix. The mean function β consists of β₀, β₁, and β₂ for the mean yield value at the farmers’ conventional fertilizer application rate, treatment effect of low fertilizer application rate, and treatment effect of high fertilizer application rate, respectively. The design matrix X consists of 3 × n elements. The first column of design matrix X consists of ones, representing the intercept term (β₀) being multiplied with each observation. When low fertilizer application rate is applied, the second column has a value of 1 for multiplication with β₁ while the third column has a value of 0 for multiplication with β₂. Conversely, when high fertilizer application rate is applied, the second column has a value of 0 while the third column has a value of 1. Therefore, β₁ and β₂ represent the change in yield value relative to β₀. Note that the spatial unit of each observation is defined as a 1.0 m × 1.0 m area to correspond with the UAV-based yield prediction unit. Accordingly, the number of observations (n) for each field varied from 263 to 2,832 according to the factors such as field size and misalignment between the UAV-based yield spatial resolution and working width of transplanters (Table 1).

The covariance function is also known as the kernel of the Gaussian process model. The exponential kernel function k between ith and jth observations is defined as follows:

(2) $\mathit{k}\left(d_{\mathit{i},\mathit{j}}|\mathit{\alpha },\mathit{\rho },\mathit{\sigma }\right)={\alpha }^2\exp \left(-\frac{1}{2}{\left(\frac{{\mathit{d}}_{\mathit{i},\mathit{j}}}{\mathit{\rho }}\right)}^2\right)+{\delta }_{\mathit{i},\mathit{j}}{\sigma }^2$(2)

where the hyperparameter α is the marginal standard deviation, the hyperparameter ρ is the length scale, d_i,j is the Euclidean distance between ith and jth observations, the hyperparameter σ is the scale of the noise, and ${\delta }_{\mathit{i},\mathit{j}}$ is the Kronecker delta function representing a value of 1 when $\mathit{i}=\mathit{j}$ otherwise the value is 0. The Euclidean distance d_i,j is calculated as follows:

(3) $d_{\mathit{i},\mathit{j}}=\sqrt{\left({\left(x_i-x_j\right)}^2+{\left(y_i-y_j\right)}^2\right)}$(3)

where x_i and x_j are the Universal Transverse Mercator (UTM) x coordinates of ith and jth observations, and y_i and y_j are the UTM y coordinates of ith and jth observations. These coordinates are based on the centroids of 1.0 m × 1.0 m yield prediction units.

The nugget-to-sill ratio in geostatistics quantifies the relative influence of short-distance variability or measurement errors (nugget) to the overall spatial variability (sill), offering insights into the strength and scale of spatial dependence in the dataset. The hyperparameters α² and σ² are equivalent to the partial sill and nugget; therefore, the nugget-to-sill ratio is evaluated. A nugget-to-sill ratio close to 0 indicates that the spatial variation is mostly structured, with a small unstructured (nugget) component. This finding suggested strong spatial autocorrelation. A nugget-to-sill ratio close to 1 indicates that the dataset has a substantial unstructured or nugget component, implying weak spatial autocorrelation. In this case, spatial observations are not strongly dependent on distance. According to the posterior distribution of the nugget-to-sill ratio, the results of the Gaussian process model can be discussed.

The Gaussian process model was established using the ‘Rstan’ package (McElreath, 2018; Stan Development Team, 2022) in R version 4.3.1 (R Core Team, 2024). Each model underwent four Markov chain Monte Carlo chains comprising 1,000 iterations each, with a 50% burn-in. This process yielded a total of 2,000 iterations for posterior inference. The convergence and stability of these chains were inspected through trace plots and R-hat values during the monitoring process (Gelman & Rubin, 1992; Gelman et al., 2013; McElreath, 2018). The trace plots visually displayed the evolution of each chain over iterations, while R-hat values served as a diagnostic measure to assess the convergence of multiple chains. A value of 1 for R-hat implies favorable convergence, while larger values indicate potential convergence issues, ill-defined prior distributions and predictors or the potential need for additional iterations. The significance of the estimated parameters was assessed at the 95% highest density interval (HDI). The source code for implementing the gaussian process regressions can be found in GitHub (https://github.com/takashi2024/StanGP_RiceStripTrial).

Net revenue calculation

Leveraging the prescribed Equations (1) and (2), the field-specific net revenue ($ ha⁻¹) associated with the yield response to varying fertilizer rates was calculated as follows:

(4) $r_{1\mathit{i}}=p_g\times \mathit{C}\times {\beta }_{1\mathit{i}}-p_f\times \left(F_c-F_l\right)$(4)

(5) $r_{2\mathit{i}}=p_g\times \mathit{C}\times {\beta }_{2\mathit{i}}-p_f\times \left(F_h-F_c\right)$(5)

where $r_{1\mathit{i}}$ and $r_{2\mathit{i}}$ represent the revenues derived from employing a lower and higher fertilizer rate, respectively, than the farmer’s rate at Bayesian iteration number i (i.e. 1, 2, … , 2000); $p_g$ is the price of fully ripened brown rice grain; $p_f$ is the price of basal fertilizer; C = 0.80 is the factor converting from rough rice to fully ripened brown rice grain; F_c is the conventional fertilizer rate used by the farmer; $F_l$ is the lower fertilizer rate; and $F_h$ is the higher fertilizer rate ($F_c$). These parameters signify the alterations in yield resulting from the adjustment in fertilizer application rates at each specific field. The prices of grain and fertilizer exhibited annual variation. Specifically, the grain prices $p_g$ in 2020, 2021, 2022, and 2023 were $1.343 kg⁻¹, $0.997 kg⁻¹, $1.076 kg⁻¹, and $1.187 kg^–1, respectively. The prices of fertilizer $p_f$ in 2020, 2021, 2022, and 2023 were $1.200 kg⁻¹, $1.172 kg⁻¹, $1.362 kg^–1, and $1.970 kg⁻¹, respectively.

Comprehensive probability map for fertilization recommendation

To holistically understand the impact of fertilizer rates on profitability while accounting for variability in fertilizer and crop prices, the number of field trials that achieved positive net revenue associated with different fertilizer rates at certain probability levels was evaluated. The analytical results should be more credible if the majority of field trials consistently demonstrate positive net revenue at higher probability levels. In the previous section, constant fertilizer and grain prices were assumed for each year. However, it’s crucial to note that this assumption is strongly post-hoc and not practical. By leveraging the variability in grain and fertilizer prices, comprehensive probability for revenue gain was generated. First, the probability of gaining positive net revenues was obtained for each field by dividing the total number of positive net revenues by the total number of Bayesian iterations for each of the lower and higher fertilizer input rates based on Equations (4)(4) $r_{1\mathit{i}}=p_g\times \mathit{C}\times {\beta }_{1\mathit{i}}-p_f\times \left(F_c-F_l\right)$(4) and (5(5) $r_{2\mathit{i}}=p_g\times \mathit{C}\times {\beta }_{2\mathit{i}}-p_f\times \left(F_h-F_c\right)$(5) ). Given the uncertainty in fertilizer and grain prices, an expansive range of grain prices $p_g$ ranging from $0.883 kg⁻¹ to $1.450 kg⁻¹ and fertilizer prices $p_f$ ranging from $1.00 kg⁻¹ to $2.35 kg⁻¹ was employed for Equations (4)(4) $r_{1\mathit{i}}=p_g\times \mathit{C}\times {\beta }_{1\mathit{i}}-p_f\times \left(F_c-F_l\right)$(4) and (5(5) $r_{2\mathit{i}}=p_g\times \mathit{C}\times {\beta }_{2\mathit{i}}-p_f\times \left(F_h-F_c\right)$(5) ). Second, the number of field trials exceeding the probability levels of 95, 90, and 80% were counted for each combination of fertilizer and grain price. The 95% probability is scientifically sound; however, it could be extremely strict for practical production level (Kyveryga et al., 2018). Thus, the probability levels were also considered at the 90 and 80% level. By doing so, 11 on-farm experimental outcomes could be comprehensively analyzed. The results were visualized using heatmaps, which provide overview of the likelihood of achieving positive net revenue under various market conditions, aiding in decision-making processes regarding fertilization in real agricultural practices. Note that this study simply presented the number of field trials achieving the positive net revenue rather than evaluating the frequency by dividing the total number of trials. This way could be more straightforward to convey the analytical concept for farmers without losing the context that ‘probability’ derived from each field trial.

Results

Spatial variation

Figure 2 shows a visual representation of the density diagram for the nugget-to-sill ratio. The majority of posterior values in most fields ranged less than 0.50. Notably, in four of eleven fields (Field 3, 4, 5, and 6), all posterior values were consistently less than 0.25, indicating robust and widespread residual spatial variation among these fields. In the context of Gaussian process modeling, which is utilized for capturing and predicting spatial dependencies, this observation underscores the suitability of applying the Gaussian process model to these fields.

Figure 2. Density plots of the nugget-to-sill ratio for 11 fields. Light green areas represent the 95% HDI. Vertical dashed lines represent the criteria classifying strong (0.00–0.25), medium (0.25–0.75), and weak autocorrelation (0.75–1.00).

Effects of fertilizer application rates on rice yield

The effects of different fertilizer rates on yield are summarized in Table 2. The results showed that all the models properly converged, as evidenced by R-hat values less than 1.10. In Fields 2 and 10, higher fertilizer application rates significantly increased yield. Moreover, in Fields 2, 4, 7, 8, 10, and 11, lower fertilizer application rates significantly decreased yield. In contrast, Fields 1 and 3 showed a significant increase in yields even with lower fertilizer application. On the other hand, Fields 5, 6, and 9 did not exhibit any significant changes in yield when the fertilizer rate was increased or decreased from the conventional application rate.

Table 2. Summary of posterior distributions of the estimated treatment impacts of fertilizer application rates on yield. The bold numbers in the columns corresponding to lower and higher fertilizer application rates indicate that the 95% HDI does not include 0.

Field ID	Farmers conventional fertilizer rate				Lower fertilizer rate than conventional rate				Higher fertilizer rate than conventional rate
	β0	2.5%	97.5%	R-hat	β1	2.5%	97.5%	R-hat	β2	2.5%	97.5%	R-hat
	t ha^−1	t ha^−1	t ha^−1		t ha^−1	t ha^−1	t ha^−1		t ha^−1	t ha^−1	t ha^−1
1	8.990	8.946	9.040	1.000	0.075	0.029	0.121	1.002	0.042	−0.004	0.089	1.000
2	9.007	8.948	9.059	1.002	−0.048	−0.076	−0.017	1.001	0.045	0.015	0.073	1.001
3	8.595	8.537	8.653	0.999	0.054	0.001	0.111	1.001	−0.024	−0.082	0.039	1.000
4	8.597	8.502	8.693	1.000	−0.106	−0.210	−0.002	1.001	0.047	−0.051	0.138	0.999
5	6.602	6.479	6.718	1.004	0.038	−0.090	0.175	1.005	0.109	−0.018	0.241	1.009
6	6.356	6.242	6.462	1.002	0.010	−0.148	0.173	1.004	0.111	−0.023	0.244	1.001
7	6.050	5.966	6.125	1.007	−0.148	−0.245	−0.043	1.008	−0.032	−0.127	0.066	1.008
8	8.646	8.404	8.782	1.006	−0.069	−0.099	−0.041	1.000	−0.002	−0.033	0.030	1.005
9	7.711	7.487	7.845	1.003	0.003	−0.015	0.021	1.000	0.001	−0.016	0.019	1.001
10	7.670	7.616	7.720	1.002	−0.054	−0.090	−0.020	1.000	0.044	0.008	0.077	1.000
11	7.441	7.373	7.498	0.999	−0.048	−0.086	−0.008	1.004	0.025	−0.024	0.072	1.002

Revenue analysis

Figure 3 shows the estimated net revenues when utilizing lower fertilizer application rates than when using the conventional rate. The visual representation elucidates several key observations. Notably, Fields 1, 3, and 9 consistently exhibited positive net revenue, emphasizing the robust and reliable financial success associated with the implementation of a lower fertilizer application rate in these fields. Even in the absence of a significant yield change in Field 9 due to alterations in fertilizer application, the recurrence of positive revenue underscores the resilience and effectiveness of this approach. Fields 2, 5, 6, 8, 10, and 11 had mostly positive net revenue, suggesting the viability of employing a lower fertilizer rate in these fields for consistently favorable financial outcomes. However, Fields 4 and 7 presented a less uniform pattern, with positive revenue observed in only a few cases, suggesting a less consistent financial benefit from the lower fertilizer rates in these specific fields.

Figure 3. Density plots of net revenue as determined by using a lower fertilizer rate than the farmer rate. Light green areas represent the 95% HDI. Year-specific grain and fertilizer prices were used for the revenue evaluation.

The estimated net revenue when utilizing higher fertilizer application rates than when using conventional rates is shown in Figure 4. Fields 8 and 9 consistently had negative net revenue, indicating a definite risk of losing profits by increasing fertilizer application rates. In fields 3, 7, 10, and 11, the probability of positive net revenue was less than 5%, further emphasizing the minimal success of employing a higher fertilizer rate for obtaining profits. Fields 1, 2, 4, and 5 showed that net revenue could be either positive or negative at the 50% HDI level. Only Field 6 had positive net revenue at more than 50%.

Figure 4. Density plots of net revenue obtained by using a higher fertilizer rate than the farmer rate. Light green areas represent the 95% HDI. Year-specific grain and fertilizer prices were used for the revenue evaluation.

Probability analysis for fertilizer recommendation

Figure 5 offers a comprehensive view of how different fertilizer application strategies can influence the probability of generating positive net revenue across field trials. In the case of the lower fertilizer rate at the 95% probability, three of eleven field trials showed positive net revenues when fertilizer price was the lowest and grain price was the highest setting. This number remained consistent even when grain price rose to the maximum value. Meanwhile, ten of eleven fields were profitable when fertilizer price was the highest and grain price was the lowest setting. Eight of eleven field trials remained profitable even when both prices were the highest settings. None of the field trials were profitable at the 95% probability with the higher fertilizer rates. At the 90 and 80% probabilities, the overall patterns of heatmaps closely resembled those at the 95% probability level. However, as the probabilities decreased, the number of field trials achieving positive net revenue increased. In 2023, fertilizer price was high and grain price was moderate, which resulted in eight of eleven field trials achieving positive net revenue at the 95% probability and nine of eleven field trials at the 90 and 80% probabilities. Intriguingly, when faced with a lower fertilizer price and higher grain price, lower fertilizer rates should be recommended consistently. Conversely, Figure 5 illustrates the challenges associated with higher fertilizer use. Throughout all years examined, no field trials demonstrated positive net revenue at the 95, 90, and 80% probabilities, underscoring the financial risks associated with higher fertilizer usage.

Figure 5. Heatmaps showing the number of field trials that achieved positive revenue gain at the 95, 90, and 80% probability levels in the case of lower fertilizer rates and higher fertilizer rates in comparison to the conventional rate used by farmers. The x-axis represents fertilizer price, while the y-axis represents grain price. Gradient color corresponds to the number of field trials. The points represent the actual fertilizer and grain prices for the corresponding years.

Discussion

The occurrence of spatial variation in crop yield posed a challenge in analyzing OFE data. A Gaussian process model was successfully applied to assess the treatment effect of fertilizer application rates for each trial while accounting for the inherent spatial yield variations in the experimental data, as evidenced by the fact that all the models converged (Table 2). A high spatial dependency of rice yield was also confirmed across the fields (Figure 2), which suggested the necessity of a spatial modeling approach for analyzing strip trials even in rice production systems. This finding is consistent with previous OFE studies on other cereal crops (Marchant et al., 2019; Tanaka, 2021). Accordingly, the effect of fertilizer application rate on rice yield could be assessed more precisely via this type of approach than via traditional statistical analysis (i.e. ordinary least squares regression or analysis of variance) for each field trial. On the other hand, high uncertainty in yield variations derived from unquantified factors (e.g. soil, weather, farmer skill, and historical crop management practices) hinders holistic data analysis from multiple field trials, as each field might have different yield-limiting factors. To overcome this issue, the present study evaluated the number of fields achieving positive profits at the different probability levels (95, 90, and 80%) using the posterior distributions derived from each Gaussian process model (Figure 5). This approach can provide valuable information on fertilization strategies while accounting for model uncertainty and variability in fertilizer and grain prices.

In this study, the results obtained with the Bayesian model enable the evaluation of net revenue associated with different fertilizer use scenarios in rice fields. These findings underscore the financial advantages of reduced fertilizer compared to conventional rates in the majority of cases (Figure 3), emphasizing the economic benefits of optimizing fertilizer usage. However, two of the eleven cases (Fields 4 and 7) showed a substantial reduction in yield with decreasing fertilizer usage (Table 2), which led to unfavorable financial outcomes (Figure 3). In contrast, in Fields 1 and 3, a lower fertilizer rate led to an increase in yield (Table 2), resulting in consistently positive net revenue across all iterations (Figure 3). However, rarely does lower fertilizer rates contribute to increasing rice yield if there is no negative lodging impact (Basak et al., 1962) or potential risk of increased disease incidence (Long et al., 2000) due to high N application above the recommended rate; such situations did not occur in the trials. This discrepancy raises concerns about potential artifacts in collected yield data or limitations of the Gaussian process model in accommodating spatial yield variability, which in turn could impact the precision of the analytical outcomes (Besag & Higdon, 1999; Katzfuss et al., 2020; Marchant et al., 2019). The study generated yield maps by employing UAV-based yield prediction model, with a moderate level of prediction accuracy (R² = 0.68). However, this way might create biases in yield data that was used for the Gaussian process models as the response variable. Provided that errors in the estimated yield are assumed to be independently and identically distributed, the noise term in Equation (2)(2) $\mathit{k}\left(d_{\mathit{i},\mathit{j}}|\mathit{\alpha },\mathit{\rho },\mathit{\sigma }\right)={\alpha }^2\exp \left(-\frac{1}{2}{\left(\frac{{\mathit{d}}_{\mathit{i},\mathit{j}}}{\mathit{\rho }}\right)}^2\right)+{\delta }_{\mathit{i},\mathit{j}}{\sigma }^2$(2) could account for these errors. Thus, the analytical outcomes should depend on the degree of yield estimation errors and the distribution of random noise. This potential problem highlighted the necessity for continued research to enhance the accuracy and reliability of yield predictions. Furthermore, Tanaka (2021) assessed the precision of various statistical approaches in detecting treatment effects from OFEs. They added assumed treatment effects and random noise to real yield map data. Similarly, further simulation research should investigate how spatial yield prediction errors could impact the precision of treatment estimates and subsequent revenue analysis.

Overall, the majority of field trials achieved positive revenue gain especially when fertilizer price was high and grain price was low (Figure 5). The result indicated the potential financial advantage of applying lower fertilizer rates rather than higher fertilizer rates. Notably, a higher fertilizer rate was not a promising strategy for obtaining positive revenue within the expected range of fertilizer and grain prices. However, the increase and reduction rates of applied fertilizer for the treatments ranged from 10 to 14%. The results indicated that a greater reduction in fertilizer application might be economically feasible, as eight of eleven field trials showed positive revenue gain at the 95% probability when fertilizer price was high and grain price was low. Japan’s strategy for sustainable food systems (Strategy Midori) includes a 30% reduction in chemical fertilizer application rates by 2050 (Ministry of Agriculture Forestry and Fisheries, 2021). Due to the continuously increasing cost of chemical fertilizer (Hellegers, 2022; Rabbi et al., 2023; Randive et al., 2021), efficiently reducing fertilizer application rates has become a critical topic for farmers. Future OFEs can inform an agronomic understanding of whether farmers can further reduce the total amount of fertilizer applied. However, farmers are typically unwilling to reduce fertilizer application rates due to unpredictable yield losses and are skeptical about the applicability of analytical outcomes from strip trials for the rest of their fields (Tanaka et al., 2023). Accordingly, some farmers may be unwilling to conduct OFEs if the fertilizer application reduction rate is more than 30% for the experimental treatments. In response to these challenges, the Bayesian data analytical framework for multiple OFEs coupled with the revenue analysis proposed in this study may help encourage farmers’ active engagement in OFEs and therefore provide recommendations on the optimal fertilizer input rate without reducing farmers’ financial profits. More importantly, as the probabilities decreased, there was an increase in the number of field trials achieving positive net revenue (Figure 5). Note that OFEs are not only for scientists but also for farmers to create and share the local knowledge (Lacoste et al., 2022). The detection of significant difference among treatments are not always be the primary focus for OFEs. Therefore, exploring how farmers respond to the outcomes based on different probability levels could be an intriguing research topic.

Conclusion

The OFE was conducted over three years in Gifu, Japan, and provided valuable insights into optimizing chemical fertilizer usage and its impact on profitability. The experimental design involved varying fertilizer rates through a strip trial approach and introduced spatial variation, challenging traditional statistical methods. The utilization of a Gaussian process model proved effective in addressing spatial intricacies and assessing treatment effects across diverse field conditions. The study highlighted the financial advantages of lower fertilizer rates than farmers’ conventional rates in the majority of cases. Furthermore, the study successfully assessed the frequency of revenue gain by counting the number of field trials achieving positive net revenue at the certain probability levels while accounting for the uncertainties in fertilizer and grain prices. The results also suggested that adopting reduced fertilizer rates could be a profitable strategy, offering valuable insights for farmers facing rising fertilizer costs. Moreover, adopting an increased fertilizer rate might not be profitable and could be equated to high-risk gambling. The graphical representations underscore the potential financial benefits of lower fertilizer rates in the face of market uncertainties. This study contributes to the ongoing debate on optimizing fertilizer usage in rice production and emphasizes the importance of integrating multiple OFE outcomes to inform decision-making in farming practices. Nevertheless, these findings relied on the UAV-based yield estimations, whose errors could potentially influence the analytical results. Further research should focus on exploring the impact of yield estimation errors on statistical robustness.

Acknowledgments

This study was supported by the JST FOREST Programme, Grant Number JPMJFR221C. We would like to thank the Japan Agricultural Cooperative Nishimino and the farmers for their permission for the field surveys.

Disclosure statement

No potential conflict of interest was reported by the author(s).