Site characterization

The study was conducted at four sites as summarized in Table 1. Deaton 6 and HF7 sites were located in the East of England (Lincolnshire) on reclaimed marshland with a shallow water table and high organic matter content. There was variation in soil physical and chemical properties across the field due to the presence of Roddons, historical features in drained marshland soils where silty clay soils follow the course of historical streams and waterways. Horse Foxhole and Buttery Hill were located in the West of England (Shropshire) on well-drained slightly stony, sandy loam soil subtended by weathered sandstone with low variation in soil nutrients across the field. Additional tuber sampling for TSD modelling (without soil analysis) was conducted at Crabtree Leasow, also located in the West of England with similar soil and weather conditions to Buttery Hill and Horse Foxhole. Thus, our study sites are representative of a broad range of soil conditions as is typical in potato production agriculture. Planting and field management was carried out by the respective farmer at each field. Consequently, land preparation was conducted similarly in all fields by ploughing at 30 cm depth followed by bed-forming at 90 cm between rows and destoning. Fertilizer was applied uniformly by broadcasting macro and micronutrients based on soil analysis as summarized in Table 1. Deaton 6 and HF7 were irrigated by drip irrigation while a hose reel irrigator was used at Horse Foxhole and Buttery Hill. All management practices were conducted uniformly across the field throughout the growing period.

Table 1. Summarized information of the study sites

1 = Number of samples collected. 2 = Main fertilizer forms were N for Nitrogen, P₂O₅, for Phosphorus, K₂O for Potassium, SO₃ for Sulphur and MgSO₄.7H₂O for Sulphur and Magnesium. 3 = Buttery Hill. 4 = Crabtree Leasow. 5 = Pentland Dell. 6 = Horse Foxhole.

Sampling design

A field survey was designed at each of the five fields using a model-based sampling approach to determine representative soil sampling locations that captured the variability in the field. The number of samples collected at all locations are as provided in Table 1. Soil sampling was done at all sites except Crabtree Leasow (thus not included in the soil analysis). Accordingly, the soil samples were collected at Buttery Hill (N = 23), HF7 (N = 24), Deaton 6 (N = 12) and Horse Foxhole (N = 23). Soil macronutrient quantities are known to correlate with organic matter content (Yang et al., Reference Yang, Luo, Lu, Schädel and Han2011), which in turn influences soil colour (Costa et al., Reference Costa, Giasson, da Silva, Coblinski and Tiecher2020). The Soil Brightness Index (SBI) as described by Mponela et al. (Reference Mponela, Snapp, Villamor, Tamene, Le and Borgemeister2020), was chosen as a substitute for spatially modelling the soil colour differences at each field and generating strata for a stratified random sampling design. At each field, the SBI was calculated at least 1 day per month for 3 months prior to crop emergence and then the average SBI was calculated. The SBI at each site was calculated using atmospherically corrected (Level-2A) satellite imagery of 10 m resolution acquired by the Sentinel-2 satellite, on manually inspected cloud-free days. Each field was delineated into three clusters by k-means clustering (k = 3) based on SBI to generate zones of relative homogeneity which formed the sampling strata. A sampling unit of 6 m by 6 m (36 m²) was chosen to cover the accuracy specification of the Garmin^TM eTrex 20 GPS receiver that was used for soil sampling. A grid of 36 m² quadrats was imposed across a rasterized SBI surface then random quadrats were drawn from each stratum.

Power analysis to determine the sample size for the survey was calculated to resolve SBI variability with a statistical power of 0.8 as recommended by Cohen (Reference Cohen1988). The effect size was estimated based on the expected within-field contrast between the k-means cluster with the lowest SBI value (dark soil) and the one with the highest SBI value (bright soil) in the sample. The effect size was therefore calculated as the difference between the mean SBI of dark soils and the combined mean of the medium and light clusters divided by the standard deviation of the entire dataset. R v4.0.2 (R Core Team, 2019) was used to calculate the sample size using the ‘pwr.t.test’ function from the ‘pwr’ package (Champely et al., Reference Champely, Ekstrom, Dalgaard, Gill, Weibelzahl, Anandkumar, Ford, Volcic and De Rosario2018). Once the sample size of each cluster was determined, sampling locations were selected by randomly selecting quadrats from the grid imposed on the SBI raster. The quadrats were georeferenced and assigned with unique identifier codes then exported as a GPx file into the GPS receiver for tracking during the soil and yield sampling. All raster analysis steps were performed using ArcGIS (ESRI, 2018).

Soil and nutrient analysis

Soil sampling was conducted before planting, but after fertilizer incorporation. Soils were sampled using a 30 cm auger. Samples were collected in triplicate from each quadrat and mixed to form a composite sample. The soil samples were then air-dried at 30°C for 72 h. After air drying, the soils were ground and sieved (<2 mm) prior to analysis. The percentages of sand, silt and clay in each soil sample were determined using the sedimentation method (Jackson et al., Reference Jackson, Farrington and Henderson1986). Hydrogen peroxide was used to oxidize organic matter, after which the particle size distribution was determined through sieving and sedimentation.

Soil samples were analysed for N, C and S using the Dumas method (Kirsten and Hesselius, Reference Kirsten and Hesselius1983). Air-dried soils (0.25 g for N and 0.15 g for C and S) were passed through a furnace at 1000°C in the presence of oxygen. The oxidized gases were then detected and measured using a thermal conductivity cell. The Olsen method (Page, Reference Page1982) was used to estimate available P in the soil. Sodium bicarbonate was used to extract P from the soil into solution and form phosphomolybdate after reaction with ammonium molybdate. The phosphomolybdate was reduced by ascorbic acid to form a blue complex whose concentration was measured spectrophotometrically at 880 nm. Concentrations of K and Mg were determined by flame photometry using ammonium nitrate as an extractant as described in (Jackson et al., Reference Jackson, Farrington and Henderson1986).

Yield data collection

At every sampling location, the number of plants within a one-metre row was counted and recorded as a measure of plant density. At harvest, all the plants in the one-metre length were carefully uprooted with a spade and the number of main stems was counted as a measure of stem density. The excavation was carried out carefully to minimize any loss of tubers. All tubers were separated from their stolon and stored for further processing.

The number of tubers at each sampling point was counted. At all locations except for Horse Foxhole, all tubers with a transversal diameter greater than or equal to 25 mm were graded using potato sizing squares into 10 mm size grades up to 65 mm, with all tubers over 65 mm in diameter placed in one bin. All tubers under 25 mm diameter were binned into one grade. After grading, the tuber number in each category was counted and its weight was determined at 0.01 g accuracy. For Horse Foxhole, the tubers were separated into typical commercial grading of 0–25 mm, 25 mm-45 mm, 45–65 mm and greater than 65 mm.

Data analysis

TSD was modelled using the Gaussian (Travis, Reference Travis1987; Aliche et al., Reference Aliche, Oortwijn, Theeuwen, Bachem, van Eck, Visser and van der Linden2019), Gamma (Aliche et al., Reference Aliche, Oortwijn, Theeuwen, Bachem, van Eck, Visser and van der Linden2019) and Weibull (Nemecek et al., 1996; Bussan et al., Reference Bussan, Mitchell, Copas and Drilias2007) distributions. The Gaussian distribution is the most widely adopted of the three distributions, with the mode equalling the mean of the distribution and its parameter estimation is simple and intuitive as described in Travis (Reference Travis1987). For the Gamma distribution, the probability density is as described in Aliche et al. (Reference Aliche, Oortwijn, Theeuwen, Bachem, van Eck, Visser and van der Linden2019), with respect to TSD in potatoes. The Weibull distribution density function was described by Nemecek et al. (Reference Nemecek, Derron, Roth and Fischlin1996). The Gamma and Weibull distributions are considered to be flexible to right-skewed data, common in potato TSD (Nemecek et al., Reference Nemecek, Derron, Roth and Fischlin1996; Bussan et al., Reference Bussan, Mitchell, Copas and Drilias2007; Aliche et al., Reference Aliche, Oortwijn, Theeuwen, Bachem, van Eck, Visser and van der Linden2019), and hence they can potentially more accurately estimate modal tuber size than the Gaussian distribution. Love and Thompson-Johns (Reference Love and Thompson-Johns1999) demonstrated that agronomic practices such as wider plant spacing can also lead to left-skewed TSD, making the Weibull distribution an ideal candidate for modelling the shape of the distribution flexibly. Figure 1 illustrates the effect of changing the shape parameter on a conceptual Weibull probability density curve. Generally, the distribution resolves to an approximately symmetrical distribution with a shape value of 3.6 and approaches right-skewness or left-skewness above or below that respectively (Lai et al., Reference Lai, Murthy, Xie and Pham2006). Thus, the Weibull distribution is ideal in this application to model both left and right-skewed data.

Fig. 1. The effect of changing the shape parameter on a conceptual Weibull probability density curve.

The MLE approach was used to obtain estimates for the mean and standard deviation of the Gaussian distribution, and the scale and shape parameters for the Weibull and Gamma distributions. The ‘fitdstrplus’ R package (Delignette-Muller and Dutang, Reference Delignette-Muller and Dutang2015) was used to obtain the parameter estimations for all three distributions from the right-censored interval potato TSD data of the current study. Accordingly, the likelihood of each parameter θ in each distribution was fitted in Eqn (1):

(1)$$\eqalign{& L( \theta ) = \mathop \prod \limits_{\,j = 1}^{N_{{\rm left}\,C}} F( {x_j^{{\rm upper}} {\rm \vert }\theta } ) \times \mathop \prod \limits_{k = 1}^{N_{{\rm right}\,C}} 1-F( {x_k^{{\rm lower}} {\rm \vert }\theta } ) \;\cr & \;\,\,\,\,\,\,\,\,\,\,\,\,\,\times \mathop \prod \limits_{m = 1}^{N_{{\rm int}\,C}} ( {F( {x_m^{{\rm upper}} {\rm \vert }\theta } ) -F( {x_j^{{\rm lower}} {\rm \vert }\theta } ) } ) } $$

Using the cumulative distributions (F), where $x_j^{{\rm upper}}$ represented the upper values defining the N ^left C left-censored observations, $x_k^{{\rm lower}}$ represented the lower values defining the N ^right C right-censored observations, $x_m^{{\rm lower}}$ and $x_m^{{\rm upper}}$ represented the intervals defining the N ^int C interval-censored observations (Delignette-Muller and Dutang, Reference Delignette-Muller and Dutang2015).

The Weibull distribution parameters were also estimated using the percentiles method (Dubey, Reference Dubey1967). The motivation for including this method was the simple closed-form procedure for calculating parameter estimates, which potentially makes it easier for non-statisticians to calculate the parameters and make an inference. The adoption of the Weibull MLE approach in TSD research is very low, with recent studies still using Gaussian mean estimates (Buhrig et al., Reference Buhrig, Thornton, Olsen, Morishita and McIntosh2015) or simply relative frequencies (Andrade et al., Reference Andrade, da Silva, Pesantes, Christensen and Zotarelli2021; van Dijk et al., Reference van Dijk, Lommen, de Vries, Kacheyo and Struik2021) to describe TSD and its changes with respect to agronomic treatments. This is attributable to non-familiarity with the approach due to the ubiquity of least-squares-based procedures and simplicity of Gaussian parameter estimates, which still provide useful models as reported in Travis (Reference Travis1987). Estimating TSD using the Weibull percentiles approach provides a simpler alternative to the poorly adopted MLE approach, contingent upon achieving a comparable empirical accuracy between the two approaches, hence their comparison in this study. Accordingly, the shape parameter ($\hat{\beta }_{{\rm Weibull}}$) of the Weibull distribution for the percentiles method was estimated by linearizing the cumulative distribution function at two different discrete diameters (e.g. 45 and 65 mm) then combining the two equations to solve for $\hat{\beta }_{{\rm Weibull}}$ as in Eqn (2):

(2)$$\hat{\beta }_{{\rm Weibull}} = ( {\ln x_i-\ln x_j} ) ^{{-}1} \times \left\{{{\rm \;ln}\left[{\displaystyle{{\ln ( {1-f( {x_j} ) } ) } \over {\ln ( {1-f( {x_i} ) } ) }}} \right]\;} \right\}$$

The scale parameter $\hat{\alpha }_{{\rm Weibull}}$ was then calculated using $\hat{\beta }_{{\rm Weibull}}$ and the cumulative density at one known quantile as in Eqn (3):

(3)$$\hat{\alpha }_{{\rm Weibull}} = \displaystyle{{x_i} \over {{[ {-\ln ( {1-f( {x_i} ) } ) } ] }^{{1 \over {{\hat{\beta }}_{{\rm Weibull}}}}}}}$$

where f(x_i..j) represented the cumulative probability of tuber number or weight at x_i or x_j tuber diameter and x_i…j where the chosen i ^th or j ^th discrete diameter of a tuber

In each sample, the estimated parameters of the four distribution-fitting approaches (Gaussian, Gamma, Weibull with MLE and Weibull with percentiles approach) were used to estimate the tuber numbers in each original discrete size grade, creating fully specified ordered discrete distributions. The estimation was done by multiplying the estimated probability of each size grade (the difference between the upper limit and lower limit cumulative probabilities of each size grade) by the observed total tuber number. The Weibull percentiles approach was fitted at the percentiles corresponding to 45 and 65 mm tuber sizes, selected because this represents the marketable range for maincrop potatoes. In this case, the percentiles were the cumulative proportions of the tuber number or weight at 45 or 65 mm, and were therefore different for every sample. The logarithm of the relative likelihood (hereinafter referred to as log relative likelihood) of the predicted discrete distribution (relative to the actual discrete distribution) was then calculated as illustrated in Lindsey (Reference Lindsey1974). The likelihood of the discrete distribution was calculated as the likelihood of a fully specified multinomial distribution in Eqn (4), where P is the probability of the j ^th category and n is the frequency of the j ^th category.

(4)$$L( P ) = \mathop \prod \limits_j P_j^{n_j} $$

The log relative likelihood was then used as a ranking index of the plausibility of the models, with the most plausible model having the highest likelihood. The log relative likelihood was calculated as in Eqn (5):

(5)$$L_{RR} = \mathop \sum \nolimits^ n_{\,j\;}\left[{\ln \left({{\'{P}}_j} \right)-\;\ln ( {P_j} ) } \right]$$

where n_j is the observed frequency of a category, $\'{P}_j$ is the probability of the predicted distribution and P _j is the likelihood of the observed distribution.

Uncertainty in the log relative likelihood was assessed using 95% confidence intervals. Fisher's information was used to estimate the variability of the estimate. The Fisher's information was derived by negating the expected value of the second derivative of the log-likelihood equation (Ly et al., Reference Ly, Marsman, Verhagen, Grasman and Wagenmakers2017), which for the multinomial distribution used was defined as in Eqn (6):

(6)$$I( {\theta_j} ) = -E\left({-\displaystyle{{n_j} \over {P_j^2 }}} \right)$$

where θ _j is the log relative likelihood of the j ^th category, n_j is the frequency of the j ^th category and P_i is the probability of the j ^th category. The reciprocal of the square of the fisher's information was used as the standard error component of the confidence interval formula.

The main purpose of fitting the theoretical continuous distributions to the discretely measured tuber size fractions is to maximize the accuracy of the estimate of the modal tuber size grade of the distribution, which is the marketable component of the production. This falls within the 45–65 mm size bands for typical maincrop varieties. Therefore, the main purpose of the modelling presented here was to determine which of the distributions-fitting approaches maximized the likelihood of this tuber size fraction, using the observed tuber numbers and weights for this size grade. The log relative likelihood estimates of the 45–65 mm size band were therefore compared for the four models. Additionally, the Root Mean Square Error (RMSE) of prediction for the frequencies or weights in the 45–65 mm size fraction were compared. The RMSE was calculated as in Eqn (7):

(7)$${\rm RMSE} = \sqrt {\displaystyle{1 \over N}\mathop \sum \limits_{i = 1}^N {( {y_{i}-\'{y}_i} )}^2} $$

where N is the number of observations, y_i is the observed value and ý_i is the predicted value. The modal tuber size of the most plausible model within the 45–65 mm size fraction was considered the best estimate of the distribution's model. For comparison purposes, the most plausible distribution fitted using the benchmark maximum likelihood approach was used to compare the modal tuber size predictions of the other models using RMSE. The RMSE and log relative likelihood were sample-specific, providing measures of the goodness of fit of the parameter estimates to the observed TSDs without generalizing uncertainties expected at out-of-sample observations. The TSD parameters are only valid for a collected sample, but within-sample hold-out methods (holding out a random size profile within the sample as a test dataset during model fitting) were considered impractical in this case as they would change the TSD of the training set. The mode of the Gaussian distribution was considered to equal the mean. The mode for the Gamma distribution was calculated in Eqn (8) and the mode for the Weibull distribution was calculated in Eqn (9), as follows:

(8)$$\widehat{{Mo}}_{{\rm Gamma}} = \hat{\gamma }_{{\rm Gamma}} \times ( {{\hat{\alpha }}_{{\rm Gamma}}-1} ) $$

(9)$$\widehat{{Mo}}_{{\rm Weibull}} = \hat{\alpha }_{{\rm Weibull}} \times \left({1-\;\displaystyle{1 \over {{\hat{\beta }}_{{\rm Weibull}}}}} \right)^{{1 \over {{\hat{\beta }}_{{\rm Weibull}}}}}$$

where $\widehat{{Mo}}_{{\rm Gamma}}$ is the estimate of the mode of the Gamma distribution, $\hat{\gamma }_{{\rm Gamma}}$ is the scale parameter of the Gamma distribution, $\hat{\alpha }_{{\rm Gamma}}$ is the shape parameter of the Gamma distribution, $\widehat{{Mo}}_{{\rm Weibull}}\;$ is the estimate of the mode of the Weibull distribution, $\hat{\propto }_{{\rm Weibull}}$ is the scale of the Weibull distribution and $\hat{\beta }_{{\rm Weibull}}$ is the shape parameter of the Weibull distribution.