# Model and auxiliary data for an accurate estimate of mean field yield

Baptiste Oger, Sébastien Roux, Gilles Le Moguédec, Bruno Tisseyre

► To cite this version:

Baptiste Oger, Sébastien Roux, Gilles Le Moguédec, Bruno Tisseyre. Model and auxiliary data for an accurate estimate of mean field yield. 13. European Conference on Precision Agriculture (ECPA 2021), Jul 2021, Budapest, Hungary. pp.669-676, 10.3920/978-90-8686-916-9_80. [hal-03833247](https://hal.inrae.fr/hal-03833247v1)

HAL Id: hal-03833247

https://hal.inrae.fr/hal-03833247v1

Submitted on 28 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Model and auxiliary data for an accurate estimate of mean field yield

B. Oger<fnref n="1" />⁻<fnref n="2" />, S. Roux<fnref n="2" />, G. Le Moguédec<fnref n="3" /> and B. Tisseyre<fnref n="1" />

baptiste.oger@supagro.fr

## Abstract

Recent research has shown how high spatial resolution auxiliary data may improve sampling for field estimates. This paper defines theoretically the variance and bias of field estimates resulting from sampling with a model-based and a mean-based inference. The factors influencing bias and variance are identified and their respective impacts quantified. The theoretical formalisation demonstrates how the characteristics of the auxiliary data sampled and their representativeness affect the mean estimation. Both model-based and mean-based approaches are then tested and compared with real and simulated data. The results highlight the differences between model-based and mean-based for estimating mean yield field values. This paper highlights the value of using a model-based approach when auxiliary data are available. In most cases results are notably better than the use of common mean-based inference.

**Keywords:** sampling, yield, vegetation index, inference, viticulture, estimator.

## Introduction

To adjust field operations (fertilization, seeding, crop protection, etc.) at the field level, growers need to know the average characteristics of their production system. They generally rely on sampling techniques, which are time-saving, to collect this information. The principle of field sampling is to estimate characteristics of a population from within-field measurements taken for a limited number of sites. As an example, sampling estimation can be done to estimate the field yield by observing relevant yield components at a limited number of measurement sites. The objective is to estimate the quantities that will be collected to fit harvesting operation and storage requirements.

Estimation methods based on field sampling mainly consist of two steps. The first one focuses on selecting the sites that should be observed. This step has been already extensively studied in the literature and precision agriculture has brought improvements by proposing methods to drive the choice of the sites according to high spatial resolution data (Tillé, 2011). The second step of estimation is based on the definition of the estimator based on samples values. The purpose of this step, on which this paper focuses, is to infer the characteristics of the whole population based on observed sample values. Surprisingly, the principles of this second step have evolved little with precision agriculture. In many studies, it is still based on classical approaches like the mean or the median of the collected samples. In recent publications, Carillo et al (2016) and Araya-Alman (2019) proposed an empirical approach to account for the spatial distribution of high resolution auxiliary data in the estimation step. These authors assumed the availability of high resolution spatial auxiliary data (i.e. vegetation index derived from multispectral remote sensing images) more or less correlated to the variable of interest to

<footnote n="1">¹ITAP, Univ. Montpellier, Montpellier SupAgro, INRAE, France</footnote>

<footnote n="2">²MISTEA, Univ. Montpellier, Montpellier SupAgro, INRAE, France</footnote>

<footnote n="3">³AMAP, Univ. Montpellier, INRAE, France</footnote>
be estimated. They proposed the use of sample values to calibrate a local model instead of calculating directly the field mean from sample values. The overall idea was to apply this model to all the field sites to estimate the variable of interest with high resolution. The mean field is then described from the estimates. Such an approach is assumed to better account for the spatial variability highlighted by the auxiliary data. This empirical approach, called model-based sampling, has been shown to be relevant to estimating grape yield from a normalized difference vegetation index (NDVI) (Carillo et al., 2016). However, the formalisation of its relevance has never been demonstrated.

This study proposes to demonstrate the potential advantages of a model-based estimator when compared to the classical arithmetic mean estimator. For each estimator, this study proposes a statistical formalism to describe expected bias and variance. This formalisation is then applied to simulated data of known characteristics to validate the relevance of the equations. Tests on real fields are performed in a second step to validate the formalisation on real cases. The examples (either real or simulated) refer to field yield estimation in viticulture with a vegetation index (NDVI) as auxiliary data.

# Theory

## Definitions

A given plot, K is defined as the set of potential within-plot sample sites, numbered from 1 to k. The spatial footprint of a sample site being defined by that of the measurement or that of the resolution of an auxiliary variable, it is assumed that there is a finite number of potential sampling sites within a field. For a given site $i \in K$, $X_i$ and $Y_i$ are respectively the auxiliary variable and the variable of interest measured at that site. The auxiliary data is assumed available at low cost and at a high spatial resolution. The auxiliary variable $X_i$ is known for every site but the variable of interest $Y_i$ is only known for a set N of sampled sites (where n is the size of the N set). It is assumed that $Y_i$ can be predicted from $X_i$ and a model. In a first approach, a basic linear model is considered (Eq. 1):

$$Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i \quad \text{with} \quad \varepsilon_i \sim N(0, \sigma^2) \text{ independent} \quad (1)$$

Let $\bar{Y}$ be the average value of Y for the whole plot:

$$\bar{Y} = \frac{1}{k} \times \sum_{i \in K} Y_i \quad (2)$$

This study assumed that measurement sites are determined (the N set is fixed). This ensures to exclude sampling-related error and to only focus on estimator related error. In the following sections, demonstrations to formalise the bias and the variance of the two considered estimators are presented; the classical mean estimator which will be called $\hat{Y}_1$ and the model-based estimator which will be called $\hat{Y}_2$. For readability reasons, the details of all the demonstrations are not presented. For more details, the reader can refer to Oger (2021) appendix available in the references.

## Mean-based inference

$\hat{Y}_1$ refers to the estimator of $\bar{Y}$ based on the mean inference. $\hat{Y}_1$ is defined in Eq. 3 as the average of the values measured $Y_i$ ($i \in N$):

$$\hat{Y}_1 = \frac{1}{n} \times \sum_{i \in N} Y_i \quad (3)$$
Main statistical properties of this estimator are given by its expected value (for the bias) and its variance (for the precision). Although the mean estimator is not based on auxiliary data, it is assumed that such data exist to characterize its properties. Its expectation (E) is then characterized in Eq. 4.

$$E[\hat{Y}_1] = \bar{Y} + \beta_1(\bar{X}_N - \bar{X}) \quad (4)$$

$\hat{Y}_1$ presents a bias which depend on $\beta_1$, the slope parameter of the linear model which relates $X_i$ to $Y_i$. For the record, $\bar{X}_N$, the auxiliary data mean value for the sample and $\bar{X}$, the auxiliary data mean value for the whole plot which is known. The difference between $\bar{X}_N$ and $\bar{X}$ models the influence of sample representativeness, showing that ideally the sample should have properties close to those of the population. Since representativeness is expressed in terms of auxiliary data, its influence depends on $\beta_1$ which characterizes the strength of the correlation between the variables $Y_i$ and $X_i$. The variance of $\hat{Y}_1$ is presented in Eq. 5.

$$Var[\hat{Y}_1] = \frac{\sigma^2}{n} \quad (5)$$

The variance can be described with $\sigma$, the residual variance of the linear model and $n$, the number of sampling sites. This result is logical, the higher the unexplained variability of the variable of interest, the more variable the estimate will be. Conversely, a larger sample size improves confidence in the estimate.

## Model-based inference

Let $\hat{Y}_i$ be the estimate for $Y_i$ based on the linear model calibrated using the sample observations. $\hat{Y}_i$ is thus based on $\hat{\beta}_0$ and $\hat{\beta}_1$, the estimated value of parameters $\beta_0$ and $\beta_1$. A second estimator, $\hat{Y}_2$, can be described from the set of values predicted by the model:

$$\hat{Y}_2 = \frac{1}{k} \times \sum_{i=1}^{k} \hat{Y}_i \quad \text{or} \quad \hat{Y}_2 = \frac{1}{k} \times \sum_{i=1}^{k} (\hat{\beta}_0 + \hat{\beta}_1 \times X_i) \quad (6)$$

As for $\hat{Y}_1$, the expectation expression is derived for $\hat{Y}_2$ in Eq. 7 (see the appendix, Oger (2021), for the details of the demonstration).

$$E[\hat{Y}_2] = \bar{Y} \quad (7)$$

Eq. 7 shows that the estimation based on a linear model is unbiased. To compare the two estimators, it is nevertheless necessary to look at the $\hat{Y}_2$ variance (Eq. 8).

$$Var[\hat{Y}_2] = \frac{\sigma^2}{n} + \frac{(\bar{X}_N - \bar{X})^2}{\sum_{i \in N} (X_i - \bar{X}_N)^2} \times \sigma^2 \quad (8)$$

As for $\hat{Y}_1$, the $\hat{Y}_2$ variance depends on $\sigma$, the residual variance of the linear model and $n$ and the number of sampling sites. Eq. 8, shows that it also depends on the sampled sites and their auxiliary data. Indeed, the minimisation of the $\hat{Y}_2$ variance requires the sample to i) be representative, so its mean $(\bar{X}_N)$ has to be as close as possible to the field mean $(\bar{X})$ and ii) present a strong dispersion around its mean $(\sum_{i \in N} (X_i - \bar{X}_N)^2)$. This aspect highlights the necessity of sampling highly dispersed values around the mean to provide the best possible estimates of the model parameters. For the record, $X_i$ values refer to the auxiliary data which is assumed to be known at all the potential sampling sites therefore $\bar{X}$ is known.
## Inference comparison

By looking at Eq. 4 and 7, $\hat{Y}_1$ appears biased while $\hat{Y}_2$ is unbiased. But at the same time, $\hat{Y}_2$ presents a higher variance than $\hat{Y}_1$ (Eq. 5 and 8). For comparison purposes theoretical Root Mean Square Error (RMSE) were computed for both estimators. The computation was performed from bias and variance as stated in Eq. 9.

$$RMSE(\hat{Y}) = \sqrt{Bias(\hat{Y})^2 + Variance(\hat{Y})} \quad (9)$$

Based on Eq. 9, RMSE of $\hat{Y}_1$ and $\hat{Y}_2$ are given respectively in Eq. 10 and Eq. 11.

$$RMSE(\hat{Y}_1) = \sqrt{\frac{\sigma^2}{n} + \beta_1^2 \times (\overline{X}_{i \in N} - \bar{X})^2} \quad (10)$$

$$RMSE(\hat{Y}_2) = \sqrt{\frac{\sigma^2}{n} + \frac{(\overline{X}_N - \bar{X})^2}{\sum_{i \in N} (X_i - \overline{X}_N)^2} \times \sigma^2} \quad (11)$$

As a consequence, inequality $RMSE(\hat{Y}_2) < RMSE(\hat{Y}_1)$ is verified if and only if:

$$\frac{\sigma^2}{\beta_1^2} < \sum_{i \in N} (X_i - \overline{X}_N)^2 \quad (12)$$

Following this result, the choice of the inference strategy must be done according to i) $\sigma^2$, the residual variance of the linear model; ii) $\beta_1$, the regression slope of the linear model used to relate $X_i$ to $Y_i$; iii) $\sum_{i \in N} (X_i - \overline{X}_N)^2$, the dispersion of the sample in terms of auxiliary data. The $\frac{\sigma^2}{\beta_1^2}$ ratio is related to the existing correlation between the estimated variable and the auxiliary data. If the regression slope is small compared to the residuals of the model, the ratio value will be high and the mean-based estimator should be preferred. Conversely, if the slope is large compared to $\sigma$ and if the sample values are dispersed around the mean, then model-based estimator should be preferred.

## Material and Method

Both estimators were applied to a yield estimation case study in viticulture. Yield estimation was performed from a common auxiliary variable, the NDVI assumed available at a high spatial resolution. The tests were carried out on real commercial plots and theoretical plots on which yield and NDVI data were simulated. The theoretical dataset was considered to verify the demonstration under controlled conditions.

### Real yield Data

Real fields used to test the method come from the experimental vineyard of INRA Pech-Rouge (43.1392; 3.1332, France). Nine vine fields representative of the south of France vineyard in terms of a training system, variety and plantation density were considered. The characteristics of each plot are shown in Table 1. For each field, the yield was measured following a regular grid of 15 m approximately. On each node of the grid, yield was measured by weighting all the fruits on 5 consecutive vines in the row. Average yield was related to the location corresponding to the central vine of each node. NDVI values were derived from a 1 m resolution multi-spectral image at veraison acquired the same
year as yield measurements. From these 1 m square image pixels, aggregation method described in (Acevedo-Opazo et al., 2008) was used to obtain 9 m square image pixels. Aggregation reduced the noise brought by canopy discontinuity and bare soil pixels. NDVI values corresponding to the grid nodes of yield measurements were extracted from the images. The final data base was a set of 313 sites with known yield and NDVI over the 9 different fields.

Table 1 presents properties of real yield data as well as numerical values involved in RMSE comparison (Eq. 12). Regarding Eq. 13, the $\frac{\sigma^2}{\beta_1^2}$ ratios are to be compared with $\sum_{i \in N}(X_i - \bar{X}_N)^2$ which depends on the sample. For NDVI data, the ratios ranged from 1 and $10^{-1}$ except for plot 22 and 88 showing that for both these fields the mean-based estimator should be preferred.

<table><caption>Table 1: Properties and numerical values for commercial data.</caption><thead><tr><th>Field</th><th>Area (ha)</th><th>Variety</th><th>β<sub>1</sub></th><th>σ</th><th>σ<sup>2</sup>/β<sub>1</sub><sup>2</sup></th><th>Correlation</th></tr></thead><tbody><tr><td>P22</td><td>1.72</td><td>Syrah</td><td>570.4</td><td>995.8</td><td>3,0</td><td>0.13</td></tr><tr><td>P63</td><td>1.33</td><td>Syrah</td><td>4511.6</td><td>672.3</td><td>2.2.10<sup>-02</sup></td><td>0.28</td></tr><tr><td>P65</td><td>0.69</td><td>Syrah</td><td>17140.0</td><td>491.2</td><td>8.2.10<sup>-04</sup></td><td>0.86</td></tr><tr><td>P76</td><td>1.14</td><td>Carignan</td><td>6501.2</td><td>617.9</td><td>9.0.10<sup>-03</sup></td><td>0.39</td></tr><tr><td>P77</td><td>1.24</td><td>Syrah</td><td>12932.3</td><td>926.0</td><td>5.1.10<sup>-03</sup></td><td>0.48</td></tr><tr><td>P80</td><td>0.54</td><td>Syrah</td><td>12002.6</td><td>690.9</td><td>3.3.10<sup>-03</sup></td><td>0.63</td></tr><tr><td>P82</td><td>1.15</td><td>Syrah</td><td>10126.3</td><td>547.8</td><td>2.9.10<sup>-03</sup></td><td>0.47</td></tr><tr><td>P88</td><td>0.85</td><td>Syrah</td><td>-1184.9</td><td>852.0</td><td>5.2.10<sup>-01</sup></td><td>-0.04</td></tr><tr><td>P104</td><td>0.81</td><td>Carignan</td><td>3758.0</td><td>1100.0</td><td>8.6.10<sup>-02</sup></td><td>0.18</td></tr></tbody></table>

## Theoretical data

A unique theoretical field was simulated based on the approach proposed by Oger et al. (2020). It corresponded to a square plot of $100m \times 100m$ with 4000 potential measurement sites (40 rows and 100 sites per row). NDVI was generated using a two-dimensional Gaussian model. Range of the semi-variogram was set to 40m and average NDVI at 0.5 with a coefficient of variation at 30%. A single sample was selected on the basis of these NDVI data according to the protocol described in the sub-section below. Multiple yield data were then simulated using a linear model as shown in Eq. 13:

$$ \text{Yield} = \beta_0 + \beta_1 \times \text{NDVI} + \varepsilon_i \quad \text{with} \quad \varepsilon_i \sim N(0, \sigma) \quad (13) $$

$\beta_0$ and $\beta_1$ were respectively set at 300 g/vine and 1400 g/vine and $\sigma$ was chosen varying from 0 to 1500 with an increment of 25. 100 yield simulations were carried out for each $\sigma$ leading to 6000 yield simulations. These particular simulations with only one NDVI, one sample and multiple yields were used to make sure to study only errors caused by the inference step and not by the site selection step, matching hypothesis made in the theoretical part.

## Sampling protocol

Samples were selected using stratified sampling approaches based on NDVI quantiles clustering as described in Carillo et al. (2016). K, the set of all sites within the plot was divided into n groups of the same size using the quantile method. The groups thus formed a partition of K. The final sample was made by selecting randomly a single sampling site
from each of the *n* group with uniform probabilities. Once the sample was selected both estimators ($\hat{Y}_1$ and $\hat{Y}_2$) were calculated. For the commercial dataset, repetitions were made with 10000 target samples per plot. Sample size was set at *n* = 8 for all experiments.

## Results

Figure 1 shows the theoretical and experimental RMSE for both models applied to simulated data. RMSEs are represented as a function of $\sigma$ values (horizontal axis). Experimental RMSEs were computed from the samples of yield values while theoretical RMSEs were derived from NDVI values of the samples based on Eq. 10 and 11. Theoretical and experimental values are very close for both inference methods (Figure 1). These results validate the demonstration leading to the formalisation of bias and variance accounting for auxiliary data under the assumption of an existing linear model between the two variables. The two RMSEs show a noteworthy difference for the lowest sigma values (from 35g/vine and below). Following Eq. 12, the theoretical sigma value from which the mean estimator starts to outperform the model estimator is equal to 1018 here. However, the difference between both RMSE for such values is small and is not clearly observable.

<figure><img src="image_2.png" /><figcaption>Figure 1: Comparison of RMSE obtained on simulated data for both estimators (mean and model-based) for different values of $\sigma$. Figure 1a presents experimental RMSE (observed errors) and Figure 1b theoretical RMSE based on simulation parameters.</figcaption></figure>

Regarding real fields, results of tests presented in Figure 2, were carried out on the nine real fields of the data set, however, for the purpose of clarity, only results obtained on four fields (22, 65, 76, and 82) are presented. Results presented were chosen to encompass the diversity of the initial 9 plots in terms of correlation coefficient and ratio between $\sigma^2$ and $\beta_1^2$. Figure 2 presents yield estimation errors (%) observed with the mean (horizontal axis) and model-based (vertical-axis) estimators. The red line ($Y = X$) separates the samples for which inference by the mean provides better estimates (below the line) from those for where model-based inference provides better estimate. For readability reasons, the 10000 points are represented by a two-dimensional density.

The difference between both inference approaches appears very clear for plot 65 with 60% of the samples giving better results with a model-based approach showing point densities more concentrated in the lower right part of the scatter plot. Plot 65 presents
indeed the highest Yield-NDVI correlation (Table 1). This result is therefore consistent with Eq. 12. For the other plots, the differences between inference methods are less significant, as shown by the near symmetry around Y=X. Surprisingly, plot 22 with no Yield-NDVI correlation shows slightly better results with model inference. Plot 82, although showing an average yield-NDVI correlation, gives slightly better results with mean inference. Regarding the other plots (Plot 63,77,80, 88 and 104) model inference gives results very similar to that of plots 77 and 82.

<figure><img src="image_3.png" /><figcaption>Figure 2: Density of yield estimation errors with a model based (vertical axis) and a mean-based inference (horizontal axis) resulting from 10,000 samples with n = 8 measurement sites. The Y=X line is the zone where the same error is observed with both inference approaches.</figcaption></figure>

## Discussion

Model-based inference gives better results overall on these examples and confirms the previous results obtained by Carillo et al. (2016). As demonstrated, the differences between both approaches tend to be rather small as $\sigma$ values increase. The results show that the model-based inference approach performs better when correlation between auxiliary data and the variable of interest is very high (when the $\frac{\sigma^2}{\beta_1^2}$ ratios). In other cases,
differences between both inference approaches remained small with a slight advantage for the model-based inference. As a result, the model-based estimator should always be considered when it is possible. Some assumptions made in this study remains questionable like: the linearity of the model relating auxiliary data with the variable of interest and the independence of the model residuals (no spatial auto-correlation). For this latter, a practical recommendation would be to space out the measurement sites based on the range of the semi-variogram of auxiliary data. Further work could lead to a better understanding of the impact of autocorrelation when predicting with a linear model. Another perspective would be to simultaneously consider the sampling step and the model-based estimation step as addressed by Brus (2000) with probability sampling.

## Conclusions

This study has helped to explain and understand why the model-based inference method, which had until now been proposed empirically, gives better results when used for field sampling. The formalisation demonstrates that gains in accuracy are strongly related to the strength of the correlation between auxiliary data and the variable of interest as well as the representativeness of the observations. The derivation of theoretical errors makes it possible to identify which approach should be favoured according to the characteristics of the plot and the sample. The model-based inference presents better results than the mean-based inference and should be recommended to drive sampling when a relevant high spatial resolution auxiliary variable is available.

## Acknowledgements

This work was supported by the French National Research Agency under the Investments for the Future Program, referred as ANR-16-CONV-0004 (#Digitag).

## References

Araya-Alman, M., Leroux, C., Acevedo-Opazo, C., Guillaume, S., Valdés-Gómez, H., Verdugo-Vásquez, N. et al. (2019). A new localized sampling method to improve grape yield estimation of the current season using yield historical data. *Precision Agriculture* 20, 445–459.

Acevedo-Opazo, C., Tisseyre, B., Guillaume, S., & Ojeda, H. (2008). The potential of high spatial resolution information to define within-vineyard zones related to vine water status. *Precision Agriculture, 9*(5), 285–302.

Brus, D. J. (2000). Using regression models in design-based estimation of spatial means of soil properties. European Journal of Soil Science, 51, 159-172

Carrillo, E., Matese, A., Rousseau, J., & Tisseyre, B. (2016). Use of multi-spectral airborne imagery to improve yield sampling in viticulture. *Precision Agriculture, 17*(1), 74-92.

Oger, B., Vismara, P., Tisseyre, B., 2020. Combining target sampling with within field route-optimization to optimise on field yield estimation in viticulture. *Precision Agriculture*. doi.org/10.1007/s11119-020-09744-0

Oger, B., 2021. Model and auxiliary data for an accurate estimate of the field mean: Appendix. https://www.researchgate.net/publication/349589524 Model and auxiliary data for an accurate estimate of the field mean

Tillé, Y. (2011). Ten years of balanced sampling with the cube method: An appraisal. Survey Methodology, 2011, vol. 37, no. 2, p. 215-226