Evaluating observed versus predicted forest biomass: R-squared, index of agreement or maximal information coefficient?

ABSTRACT

The accurate prediction of forest above-ground biomass is nowadays key to implementing climate change mitigation policies, such as reducing emissions from deforestation and forest degradation. In this context, the coefficient of determination ($R^2$) is widely used as a means of evaluating the proportion of variance in the dependent variable explained by a model. However, the validity of $R^2$ for comparing observed versus predicted values has been challenged in the presence of bias, for instance in remote sensing predictions of forest biomass. We tested suitable alternatives, e.g. the index of agreement ($d$) and the maximal information coefficient ($MIC$). Our results show that $d$ renders systematically higher values than $R^2$, and may easily lead to regarding as reliable models which included an unrealistic amount of predictors. Results seemed better for $MIC$, although $MIC$ favoured local clustering of predictions, whether or not they corresponded to the observations. Moreover, $R^2$ was more sensitive to the use of cross-validation than $d$ or $MIC$, and more robust against overfitted models. Therefore, we discourage the use of statistical measures alternative to $R^2$ for evaluating model predictions versus observed values, at least in the context of assessing the reliability of modelled biomass predictions using remote sensing. For those who consider $d$ to be conceptually superior to $R^2$, we suggest using its square $d^2$, in order to be more analogous to $R^2$ and hence facilitate comparison across studies.

KEYWORDS:

• Model assessment
• overfitting
• biomass
• LIDAR

Introduction

Obtaining accurate predictions of forest above-ground biomass ($AGB$) is nowadays key to estimating total carbon stock and implementing climate change mitigation policies, such as reducing emissions from deforestation and forest degradation (REDD) (Agrawal, Nepstad, & Chhatre, 2011; Eggleston, Buendia, Miwa, Ngara, & Tanabe, 2006; UNFCCC, 2014). To facilitate REDD implementation, tree allometry models are currently being developed in order to avoid the need for destructive sampling (Basuki, van Laake, Skidmore, & Hussin, 2009; Chave et al., 2005, 2014; Cuny et al., 2015; Iais et al., 2011; Krainovic, Almeida, & Sampaio, 2017; Marshall et al., 2012; Sawadogo et al., 2010). Methods combining field measurements and remote sensing technologies can provide robust and cost-effective means for reliably predicting $AGB$ from large areas of forests (Asner, 2011; Clark & Kellner, 2012; Goetz & Dubayah, 2014). For this reason, earth observation technologies using remote sensors are also nowadays extensively used to obtain predictions over large areas (Adnan et al., 2019; Asner & Mascaro, 2014; Bottalico et al., 2017; Hansen et al., 2013; Næsset, 2004; Tyukavina et al., 2017). Methods developed for $AGB$ estimation made use of the predictive power of spectral sensors (e.g. Franco-Lopez, Ek, & Bauer, 2001), LIDAR (e.g. Coomes et al., 2017; d’Oliveira, Reutebuch, McGaughey, & Andersen, 2012; Domingo, Lamelas, Montealegre, García-Martín, & de la Riva, 2018; García, Riaño, Chuvieco, & Danson, 2010; Montealegre, Lamelas-Gracia, García-Martín, de la Riva-Fernández, & Escribano-Bernal, 2017) or combinations of both (e.g. Asner, 2009; Bright, Hicke, & Hudak, 2012; Egberth et al., 2017; Estornell, Ruiz, Velázquez-Martí, & Hermosilla, 2012; Hernando et al., 2019). Zolkos et al., (2013) compiled a comprehensive revision of studies obtaining remote sensing predictions of $AGB$ by calibrating regional models with field plot information.

The correct evaluation of $AGB$ models is key to REDD, as errors may propagate when $AGB$ estimations are upscaled (Chave et al., 2004; Molto, Rossi, & Blanc, 2013; Valbuena et al., 2016). However, there is at the moment a critical lack of consensus on good practices for predicting forest $AGB$, both in the determination of allometric biomass equations (Sileshi, 2014; Temesgen, Affleck, Poudel, Gray, & Sessions, 2015) and in the evaluation of $AGB$ predictions using remote sensing methods (Loew et al., 2017; Valbuena et al., 2017a). The large variety of modelling methods and approaches for evaluating their accuracy complicates meta-analyses comparing various approaches (Loew et al., 2017; Zolkos et al., 2013). In a recent article (Valbuena et al., 2017a), we focused on the evaluation of final $AGB$ predictions from remote sensing with measures of accuracy and precision. We recommended the use of Piñeiro, Perelman, Guerschman, and Paruelo (2008) hypothesis test, considered the usability of Theil’s (1958) partial inequality coefficients, and suggested means for preventing overfitting to the sample (Valbuena et al., 2017a). In this paper, we focus on the evaluation of agreement between the $AGB$ observed and the $AGB$ predicted by the model.

The agreement between the observed and predicted values is an important indicator in model evaluation, and the best approach for its assessment is a subject of current discussion (Duveiller, Fasbender, & Meroni, 2016; Loew et al., 2017; Simon & Tibshirani, 2011; Willmott, Robeson, & Matsuura, 2012). The coefficient of determination ($R^2$) is widely used as the most suitable statistic for describing the agreement of model predictions with those observed empirically. $R^2$ expresses the proportion of variance in the dependent variable explained by a model. However, the use of $R^2$ in performance evaluation for predictive methods has been criticized for not being necessarily related to the accuracy of predictions (e.g. Fox, 1981; Paruelo, Jobbágy, Sala, Lauenroth, & Burke, 1998; Willmott, 1982). Alternatives have been proposed, such as the index of agreement ($d$), which can be interpreted as the relative prediction error (Willmott, 1981; Willmott & Wicks, 1980). Ever since, it has been popular in a variety of fields such as hydrology (Legates & McCabe, 1999; López-Moreno, Latron, & Lehmann, 2010), agriculture (Nendel et al., 2013), neurology (Ganpule et al., 2017), meteorology (Aschonitis et al., 2017; Morsy, El-Sayed, & Ouda, 2016), forest ecology (Ibrom et al., 2007; Ward, Bell, Clark, & Ram, 2013), or climate change (Bring & Destouni, 2014; Gaitan, Hsieh, & Cannon, 2014; Oyler, Ballantyne, Jencso, Sweet, & Running, 2015). The use of Willmott’s $d$ in remote sensing literature has, however, been marginal (Almeida et al., 2016; García et al., 2010; Grzegozewski, Johann, Uribe-Opazo, Mercante, & Coutinho, 2016; Wachholz de Souza, Mercante, Johann, Camargo Lamparelli, & Uribe-Opazo, 2015; Yebra & Chuvieco, 2009). Since many of these authors recommended the use of $d$ above $R^2$, this research was devoted to study the convenience of using $d$ for reporting the agreement between remote sensing $AGB$ predictions and their observed values. In addition to $d$, we postulated that using the recently developed maximal information coefficient ($MIC$) (Reshef et al., 2011) as a means for evaluating observed versus predicted values may be advantageous. $MIC$ was conceived as a non-parametric alternative to correlation, in the context of evaluating relationships between explanatory and response variables (Speed, 2011). Despite being very recent, it has already been found useful in a wide variety of contexts, such as microbiology (Thomas, Bordron, Eveillard, & Michel, 2017), medicine (Vallières et al., 2017), big data analysis (Chen & Yang, 2016), and also remote sensing (Görgens, Valbuena, & Rodríguez, 2017). Although it has been suggested that the robustness of this statistic against noise should be studied in terms of measuring predictive accuracy of models (Loew et al., 2017; Murrell, Murrell, & Murrell, 2014), to our knowledge no previous study has considered the suitability of $MIC$ in the context of evaluating the agreement between modelled predictions of $AGB$ and their corresponding observed values.

In this article, we consider the adequacy of using alternatives to the $R^2$ for evaluating the agreement between observed and predicted values. The alternatives considered were $d$ and $MIC$, and the study was conceived as a further aspect of prediction accuracy to be revised in the context of $AGB$ predictions from remote sensing, additional to those outlined in Valbuena et al. (2017a). We, therefore, followed up the same experimental design to allow direct comparison and contrasting against additional measures of accuracy: results are presented for unviable modelling approaches alongside correct ones, in order to highlight whether the alternative measures of agreement reliably make a difference among them.

Material and methods

Field and remote sensing datasets

The study employed data from $n=37$ field plots composed of two concentric circles of radii 10 m and 20 m (Valbuena et al., 2013b), which were obtained in the summer 2006 at the Pinus sylvestris-dominated forests of Valsaín (Spain, approx. lat.: 41°04ʹ N, lon.: 4°09ʹ W; 1.3–1.5 km a.s.l.). Within the inner circle, all trees were measured, including seedlings and saplings, whereas the outer circle only were measured trees with diameter at breast height (dbh, cm) above 10 cm. Plot expansion factors were used to expand the forest information sampled within the inner plot to the outer plot (Valbuena, Packalen, Mehtätalo, Garcia-Abril, & Maltamo, 2013c). Treetop heights (h, m) were determined for every individual tree using a Vertex III Hypsometer (Haglof, Sweden). Plot centres were staked out using a HiPer-Pro (Topcon, California) receiver set at 1–2 m above the ground for differentially corrected global navigation satellite systems (GNSS) positioning (Valbuena, Mauro, Rodriguez-Solano, & Manzanera, 2012). Tree allometry (Montero, Ruiz-Peinado, & Muñoz, 2005) was used to calculate plot $AGB$ (Mg·ha⁻¹) from these field measurements. While Ruiz-Peinado, Del Rio, and Montero (2011) improved Montero et al.’s (2005) models by including a property of additivity among biomass components, we considered such property would not propagate to the remote sensing predictions with the methodology employed in this article (Hernando et al., 2019). The remote sensing survey was carried out in September 2006 using a LIDAR system (ALS50-II from Leica Geosystems, Switzerland) and a multispectral camera (DMC from Zeiss-Intergraph, Germany). The LIDAR pulse density was 1.15 pulses·m⁻², and approximate footprint diameter was 0.5 m at the ground. Using Terrascan (Terrasolid, Finland) returns were classified and those at the ground were interpolated into a digital terrain model (Valbuena, Mauro, Arjonilla, & Manzanera, 2011). DMC original multispectral bands had an approximate spatial resolution of 60 cm. The fusion of both remote sensors’ information was done with an assured perfect fit via a back-projection data fusion algorithm (Valbuena et al., 2013a, 2011). The correspondence between the field and remote sensing data was guaranteed with survey-grade and differentially corrected GPS positioning (Valbuena et al., 2012). Remote sensing predictors were computed from statistical descriptors of the distributions of heights from the LIDAR data, and normalised difference vegetation index (NDVI) values from the multispectral sensor (Manzanera et al., 2016). More details about the data employed and processing steps can be further scrutinized from Valbuena et al. (2017b), Appendix A, and the above-listed references. Readers interested in comparable studies in Mediterranean pine forest ecosystems in Spain may refer to a recent review in Gómez et al. (2019), or individual studies using LIDAR (Adnan et al., 2019; Bottalico et al., 2017; Estornell, Ruiz, Velázquez-Martí, & Hermosilla, 2011; García et al., 2010; Gonzalez-Ferreiro, Dieguez-Aranda, & Miranda, 2012; González-Olabarria, Rodríguez, Fernández-Landa, & Mola-Yudego, 2012; Montealegre, Lamelas, de la Riva, García-Martín, & Escribano, 2016; Montealegre et al., 2017; Valbuena et al., 2013c), in combination with multispectral sensors (Estornell et al., 2012; Hernando et al., 2019; Manzanera et al., 2016; Valbuena et al., 2013a, 2017b, 2011), or comparing methods for selection of predictor variables (Garcia-Gutierrez et al., 2014; Valbuena et al., 2017a) or prediction methods (García-Gutiérrez, Martínez-Álvarez, Troncoso, & Riquelme, 2015; Guerra-Hernández et al., 2016; Domingo, Lamelas-Gracia, Montealegre-Gracia, & de la Riva-Fernández, 2017; Domingo et al., 2018; Valbuena et al., 2016).

Modelling alternatives compared

The values of $R^2$, $d$ and $MIC$ calculated in this study were obtained for the same modelling alternatives reported in Valbuena et al. (2017a), in order to allow direct comparisons according to the conclusions reached in that study. These consisted in predictions of plot-level $AGB$ using parametric models and a non-parametric method: best-subset, step-wise and most similar neighbours, all of them including ‘restricted’ versions according to Piñeiro et al. (2008) and Valbuena et al. (2017a). All computations were carried out using the R statistical environment (version 3.3.1; R Development Core Team, 2016). The parametric models were power models adjusted in their linear form, i.e. using a natural logarithm transform of the response variable (e.g. Asner & Mascaro, 2014; Hudak et al., 2006; Næsset, 2002), which were bias-corrected when transformed back to the original scale (Baskerville, 1972; Sprugel, 1983). A first alternative, here forth denominated “best-subset” (Hudak et al., 2006; Miller, 1984), selected model predictors via exhaustive evaluation of all possible combinations using package “leaps” of R (Lumley & Miller, 2009), and minimization of Mallows’ Cp (Mallows, 1973) as predictor selection criteria. Next, there was a modification of this alternative here forth denominated “best-subset restricted overfitting” (Valbuena et al., 2017a). It consisted in further constraining the “best-subset” result by restricting the inflation of sum of squares to a 10% in the cross-validation, to avoid models overfitted to the sample (Ehrenberg, 1982; Lipovetsky, 2013; Weisberg, 1985), plus positive results in Piñeiro et al.’s (2008) hypothesis test to the observed versus predicted fit. The next alternative employed the same power model but variables were selected according to a step-wise procedure (Næsset, 2002; Weisberg, 1985) and the corrected Akaike Information Criterion (AIC) (Burnham & Anderson, 2002; Sugiura, 1978), as implemented in the function “stepAIC” of R (Venables & Ripley, 2002). We denominate this method “step-wise”, and also presented a “step-wise restricted overfitting” modification using the hypothesis test and the limiting criterion of 10% inflation in sums of squares (Valbuena et al., 2017a).

The other alternative employed was the non-parametric prediction. This was based on the most similar neighbour (MSN) method, one type of machine learning approach among the so-called group of nearest neighbour methods (Franco-Lopez et al., 2001; McInerney, Suárez, Valbuena, & Nieuwenhuis, 2010; McRoberts, Nelson, & Wendt, 2002). MSN predictions were carried out using the “yaImpute” package of R (Crookston & Finley, 2007). We set the algorithm to use inverse distance weighting averaging of three neighbours, based on previous experience (Almeida et al., 2016; Eskelson et al., 2009; Valbuena, Vauhkonen, Packalén, Pitkanen, & Maltamo, 2014). In the case of MSN, we employed canonical regression analysis (Cohen, Maiersperger, Gower, & Turner, 2003; Manzanera et al., 2016) to recursively select the predictors. The final selection criterion was the root mean squared error minimization with the same above-mentioned limiting restrictions of Piñeiro et al.’s (2008) hypothesis tests and avoiding overfitting (Valbuena et al., 2017a). Additionally, we calculated results also for a wide range of number of predictors $p=1\dots 30$, with the intention to emphasize whether the chosen statistics of agreement between observed and predicted would show any differences for unrealistically low $n/p$ ratios.

Evaluating the agreement between observed and predicted $AGB$

For all the modelling alternatives considered, we computed predictions using two options: first, (1) the model fit employing the entire dataset of $i=1\dots n$ forest plots used to train the model (i.e. model residuals without external validation), and also (2) a leave-one-out cross-validation removing each case $i$ from the training data as a prior step to the whole modelling process. In results, these two options are denoted using superscripts (or subscripts) $fit$ and $cv$, respectively. Thus, the training dataset used was the vector of measured $AGB$ values at plot-level $O=ob{s}_i$ where $ob{s}_i=ob{s}_1,\dots ,ob{s}_n$, from which we obtained the vector of modelled predictions $P=pr{e}_i$ where $pr{e}_i=pr{e}_1,\dots ,pr{e}_n$. These predictions where either those fitted in model training $\left(pr{e}_i^{fit}\right)$, or those calculated using cross-validation $\left(pr{e}_i^{cv}\right)$, since both were employed to compute the range of statistical measures considered for describing the agreement between observed and predicted. Therefore, $pr{e}_i$ denotes either $pr{e}_i^{fit}$ or $pr{e}_i^{cv}$ in the following equations detailing the measurements of agreement considered:

(1) The coefficient of determination, which shows the ratio of the sum of squared residuals to the total sum of squares:

(1) $\begin{array}{r}{R}^2=1-\left[{\sum }_{i=1}^n{\left(pr{e}_i-ob{s}_i\right)}^2/{\sum }_{i=1}^n{\left(ob{s}_i-\overline{obs}\right)}^2\right]\end{array}$(1)

(2) Willmott’s index of agreement, which was originally suggested to be (Willmott & Wicks, 1980):

(2) $\begin{array}{rl}& d=1-\left[{\sum }_{i=1}^n{\left(pr{e}_i-ob{s}_i\right)}^2\right./\\ & \left.{\sum }_{i=1}^n{\left(\left|pr{e}_i-\overline{obs}\right|+\left|ob{s}_i-\overline{obs}\right|\right)}^2\right]\end{array}$(2)

It was later suggested, however, that squaring the residuals may be an inconvenient approach (Willmott et al., 1985), and therefore they suggested the following refinement:

(3) $\begin{array}{rl}& d_1=1-\left[{\sum }_{i=1}^n\left|pr{e}_i-ob{s}_i\right|/\right.\\ & \left.{\sum }_{i=1}^n\left(\left|pr{e}_i-\overline{obs}\right|+\left|ob{s}_i-\overline{obs}\right|\right)\right]\end{array}$(3)

And also another further modification was suggested more recently (Willmott et al., 2012):

(4) $d_r=1-\left[{\sum }_{i=1}^n\left|pr{e}_i-ob{s}_i\right|/2{\sum }_{i=1}^n\left|ob{s}_i-\overline{obs}\right|\right]$(4)

The full version of the refined index of agreement $d_r$ consists of inverting the fraction and subtracting one from it (Willmott et al., 2012: Equation (5)(5) $\begin{array}{r}MI(P,O)=\sum _{P,O}p(P,O)\cdot log[p(P,O)/p(P)p(O)]\end{array}.$(5) ), in order to accommodate poorly performing models. However, this was not necessary in our case since the agreement between observed and predicted was higher for all the modelling alternatives considered.

(3) The maximal information coefficient, which measures the entropy of the relationship (Reshef et al., 2011; Speed, 2011). It is based on the naïve mutual information $MI\left(P,O\right)$ (Linfoot, 1957). The computation of $MIC$ involves binning $P$ and $O$, calculating relative abundances $p\left(P,O\right)$ at each grid of the resulting cell, and consider the overall uncertainty in their relationship using Shannon’s (1948) entropy:

(5) $\begin{array}{r}MI(P,O)=\sum _{P,O}p(P,O)\cdot log[p(P,O)/p(P)p(O)]\end{array}.$(5)

In Reshef et al.’s (2011) algorithm, $MIC$ results from the bin size that maximizes $MI\left(P,O\right)$:

(6) $MIC=\underset{P,O_{total}<B}{max}\left\{MI\left(P,O\right)/log\left[min\left(P,O\right)\right]\right\}$(6)

We employed the R implementation for $MIC$ computation available in package “minerva” (Filosi, Visintainer, & Albanese, 2014).

While the calculation of $d$ resembles that of a sample Pearson’s correlation coefficient ($r$), $MIC$ is customarily interpreted as a non-parametric version of a correlation coefficient (Görgens et al., 2017). Also, Simon and Tibshirani (2011) recommended the use of Brownian distance correlation (Székely & Rizzo, 2009) above $MIC$. For these reasons, we also calculated from our dataset the following statistics for additional discussion on the most convenient approach for assessing the agreement between observed versus predicted in model evaluation:

(4) The coefficient of correlation between observed and predicted values, calculated from Pearson’s product moment correlation:

(7) $\begin{array}{rl}& r=\left[{\sum }_{i=1}^n\left(pr{e}_i-\overline{pre}\right)\left(ob{s}_i-\overline{obs}\right)\right]/\\ & \left[\sqrt{{\sum }_{i=1}^n{\left(pr{e}_i-\overline{pre}\right)}^2}\sqrt{{\sum }_{i=1}^n{\left(ob{s}_i-\overline{obs}\right)}^2}\right]\end{array}$(7)

(5) The distance correlation ($dCor$) is an approach to evaluating the relationship between two variables based on their energy (Székely & Rizzo, 2017), i.e., in this case, the distances between $\left(P,O\right)$ data in the metric space. We employed the distance correlation statistic implemented in the package “energy” of R (Rizzo & Szekely, 2017).

All the alternatives for evaluating the degree of agreement between observed and predicted were calculated both from the model fit predictions $\left(P^{fit}\right)$ and the cross-validated versions $\left(P^{cv}\right)$: $R_{fit}^2$, $R_{cv}^2$ (Equation (1)(1) $\begin{array}{r}{R}^2=1-\left[{\sum }_{i=1}^n{\left(pr{e}_i-ob{s}_i\right)}^2/{\sum }_{i=1}^n{\left(ob{s}_i-\overline{obs}\right)}^2\right]\end{array}$(1) ), $d^{fit}$, $d^{cv}$ (Equation (2)(2) $\begin{array}{rl}& d=1-\left[{\sum }_{i=1}^n{\left(pr{e}_i-ob{s}_i\right)}^2\right./\\ & \left.{\sum }_{i=1}^n{\left(\left|pr{e}_i-\overline{obs}\right|+\left|ob{s}_i-\overline{obs}\right|\right)}^2\right]\end{array}$(2) ), and $MI{C}^{fit}$, $MI{C}^{cv}$ (Equation (5)(5) $\begin{array}{r}MI(P,O)=\sum _{P,O}p(P,O)\cdot log[p(P,O)/p(P)p(O)]\end{array}.$(5) ); plus the refined versions of Willmott’s agreement $d_1^{fit}$, $d_1^{cv}$ (Equation (3)(3) $\begin{array}{rl}& d_1=1-\left[{\sum }_{i=1}^n\left|pr{e}_i-ob{s}_i\right|/\right.\\ & \left.{\sum }_{i=1}^n\left(\left|pr{e}_i-\overline{obs}\right|+\left|ob{s}_i-\overline{obs}\right|\right)\right]\end{array}$(3) ), and $d_r^{fit}$, $d_r^{cv}$ (Equation (4)(4) $d_r=1-\left[{\sum }_{i=1}^n\left|pr{e}_i-ob{s}_i\right|/2{\sum }_{i=1}^n\left|ob{s}_i-\overline{obs}\right|\right]$(4) ), and the additional statistics of correlation $r^{fit}$, $r^{cv}$, $dCo{r}^{fit}$and $dCo{r}^{cv}$. The relative merits of each of the proposed statistical measures of agreement between predicted and observed were evaluated by analysing the results provided by the different alternative prediction methods to the same dataset. We also compared results obtained while increasing the number of predictors in MSN. Taking into account that many of the modelling alternatives have been already ruled out as unreliable by additional statistical criteria (Valbuena et al., 2017a), our objective was to observe whether any measures of the agreement would reveal similar conclusions or otherwise conceal the unreliability of predictions. Spearman’s rank correlation coefficient ($\rho$) was employed to assess redundancy among these measures since it tests whether two methods would rank alternatives in a similar manner.

Results

All the statistical measures of agreement between observed and predicted considered rendered most modelling alternatives reliable. Table 1 compares the results from different prediction method alternatives. The values of all measures have been expressed in per-cent units, with the intention to use them for interpreting the proportion of variance in the observed values that is explained by the model predictions, as $R^2$ is interpreted. It is worth noting the high values obtained by the unrestricted version of the step-wise selection, which unrealistically resulted in an over-parameterised model with $p=23$ predictors in spite of the use of AIC. Figure 1 shows bar plots comparing the results for parametric power models, according to whether they included or lacked the overfitting restrictions to variable selection. The versions which restricted the overfitting were more realistic. Willmott’s $d$ showed less contrast among methods than $R^2$ or $MIC$. It also showed lesser differences between values obtained using the whole dataset and cross-validated results.

Table 1. Comparison of diagnoses for different modelling methods and variable selection alternatives to obtain above-ground biomass ($AGB$, Mg·ha⁻¹) predictions

		Best-subset	Best-subset restricted overfitting	Step-wise	Step-wise restricted overfitting	MSN restricted overfitting
Number of predictors* (p)		8	2	23	2	5
Pre./Obs. agreement	$R_{fit}^2$ (%)	93.1	77.0	98.0	78.9	69.4
	$R_{cv}^2$ (%)	88.9	73.4	83.6	75.7	76.3
	$d^{fit}$ (%)	96.6	88.6	99.0	89.5	77.7
	$d^{cv}$ (%)	94.5	86.9	91.7	88.0	81.0
	$MI{C}^{fit}$ (%)	84.2	63.9	99.9	91.2	74.4
	$MI{C}^{cv}$ (%)	75.5	61.0	88.8	81.3	66.4

MSN: most similar neighbour. $R_{fit}^2$: residual coefficient of determination (Equation (1)(1) $\begin{array}{r}{R}^2=1-\left[{\sum }_{i=1}^n{\left(pr{e}_i-ob{s}_i\right)}^2/{\sum }_{i=1}^n{\left(ob{s}_i-\overline{obs}\right)}^2\right]\end{array}$(1) ). $R_{cv}^2$: cross-validated coefficient of determination (Equation (1)(1) $\begin{array}{r}{R}^2=1-\left[{\sum }_{i=1}^n{\left(pr{e}_i-ob{s}_i\right)}^2/{\sum }_{i=1}^n{\left(ob{s}_i-\overline{obs}\right)}^2\right]\end{array}$(1) ). $d^{fit}$: residual Willmott’s index of agreement (Equation (2)(2) $\begin{array}{rl}& d=1-\left[{\sum }_{i=1}^n{\left(pr{e}_i-ob{s}_i\right)}^2\right./\\ & \left.{\sum }_{i=1}^n{\left(\left|pr{e}_i-\overline{obs}\right|+\left|ob{s}_i-\overline{obs}\right|\right)}^2\right]\end{array}$(2) ). $d^{cv}$: cross-validated Willmott’s index of agreement (Equation (2)(2) $\begin{array}{rl}& d=1-\left[{\sum }_{i=1}^n{\left(pr{e}_i-ob{s}_i\right)}^2\right./\\ & \left.{\sum }_{i=1}^n{\left(\left|pr{e}_i-\overline{obs}\right|+\left|ob{s}_i-\overline{obs}\right|\right)}^2\right]\end{array}$(2) ). $MI{C}^{fit}$: residual maximal information coefficient (Equation (6)(6) $MIC=\underset{P,O_{total}<B}{max}\left\{MI\left(P,O\right)/log\left[min\left(P,O\right)\right]\right\}$(6) ). $d^{cv}$: cross-validated maximal information coefficient (Equation (6)(6) $MIC=\underset{P,O_{total}<B}{max}\left\{MI\left(P,O\right)/log\left[min\left(P,O\right)\right]\right\}$(6) ). Coefficients have been multiplied by 100 to yield percentage units.

*The actual predictors selected with each method are detailed in Appendix A.

Figure 1. Comparison of results for alternative statistics of agreement – Willmott’s (1981) index of agreement ($d$), maximal information coefficient ($MIC$; Reshef et al., 2011) and coefficient of determination (R²) – for power models with various variable selection methods, including versus lacking the overfitting restrictions (Valbuena et al., 2017a)

The effects of increasing $p$ in MSN predictions are expressed in Figure 2. The cross-validated coefficients of agreement ranged $R_{cv}^2=63.7-87.1\mathrm{\%}$ for $p=5-28$, however dropping as low as $R_{cv}^2=33.1\mathrm{\%}$ for $p=30$ or $R_{cv}^2=38.7\mathrm{\%}$ for $p=2$. Results in Figure 2 suggest that for small $p$ there may be little divergence whether coefficients of determination are obtained from the model fit $R_{fit}^2$ or cross-validated values $R_{cv}^2$. Willmott’s index of agreement $d$ showed similar patterns as $R_{cv}^2$ (Figure 2), and they were highly correlated $\rho \left(d^{cv},R_{cv}^2\right)=0.99$. This was not the case for $MI{C}^{cv}$ $(\rho \left(MI{C}^{cv},R_{cv}^2\right)=0.65)$, which was, therefore, less redundant than $d$. The difference between non-cross-validated and cross-validated results was, in general, more pronounced for $R^2$ or $MIC$ than for $d$, for the parametric models (Figure 1) as well as for the various MSN predictions (Figure 2).

Figure 2. Evolution of measures of agreement between observed and predicted by most similar neighbour (MSN) method for increasing the number of predictors ($p$). Dashed lines show values obtained using the whole training dataset, whereas solid lines were yielded by cross-validated predictions

In light of the results obtained in Figure 2, we deduced that $d^{cv}$ was systematically higher than $R_{cv}^2$, ranging as much as $d=59.2-93.8\mathrm{\%}$ for all $p=1-30$. Hence, the proportion of explained variance may simply seem larger when using $d$ instead of $R^2$. For this reason, we further investigated the convenience of using the alternative formulations for Willmott’s index of agreement – $d_1$ and $d_r$ –, which are all compared in Figure 3. Results were again either just lower or higher, but still redundant to the information already provided by $R^2$, since $\rho \left(d_1^{cv},R_{cv}^2\right)=0.99$ and $\rho \left(d_r^{cv},R_{cv}^2\right)=0.98$. Moreover, the latest refined version of Willmott’s index – $d_r^{cv}$ –, showed less contrast among alternatives than $R_{cv}^2$, i.e. with smoother fluctuations along increasing $p$, which rendered it less useful to discriminate the reliability of predictions. Also, when fitting the whole training dataset with very high $p$ these modified versions – $d_1^{fit}$ and $d_r^{fit}$ – obtained very low values, which seems unfair in light of the results obtained by $R_{fit}^2$, $d^{fit}$ and $MI{C}^{fit}$ for $p>22$ in Table 1.

Figure 3. Comparison of the different versions of Willmott’s index of agreement between observed and predicted (Willmott et al., 1985, 2012; Willmott & Wicks, 1980) against the coefficient of determination $R^2$, when the most similar neighbour (MSN) method increased the number of predictors ($p$). Dashed lines show values obtained using the whole training dataset, whereas solid lines were yielded by cross-validated predictions

Discussion

In studies carrying out remote sensing $AGB$ predictions, the most popular statistic employed to evaluate the agreement between observed and predicted is undoubtedly $R^2$ (e.g. Bright et al., 2012; Chen & Zhu, 2013; d’Oliveira et al., 2012; Erdody & Moskal, 2010; Hudak et al., 2006; Latifi et al., 2015; McInerney et al., 2010; Næsset, 2002; Straub, Tian, Seitz, & Reinartz, 2013; Valbuena et al., 2014; Wing et al., 2012; Zhao, Popescu, & Nelson, 2009). It is, however, well known that in presence of bias, $R^2$ may be high despite a poor correspondence between observed and predicted (see, e.g., Valbuena et al., 2014; White, Coops, & Scott, 2000). For this reason, it has been argued that $R^2$ may not necessarily be related to the accuracy of predictions (Paruelo et al., 1998), and hence Willmott’s (1981) $d$ has been suggested as a more appropriate alternative to evaluate remote sensing predictions (Almeida et al., 2016; García et al., 2010; Yebra & Chuvieco, 2009), in line with discussions held in the broader realm of statistical modelling (Fox, 1981; Willmott, 1982). In our results, however, we observed an underperformance of $d$ as a statistic truly expressing the degree of agreement between observed and predicted. Even having been obtained after cross-validation, values of $d$ summarized in Table 1 and Figure 1 were quite large in all cases, regardless of results in other measures of model performance (see Valbuena et al., 2017a). On the other hand, $R_{cv}^2$ obtained sounder results which were more likely to discriminate reliable from unreliable alternatives, whereas Figure 2 shows that $d$ may simply yield values systematically larger than $R_{cv}^2$. The modifications to $d$ did not necessarily overcome this problem, since Willmott et al.’s (1985) $d_1$ was instead systematically lower than $R^2$, and the more recent Willmott et al.’s (2012) $d_r$ showed similar but smoother patterns of variation (Figure 3). As a conclusion, we suggest that $R^2$ may be a better alternative than $d$, in light of our results. Although Willmott (1981) correctly affirmed that $R^2$ may show large values if predictions are evenly distributed through the range but yet biased, we suggest that $R_{cv}^2$ may simply be used along with a hypothesis test ensuring the absence of bias, such as those proposed by Piñeiro et al. (2008).

Willmott’s (1981, 1982) put forward a number of arguments in favour of using $d$, and many of them point out that $d$ may be conceptually superior to $R^2$ for describing the agreement in the observed versus predicted relationship for model evaluation (García et al., 2010; Yebra & Chuvieco, 2009). One interesting remark arises when comparing the calculation of Willmott’s (1981) index $d$ of agreement with that for Pearson’s sample correlation coefficient $r$. It is noteworthy to mention that, although $R^2$ expresses the proportion of variance in the dependent variable explained by a model while $r$ the dependence between two variables, $R^2$ and $r$ are closely related to one another. In particular, $R^2$ is equivalent to the square of Pearson’s correlation ($r^2$) when using a simple linear regression with intercept fit. Although this equivalence is not valid for MSN, we explored the possibility that $R^2$ and $r^2$ may be fairly close by analogy. Figure 4 shows that was indeed the case, as it was also for the square of Willmott’s index of agreement ($d^2$). Thus, we argue that studies with a particular interest in using Willmott’s (1981, 1982) suggestions for evaluating model predictions (Almeida et al., 2016; García et al., 2010; Grzegozewski et al., 2016; Wachholz de Souza et al., 2015; Yebra & Chuvieco, 2009) should report $d^2$ instead, in order to make it more comparable to more widespread studies using $R^2$ (e.g. Næsset, 2002, 2004; Hudak et al., 2006; Zhao et al., 2009; Erdody & Moskal, 2010; McInerney et al., 2010; Bright et al., 2012; d’Oliveira et al., 2012; Wing et al., 2012; Chen & Zhu, 2013; Straub et al., 2013; Asner & Mascaro, 2014; Valbuena et al., 2014; Latifi et al., 2015; Zolkos et al., 2013).

Figure 4. Comparison of original versions of Willmott’s index of agreement $d$ and other measures of correlation (dashed lines), and their squared versions (solid lines) against the coefficient of determination $R^2$

The case of using maximal information coefficient (Reshef et al., 2011) could also be supported by analogy to the relationship between $R^2$ and $r^2$, since $MIC$ is usually regarded as a non-parametric version of correlation (Chen & Yang, 2016; Görgens et al., 2017; Thomas et al., 2017; Vallières et al., 2017). Results for $MIC$ in Table 1 and Figure 2 can be best interpreted by observing the scatterplots of observed versus cross-validated predictions published in Valbuena et al. (2017a: Figures 2–3). The remarkably higher values of $MIC$ with respect to $R^2$ for, e.g., $p=8$ or $p=10$ (Figure 2), can be explained by the fact that these scatterplots show clustered cases in the observed versus predicted metric space (see Valbuena et al., 2017a: Figure 2). In that sense, we deduct that $MIC$ simply favours local clustering of predictions, whether or not they correspond to the observations. We, therefore, discourage the use of $MIC$ as a substitute of $R^2$ for model evaluation. Since Simon and Tibshirani (2011) recommended the use of distance correlation $dCor$ (Székely & Rizzo, 2017) as an alternative to $MIC$, we also included it in Figure 4. The analogy of squaring $d$ does not work for $MIC$, which simply showed unrealistically low values when squared. The square of distance correlation ($dCo{r}^2$), however, also obtained values which were fairly similar to those for $R^2$ (Figure 4).

Another unexpected result obtained in Figure 2 was the fact that differences between $R_{fit}^2$ and $R_{cv}^2$ only became apparent for large $p$, even for such a small $n$ as it was employed in our study case. This result suggests that, in spite of the popularity of cross-validation in the assessment of forest $AGB$ using remote sensing (e.g. Franco-Lopez et al., 2001; García et al., 2010; Hudak et al., 2006; Hudak, Crookston, Evans, Hall, & Falkowski, 2008; Latifi et al., 2015; McInerney et al., 2010; McRoberts et al., 2002; Næsset, 2002; Packalén & Maltamo, 2008; Valbuena et al., 2013a; Wing et al., 2012), reports on results obtained by leave-one-out cross-validation could easily converge to those yielded by the model fit (i.e. using the same dataset for training and validation). This may put into question the need for cross-validation at all since it makes little difference for those models with more realistic number of predictors $\left(p\le 3-5\right)$ (Figure 2). Although cross-validation may add robustness to an analysis with small sample sizes and few highly influential cases, our results show little difference compared to using the entire training dataset (Valbuena et al., 2017a). Future research should be devoted to further investigating the reliability of cross-validation, accounting for trade-offs between model generality and statistical power (Cohen et al., 2003) faced when pondering between using cross-validation or separating part of the available field information for independent model validation. It is useful for practical forest inventory to mention that models restricted according to Valbuena et al. (2017a), despite using less predictor variables, obtained similar performances in light of all the accuracy assessment figures analysed.

Conclusions

These are the conclusions that can be drawn from the results obtained in this study and their subsequent discussion. First, we observed an underperformance of Willmott’s (1981) $d$ in comparison to $R^2$, since the latter was more sensitive to actual model performance. Second, the refinements later introduced to that same index (Willmott et al., 1985, 2012) did not solve the observed shortcomings either. Third, we wish to put forward a recommendation to use a squared version of Willmott’s index $d^2$ for those who prefer it above $R^2$, to allow better comparison with studies using $R^2$. Fourth, the maximal information coefficient $MIC$ is not at all suitable for comparing relationships between observed and predicted in model evaluation. And finally, for a small number of predictors, the difference between using cross-validation or not using it at all can be negligible, which shows how useful can be its comparison for evaluating overfitting. More research would be needed to explore the usefulness of different cross-validation alternatives in the specific topic of remote sensing predictions of forest $AGB$.

Acknowledgments

This work was partially supported by the Spanish Directorate General for Scientific and Technical Research under Grant CGL2013-46387-C2-2-R. We also thank the Valsain Forest Centre, of the National Park Body (Spain), and Prof. David A. Coomes (University of Cambridge) for their valuable help. Dr Valbuena’s work is supported by an EU Horizon 2020 Marie Sklodowska-Curie Action entitled “Classification of forest structural types with LIDAR remote sensing applied to study tree size-density scaling theories” (LORENZLIDAR-658180). Danilo Almeida acknowledges support from São Paulo Research Foundation (FAPESP) (grant 2018/21338-3).

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the Spanish Directorate General for Scientific and Technical Research [Grant CGL2013-46387-C2-2-R]; EU Horizon 2020 Marie Sklodowska-Curie Action [LORENZLIDAR-658180]; São Paulo Research Foundation (FAPESP) [grant 2018/21338-3].

Appendices

Appendix A. Details on the predictors selected by each of these methods. These models were the same detailed in Valbuena et al. (2017a)

Best-subset	Best-subset restricted overfitting	Step-wise	Step-wise restricted overfitting	MSN restricted overfitting
H.GC	H.mean	H.GC	H.P75	H.L2
H.L3	NDVI.P75	H.CV	H.L3	H.L3
H.Mean		NDVI.GC		H.AAD
H.P95		H.L2		H.Mean
H.IQR		H.L4		H.GC
H.P40		H.L3
NDVI.Var		H.AAD
NDVI.AAD		H.Mean
		H.SD
		H.Kurt
		H.P90
		H.P95
		H.P80
		H.P75
		NDVI.SD
		NDVI.AAD
		NDVI.L1
		NDVI.L3
		NDVI.P25
		NDVI.P05
		NDVI.P75
		NDVI.P70
		NDVI.Var

Descriptors of predictor variables are detailed in Manzanera et al. (2016) and Valbuena et al (2017b). Names are composed of a prefix and a suffice separated by a dot (.). Prefices refer to the attribute characteristic used: height from LIDAR sensor (H) or normalized vegetation index from multispectral sensor (NDVI). Suffices refer to the statistical descriptor of the distribution of that attribute characteristic: mean (Mean), variance (Var), standard deviation (SD), average absolute deviation (AAD), percentiles (P05, P25, P40, P70, P75, P80, P90, P95), L-moments (L1, L2, L3, L4), Gini coefficient (GC), inter-quartile range (IQR).