{'fraction_non_alphanumeric': 0.030616150019135095, 'fraction_numerical': 0.013394565633371604, 'mean_word_length': 4.989875835721108, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Assessment of research frameworks for on-farm experimentation through a simulation study of wheat yield in Japan

On-farm experiments can provide farmers with information on more efficient crop management in their own fields. Recent developments in precision agricultural technologies, such as yield monitoring and variable-rate application technology, allow farmers to implement on-farm experiments easily. Research frameworks including the experimental design, data collection and preprocessing steps, and the statistical analysis method strongly influences the precision of the experiment. Conventional statistical approaches (e.g., ordinary least squares regression) may not be appropriate for on-farm experiments because they are not capable of accurately accounting for the underlying spatial variations in a particular response variable (e.g., yield data). We explored the effect of sensor types, data preprocessing, experimental designs, and statistical approaches on the type I error rates and estimation accuracy through a simulation study hypothesized to conduct on-farm experiments in 3 wheat fields in Japan. Isotropic and anisotropic spatial linear mixed models were established for comparison with ordinary least squares regression models. The repeated designs were not sufficient to reduce both the risk of a type I error and the estimation bias. A combination of a repeated design and an anisotropic model is required to improve the precision of the experiments. Although the anisotropic model sometimes did not contribute to reducing the bias efficiently, the results of the anisotropic model showed large standard errors, especially when the estimates had large biases. This finding highlights an advantage of anisotropic models since they enable experimenters to consider the reliability of the estimates when needed.


Introduction
On-farm experimentation is a means of farmer-centric research and extension that examines the effect of crop management (e.g., fertilizer application, irrigation, and pest control) and variety selection on crop productivity in farmers' own fields [1,2]. Since the last century, agricultural experiments have been primarily performed by researchers in experimental fields under highly controlled conditions to ensure the accuracy of the estimated treatment effects. This typical research approach has contributed to improving our understanding of crop physiology and to developing agronomic practices, but the research results are not straightforwardly tailored to entire fields or regions. Thus, farmers and crop advisors have learned how to adjust crop management techniques by trial and error on farms [3]. Recent developments in precision agricultural technologies, such as yield monitoring for combine harvesters and variable-rate application technology (e.g., fertilizers, seeds, and herbicides), allow farmers to run on-farm experiments easily on their own farms [4,5,6]. On-farm experimentation has been gaining popularity by providing information on the best crop management techniques for specific regions, farmers, and fields [2]. Moreover, new challenges for the implementation of on-farm experiments have arisen.To evaluate the precise effect of genotype and crop management on crop productivity, arXiv:2004.12741v1 [stat.AP] 27 Apr 2020 conventional small-plot experiments have been widely applied in agricultural research. These conventional small-plot experiments depend on the combination of three basic principles of experimental design (randomization, replication, and local control) and statistical approaches such as analysis of variance (ANOVA), which were established by Fisher [7]. This fundamental agronomic research framework effectively separates the spatial variations and measurement errors from the observed data (e.g., the grain yield) to detect the significance of the treatment effect [1]. Widely used conventional statistical approaches, including ANOVA or ordinary least squares (OLS) regression, depend on the assumption that errors are independent. However, soil properties and crop yield are not spatially distributed at random, and similar values are observed near each other, which is called spatial autocorrelation [8]. Spatial autocorrelation in a response variable (e.g., the crop yield) violates a conventional statistical assumption, which leads to unreliable inferences (e.g., overestimation of the treatment effects) [9,10]. Thus, the conventional approach is not directly applicable to on-farm experiments, which generally have larger areas and simpler experimental arrangements.
To account for the underlying spatial structures, model-based geostatistics have been developed in disciplines related to agricultural science [11]; thus, geostatistical approaches have been applied to analyze on-farm data. For instance, yield data were derived from chessboard or repeated strip trials were kriged (interpolated) for the other treatment plots, and a yield response model was successfully established with a regression model [12,13]. However, it does not involve straightforward estimation with multiple kriging processes; it requires repeated and complicated experimental designs; thus, farmers cannot accept this approach easily. On-farm experimental designs must consider technical, agronomic, and economic restrictions, but they also have to be statistically robust to provide reliable estimates of treatment effects [10]. A simulation study examined the effect of spatial structures, experimental designs, and estimation methods on type I error rates and bias of treatment estimates by using geostatistical approaches [14]. This study indicated that higher autocorrelation significantly increased the type I error rates, while a spatial linear mixed model reduced them regardless of the experimental design, and more randomized and repeated experimental designs (e.g., split-planter, strip trials, and chessboard designs) increased the accuracies of the treatment estimates. Furthermore, Marchant et al. (2019) [6] demonstrated that a spatial linear mixed model representing anisotropic spatial variations could successfully evaluate the treatment effect even in simple on-farm experimental designs (e.g., strip trials). Appropriate experimental designs and statistical approaches may vary according to the spatial variability and farmers' available machinery, and a tradeoff between the simplicity of the experimental design and the desired precision of the outcome should be considered when conducting on-farm experiments [14]. Although there may not be a versatile framework for the implementation of on-farm experiments, the accumulation of more knowledge regarding the relationships between experimental designs and statistical approaches should be examined in a variety of targets.
Crop performance can be measured by several types of sensors. Yield monitoring by a combine harvester is one of the most frequently used measurement methods for on-farm experimentation. However, the yield monitor data inherently include artifacts and noise not related to the actual yield variation; thus, filters for removing potential artifacts have been advised [6]. Data preprocessing is very important for obtaining reliable results. Moreover, a large number of spatial data points are generated by sensor measurements, which is becoming a challenging task because it increases the computational time for spatial analyses. Remote sensing technologies such as unmanned aerial vehicles (UAVs) and satellite imageries provide valuable spatial datasets for on-farm experimentation, and they often generate substantial spatial data points. To achieve a feasible computational time and the desired precision of the outcome, a data preprocessing framework (e.g., smoothing or reducing the number of spatial data points) should be considered.
The objectives of this study were to assess the effects of different sensor types, data preprocessing, experimental designs, and statistical approaches on the type I error rates and estimation accuracy through a simulation study of on-farm experiments on wheat production in Japan. We further examined the inference framework for experimenters from the perspective of model uncertainty. The predicted yield data were derived from remotely sensed imageries and commercial yield monitors for combine harvesters. Several data preprocessing methods (e.g., the grid size for in the averaging yield or a reduction in the yield data points) and hypothetical experimental designs were assumed to have been applied to those datasets. A state-of-the-art spatial anisotropic model was developed to account for spatial autocorrelation and to reduce the estimation bias, and it was compared with OLS regression and standard spatial isotropic models as traditional statistical approaches.

Yield data and experimental designs
We explored experimental designs and statistical approaches suitable for on-farm experiments in Japan through a simulation study of winter wheat yield. Three fields were used for the simulation study (Table 1). Currently, yield monitors on commercial combine harvesters are still not prevalent in Japan. Thus, crop yield maps of two upland fields converted from paddy fields (Fields 1 and 2) were derived from a previous study performed by Zhou  These fields were remotely sensed by a commercial multispectral camera (Sequoia+, Parrot, France) mounted on an UAV at the grain-filling stage in 2018 and 2019. Briefly, winter wheat yield was predicted using a linear regression model with a predictor of enhanced vegetation index 2 derived from the imageries. The ground sample distance of those imageries was 0.06 m pixel −1 . As noted above, yield monitor data from Japanese farmers' fields are not commonly available, so yield monitor data were collected in 2016 at a demonstration farm (Field 3) of New Holland HFT Japan, Inc. by a commercial combine harvester CX8.70 (New Holland, Belgium) (Fig. 3). The fields were located in Tomakomai, Hokkaido, Japan (42 • 45' N, 141 • 44' E). The entire field size was 8.45 ha, and the area used for the simulation was 100 × 200 m (2.00 ha). The header width was approximately 5 m: a typical size for a large-scale farm in Hokkaido, Japan. The yield monitor on the combine harvester recorded the crop yield along rows at an interval of approximately 1.3 m.
Two fertilizer treatments (control and treatment plots) were assumed to have been applied in the fields. The experimental design might be very important to reduce the risk of type I error and to evaluate the treatment effect precisely. From a practical viewpoint, we assumed that a tractor equipped with an 18-m working width broadcaster was used for the on-farm experiments, which is equivalent to 1.5 round trips along the long-side direction for Fields 1 and 2 and is equivalent to 2.5 round trips for Field 3. All fields were used for the experiments, but the edges (6 m from long-side edges and 15 m from the narrow-side edges) were excluded from the analyses to avoid edge effects in Fields 1 and 2. Four experimental designs were simulated in Field 1 ( Fig. 1) and Field 2 ( Fig. 2), and 5 experimental designs were simulated in Field 3 (Fig. 3). A simple strip trial (D1) was the easiest and most practical experimental design for farmers. A simple split-plot trial (D2) was established by splitting the experimental plots perpendicular to the farming operations, and it might require more complicated manual operation than D1. Just for Field 3, we further considered repeated strip trials (D3) because the field size was sufficiently large to implement it. A combination of strip(s) and split-plot trials (D3 for Field 1 and 2, D4 for Field 3) was established. A more repeated systematic design was established (D4 for Fields 1 and 2, D5 for Field 3), which could not be implemented without VRA technology. Note that our analysis was performed through a simulation study, but the dataset was based on real collected data, which allows us to examine the effect of anisotropy and unpredictable variations in the actual fields.   The border points indicate the area that was not used for the analysis due to potential edge effects.

Statistical analysis
OLS regression is based on the assumption that the errors are independent. In assessing the significance of the treatment effect on crop yield in on-farm experiments, spatial autocorrelation should be considered because the OLS estimator increases the risk of type I error [9,10]. Thus, a spatial linear mixed model was used to evaluate the effects of hypothetical treatment on wheat yield. The spatial linear mixed model is written as where y is a vector of length n of the response variable, n is the number of measurements, X is the n × p fixed design matrix with the values of the vector of size p, which is the number of independent variables, β is a vector of length p of the fixed-effects coefficients, and is a vector of length n of the random effects with covariance matrix V . The random effects are assumed to be spatially correlated, and the exponential function was used for the covariance estimation. The covariance function is written as where c 0 is the nugget variance, c 1 is the sill variance, h is the distance between the two measurements, and a is the distance parameter. The variogram is written as The above model is a geometrically isotropic model, which has the same parameters for the covariance functions in all directions. Marchant et al. (2019) [6] reported that yield monitor data showed strong anisotropy between the direction of the combine harvester's rows and perpendicular to the direction of the rows; thus, a product-sum covariance model [15] was used to model the variation along each direction. Furthermore, the directions of farming operations (e.g., sowing and fertilizer application) showed different spatial variabilities in the crop yield surveyed by manual harvesting, and a geometrically anisotropic model for stretching the range parameter according to certain directions was not sufficient to fit the variograms because the sill variance was unlikely to be identical for each direction, as suggested by Zhou et al. (under review). To establish an anisotropic model, direction-specific covariance functions should be parameterized. The sum-metric model presented by Bilonick [16] was used in this study. The sum-metric model is written as where h x is the lag perpendicular to the farming operation (along rows), h y is the lag in the direction of travel (within rows), and h xy is the lag obtained by introducing a geometric anisotropy ratio α. The sum-metric model has been used for fitting space-time variograms previously [17], and its advantage is that the combination of 2-directional static components and 1 dynamic component of covariance function is easily interpretable in physical sense [18]. Thus, the isotropic model has 3 parameters to be estimated, and the anisotropic model has 8 parameters to be estimated. For the estimation of these random effects parameters, a restricted maximum likelihood (REML) estimator was used. The advantage of the REML estimator is that it does not depend on the unknown fixed effects; therefore, the estimates are less biased than maximum likelihood estimates [19]. Consequently, the REML estimator calculated the fixed-effect coefficient β and its standard error. The statistical significance of the experimental treatment was assessed by z statistics, and two-sided p-values were computed. The preferred model was evaluated based on the lower values of the Akaike information criterion (AIC) [20] between the geometrically isotropic and anisotropic models. Most of the available applications for model-based geostatistics, such as the geoR package [21] implemented in the R environment [22], can only fit isotropic variograms based on the REML estimator. The R package gstat [23] is available for fitting space-time variograms, but the REML estimator is not implemented. Therefore, it is necessary to develop a user-friendly application that can fit an anisotropic model based on the REML estimator (Zhou et al. under review). All of the computations were conducted using MATLAB [24], and the documented MATLAB source code is available at GitHub (https://github.com/takashit754/geostat). Finally, we compared the performances of 3 models: OLS regression models, spatial isotropic linear mixed models, and spatial anisotropic linear mixed models. For the spatial linear mixed models, we further evaluated the residual experimental and fitted variograms to check the underlying spatial variations.
It is not feasible to use raw yield data in linear mixed models due to the high computational cost. Averaged values in a certain grid size are generally used for mapping high-resolution spatial data; however, the grid size significantly affects the within-field spatial structures [25]. Thus, three grid sizes of 1.0 × 1.0, 2.5 × 2.5, and 5.0 × 5.0 m were used to generate averaged yield data points derived from the UAV-based imageries in Fields 1 and 2. Furthermore, we preprocessed the yield monitor data in Field 3 according to the method proposed Marchant et al. (2019) [6] with modifications. Briefly, the coordinates were rotated as the combine harvester traveled in the y direction. The x coordinates were adjusted across each row in an exact straight line. Small variations in the coordinates and noise in the yield monitor data located near each other may prevent the estimation of the model parameters. A square lattice of points was preferred to fit the anisotropic model, as it could accommodate the computational cost and improve the precision of the maximum likelihood calculations. Therefore, intervals of 2.5 and 5.0 m in the y direction were used for averaging the yield monitor data points. For the interval of 2.5 m, each averaged yield data point contained 1-2 raw yield data points, while for the interval of 5.0 m, each averaged yield data point contained 3-4 raw yield data points.
The simulated type I error rates were presented as p-values of the experimental treatment on the assumption that the treatment population parameter was zero. To assess the effect of the grid size in the averaging, experimental designs, and statistical models on the simulated type I error rates, three-way analysis of variance (ANOVA) based on type III sums of squares was performed.
To assess the estimated accuracy, randomly generated numbers from a Gaussian distribution (mean±SD: 0.30±0.10 t ha −1 ) were added to the yield data for each point in the treatment plots. Then, the bias was estimated by computing the difference between the fixed-effect coefficient β from each statistical model and the population parameter yielded by a Gaussian random number generator (approximately 0.3 t ha −1 ). In addition, 95% confidence intervals were calculated using the standard error estimated from each model. To assess the effect of the grid size, experimental designs, and statistical models on the absolute bias, three-way ANOVA based on type III sums of squares was performed.

Simulated type I error rates
The effects of the experimental design, grid size, and model on the simulated type I error rates were evaluated by three-way ANOVA ( Table 2). There were significant effects of the experimental design on the simulated type I error rates in Fields 2 and 3. The simulated type I error rates for each field are shown in Tables 3 and 4, respectively. The simulated type I error rates were greater than the significance level (>0.05) in the simplest design (D1) for Field 1 (Table 3). Moreover, they were greater than the significance level (>0.05) in the most repeated and complicated design (D5) in Fields 2 and 3 (Tables 3 and 4, respectively). Although a smaller grid size tended to result in lower simulated type I error rates for all the fields, the grid size significantly affected the simulated type I error rates only for Field 2 ( Table 2). Overall, lower simulated type I error rates were frequently observed in the OLS model than in the other models. For the anisotropic model, there were significant simulated type I error rates only in D2 and D1 for Field 1 and Field 2, respectively. There were significant effects of model selection for Fields 1 and 2 ( Table 2). There was a significant interaction between the experimental design and the model in Field 2 only.

The precision of estimates
The effects of experimental design, grid size, and model on the bias were evaluated by three-way ANOVA ( Table 5). The experimental design had a significant effect on the bias for all the fields. Model selection had a significant effect on the bias for Field 3 only. There were significant interactions between the experimental design and model for all the fields. The bias and 95% confidence intervals for the 2.5-m grid are shown in Fig. 4 since the grid size had no significant effect on the bias ( Table 5). The closer the bias is to zero, the closer the estimates are to the actual treatment values yielded by the Gaussian random number generator (approximately 0.3 t ha −1 ). Overall, the OLS model had narrower 95% confidence intervals, and the 95% confidence intervals did not contain zero more frequently than the other models. Furthermore, the 95% confidence intervals of the OLS model did not contain zero even in the most repeated and complicated design (D5) of Field 1. The bias were greatly reduced by using an anisotropic model in Field 3.

Variograms
The fitted 2-directional residual variograms are shown in Fig. 5. In Field 1, the sill variance in the fitted variogram was 0.011 (t ha −1 ) 2 in the direction of farming operations (y), while it was 0.004 (t ha −1 ) 2 in the direction perpendicular to the farming operations (x). In Field 2, the sill variance in the fitted variogram was approximately 0.003 (t ha −1 ) 2 in both directions. However, the range parameter that reaches the sill variance at the 95% level was 129.1 m in the direction of farming operations (y), while it was much larger than the maximum lag (1382 m) in the direction perpendicular to the farming operations (x). Consequently, there were likely only nugget effects in the direction perpendicular to the farming operations (x) in Fields 1 and 2. In contrast, in Field 3, the sill variance in the direction perpendicular to the farming operations (x) was 7 times larger (0.028 (t ha −1 ) 2 than in the direction of the farming operations (y) (0.196 (t ha −1 ) 2 . Overall, there was more or less anisotropy for all fields, and the largest differences in the sill variance between the directions were detected in Field 3.

Discussion
The sensor types, data preprocessing, experimental designs, statistical approaches, and within-field spatial structures affect the precision in on-farm experiments [14,6]; thus, there are many possibilities for the best experimental design and statistical approaches in on-farm experiments [2]. We explored those effects on the type I error rates and estimated accuracy through a simulation study that generated hypothetical treatments on real wheat yield data in Japanese fields. It is not feasible to establish a complete and versatile framework for on-farm experiments, and it is necessary to improve the experimental setup method according to the conditions. Several important implications for experimenters can be drawn from our results.
In Fields 2 and 3, the experimental designs significantly affected the simulated type I error rates (Table 3), and the most repeated and complicated design (D5) successfully accommodated the simulated type I error rates ( Table 4). The results agreed with a previous study, which reported that designs with fewer replications and larger experimental units tended to increase the risk of type I error [14]. However, in Field 1, there were no significant effects of the experimental design on the simulated type I error rates. Furthermore, the simulated type I error rate was relatively higher in the simplest design (D1) than in the other repeated designs in Field 1 (Table 4). These results indicated that repeated designs are not sufficient to avoid the risk of type I error. Moreover, the model was also a significant factor affecting the simulated type I error rates in Fields 1 and 2. Although the simulated type I error rates were not significantly affected by the model in Field 3 (p=0.059), they were always greater than the significance level (>0.05) in the anisotropic model. Similarly, the experimental designs significantly affected the bias for all the fields (Table 5), and the repeated designs (D3 and D4 for Field 2 and D4 and D5 for Field 3) greatly reduced the bias (Fig. 4). Consequently, the combination of a repeated design and an anisotropic model might be a solution for avoiding the risk of type I error and reducing the estimation bias. However, between the experimental design and model, there was a significant interaction on the simulated type I error rates for Field 2 and significant interactions on the bias for all the fields. These interactions indicated that the best model could be varied according to the experimental design and field conditions. Therefore, this recommendation is still restrictive to obtaining robust outcomes through on-farm experiments.
The fitted residual variograms showed strong anisotropy with different sill variances between the directions for yield monitor data (Field 3) rather than the remotely sensed yield data (Fields 1 and 2) (Fig. 5). The bias was greatly reduced by using an anisotropic model in D2 for Field 3. Furthermore, the anisotropic model had 95% confidence intervals that contained zero for four designs (D1, D3, D4 and D5), although OLS and isotropic models had 95% confidence intervals that contained zero only for two designs (D4 and D5). These results are in agreement with the finding from Marchant et al. (2019) [6], who demonstrated that isotropic models are not appropriate to account for spatial autocorrelation; thus, they applied a product-sum model [15], which is more complex than the standard isotropic model. We also demonstrated that the sum-metric model [16], as the state-of-the-art model, separated the underlying spatial variation from the yield monitor data to evaluate the treatment effects more accurately than the isotropic model. Thus, an anisotropic model is necessary for on-farm experiments that use a yield monitor on a combine harvester. In this study, we tested the exponential function in all directions, although there are many available covariance functions, such as the Matérn function. The selection of the covariance function should be considered in future studies.
For the remotely sensed data, the anisotropic model successfully reduced the risk of type I error (Table 3), although it did not contribute to reducing the bias (Table 5). There was also more or less anisotropy in the remotely sensed data (Fields 1 and 2) (Fig. 5). Furthermore, the AIC always showed smaller values in the anisotropic model than in the isotropic model (data not shown). Thus, an anisotropic model is recommended for the analysis of on-farm data even derived from remotely sensed data. Moreover, this anisotropy may be attributed not only to the direction of farming operations but also to the field size or shape. Typical Japanese paddy fields are rectangular and small (e.g., less than 1 ha), so a narrow side may not be sufficiently large to exhibit developed spatial structures. For instance, in Fields 1 and 2, the maximum lag on the narrow side was only 38 m, while it was 100 m in Field 3. There was apparently no sill variance and only nugget effects along the narrow side. Therefore, these two small-scale paddy fields may inherently contain anisotropy. It is noteworthy that these findings may be specific to small-scale fields in Asian countries. Previous on-farm experiments have been carried out using commercial yield monitors on combine harvesters or UAV-based remote sensing in large-scale fields ranging from 8 to 16 ha in the UK [6] and from 10 to 100 ha in the US, South America, and South Africa [26]. Further research is still required to confirm the factors underlying this anisotropy and to explore the best research framework to implement on-farm experiments appropriate for Asian countries with relatively small fields.
A smaller grid size tended to decrease the simulated type I error rates, although it was significant in Field 2 only ( Table 2). The grid size did not have a significant effect on the bias ( Table 4). Considering that all range parameters were greater than the maximum grid size of this study (5 m) and that a small grid size increases the computational time, it is not necessary to use a detailed yield map such as a 1.0-m grid. Furthermore, downscaled coarse spatial data can reduce the inconsistencies in the model assumptions to fit the spatial variability better than raw fine spatial data [27]. However, one should be aware that a larger grid size would ignore the small-scale variations in the raw data [25], which may influence the estimations of the nugget and sill variance and further treatment effects.
To reiterate, an anisotropic model is an incomplete solution but provides more robust outcomes than traditional statistical approaches. An anisotropic model is advantageous since it covered larger standard errors when the estimates had large biases (Fig. 4). Moreover, the 95% confidence intervals of the OLS and isotropic models were generally narrow, and the hypothetical treatment effects fell outside of them. To make a reliable decision according to the results of on-farm experiments, experimenters should keep in mind that estimates may not always be precise, but they can consider how much it could vary, as this variability would result in an adverse scenario. Farmers are not interested in whether the treatment is significant at the 0.05 level but rather in whether there will be a return on investment [1]. Therefore, as previous studies determined [28,10], it is necessary to examine how much possible marginal profits could vary by using real on-farm data in further studies.

Conclusions
The outcomes of on-farm experiments can support farmers' decision-making processes, while inappropriate procedures would result in incorrect interpretations. The repeated experimental designs examined here can partially contribute to reducing the risk of type I error and bias, but they were not sufficient to separate the underlying spatial variations from the yield data to evaluate the treatment effects accurately. Although we have many choices for the experimental design and statistical approaches, the combination of repeated designs and anisotropic models provides more reliable outcomes than the other methods to avoid issues arising from on-farm experiments. The results of the anisotropic model showed large standard errors, especially when the estimates had large biases. Considering that the aim of on-farm experiments is to provide farmers with information on economic feasibility, these statistical characteristics of anisotropic models are advantageous, as experimenters have opportunities to infer the analytical results conservatively. To examine the effect of these statistical characteristics on farmers' decision-making process, economic analysis is needed by using real on-farm data in the future.