Earth Observation Applications to Agricultural and Environmental Statistics

Abstract

Pixels classifiers do not allow assessment of uncertainty due to classification errors. This is a major limitation to the direct use of crop type maps (CTM) as agricultural or environmental statistics. Here, we include the results of our recent research to overcome this limitation, using statistical models that integrate ground and CTM data to evaluate any source of uncertainty, including classification errors. We use linear models for field data that can be treated as values of a continuous variable, such as crop acreage observed at parcel level, and multinomial logit models for field data that can be treated as values of a categorical variable, such as the type of crop observed at the pixel level. Among the tasks required to produce agricultural and environmental statistics, we focus on three for which CTM are especially useful: improving the effectiveness of currently used methods, increasing the granularity of statistics by producing small area estimations, and optimizing the sample design. A prototype for integrating ground and CTM data to achieve each of them is included, together with the results of its application to data provided by two National Statistical Offices.

1. Introduction

“End hunger, achieve food security, and promote sustainable agriculture” is the second (the first is “end poverty”) Sustainable Goal in the United National 2030 Development Agenda. To implement policies and practices in pursuit of these goals, and to monitor and control such policies, detailed information about land uses, crops acreage, and crops yield is necessary.

International programs to improve agricultural and environmental statistics have been recently launched by the World Bank and the Food and Agriculture Organization (FAO) of the United Nations. Looking for synergies with these programs, and with the aim of facilitating the use of remote sensing (RS) information at the National Statistical Offices (NSOs) supporting agricultural and environmental statistics, the European Space Agency (ESA) launched the Sentinels for Agricultural Statistics (Sen4Stat) project.

In this chapter, along with well-established results in the RS literature (the statistical classifiers for crop type maps (CTM) generation), we include the main results of our recent research in the applications of CTM to agricultural and environmental statistics, carried out within the framework of Sen4Stat and partially published in the Statistical Journal of the IAOS [1]. To date, pixel classifiers do not allow the assessment of uncertainty due to classification errors. This is a major limitation to directly using CTM data as statistical estimates, and here, we focus on how to use statistical models to integrate CTM with ground data to produce efficient agricultural and environmental statistics, taking into account any source of uncertainty, including the classification errors.

In large areas, where the sample size for collecting ground data is large enough, we follow an approach based on designing a probabilistic scheme to select the sample and on statistical models to integrate the ground data with CTM data. If ground data can be treated as values of a continuous variable, such as crop acreage observed at the parcel level, we use linear models. If ground data can be treated as values of a categorical variable, such as the type of crop observed at the pixel level, we use multinomial logit models. This design-based approach is robust in the sense that the inference is based not on the model, but on the sampling design, and the estimates are design-consistent. In small areas, where the sample size is small, the design-based approach is not accurate enough for most applications, and the inference is based on the model.

We focus on three of the tasks necessary to produce agricultural and environmental statistics, for which CTM is especially useful: improving the efficiency of currently used methods, increasing the granularity of statistics, and optimizing sample design. A ground and CTM data integration prototype to conduct each of these tasks is included, along with examples of its application to data provided by National Statistical Offices.

2. Pixels classification error

Many RS applications involve the classification of the $i=1,2,\cdots ,N$ population pixels, using multispectral ($L$ bands) reflectance data, ${\mathbf{x}}_i=\underset{1\le l\le L}{\mathit{row}}\left(x_{\mathit{li}}\right)$ [2]. In agricultural and environmental statistics, the product of pixels’ classification most often used is a CTM, where a class is associated with each one of the type of crops or land uses, $k=1,2,\cdots ,K$.

However, pixels classifiers are subject to error, and as a result, CTM are subject to uncertainty. In this section, we will show that pixel classifiers do not allow the assessment of this uncertainty. This is a major limitation to use directly CTM as agricultural or environmental statistics, because in addition to estimates of the study population characteristics, measures of the uncertainty of these estimates are required.

2.1 Classifiers

Let $i\in k$ denote the event “pixel $i$ belongs to class $k$.” Except for the few pixels included in the sample, there is uncertainty about the occurrence or not of this event, and to model this uncertainty, as well as the instrument measurement errors, we consider the joint probability function, $f\left(i\in k,{\mathbf{x}}_i\right)$, of the two events “pixel $i$ belongs to class $k$ and its spectral measure is ${\mathbf{x}}_i$.” This function is used to define the decision rule: “Classify $i$ in the class $k$ for which the probability $f\left(i\in k,{\mathbf{x}}_i\right)$ is maximum: $f\left(i\in k,{\mathbf{x}}_i\right)=\underset{k=1,2,\cdots ,K}{max}f\left(i\in k,{\mathbf{x}}_i\right)$.” Three main approaches proposed in the literature to implement this rule are (i) maximum likelihood, (ii) maximum entropy, and (iii) multinomial regression.

The maximum likelihood (ML) approach [2] assumes that ${\mathbf{x}}_i$ is distributed, within a given class $k$, according to a multivariate normal, $\left.{\mathbf{x}}_i\right|i\in k\to N_L\left({\mathbf{m}}_k,{\mathbf{V}}_k\right)$, of mean ${\mathbf{m}}_k=E\left(\left.{\mathbf{x}}_i\right|i\in k\right)$, and covariance matrix ${\mathbf{V}}_k=\mathit{Var}\left(\left.{\mathbf{x}}_i\right|i\in k\right)$. Under this assumption, the decision rule $f\left(i\in k,{\mathbf{x}}_i\right)=\underset{k=1,2,\cdots ,K}{max}f\left(i\in k,{\mathbf{x}}_i\right)$ is equivalent to classify $i$ in the class $k$ if $g_k\left({\mathbf{x}}_i\right)=\underset{k=1,2,\cdots ,K}{max}{g}_k\left({\mathbf{x}}_i\right)$, where $\left\{g_k\left({\mathbf{x}}_i\right)=logf\left(i\in k\right)-\frac{1}{2}log\left|det{\mathbf{V}}_k\right|-\left[\frac{1}{2}\left({\mathbf{x}}_i-{\mathbf{\mu }}_k\right){\mathbf{V}}_k^{-1}{\left({\mathbf{x}}_i-{\mathbf{\mu }}_k\right)}^T\right];i=1,2,\cdots ,N;k=1,2,\cdots ,K\right\}$ (see Appendix 1). ${\mathbf{m}}_k$ and ${\mathbf{V}}_k$ are unknown but can be estimated using a simple sample of $n_k$ pixels and the estimators ${\widehat{\mathbf{m}}}_k=\underset{1\le l\le L}{\mathit{col}}\left({\widehat{m}}_{\mathit{lk}}\right)=\underset{1\le l\le L}{\mathit{col}}\left(\frac{1}{n_k}{\sum }_{i=1}^{n_k}{x}_{\mathit{lki}}\right)$ and ${\widehat{\mathbf{V}}}_k=\underset{1\le l\le L}{\mathit{col}}\underset{1\le l^{\prime }\le L}{\mathit{row}}\left(\frac{1}{n_k-1}{\sum }_{i=1}^{n_k}\left(x_{\mathit{lki}}-{\widehat{m}}_{\mathit{lk}}\right)\left(x_{l^{\prime }\mathit{ki}}-{\widehat{m}}_{l^{\prime }k}\right)\right)$.

Two disadvantages of this approach are that normality is a very restrictive assumption for reflectance data and that it does not allow a direct computation of class probability, ${\mu }_{\mathit{ik}}=P\left(\left.i\in k\right|{\mathbf{x}}_i\right)$ (given ${\mathbf{x}}_i$, ${\mu }_{\mathit{ik}}$ is the probability that crop $k$ covers pixel $i$).

The maximum entropy (ME) approach has the advantage that it does not require any assumptions about statistical distribution and directly provides estimates of the probability of each class. The entropy is an uncertainty measure: for the $i^{\mathit{th}}$ pixel, the quantity of uncertainty in ${\mathbf{\mu }}_i=\underset{1\le k\le K}{\mathit{col}}\left({\mu }_{\mathit{ik}}\right)$ is $H\left({\mathbf{\mu }}_i\right)=-{\sum }_{k=1}^K{\mu }_{\mathit{ik}}ln{\mu }_{\mathit{ik}}$. To estimate $\mathbf{\mu }=\underset{1\le i\le N}{\mathit{col}}\underset{1\le k\le K}{\mathit{col}}\left({\mu }_{\mathit{ik}}\right)$, we use crops data observed in a sample of $n$ pixels, $\mathbf{y}=\underset{1\le i\le n}{\mathit{col}}\underset{1\le k\le K}{\mathit{col}}\left(y_{\mathit{ik}}\right)$ where $y_{\mathit{ik}}=1$ if crop $k$ covers pixel $i$ and $y_{\mathit{ik}}=0$ otherwise, together with RS data $\mathbf{X}=\underset{1\le i\le n}{\mathit{col}}\underset{1\le l\le L}{\mathit{row}}\left(x_{\mathit{li}}\right)$. The maximum entropy principle is to choose the class probability values, ${\mu }_{\mathit{ik}}$, compatible with the observed data, $\mathbf{y}$, that makes $H\left(\mathbf{\mu }\right)=-{\mathbf{\mu }}^{T}ln\mathbf{\mu }=-{\sum }_{i=1}^n{\sum }_{k=1}^K{\mu }_{\mathit{ik}}ln{\mu }_{\mathit{ik}}$ maximum.

The condition of compatibility between $\mathbf{\mu }$ and $\mathbf{y}$ is specified by the model $\mathbf{y}=\mathbf{\mu }+\mathbf{\epsilon }$, where $\mathbf{\epsilon }=\underset{1\le i\le n}{\mathit{col}}\underset{1\le k\le K}{\mathit{col}}\left({\epsilon }_{\mathit{ik}}\right)$ and ${\epsilon }_{\mathit{ik}}$ is an unobserved random perturbation. The RS data are introduced by transforming this model into a set of consistency restrictions: $\left({\mathbf{I}}_K\otimes {\mathbf{X}}^T\right)\mathbf{y}=\left({\mathbf{I}}_K\otimes {\mathbf{X}}^T\right)\mathbf{\mu }+\left({\mathbf{I}}_K\otimes {\mathbf{X}}^T\right)\mathbf{\epsilon }$, to which we impose the condition of uncorrelation between spectral measures and random perturbations, $\left({\mathbf{I}}_K\otimes {\mathbf{X}}^T\right)\mathbf{\epsilon }=\mathbf{0}$, to reduce the consistency restrictions to $\left({\mathbf{I}}_K\otimes {\mathbf{X}}^T\right)\mathbf{y}=\left({\mathbf{I}}_K\otimes {\mathbf{X}}^T\right)\mathbf{\mu }$.

The number $\mathit{nK}$ of unknown class probabilities, $\mathbf{\mu }=\underset{1\le i\le n}{\mathit{col}}\underset{1\le k\le K}{\mathit{col}}\left({\mu }_{\mathit{ik}}\right)$, is higher than the number $\mathit{LK}$ of consistency restrictions, and as a result, there is not a unique solution for $\mathbf{\mu }$. Among the countless solutions, we choose the one making $H\left(\mathbf{\mu }\right)$ maximum, subject to the consistency restrictions and the adding-up normalization conditions, $\underset{n}{\mathit{diag}}\left({\mathbf{1}}_K^T\right)\mathbf{\mu }={\mathbf{1}}_n$, where ${\mathbf{1}}_K$ is a $\left(K\times 1\right)$ vector of ones, and ${\mathbf{1}}_n$ is a $\left(n\times 1\right)$ vector of ones. This solution is found by solving the Lagrange function $\mathrm{\Psi }=-{\mathbf{\mu }}^{T}ln\mathbf{\mu }+{\mathrm{\beta }}^T\left[\left({\mathbf{I}}_K\otimes {\mathbf{X}}^T\right)\mathbf{\mu }-\left({\mathbf{I}}_K\otimes {\mathbf{X}}^T\right)\mathbf{y}\right]+{\mathrm{\lambda }}^T\left({\mathbf{1}}_n-\underset{n}{\mathit{diag}}\left(1_K^T\right)\mathbf{\mu }\right)$, where $\mathbf{\beta }=\underset{1\le k\le K}{\mathit{col}}\left({\mathbf{\beta }}_k\right)=\underset{1\le k\le K}{\mathit{col}}\underset{1\le l\le L}{\mathit{col}}\left({\beta }_{\mathit{kl}}\right)$ and $\mathrm{\lambda }=\underset{1\le i\le n}{\mathit{col}}\left({\lambda }_i\right)$.

The multinomial (MNL) approach is similar to the maximum entropy. In fact, the entropy $H\left(\mathbf{\mu }\right)$ can be derived from the multinomial distribution (see Appendix 1). Both approaches directly provide coincident class probabilities estimates, but the former requires to specify a function relating these probabilities with the RS data, and we consider the logit function, ${\mu }_{\mathit{ik}}=\frac{e^{{\sum }_{l=1}^L{\beta }_{\mathit{kl}}{x}_{\mathit{li}}}}{1+{\sum }_{k=1}^{K-1}{e}^{{\sum }_{l=1}^L{\beta }_{\mathit{kl}}{x}_{\mathit{li}}}}$, where the coefficients $\mathbf{\beta }=\underset{1\le k\le K}{\mathit{col}}\left({\mathbf{\beta }}_k\right)=\underset{1\le k\le K}{\mathit{col}}\underset{1\le l\le L}{\mathit{col}}\left({\beta }_{\mathit{kl}}\right)$ coincide with the Lagrange parameters linked to the consistency restrictions in the ME approach. Consequently, the class probability estimates ${\widehat{\mu }}_{\mathit{ik}}$ based on ME coincide with those based on MNL [3]: for any $i=1,2,\cdots ,N$ it is ${\widehat{\mu }}_{\mathit{ik}}=\frac{e^{{\sum }_{l=1}^L{\widehat{\beta }}_{\mathit{kl}}{x}_{\mathit{li}}}}{1+{\sum }_{k=1}^{K-1}{e}^{{\sum }_{l=1}^L{\widehat{\beta }}_{\mathit{kl}}{x}_{\mathit{li}}}}$ for $k=1,2,\cdots ,K-1$ and ${\widehat{\mu }}_{\mathit{iK}}=\frac{1}{1+{\sum }_{k=1}^{K-1}{e}^{{\sum }_{l=1}^L{\widehat{\beta }}_{\mathit{kl}}{x}_{\mathit{li}}}}$, where $\widehat{\mathbf{\beta }}=\underset{1\le k\le K}{\mathit{col}}\underset{1\le l\le L}{\mathit{col}}\left({\widehat{\beta }}_{\mathit{kl}}\right)$ is the $\mathbf{\beta }$ estimate.

Although MNL has the disadvantage over ME that it requires specifying a functional form to link ${\mu }_{\mathit{ik}}$ and ${\mathbf{x}}_i$, it has the advantage over ML and ME that it allows to calculate the accuracy of the estimates by integrating ground and RS data (even in the case where the classifier used is ML or ME), taking into account the uncertainty due to classification errors.

2.2 Classification errors

The decision rule is subject to errors, and as a result, the classification is subject to uncertainty. Let $i\subset k$ denote the decision to “include pixel $i$ in class $k$.” This decision is subject to two types of error: the type one, called false negative and denoted by $i\not\subset k\left|i\in k\right.$, consists of excluding from a given class $k$ a pixel $i$ that is actually of that class, and the error of type two, called false positive and denoted by $i\subset k\left|i\notin k\right.$, consists of including a pixel $i$ in a given class that is not of that class.

The computation of the classification error probabilities $P\left(i\subset k\left|i\notin k\right.\right)$ and $P\left(i\not\subset k\left|i\in k\right.\right)$ is almost never possible, mainly due to the high dimension $L$ of RS data [2]. In effect, $P\left(i\subset k\left|i\notin k\right.\right)=P\left(i\subset k\right)\left[1-P\left(\left.i\in k\right|i\subset k\right)\right]/\left[1-P\left(i\in k\right)\right]$ and $i\subset k\iff {\widehat{\mathbf{\beta }}}_k-{\widehat{\mathbf{\beta }}}_{k^{\prime }}>0;k\ne k^{\prime }$ so that $P\left(i\subset k\right)=P\left({\widehat{\mathbf{\beta }}}_k-{\widehat{\mathbf{\beta }}}_{k^{\prime }}>0\right);k\ne k^{\prime }$. In the MNL approach, the ${\widehat{\mathbf{\beta }}}_k-{\widehat{\mathbf{\beta }}}_{k^{\prime }}$ distribution is asymptotically normal L-dimensional, and the computation of $P\left({\widehat{\mathbf{\beta }}}_k-{\widehat{\mathbf{\beta }}}_{k^{\prime }}>0\right);k\ne k^{\prime }$ is complex due to the high $L$ values.

Instead of $P\left(i\subset k\left|i\notin k\right.\right)$ and $P\left(i\not\subset k\left|i\in k\right.\right)$, in practice, several ratios are calculated using the sample data (or, more often, by dividing the sample data into two groups, one for training and the other for testing the classifier), but this does not allow us to evaluate the uncertainty due to classification errors. For instance, if we denote by $n_{i\in k}$ the number of sample pixels that actually are of class $k$, and by $n_{i\in k\cap i\subset k}$ those among $n_{i\in k}$ that are correctly classify in $k$, then $n_{i\in k\cap i\subset k}/n_{i\in k}$ is an estimate of the so-called producer’s accuracy for class $k$. In the same way, if we denote by $n_{i\subset k}$ the number of sample pixels included in the class $k$, then $n_{i\in k\cap i\subset k}/n_{i\subset k}$ is an estimate of the so-called user’s accuracy for class $k$. An estimate of the so-called overall accuracy is ${\sum }_{k=1}^K{n}_{i\in k\cap i\subset k}/{\sum }_{k=1}^K{n}_{i\in k}$.

Accounting for the uncertainty due to classification errors is key, and that is why the data resulting from classification cannot be used directly as agricultural or environmental statistics. To take into account any uncertainties, including those due to classification errors, the CTM are integrated with the ground data, using statistical models in the way we show in the following sections.

3. Design-based approach: large area estimation

The approach to producing official statistics is based on a probabilistic scheme designed to select the sample from which to collect the required ground data. This scheme allows us to assign a probability of inclusion in the sample to each population unit. To estimate the characteristics of the study population, we use design-consistent estimators whose distribution is tightly concentrated around the true characteristics value, if the sample size is large enough. The uncertainty of these estimates is assessed using the statistical distribution defined by the inclusion probabilities.

The sample is selected using an area frame or a list frame [4, 5, 6]. An example of an area frame is that used by Spain: the national topographic map divided into sampling units, each of which is a square of 700 m on each side (49 hectares). The inclusion probabilities is equal to $n/N$, where $n$ and $N$ are the number of sampling units in the sample and in the population, respectively [6]. France and the United States also use area frames but using points as sampling units, rather than polygons. France to estimate crop acreage [7], and the United States to collect environmental data and develop a national inventory of natural resources [8, 9].

An example of list frame is that used by Senegal: it is based on the 2013 population census [10]. The sampling unit is a household with agricultural activities. The sampling scheme consists of two stages. In the first stage, a sample of Enumeration Areas (EA) is selected with replacement and probabilities proportional to the number of households in the EA. In the second stage, from each selected EA, a sample of households is selected without replacement and with equal probability among households with agricultural activities. To integrate the data collected on the household sample with the CTM data, a georeferenced pixel is selected in each parcel within the two-stage sample.

To assess uncertainties, and particularly those due to classification errors, we use statistical models. The choice of a model depends on ground data, and we consider linear models, for ground data collected at the parcel level, which can be treated as values of a continuous variable, and multinomial logit models for ground data collected at the pixel level, which can be treated as values of a categorical variable.

3.1 Ground data collected at the parcel level

Let the sampling unit, $i$, be a parcel or a cluster of parcels (say a polygon, as in Spain), and let $i=1,2,\cdots ,N$ be the set of sampling units of the population. We consider a linear model $y_i={\mathbf{x}}_i\mathbf{\beta }+{\epsilon }_i$ that relates ground data observed on the study crop in the $i^{\mathit{th}}$ sampling unit, $y_i$, with the CTM data, ${\mathbf{x}}_i$, aggregated at the $i^{\mathit{th}}$ level. For a given crop, ${\mathbf{x}}_i$ is a row vector with only two components; the first is the constant 1, $x_{1i}=1$, and the second, $x_{2i}=x_i$, is the number of pixels classified in the CTM as belonging to the study crop in $i$: for convenience, we will denote ${\mathbf{x}}_i$ in the general form ${\mathbf{x}}_i=\underset{1\le l\le L}{\mathit{row}}\left(x_{\mathit{li}}\right)$ where $L=2$, $x_{1i}=1$ and $x_{2i}=x_i$. The parameters vector $\mathbf{\beta }$ is unknown and must be estimated.

In sampling units not included in the sample, only ${\mathbf{x}}_i$ is observed. Given ${\mathbf{x}}_i$, the expected value of the unknown $y_i$ is modeled by ${\mathit{Ey}}_i={\mathbf{x}}_i\mathbf{\beta }$, and the deviation of $y_i$ from ${\mathbf{x}}_i\mathbf{\beta }$ is modeled by ${\epsilon }_i$. This deviation is a source of uncertainty, mainly due to classification errors, and we will show how to assess it using the ${\epsilon }_i$ sampling variance.

The survey variable total in the population is $y_N={\sum }_{i=1}^N{y}_i$. As a result of including the component $x_{1i}=1$ in the model, $y_N$ can be written as $y_N={\mathbf{x}}_N{\mathrm{B}}_N$, where ${\mathbf{x}}_N={\sum }_{i=1}^N{\mathbf{x}}_i=\left[N{x}_N\right]$, and $x_N={\sum }_{i=1}^N{x}_i$ is the total number of pixels classified in the CTM as belonging to the study crop, and ${\mathbf{B}}_N={\left({\mathbf{X}}_N^T{\mathbf{X}}_N\right)}^{-1}{\mathbf{X}}_N^T{\mathbf{y}}_N$ designates the vector of population regression coefficients ${\mathbf{B}}_N=\underset{1\le l\le L}{\mathit{col}}\left(B_l\right)$ (with $L=2$), where ${\mathbf{X}}_N=\underset{1\le i\le N}{\mathit{col}}\underset{1\le l\le L}{\mathit{row}}\left(x_{\mathit{li}}\right)$ (with $L=2$) and ${\mathbf{y}}_N=\underset{1\le i\le N}{\mathit{col}}\left(y_i\right)$ [11]. Since ${\mathbf{x}}_i$ is known for every population unit $i=1,2,\cdots ,N$, then ${\mathbf{x}}_N$ and ${\mathbf{X}}_N$ are known. However, ${\mathbf{y}}_N$ is unknown and so is ${\mathbf{B}}_N$.

To estimate ${\mathbf{B}}_N$, we use the data collected in the sample $n$, selected from the population with inclusion probabilities $\left\{{\pi }_i;i=1,2,\cdots ,N\right\}$. The design-based estimator ${\widehat{\mathbf{B}}}_{\pi }={\left({\sum }_{i=1}^n\frac{{\mathbf{x}}_i^T{\mathbf{x}}_i}{{\pi }_i}\right)}^{-1}{\sum }_{i=1}^n\frac{{\mathbf{x}}_i^T{y}_i}{{\pi }_i}$ is design-consistent for ${\mathbf{B}}_N$, and as a result, the synthetic or projective estimator ${\widehat{y}}_N={\mathbf{x}}_N{\widehat{\mathbf{B}}}_{\pi }$ is design-consistent for $y_N$ [11].

The error of ${\widehat{y}}_N$ as an estimator of $y_N$ is ${\widehat{y}}_N-y_N={\mathbf{x}}_N\left({\widehat{\mathbf{B}}}_{\pi }-{\mathbf{B}}_N\right)$, and its sampling variance is $V\left({\widehat{y}}_N-y_N\right)={\mathbf{x}}_{N}V\left({\widehat{\mathbf{B}}}_{\pi }-{\mathbf{B}}_N\right){\mathbf{x}}_N^T$, where $V\left({\widehat{\mathbf{B}}}_{\pi }-{\mathbf{B}}_N\right)={\left({\sum }_{i=1}^N{\mathbf{x}}_i^T{\mathbf{x}}_i\right)}^{-1}V{\sum }_{i=1}^n\frac{{\mathbf{x}}_i^T{\mathbf{\epsilon }}_{\mathit{\text{iN}}}}{{\pi }_i}{\left({\sum }_{i=1}^N{\mathbf{x}}_i^T{\mathbf{x}}_i\right)}^{-1}$ and ${\epsilon }_{\mathit{\text{iN}}}=y_i-{\mathbf{x}}_i{\mathbf{B}}_N$. The sampling variance is a measure of uncertainty that depends both on the deviation ${\epsilon }_{\mathit{\text{iN}}}$ between the true value of $y_i$ and its expected value based on ${\mathbf{x}}_i$ and on the inclusion probabilities, ${\pi }_i$. A design-consistent estimator of the sampling variance is $\widehat{V}\left({\widehat{y}}_N-y_N\right)=V\left({\sum }_{i=1}^n\frac{{\widehat{\epsilon }}_{\mathit{\text{iN}}}}{{\pi }_i}\right)$, where $V\left(.\right)$ is the design-based variance and ${\widehat{\epsilon }}_{\mathit{\text{iN}}}=y_i-{\mathbf{x}}_i{\widehat{\mathbf{B}}}_{\pi }$.

The asymptotic distribution of ${\widehat{y}}_N$ is normal, ${\left[V\left({\widehat{y}}_N-y_N\right)\right]}^{-1/2}\left({\widehat{y}}_N-y_N\right)\to N\left(0,1\right)$, and can be used for constructing uncertainty measures in terms of probability for $y_N$, such as confidence intervals.

To estimate the total of the survey variable, $y_{N_R}={\sum }_{i=1}^{N_R}{y}_i$, in large areas that are part of the population, such as a region or a province $R$ with $N_R$ sampling units within a country, we use the sample of size $n$ selected from the population with inclusion probabilities $\left\{{\pi }_i;i=1,2,\cdots ,N\right\}$ and the estimator ${\widehat{y}}_{N_R}={\mathbf{x}}_{N_R}{\widehat{\mathbf{B}}}_{\pi }+\frac{N_R}{{\widehat{N}}_R}{\sum }_{i=1}^{n_R}\frac{y_i-{\mathbf{x}}_i{\widehat{\mathbf{B}}}_{\pi }}{{\pi }_i}$, where $n_R$ is number of sampling units included in the sample $n$ that fall in R, and ${\widehat{N}}_R={\sum }_{i=1}^n\frac{I_i}{{\pi }_i}$ is an estimator of $N_R$, where $I_i=1$ if $i$ is in R and $I_i=0$ otherwise. A measure of the ${\widehat{y}}_{N_R}$ uncertainty is its sampling variance $V\left({\widehat{y}}_{N_R}-y_{N_R}\right)=N_R^{2}V\frac{1}{{\widehat{N}}_R}\left({\sum }_{i=1}^{n_R}\frac{{\epsilon }_{\mathit{\text{iN}}}}{{\pi }_i}\right)$, and a design-based estimator of $V\left({\widehat{y}}_{N_R}-y_{N_R}\right)$ is $\widehat{V}\left({\widehat{y}}_{N_R}-y_{N_R}\right)=N_R^{2}V\frac{1}{{\widehat{N}}_R}\left(\sum _{i=1}^{n_R}\frac{{\widehat{\epsilon }}_{\mathit{\text{iN}}}}{{\pi }_i}\right)$.

3.1.1 Crop type map relative efficiency

The relative efficiency of CTM data is ${\mathit{RE}}_{\mathit{CTM}}=V_G{\left(V_{G+\mathit{CTM}}\right)}^{-1}$, where $V_G$ is the current sampling variance based on ground data alone, and $V_{G+\mathit{CTM}}$ is the sampling variance based on the integration of ground and CTM data. If ${\mathit{RE}}_{\mathit{CTM}}>1$, then the current sampling variance may be reduced by the amount $V_G-V_{G+\mathit{CTM}}=\left(1-{\mathit{RE}}_{\mathit{CTM}}^{-1}\right)V_G$, without increasing the sample size. As a result, the ground sample size can be reduced to $n_{G+\mathit{CTM}}={\mathit{RE}}_{\mathit{CTM}}^{-1}{n}_G$ sampling units, without loss of accuracy. $n_G$ is the current sample size for ground data collection.

For ${\pi }_i=n/N$, the sampling variance based on ground data alone is $V_G=V{\sum }_{i=1}^n\frac{y_i}{n/N}=N^2\left(1-n/N\right)1/n\left(n-1\right){\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2$ and the sampling variance using both ground and CTM data is $V_{G+\mathit{CTM}}=V{\sum }_{i=1}^n\frac{y_i-{\mathbf{x}}_i{\widehat{\mathbf{B}}}_{\pi }}{n/N}=N^2\left(1-n/N\right)1/n\left(n-2\right){\sum }_{i=1}^n{\left(y_i-{\mathbf{x}}_i{\widehat{\mathbf{B}}}_{\pi }\right)}^2$. As a result, ${\mathit{RE}}_{\mathit{CTM}}\simeq {\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2{\left({\sum }_{i=1}^n{\left(y_i-{\mathbf{x}}_i{\widehat{\mathbf{B}}}_{\pi }\right)}^2\right)}^{-1}$ and $n_{G+\mathit{CTM}}=\left({\sum }_{i=1}^n{\left(y_i-{\mathbf{x}}_i{\widehat{\mathbf{B}}}_{\pi }\right)}^2\right){\left({\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2\right)}^{-1}{n}_G$.

3.1.2 Example

Our data source is Sen4Stat. The goal of this project was focused on developing an open-source system for generating CTM from Sentinel images, using data provided by the NSOs (https://www.esa-sen4stat.org/). We use field data on crops collected at the parcel level by the NSO of Spain, to illustrate the integration of continuous ground data and CTM data using linear models.

3.1.2.1 Crop acreage estimates

We consider a large area of size 100 × 100 km in Castilla y León (Spain), with X coordinates (in meters): (300,000,400,000) and Y coordinates (4,600,000,4,700,000). The sampling unit is a square with a side of 700 m, and the sample size is 419 sampling units. The study crop is barley, and the ground data were observed at the parcel level. However, data were aggregated at the sampling unit level for calculations so that $y_i$ represents the ground data on the barley acreage in sampling unit $i$ and in ${\mathbf{x}}_i=\left[1x_i\right]$, $x_i$ is the number of pixels classified in CTM as belonging to barley in this same sampling unit $i$.

Model parameter estimates are ${\widehat{\mathbf{B}}}_{\pi }={\left[0.0460.093\right]}^T$, and barley acreage estimates based both on ground data alone and on the integration of ground and CTM data are found in Table 1. As a result of this integration, the accuracy of the estimates improves considerably: the sampling error (or coefficient of variation, which is the ratio between the square root of the sampling variance and the acreage estimate in %) and the wide of the confidence interval (within whose limits the true value of the barley acreage lies, with a high probability (95%) were reduced by half.

Data	Acreage (hectare)	Uncertainty			Relative efficiency of CTM data
		95% Confidence interval (hectare)		Sampling error (CV%)
		Limits	Wide
Ground	236165.4	Lw:15951.7	40427.3	4.37	—
		Up: 56379
Ground + CTM	228550.1	Lw:19699.8	17700.5	1.98	5.2
		Up:37400.3

Table 1.

Barley acreage estimates in a 100 × 100 km area of Castilla y León, 2018.

The CTM relative efficiency, ${\mathit{RE}}_{\mathit{CTM}}$, was high: using CTM data, the current ground sample size could be reduced to less than one-fifth without loss of accuracy.

3.1.2.2 Crop yield estimates

Let $y_{\mathit{ij}}$ be the yield of barley observed in parcel $j$ of sampling unit $i$, and let ${\mathbf{x}}_{\mathit{ij}}=\underset{1\le l\le L}{\mathit{row}}\left(x_{\mathit{ijl}}\right)$ be a $1\times L$ row vector with $L=2$, $x_{\mathit{ij}1}=1$ in the first position, and a normalized difference vegetation index (NDVI) in the second position, . The estimates based both on ground data alone and on the integration of ground and NDVI data are found in Table 2.

Data	Yield (kg/hectare)	Uncertainty			Relative efficiency of RS data
		95% Confidence interval (kg/hectare)		Sampling error (CV%)
		Limits	Wide
Ground	4213.617	Lw: 4093.97	239.293	1.45	—
		Up: 4333.263
Ground + RS	4155.688	Lw: 4033.231	244.915	1.50	0.95
		Up: 4278.146

Table 2.

Barley yield estimates in a 100 × 100 km area of Castilla y León, 2018. NDVI.

Although the vegetation index contribution to explain crop yield is statistically significant, its contribution to improve yield estimator accuracy is low: the sampling error is of the same order of magnitude as if ground data alone were used. The relative efficiency is nearly 1, so the current sample size could barely be reduced without loss of accuracy. Additional research is required to make production estimates more reliable through improved yield estimation from RS data.

3.1.2.3 Crop production estimates

Table 3 shows the estimate production of barley, calculated as the product of the crop acreage estimate and yield estimate.

Data	Production (tons) (1000 kg)	Uncertainty
		95% Confidence interval (tons: 1000 kg)		Sampling error (CV%)
		Limits	Wide
Ground	987,841	Lw: 903640	168,402	4.35
		Up: 1072042
Ground + RS	965,733	Lw: 915194	101,078	2.67
		Up: 1016272

Table 3.

Barley production estimates in a 100 × 100 km area of Castilla y León, 2018.

Using RS data, the production estimator accuracy increased considerably: the estimation error decreased by half, even if the RS data is not reliable for crop yield.

3.1.2.4 Crop acreage estimates at the provincial level

We estimate barley acreage in that part of the provinces within the 100 × 100 km area of Castilla y León (Spain). Results are in Table 4.

Province	Ground data alone		Ground & CTM data		Relative efficiency of CTM data
	Acreage (has.)	Sampling error (CV%)	Acreage (has.)	Sampling error (CV%)
León	6853.2	24.15	6834.5	16.20	2.3
Palencia	88602.0	7.36	90535.3	3.33	4.7
Valladolid	128209.5	5.57	119707.4	2.66	5.1
Zamora	12324.2	17.71	10948.2	8.16	6.4
Total area	235989.1	4.37	228028.5	1.98	5.2

Table 4.

Barley acreage estimates at provincial level.

The estimate accuracy at provincial level is less than that at the national level, but the CTM relative efficiency is similar. Where the sampling variance using ground data alone is high (León), its integration with CTM data reduces the sampling error below the standard limits (20%) considered acceptable in official statistics.

3.2 Ground data collected at the pixel level

We use pixel $i$ as sampling unit, in the same way that a point is used, that is, treating the ground data observed at the pixel level as values of a categorical variable. To each pixel $i=1,2,\cdots ,N$ of the study area, we associate a survey vector ${\mathbf{y}}_i=\underset{1\le k\le K}{\mathit{col}}\left(y_{\mathit{ik}}\right)$, where $y_{\mathit{ik}}=1$ if crop $k$ covers pixel $i$ and $y_{\mathit{ik}}=0$ otherwise. $K$ is the total number of crops so that the constraint ${\sum }_{k=1}^K{y}_{\mathit{ik}}=1$ holds.

The total of this survey vector in the population is ${\mathbf{y}}_N={\sum }_{i=1}^N{\mathbf{y}}_i=\underset{1\le k\le K}{\mathit{col}}\left({\sum }_{i=1}^N{y}_{\mathit{ik}}\right)=\underset{1\le k\le K}{\mathit{col}}\left(y_{\mathit{Nk}}\right)$, where $y_{\mathit{Nk}}$ is the total number of pixels covered by crop $k$. It is assumed that a pixel is covered by a single crop (or that a class of mixed pixels is included) so that the area covered by $k$ is $A_k={\mathit{ay}}_{\mathit{Nk}}$, where $a$ is the area of the terrain represented by a pixel.

To estimate ${\mathbf{y}}_N$, we follow a design-based approach using ground data collected in a sample of $n$ pixels, selected from the population with inclusion probabilities $\left\{{\pi }_i;i=1,2,\cdots ,N\right\}$. To integrate these sample data with CTM data, we use a multinomial logit model [12, 13, 14]. We assume that the survey vector ${\mathbf{y}}_i=\underset{1\le k\le K}{\mathit{col}}\left(y_{\mathit{ik}}\right)$ follows a multinomial distribution $\mathit{MN}\left(1,{\mathbf{\mu }}_i\right)$, where ${\mathbf{\mu }}_i=\underset{1\le k\le K}{\mathit{col}}\left({\mu }_{\mathit{ik}}\right)$ and ${\mu }_{\mathit{ik}}$ is the probability that crop $k=1,2,\cdots ,K$ covers pixel $i$ (probability of $y_{\mathit{ik}}=1$ and $y_{i{k}^{\prime }}=0;\forall k^{\prime }\ne k$), with the constraint ${\sum }_{k=1}^K{\mu }_{\mathit{ik}}=1$.

To link ${\mu }_{\mathit{ik}}$ with CTM data, we use a logit model ${\mu }_{\mathit{ik}}=\frac{exp\left({\mathbf{x}}_i{\mathbf{\beta }}_k\right)}{{\sum }_{k=1}^{K}exp\left({\mathbf{x}}_i{\mathbf{\beta }}_k\right)}$, where ${\mathbf{x}}_i=\underset{1\le l\le L}{\mathit{row}}\left(x_{\mathit{il}}\right)$ is an indicator vector of the class where pixel $i$ is located within the CTM ($x_{\mathit{il}}=1$ if pixel $i$ is in the class $l$ and $x_{\mathit{il}}=0$ otherwise), and ${\mathbf{\beta }}_k=\underset{1\le l\le L}{\mathit{col}}\left({\beta }_{\mathit{kl}}\right)$ is an unknown parameter vector. To estimate ${\mu }_{\mathit{ik}}$, it is sufficient to obtain estimates of ${\mathbf{\beta }}_k$ for $k=1,2,\cdots ,K-1$ [12]. The probability of crop $K$ is ${\mu }_K=1-\sum _{k=1}^{K-1}{\mu }_{\mathit{ik}}$ and can be estimated using the estimates ${\widehat{\mathbf{\beta }}}_{\mathit{\pi k}}$ for $k=1,2,\cdots ,K-1$.

The design-based parameter estimator ${\widehat{\mathbf{B}}}_{\pi }=\underset{1\le k\le K-1}{\mathit{col}}\left({\widehat{\mathbf{B}}}_{\mathit{\pi k}}\right)=\underset{1\le k\le K-1}{\mathit{col}}\underset{1\le l\le L}{\mathit{col}}\left({\widehat{B}}_{\mathit{\pi kl}}\right)$ is design-consistent and can be found iteratively using ${\widehat{\mathbf{B}}}_{\pi }^{\left(m+1\right)}={\widehat{\mathbf{B}}}_{\pi }^{\left(m\right)}+{\left.{\left[{\left.{\sum }_{i=1}^n\frac{\mathbf{H}\left(y_i,\mathbf{\beta }\right)}{{\pi }_i}\right|}_{\mathbf{B}={\widehat{\mathbf{B}}}_{\pi }^{\left(m\right)}}\right]}^{-1}{\sum }_{i=1}^n\frac{\mathbf{b}\left(y_i,\mathbf{\beta }\right)}{{\pi }_i}\right|}_{\mathbf{B}={\widehat{\mathbf{B}}}_{\pi }^{\left(m\right)}}$ [11], where $\mathbf{b}\left(y_i,\mathbf{\beta }\right)=\underset{1\le k\le K-1}{\mathit{col}}\underset{1\le l\le L}{\mathit{col}}\left[\left(y_{\mathit{ik}}-{\mu }_{\mathit{ik}}\right)x_{\mathit{il}}\right]$ and $\mathbf{H}\left(y_i,\mathbf{\beta }\right)=\left[\underset{1\le k\le K-1}{\mathit{col}}\underset{1\le k^{\prime }\le K-1}{\mathit{row}}\left({\delta }_{k{k}^{\prime }}{\mu }_{\mathit{ik}}-{\mu }_{\mathit{ik}}{\mu }_{i{k}^{\prime }}\right)\right]\otimes {\mathbf{x}}_i^T{\mathbf{x}}_i$, with ${\delta }_{k{k}^{\prime }}=1$ if $k=k^{\prime }$ and ${\delta }_{k{k}^{\prime }}=0$ otherwise.

A design-consistent estimator of ${\mathbf{y}}_N=\underset{1\le k\le K}{\mathit{col}}\left(y_{\mathit{Nk}}\right)$ is ${\widehat{\mathbf{y}}}_N={\sum }_{i=1}^N{\widehat{\mathbf{\mu }}}_i+\frac{N}{{\widehat{N}}_p}{\sum }_{i=1}^n\frac{{\mathbf{y}}_i-{\widehat{\mathbf{\mu }}}_i}{{\pi }_i}=\underset{1\le k\le K-1}{\mathit{col}}\left({\widehat{y}}_{\mathit{Nk}}\right)$, where ${\widehat{\mathbf{\mu }}}_i=\underset{1\le k\le K-1}{\mathit{col}}\left({\widehat{\mu }}_{\mathit{ik}}\right)$, ${\widehat{\mu }}_{\mathit{ik}}=\frac{exp\left({\mathbf{x}}_i{\widehat{\mathbf{B}}}_{\mathit{\pi k}}\right)}{1+{\sum }_{k=1}^{K-1}exp\left({\mathbf{x}}_i{\widehat{\mathbf{B}}}_{\mathit{\pi k}}\right)}$ for $k=1,2,\cdots ,K-1$, ${\widehat{\mu }}_{\mathit{iK}}=1-{\sum }_{k=1}^{K-1}{\widehat{\mu }}_{\mathit{ik}}=\frac{1}{1+{\sum }_{k=1}^{K-1}exp\left({\mathbf{x}}_i{\widehat{\mathbf{B}}}_{\mathit{\pi k}}\right)}$, ${\widehat{N}}_p={\sum }_{i=1}^n\frac{1}{{\pi }_i}$, and ${\widehat{y}}_{\mathit{Nk}}={\sum }_{i=1}^N{\widehat{\mu }}_{\mathit{ki}}+\frac{N}{{\widehat{N}}_p}{\sum }_{i=1}^n\frac{y_{\mathit{ki}}-{\widehat{\mu }}_{\mathit{ki}}}{{\pi }_i}$ is the estimator of the total number of pixels covered by crop $k$ for $k=1,2,\cdots ,K-1$; for crop $K$, it is ${\widehat{y}}_{\mathit{NK}}=N-{\sum }_{k=1}^{K-1}{\widehat{y}}_{\mathit{Nk}}$. The sampling covariance matrix of ${\widehat{\mathbf{y}}}_N$ is $\mathbf{V}{\widehat{\mathbf{y}}}_N=N^2\mathbf{V}\frac{1}{{\widehat{N}}_p}{\sum }_{i=1}^n\frac{{\mathbf{y}}_i-{\widehat{\mathbf{\mu }}}_i}{{\pi }_i}=\underset{1\le k\le K-1}{\mathit{col}}\underset{1\le k^{\prime }\le K-1}{\mathit{row}}\left(N^2\mathit{Cov}\left({\widehat{y}}_{\mathit{rk}},{\widehat{y}}_{r{k}^{\prime }}\right)\right)$, where ${\widehat{y}}_{\mathit{rk}}=\frac{1}{{\widehat{N}}_p}{\sum }_{i=1}^n\frac{y_{\mathit{ki}}-{\mu }_{\mathit{ki}}}{{\pi }_i}$ is a function of ${\widehat{N}}_{\mathit{pk}}={\sum }_{i=1}^n\frac{y_{\mathit{ki}}}{{\pi }_i}$ (because ${\widehat{N}}_p={\sum }_{k=1}^K{\widehat{N}}_{\mathit{pk}}$), the estimators of the total number of pixels covered by crop $k=1,2,\cdots ,K$ based on ground data alone, and of ${\widehat{r}}_{\mathit{kN}}={\sum }_{i=1}^n\frac{y_{\mathit{ki}}-{\widehat{\mu }}_{\mathit{ki}}}{{\pi }_i}$, the estimator of the total of the deviations between $y_{\mathit{ki}}$ and ${\mu }_{\mathit{ki}}$ in the population.

The sampling variance is an uncertainty measure that depends both on the deviation between the true, ${\mathbf{y}}_i$, value and its expected value estimate, ${\widehat{\mathbf{\mu }}}_i$, and on the sample design, ${\pi }_i$. We focus on the sampling variance $V{\widehat{y}}_{\mathit{Nk}}$ for $k=1,2,\cdots ,K-1$, which are the elements along the diagonal of $\mathbf{V}{\widehat{\mathbf{y}}}_N$. Let ${\widehat{\mathbf{G}}}_k=\underset{1\le j\le K+1}{\mathit{row}\left(g_j\right)}=\left[{\widehat{r}}_{\mathit{kN}}{\widehat{N}}_{p1}{\widehat{N}}_{p2}\cdots {\widehat{N}}_{\mathit{pK}}\right]$ be a row vector whose $K+1$ components are the estimators on which ${\widehat{y}}_{\mathit{rk}}$ depends. The ${\widehat{y}}_{\mathit{Nk}}$ sampling variance is $V{\widehat{y}}_{\mathit{Nk}}=\underset{1\le j\le K+1}{\mathit{row}}\left(\frac{\partial {\widehat{y}}_{\mathit{rk}}}{\partial g_j}\right)\mathbf{V}{\widehat{\mathbf{G}}}_k\underset{1\le j\le K+1}{\mathit{col}}\left(\frac{\partial {\widehat{y}}_{\mathit{rk}}}{\partial g_j}\right)$, where $\mathbf{V}{\widehat{\mathbf{G}}}_k=\underset{1\le j\le K+1}{\mathit{col}}\underset{1\le j^{\prime }\le K+1}{\mathit{row}}\left(\mathit{Cov}\left(g_j,g_{j^{\prime }}\right)\right)$ is the design-based covariance matrix of ${\widehat{\mathbf{G}}}_k$. The sampling covariance matrix of ${\widehat{\mathbf{y}}}_N$ is estimated by replacing ${\mathbf{\mu }}_i$ with ${\widehat{\mathbf{\mu }}}_i$ in$\mathbf{V}{\widehat{\mathbf{y}}}_N$.

3.2.1 Crop type map relative efficiency

The sampling variance of the $y_{\mathit{Nk}}$ estimator based on ground data alone is $V_G=N^{2}V\frac{1}{{\widehat{N}}_p}{\sum }_{i=1}^n\frac{y_{\mathit{ki}}}{{\pi }_i}$ and that based on its integration with CTM data is $V_{G+\mathit{CTM}}=N^{2}V\frac{1}{{\widehat{N}}_p}{\sum }_{i=1}^n\frac{y_{\mathit{ki}}-{\mu }_{\mathit{ki}}}{{\pi }_i}$. The relative efficiency of CTM data for estimating the total number of pixels covered by crop $k$ is ${\mathit{RE}}_{\mathit{\text{CTMk}}}=V\frac{1}{{\widehat{N}}_p}{\sum }_{i=1}^n\frac{y_{\mathit{ki}}}{{\pi }_i}\times {\left(V\frac{1}{{\widehat{N}}_p}{\sum }_{i=1}^n\frac{y_{\mathit{ki}}-{\mu }_{\mathit{ki}}}{{\pi }_i}\right)}^{-1}$.

3.2.2 Example: crop acreage estimates

We use field data collected at the pixel level by the NSO of Senegal, to illustrate the integration of categorical ground data and CTM data using multinomial logit models. The crop covering each pixel included in the sample selected in the Nioro Department was observed in 2021. In the CTM, pixels covering agricultural areas were classified into four crop types: (i) maize, (ii) millet, (iii) groundnut, and (iv) other crops, and the data are encode as ${\mathbf{x}}_i=\left[1000\right]$ for any pixel $i$ of the class maize, ${\mathbf{x}}_i=\left[0100\right]$ for any pixel $i$ of the class millet, ${\mathbf{x}}_i=\left[0010\right]$ for any pixel $i$ of the class groundnut, and ${\mathbf{x}}_i=\left[0001\right]$ for any pixel $i$ of the class other crops.

We focus on estimating the acreage of the two main crops: millet and groundnut. The remaining crops (mainly maize) are included in a third category called other, whose estimated probability is one minus the estimated probability of millet and groundnut. Model parameter estimates are in Table 5.

Crop	Crop type map
	Maize	Millet	Groundnut	Other crops
Millet	−0.010655	1.47958334	0.88931346	8.76504093
Groundnut	−0.659204	−0.59871700	2.89873948	1.21898514

Table 5.

Model parameter estimates (${\widehat{\mathrm{B}}}_{\mathit{\text{πmillet}}}$ and ${\widehat{\mathrm{B}}}_{\mathit{\text{πgroundnut}}}$).

Crop acreage estimates of millet and groundnut based on ground and CTM data are in Table 6.

Crop	Acreage (Hectare)	Uncertainty
			Sampling error (CV%)	Limits of 95% Confidence interval
				Lower	Upper	Wide
Millet	89,215		4.11	81978.88	96330.4	14351.52
Groundnut	78,815		3.71	73089.15	84550.98	11461.82

Table 6.

Crop acreage estimates, Nioro (Senegal) 2021.

The CTM data efficiency estimate is in Table 7.

Crop type	Standard errors of proportion estimators		Relative efficiency of CTM data
	Ground data alone	Ground and CTM data
Millet	3.37	1.90	3.13
Groundnut	3.34	1.52	4.80

Table 7.

Efficiency of CTM data for crop acreage estimation, Nioro (Senegal).

The integration of ground data observed at the pixel level and CTM data improves the accuracy of the estimators based on ground data alone: the standard error is cut by a half, without loss of design-based consistency. As a result, the size of the current sample used for collecting ground data could be reduced to less than a third, without loss of accuracy.

Table 8 shows crop acreage estimates in the districts Medina Sabakh, Paoskoto, and Wack Ngouna.

District	Millet		Groundnut
	Acreage (Has.)	Sampling error (CV%)	Acreage (has.)	Sampling error (CV%)
Medina Sabakh	20067.21	8.6	19765.36	7.3
Paoskoto	38316.02	5.3	35018.23	4.0
Wack Ngouna	30831.77	11.9	24030.71	10.7
Total Nioro	89215.00	4.11	78815.00	3.71

Table 8.

Crop acreage estimates at the district level, Nioro (Senegal).

The estimate accuracy at the district level is lower than at the department level. However, the sampling error is within the standard limit (CV < 20%) for labeling as official statistics.

4. Model-based approach: small area estimation

In small areas, such a municipality, the sample size is small, and as a result, the design-based approach is not accurate enough for most applications. In these areas, inference is based on the model, rather than the sampling design, and the estimates are not robust, but model-dependent.

4.1 Ground data collected at the parcel level

Let the sampling unit, $i$, be a parcel or a cluster of parcels (say a polygon, as in Spain), and let $y_{\mathit{di}}$ be the value of the survey variable in a sampling unit $i=1,2,\cdots ,N_d$ of the small area $d=1,2,\cdots ,D$. Let $\left\{\left(y_{\mathit{di}},{\mathbf{x}}_{\mathit{di}}\right),i=1,2,\cdots ,n_d\right\}$ be the ground data set, $y_{\mathit{di}}$, and the CTM data set, ${\mathbf{x}}_{\mathit{di}}=\left[1x_{d2i}\right]$, in the $n_d$ sampling units from the whole sample $n$ that fall in $d$. Here, $x_{d2i}$ is the number of pixels classified in CTM as belonging to the study crop within the sampling unit $i$ of the small area $d$.

To estimate $y_{N_d}={\sum }_{i=1}^{N_d}{y}_{\mathit{di}}$, the total of the survey variable in $d$, we use CTM data and the linear mixed model $y_{\mathit{di}}={\mathbf{x}}_{\mathit{di}}\mathbf{\beta }+u_d+{\epsilon }_{\mathit{di}}$, where $\left(u_d,{\epsilon }_{\mathit{di}}\right)$ are zero-mean independent random variables of variances $\left({\sigma }_u^2,{\sigma }_e^2\right)$.

The model-based estimator of $y_{N_d}$ is ${\widehat{y}}_{N_d}=N_d\left[\left(1-{\widehat{\gamma }}_d\right){\overline{\mathbf{x}}}_{N_d}\widehat{\mathbf{\beta }}+{\widehat{\gamma }}_d\left({\overline{y}}_{n_d}+\left({\overline{\mathbf{x}}}_{N_d}-{\overline{\mathbf{x}}}_{n_d}\right)\widehat{\mathbf{\beta }}\right)\right]$, where $\widehat{\mathbf{\beta }}={\left({\mathbf{X}}^T{\widehat{\mathbf{V}}}^{-1}\mathbf{X}\right)}^{-1}{\mathbf{X}}^T{\widehat{\mathbf{V}}}^{-1}\mathbf{y}$ is the generalized least-square estimator of $\mathbf{\beta }$, with $\mathbf{y}=\underset{1\le d\le D}{\mathit{col}}\underset{1\le i\le n_d}{\mathit{col}}\left(y_{\mathit{di}}\right)$, $\mathbf{X}=\underset{1\le d\le D}{\mathit{col}}\underset{1\le i\le n_d}{\mathit{col}}\left({\mathbf{x}}_{\mathit{di}}\right)$, ${\widehat{\mathbf{V}}}^{-1}=\underset{1\le d\le D}{\mathit{diag}}\left({\widehat{\mathbf{V}}}_d^{-1}\right)$ and ${\widehat{\mathbf{V}}}_d^{-1}=\frac{1}{{\widehat{\sigma }}_e^2}{\mathbf{I}}_{n_d}-\frac{{\widehat{\gamma }}_d}{n_d{\widehat{\sigma }}_e^2}{\mathbf{1}}_{n_d}{\mathbf{1}}_{n_d}^T$, where ${\widehat{\gamma }}_d=\frac{{\widehat{\sigma }}_u^2}{{\widehat{\sigma }}_u^2+{\widehat{\sigma }}_e^2/n_d}$. Here, ${\widehat{\sigma }}_e^2=\frac{{\mathbf{e}}^T\mathbf{e}}{n-D-1}$ and ${\widehat{\sigma }}_u^2=\frac{{\mathrm{u}}^T\mathrm{u}‐\left(n-2\right){\widehat{\sigma }}_e^2}{n_{\ast }}$, with $n_{\ast }=\sum _{d=1}^D{n}_d\left[1-n_d{\overline{\mathbf{x}}}_d{\left({\mathbf{X}}^T\mathbf{X}\right)}^{-1}{\overline{\mathbf{x}}}_d^T\right]$, are unbiased estimators of the variance components, where ${\mathbf{e}}^T\mathbf{e}$ is the sum of squared residuals in the model fitted by ordinary least squares, taking as fixed the small-area effect $u_d.{\mathbf{u}}^T\mathbf{u}$ is the sum of squared residuals in the model fitted by ordinary least squares, with $u_d=0$. ${\overline{\mathbf{x}}}_{N_d}$ and ${\overline{\mathbf{x}}}_{n_d}$ are the population and sample means of the vector ${\mathbf{x}}_{\mathit{di}}$, respectively.

An unbiased estimator of the total mean-square error estimator $\mathit{MSE}\left({\widehat{y}}_{N_d}\right)$ is [15]: $\widehat{M}\mathit{SE}\left({\widehat{y}}_{N_d}\right)=N_d^2{\left(1-\frac{n_d}{N_d}\right)}^2\left[h_{1d}\left({\widehat{\sigma }}_v^2,{\widehat{\sigma }}_e^2\right)+h_{2d}\left({\widehat{\sigma }}_v^2,{\widehat{\sigma }}_e^2\right)+2h_{3d}\left({\widehat{\sigma }}_v^2,{\widehat{\sigma }}_e^2\right)\right]$, where:

$$\begin{array}{l}{h}_{\mathit{1}d}\left({\widehat{\sigma }}_v^{\mathit{2}},{\widehat{\sigma }}_e^{\mathit{2}}\right)={\widehat{\gamma }}_d\left(\frac{{\widehat{\sigma }}_e^{\mathit{2}}}{n_d}\right)+{\left(1-\frac{n_d}{N_d}\right)}^2\frac{\left(N_d-n_d\right)}{N_d^2}\\ h_{\mathit{2}d}\left({\widehat{\sigma }}_v^{\mathit{2}},{\widehat{\sigma }}_e^{\mathit{2}}\right)={\widehat{\sigma }}_e^{\mathit{2}}\left({\overline{\mathbf{X}}}_{N_d-n_d}-{\widehat{\gamma }}_d{\overline{\mathbf{x}}}_d\right){\mathbf{A}}^{-1}{\left({\overline{\mathbf{X}}}_{N_d-n_d}-{\widehat{\gamma }}_d{\overline{\mathbf{x}}}_d\right)}^T\\ h_{3d}\left({\widehat{\sigma }}_v^{\mathit{2}},{\widehat{\sigma }}_e^{\mathit{2}}\right)=\frac{1}{n_d^2}\frac{1}{{\left({\widehat{\sigma }}_v^{\mathit{2}}+\frac{{\widehat{\sigma }}_e^{\mathit{2}}}{n_d}\right)}^3}\left[{\left({\widehat{\sigma }}_e^{\mathit{2}}\right)}^{2}V{\widehat{\sigma }}_v^{\mathit{2}}+{\left({\widehat{\sigma }}_v^{\mathit{2}}\right)}^{2}V{\widehat{\sigma }}_e^{\mathit{2}}-2{\widehat{\sigma }}_e^{\mathit{2}}{\widehat{\sigma }}_v^{\mathit{2}}\mathit{Cov}\left({\widehat{\sigma }}_e^{\mathit{2}},{\widehat{\sigma }}_v^{\mathit{2}}\right)\right]\end{array}$$

(where)

$$\begin{array}{l}\mathbf{A}=\sum _{d=1}^D\left(\sum _{i=1}^{n_d}{\mathbf{x}}_{\mathit{di}}^T{\mathbf{x}}_{\mathit{di}}-{\widehat{\gamma }}_d{n}_d{\overline{\mathbf{x}}}_{n_d}^T{\overline{\mathbf{x}}}_{n_d}\right)\\ V\left({\widehat{\sigma }}_v^{\mathit{2}}\right)=\frac{2}{n_{\ast }^2}\left[\frac{1}{n-D-1}\left(D-1\right)\left(n-2\right){\left({\widehat{\sigma }}_e^{\mathit{2}}\right)}^2+2n_{\ast }{\widehat{\sigma }}_e^{\mathit{2}}{\widehat{\sigma }}_v^{\mathit{2}}+n_{\ast \ast }{\left({\widehat{\sigma }}_v^{\mathit{2}}\right)}^2\right]\\ V\left({\widehat{\sigma }}_e^{\mathit{2}}\right)=\frac{2{\left({\widehat{\sigma }}_e^{\mathit{2}}\right)}^2}{n-D-1}\\ \mathit{Cov}\left({\widehat{\sigma }}_e^{\mathit{2}},{\widehat{\sigma }}_v^{\mathit{2}}\right)=-\frac{1}{n_{\ast }}\left(D-1\right)V{\widehat{\sigma }}_e^{\mathit{2}}\end{array}$$

Here, $n_{\ast \ast }={\sum }_{d=1}^D{n}_d^2\left(1-n_d{\overline{\mathbf{x}}}_{n_d}{\mathbf{A}}_1^{-1}{\overline{\mathbf{x}}}_{n_d}^T\right)+\mathit{tr}{\left({\mathbf{A}}_1^{-1}{\sum }_{d=1}^D{n}_d^2{\overline{\mathbf{x}}}_{n_d}^T{\overline{\mathbf{x}}}_{n_d}\right)}^2$. Because ${\mathbf{A}}_1={\sum }_{d=1}^D{\sum }_{i=1}^{n_d}{\mathbf{x}}_{\mathit{di}}^T{\mathbf{x}}_{\mathit{di}}$, $n_{\ast \ast }$ may be simplified to $n_{\ast \ast }={\sum }_{d=1}^D{n}_d^2\left[1-{\overline{\mathbf{x}}}_{n_d}{\left({\mathbf{X}}^T\mathbf{X}\right)}^{‐1}{\overline{\mathbf{x}}}_{n_d}^T\right]=n_{\ast }-n+{\sum }_{d=1}^D{n}_d^2$.

4.1.1 Example

Table 9 shows barley acreage estimates at the municipality level in that part of Zamora province within the 100 × 100 km area of Castilla y León (Spain).

Municipality	Sample size	Estimates
		Acreage (Hectares)	Sampling error (CV%)
Belver de los Montes	1	212.96	29.1
Castroverde	3	2914.22	8.0
Pinilla de Toro	4	963.30	10.0
Quintanilla del Monte	1	466.65	20.3
Toro	1	615.91	14.0
Vezdemarbán	3	1358.22	12.6
Villalpando	2	560.05	39.1
Villamayor de Campos	1	1056.23	11.1
Villanueva del Campo	1	784.03	13.2
Villar de Fallaves	1	844.16	11.0
Villardondiego	1	516.40	11.5
Villavendimio	1	656.07	10.4
Total Zamora	20	10948.2	8.16

Table 9.

Barley acreage estimates at municipality level.

Even in municipalities where the sample size is as small as a single sampling unit, CTM data contribute to reducing sampling error below the standard limit (20%) considered to label estimates as official statistics.

4.2 Ground data collected at the pixel level

We consider a population of $i=1,2,\cdots ,N$ pixels grouped in a set of $d=1,2,\cdots ,D$ small areas (municipalities) such that $N={\sum }_{d=1}^D{N}_d$, where $N_d$ is the number of pixels in the small area $d$. Let ${\mathbf{y}}_{\mathit{di}}=\underset{1\le k\le K}{\mathit{col}}\left(y_{\mathit{dik}}\right)$ be the survey vector where $y_{\mathit{dik}}=1$ when crop $k$ covers pixel $i$ and $y_{\mathit{dik}}=0$ otherwise.

To estimate ${\mathbf{y}}_{N_d}={\sum }_{i=1}^{N_d}{\mathbf{y}}_{\mathit{di}}=\underset{1\le k\le K}{\mathit{col}}\left(y_{N_{d}k}\right)$, the total of the survey vector over a small area $d$, where $y_{N_{d}k}={\sum }_{i=1}^{N_d}{y}_{\mathit{dik}}$ is the number of pixels of $d$ covered by crop $k$, we use CTM data, and we assume that ${\mathbf{y}}_{\mathit{di}}$ follows a multinomial distribution $\mathit{MN}\left(1,{\mathbf{\mu }}_{\mathit{di}}\right)$, where ${\mathbf{\mu }}_{\mathit{di}}=\underset{1\le k\le K}{\mathit{col}}\left({\mu }_{\mathit{dik}}\right)$ and ${\mu }_{\mathit{dik}}$ is the probability of $y_{\mathit{dik}}=1$ and $y_{\mathit{di}{k}^{\prime }}=0;\forall k^{\prime }\ne k$, with the constraint ${\sum }_{k=1}^K{\mu }_{\mathit{dik}}=1$. To link ${\mu }_{\mathit{dik}}$ with CTM data, we use a logit mixed model ${\mu }_{\mathit{dik}}=\frac{exp\left({\mathbf{x}}_{\mathit{di}}{\mathbf{\beta }}_k+v_{\mathit{dk}}\right)}{1+{\sum }_{k=1}^{K-1}exp\left({\mathbf{x}}_{\mathit{di}}{\mathbf{\beta }}_k+v_{\mathit{dk}}\right)}$, where ${\mu }_{\mathit{dik}}$ is the probability that crop $k$ covers $i$ of $d$, ${\mathbf{x}}_{\mathit{di}}=\underset{1\le l\le L}{\mathit{row}}\left(x_{\mathit{dil}}\right)$ is an indicator vector whose components are $x_{\mathit{dil}}=1$ if pixel $i$ of $d$ is classified in the crop-type class $l=1,2,\cdots ,L$ and $x_{\mathit{dil}}=0$ otherwise, and ${\mathbf{x}}_{\mathit{di}}{\mathbf{\beta }}_k+v_{\mathit{dk}}={\eta }_{\mathit{dik}}$ is a linear mixed predictor consisting of a fixed ${\mathbf{x}}_{\mathit{di}}{\mathbf{\beta }}_k$ and a random component, $v_{\mathit{dk}}$, associated with the small areas and specifics for each crop $k$.

Here, ${\mathbf{\beta }}_k=\underset{1\le l\le L}{\mathit{col}}\left({\beta }_{\mathit{kl}}\right)$ is a vector of parameters measuring the effects of the CTM predictors ${\mathbf{x}}_{\mathit{di}}$ on ${\mathbf{\mu }}_{\mathit{di}}$, for example, for a pixel of the crop-type class $l=1,2,\cdots ,L$ we have ${\mu }_{\mathit{dik}}=\frac{exp\left({\beta }_{\mathit{kl}}+v_{\mathit{dk}}\right)}{1+{\sum }_{k=1}^{K-1}exp\left({\beta }_{\mathit{kl}}+v_{\mathit{dk}}\right)}$ for crops $k=1,2,\cdots ,K-1$ and ${\mu }_{\mathit{diK}}=\frac{1}{1+{\sum }_{k=1}^{K-1}exp\left({\beta }_{\mathit{kl}}+v_{\mathit{dk}}\right)}$ for crop $K$. The random effect, $v_{\mathit{dk}}$, accounts for the heterogeneity of the survey vector ${\mathbf{y}}_{\mathit{di}}=\underset{1\le k\le K}{\mathit{col}}\left(y_{\mathit{dik}}\right)$ among small areas, not explained by CTM data. It is assumed that $v_{\mathit{dk}}$ is a realization of a zero-mean random variable whose small-area variance is ${\mathit{Vv}}_{\mathit{dk}}={\sigma }_k^2$ for crop $k$ and whose cross-crop covariance is $\mathit{Cov}\left(v_{\mathit{dk}},v_{d{k}^{\prime }}\right)={\sigma }_{k{k}^{\prime }}^2$ so that the vector of small-area random effects $\mathbf{V}=\underset{1\le d\le D}{\mathit{col}}\underset{1\le k\le K-1}{\mathit{col}}\left(v_{\mathit{dk}}\right)$ has mean $E\mathbf{V}=\mathbf{0}$ and covariance matrix $\mathbf{\Sigma }=\mathbf{\theta }\otimes {\mathbf{I}}_D$, where $\mathbf{\theta }=\underset{1\le k\le K-1}{\mathit{row}}\underset{1\le k^{\prime }\le K-1}{\mathit{col}}\left({\sigma }_{k{k}^{\prime }}^2\right)$ and ${\sigma }_{k{k}^{\prime }}^2={\sigma }_k^2$ if $k=k^{\prime }$.

To estimate ${\mathbf{y}}_{N_d}$, we observe ${\mathbf{y}}_{\mathit{di}}$ only in the $n_d$ pixels that fall within the small area $d$, that is, a small part of the full sample of size $n={\sum }_{d=1}^D{n}_d$; the values of ${\mathbf{y}}_{\mathit{di}}$ in the remaining $N_d-n_d$ pixels are unknown, and we want to predict them using the model. We consider ${\mathbf{y}}_{N_d}$ decomposed into two ${\mathbf{y}}_{\mathit{di}}$ aggregates ${\mathbf{y}}_{N_d}={\sum }_{i_s=1}^{n_d}{\mathbf{y}}_{{\mathit{di}}_s}+{\sum }_{i_r=1}^{N_d-n_d}{\mathbf{y}}_{{\mathit{di}}_r}$: one known, ${\sum }_{i_s=1}^{n_d}{\mathbf{y}}_{{\mathit{di}}_s}$, and the other, ${\sum }_{i_r=1}^{N_d-n_d}{\mathbf{y}}_{{\mathit{di}}_r}$, unknown. We use ${\widehat{\mathbf{y}}}_{N_d}={\sum }_{i_s=1}^{n_d}{\mathbf{y}}_{{\mathit{di}}_s}+{\sum }_{i_r=1}^{N_d-n_d}{\widehat{\mathbf{\mu }}}_{{\mathit{di}}_r}$ as a predictor of ${\mathbf{y}}_{N_d}$, where ${\widehat{\mathbf{\mu }}}_{{\mathit{di}}_r}=\underset{1\le k\le K-1}{\mathit{col}}\left({\widehat{\mu }}_{{\mathit{di}}_{r}k}\right)$ with ${\widehat{\mu }}_{{\mathit{di}}_{r}k}=\frac{exp\left({\mathbf{x}}_{{\mathit{di}}_r}{\widehat{\mathbf{\beta }}}_k+{\widehat{v}}_{\mathit{dk}}\right)}{1+{\sum }_{k=1}^{K-1}exp\left({\mathbf{x}}_{{\mathit{di}}_r}{\widehat{\mathbf{\beta }}}_k+{\widehat{v}}_{\mathit{dk}}\right)}$ for crop $k=1,2,\cdots ,K-1$, and ${\widehat{\mu }}_{{\mathit{di}}_{r}K}=\frac{1}{1+{\sum }_{k=1}^{K-1}exp\left({\mathbf{x}}_{{\mathit{di}}_r}{\widehat{\mathbf{\beta }}}_k+{\widehat{v}}_{\mathit{dk}}\right)}$ for crop $K$. Here, ${\widehat{\mathbf{\beta }}}_k$ is the ${\mathbf{\beta }}_k$ estimate and ${\widehat{v}}_{\mathit{dk}}$ is the predictor of the small-area random effect $v_{\mathit{dk}}$ in $\mathbf{v}=\underset{1\le d\le D}{\mathit{col}}\underset{1\le k\le K-1}{\mathit{col}}\left(v_{\mathit{dk}}\right)$. To estimate both the fixed-effect $\mathbf{\beta }=\underset{1\le k\le K-1}{\mathit{col}}\left({\mathbf{\beta }}_k\right)$ and random-effect $\mathbf{\theta }=\underset{1\le k\le K-1}{\mathit{row}}\underset{1\le k^{\prime }\le K-1}{\mathit{col}}\left({\sigma }_{k{k}^{\prime }}^2\right)$ model parameters, we use maximum-likelihood estimators and the big sample of size $n={\sum }_{d=1}^D{n}_d$. As a result, both predictors ${\widehat{\mathbf{\mu }}}_{{\mathit{di}}_r}$ and ${\widehat{\mathbf{y}}}_{N_d}$ are asymptotically unbiased.

The mean squared error of the predictor ${\widehat{\mathbf{y}}}_{N_d}$ of ${\mathbf{y}}_{N_d}$ is $\mathit{MSE}\left({\widehat{\mathbf{y}}}_{N_d}\right)=E\left({\widehat{\mathbf{y}}}_{N_d}-{\mathbf{y}}_{N_d}\right){\left({\widehat{\mathbf{y}}}_{N_d}-{\mathbf{y}}_{N_d}\right)}^T=E{\left[{\sum }_{i_r=1}^{N_d-n_d}\left({\widehat{\mathbf{\mu }}}_{{\mathit{di}}_r}-{\mathbf{y}}_{{\mathit{di}}_r}\right)\right]}^2$, where ${\sum }_{i_r=1}^{N_d-n_d}\left({\widehat{\mathbf{\mu }}}_{{\mathit{di}}_r}-{\mathbf{y}}_{{\mathit{di}}_r}\right)={\sum }_{i_r=1}^{N_d-n_d}\left({\widehat{\mathbf{\mu }}}_{{\mathit{di}}_r}-{\mathbf{\mu }}_{{\mathit{di}}_r}\right)\left({\mathbf{y}}_{{\mathit{di}}_r}-{\mathbf{\mu }}_{{\mathit{di}}_r}\right)={\widehat{\mathbf{\tau }}}_{\mathit{dr}}-{\mathbf{\tau }}_{\mathit{dr}}-{\mathbf{\epsilon }}_{\mathit{dr}}$. Here, ${\mathbf{\tau }}_{\mathit{dr}}={\sum }_{i_r=1}^{N_d-n_d}E{\mathbf{y}}_{{\mathit{di}}_r}={\sum }_{i_r=1}^{N_d-n_d}{\widehat{\mathbf{\mu }}}_{{\mathit{di}}_r}$ is the expected value of the survey vector in the $N_d-n_d$ non-observed pixels, ${\widehat{\mathbf{\tau }}}_{\mathit{dr}}={\sum }_{i_r=1}^{N_d-n_d}{\widehat{\mathbf{\mu }}}_{{\mathit{di}}_r}$ is the predictor of ${\mathbf{\tau }}_{\mathit{dr}}$, and ${\mathbf{\epsilon }}_{\mathit{dr}}={\sum }_{i_r=1}^{N_d-n_d}{\mathbf{\epsilon }}_{{\mathit{di}}_r}$ with ${\mathbf{\epsilon }}_{{\mathit{di}}_r}=\left({\mathbf{y}}_{{\mathit{di}}_r}-{\mathbf{\mu }}_{{\mathit{di}}_r}\right)$ is the unobservable error vector.

Then, $\mathit{MSE}\left({\widehat{\mathbf{y}}}_{N_d}\right)=E{\left[{\widehat{\mathbf{\tau }}}_{\mathit{dtr}}-{\mathbf{\tau }}_{\mathit{dtr}}-{\mathbf{\epsilon }}_{\mathit{dtr}}\right]}^2=\mathit{MSE}\left({\widehat{\mathbf{\tau }}}_{\mathit{dr}}\right)+E{\mathbf{\epsilon }}_{\mathit{dr}}{\mathbf{\epsilon }}_{\mathit{dr}}^T-E\left({\widehat{\mathbf{\tau }}}_{\mathit{dr}}-{\mathbf{\tau }}_{\mathit{dr}}\right){\mathbf{\epsilon }}_{\mathit{dr}}^T-E{\mathbf{\epsilon }}_{\mathit{dr}}{\left({\widehat{\mathbf{\tau }}}_{\mathit{dr}}-{\mathbf{\tau }}_{\mathit{dr}}\right)}^T$. If we assumed that ${\mathbf{\epsilon }}_{\mathit{dr}}$ and ${\widehat{\mathbf{\tau }}}_{\mathit{dr}}$ are independent, then the covariance term is null, $E{\mathbf{\epsilon }}_{\mathit{dr}}{\left({\widehat{\mathbf{\tau }}}_{\mathit{dr}}-{\mathbf{\tau }}_{\mathit{dr}}\right)}^T=0$, and the mean square error of ${\widehat{\mathbf{y}}}_{N_d}$ reduces to $\mathit{MSE}\left({\widehat{\mathbf{y}}}_{N_d}\right)=\mathit{MSE}\left({\widehat{\mathbf{\tau }}}_{\mathit{dr}}\right)+E{\mathbf{\epsilon }}_{\mathit{dr}}{\mathbf{\epsilon }}_{\mathit{dr}}^T$, where $E{\mathbf{\epsilon }}_{\mathit{dr}}{\mathbf{\epsilon }}_{\mathit{dr}}^T\simeq E\left(\left.{\mathbf{\epsilon }}_{\mathit{dr}}{\mathbf{\epsilon }}_{\mathit{dr}}^T\right|{\mathbf{V}}_d\right)={\sum }_{i_r=1}^{N_d-n_d}{\mathbf{Vy}}_{{\mathit{di}}_r}$, where ${\mathbf{Vy}}_{{\mathit{di}}_r}$ is the ${\mathbf{y}}_{{\mathit{di}}_r}$ covariance matrix.

Appendix 2 for more details on parameter estimations, prediction of random effects, and mean squared error estimations.

4.2.1 Example

We use the field data described in Section 3.2.2. Here, we encode the CTM data as follow: ${\mathbf{x}}_i=\left[100\right]$ for a pixel $i$ in the crop-type class maize and other crops, ${\mathbf{x}}_i=\left[010\right]$ for crop-type class millet, and ${\mathbf{x}}_i=\left[001\right]$ for crop-type class groundnut. We assume that random effects have the same variance ${\sigma }^2$ and are independent among crops: $\mathit{Cov}\left(v_{\mathit{dk}},v_{d{k}^{\prime }}\right)=0$ for $k\ne k^{\prime }$ and $k=1,2,\cdots ,K-1$. The model parameter estimates are in Table 10. The values of the z-Wald statistics, for ${\mathbf{\beta }}_k$ parameters, and $\lambda$ likelihood ratio, for the variance of the random effects ${\sigma }^2$, are shown in parentheses.

Crop $\left(k\right)$	${\widehat{\mathrm{\beta }}}_k=\underset{1\le l\le L}{\mathit{col}}\left({\widehat{\mathrm{\beta }}}_{\mathit{kl}}\right)$ for crop-type class $\left(l\right)$			Variance of the small area effect $\left({\widehat{\sigma }}^2\right)$
	Maize and other	Millet	Groundnut
Millet	0.109 (z = 0.240)	1.429 (z = 6.236)	0.559 (z = 1.608)	0.028 ($\lambda$ = 24.636)
Groundnut	−0.662 (z = −1.153)	−0.587 (z = −1.711)	2.522 (z = 8.494)

Table 10.

Model parameter estimates.

According to the z-Wald statistics, the parameters associated with the crop-type classes millet and groundnut are significant for both millet and groundnut (at level < 0.1). The parameters associated with the crop-type class maize and other are not significant (at the 10% level) for both millet and groundnut. According to the likelihood ratio test, the random small-area effect is also significant (at level < 0.01).

Crop-acreage estimates for millet and groundnut in four communities in Senegal’s Nioro region are in Table 11.

Municipality	Population (pixels)	Sample (pixels)	Crop acreage (hectares) (Error: RRMSE%)
			Millet	Groundnut
Paos Koto	1,163,948	6	4500 (26.4)	1957 (56.5)
Taiba Niassene	1,163,096	14	4902 (34.3)	2000 (74.5)
Keur Madongo	302,088	2	1281 (24.9)	577 (45.1)
Wack Ngouna	2,766,043	12	9280 (25.0)	4167 (47.1)

Table 11.

Crop acreage estimation at municipality level, Nioro (Senegal).

The integration of ground and CTM data makes possible to evaluate the estimates accuracy using the relative root mean squared error, which is the ratio between root mean squared error and crop acreage. It is remarkable that in small areas where the sampling fraction is almost null, this integration reduces the estimation error to levels close to the limits that allow the estimates to be labeled as official: for millet, the estimation error is close to 20%, and for groundnut, the root mean squared error is similar but the acreage is about half the acreage of millet so that, as a result, the relative root mean squared error is about twice than for millet.

5. Optimizing the sample design

The objective is to minimize the sampling variance, that is, maximizing the estimates accuracy, without increasing the cost. The sampling variance depends on the spatial correlation structure of the survey variable, $y_i$. The correlation ${\rho }_y\left(\mathit{\text{dist}}\left(s_i,s_{i^{\prime }}\right)\right)$ between two observations $\left(y_i,y_{i^{\prime }}\right)$ at the points of coordinates $s_i$ and $s_{i^{\prime }}$ is a positive and decreasing function of the distance, $d=\mathit{\text{dist}}\left(s_i,s_{i^{\prime }}\right)$: it decreases when the distance between these points increases. Details on how to use a land-use map to estimate the spatial correlation structure of crops data, using variogram and correlogram models [16], and on how to use this structure for sample design are in [17]. Here, we propose to use a CTM to estimate the parameters of the correlogram model. We use the estimated correlogram to evaluate the (anticipated) sampling variance [18] as a function of the design variables (sampling unit size and sample size) [19, 20].

We consider a simple case: to estimate the total of the survey variable in the population, $y_N={\sum }_{i=1}^N{y}_i$, using a simple random sample (SRS) scheme (equal probability of inclusion, ${\pi }_i=n/N$) and the estimator ${\widehat{y}}_{N_{\mathit{SRS}}}=N/n{\sum }_{i=1}^n{y}_i$. The sampling variance is $V{\widehat{y}}_{N_{\mathit{SRS}}}=N^2\left(1-n/N\right)S^2/n$, where and ${\overline{y}}_N=y_N/N$ are the variance and the mean of the survey variable in the population, respectively.

Pixels are arranged in rows and a columns, and we consider the sampling unit as a square block of $n_0$ pixels and identify each sampling unit $i=1,2,\cdots ,N$ by the row and column of the pixel in the lower left corner of the square block. The distance $d=\mathit{\text{dist}}\left(s_i,s_{i^{\prime }}\right)$ can be expressed as $d=\sqrt{u^2+v^2}$, where $u$ is the number of sampling units between $s_i$ and $s_{i^{\prime }}$ in the row direction, and $v$ is the number of sampling units between $s_i$ and $s_{i^{\prime }}$ in the column direction. The two most often used models to assess the structure of ${\rho }_y\left(d\right)$ are the exponential, $\rho \left(\left.u,v\right|a,\tau \right)=\left(1-\tau \right)e^{-d/a}$, and the spherical $\rho \left(\left.u,v\right|a,\tau \right)=\left(1-\tau \right)e^{-d/a}\left[1-\frac{3}{2}\frac{d}{a}+\frac{d^3}{2a^3}\right]$ if $d\le a$ and $\rho \left(\left.u,v\right|a,\tau \right)=0$ if $d>a$.

The model parameters are the range rate $a$ and ratio $\tau ={\tau }_0/\left({\tau }_0+{\tau }_d\right)$. Here, ${\tau }_0$ is the nugget effect, that is, the variation at or near the origin (independent of distance); ${\tau }_d$ is the partial sill (a function of distance d between sampling points); and $\left({\tau }_0+{\tau }_d\right)$ is the sill, that is, the maximum variation far from the origin.

The anticipated variance is the expected value (based on the correlogram model) of the sampling variance $\mathit{AV}{\widehat{y}}_{N_{\mathit{SRS}}}=\mathit{EV}{\widehat{y}}_{N_{\mathit{SRS}}}=N^2\left(1-n/N\right){\mathit{ES}}^2/n$, where ${\mathit{ES}}^2={\sigma }^2\mathrm{\Psi }\left(\left.N,n_0\right|a,\tau \right)$, ${\sigma }^2$ is the variance of $y_i$, and $\mathrm{\Psi }\left(N,n_0\left|a,\tau \right.\right)=n_0\left({\mathit{Nn}}_0-1\right)/\left(N-1\right)\left[1-\mathrm{\Phi }\left({\mathit{Nn}}_0\left|a,\tau \right.\right)\right]-{\mathit{Nn}}_0\left(n_0-1\right)/\left(N-1\right)\left[1-\mathrm{\Phi }\left(n_0\left|a,\tau \right.\right)\right]$. Here, $\mathrm{\Phi }\left({\mathit{Nn}}_0\left|a,\tau \right.\right)$ is the average correlation between pairs of pixels over the $C_{{\mathit{Nn}}_0}^2$ pairs in the population, and $\mathrm{\Phi }\left(n_0\left|a,\tau \right.\right)$ is the average correlation between pairs of pixels over the $C_{n_0}^2$ pairs in a sampling unit.

The sample design is optimized to find the sampling unit size, $n_0$, and the sample size, $n$, that minimize the sampling variance, $\underset{\left\{\left(n,n_0\right)\right\}}{min}\mathit{AV}{\widehat{y}}_{N_{\mathit{SRS}}}=\underset{\left\{\left(n,n_0\right)\right\}}{min}{N}^2\left(1-n/N\right){\sigma }^2/n\mathrm{\Psi }\left(N,n_0\left|a,\tau \right.\right)$ subject to the cost constraint, $C_0+C_w{\mathit{nn}}_0+C_k\sqrt{\mathit{\text{An}}}\le C$. Here, $C$ is the total budget, $C_0$ is a fixed cost independent of the sample size, $C_w$ is the interviewing cost, $C_k$ is a cost per distance unit, and A is the surveyed area. The solution to this problem is the optimum sampling unit size $n_0$ and optimum sample size $n$. In addition to the budget, $C$, this optimum solution is conditioned to the correlogram model parameters $\left(a,\tau \right)$.

5.1 Example

For illustrations, we use the case study in [20]. The data source in this case study was a land-use map, but we establish a parallelism between the map and the satellite image to show how the proposed approach should be applied when a CTM is used, instead of a land-use map. The map is digitized, but the minimum level of disaggregation is a square area of side 1 km so that the minimum sampling unit size in the CTM would be a square block of side 100 pixels of 10 m and the survey variable, $y_i$, would be the number of pixels of the study crop type within sampling unit $i$.

The studied region was an area of 200 km by 200 km in Castilla y León, Spain, that is, 20,000 by 20,000 pixels of 10 m and $N=40000$ sampling units of 10,000 pixels. We computed the empirical variogram for nonirrigated herbaceous crops using the moment estimator [17], $2\widehat{\gamma }\left(d\right)=\frac{1}{\left|N\left(d\right)\right|}\sum _{N\left(d\right)}{\left(y_i-y_{i^{\prime }}\right)}^2$, where $N\left(d\right)=\left\{\left(s_i,s_{i^{\prime }}\right)\left|\mathit{\text{dist}}\left(\left(s_i,s_{i^{\prime }}\right)=d\right)\right.\right\}$. A spherical model was fitted to the empirical variogram and found as parameter estimates, ${\widehat{\tau }}_0=728.96,{\widehat{\tau }}_d=539.60$, and $\widehat{a}=46.17$. Using these parameter estimates, we compute the estimate $\mathrm{\Phi }\left(N\left|\widehat{a},\widehat{\tau }\right.\right)=0.032$. The parameter ${\sigma }^2$ is estimated by ${\widehat{\sigma }}^2=\frac{S^2}{\left[1-\mathrm{\Phi }\left(N\left|\widehat{a},\widehat{\tau }\right.\right)\right]}$, where $S^2=\frac{1}{1-N}{\sum }_{i=1}^N{\left(y_i-\overline{y}\right)}^2$ and $\overline{y}=\frac{1}{N}{\sum }_{i=1}^N{y}_i$.

In addition to the sampling unit of minimum size $\left(n_0=1\right)$, we considered a set of sampling unit sizes ranging from $n_0=2$ to $n_0=10$ times the sampling unit of minimum size $\left(n_0=1\right)$. To estimate $\mathrm{\Phi }\left(n_0\left|a,\tau \right.\right)$, it is assumed that the sampling unit of size $n_0$ is a square lattice of side $\sqrt{n_0}$ and that the distance between a pair of its components is $d=\sqrt{2n_0}$ so that $\mathrm{\Phi }\left(n_0\left|\widehat{a},\widehat{\tau }\right.\right)=\left(1-\widehat{\tau }\right)\left[1-\frac{3}{2}\frac{\sqrt{2n_0}}{\widehat{a}}+\frac{{\left(\sqrt{2n_0}\right)}^3}{2{\widehat{a}}^3}\right]$.

The optimization problem is

$$\begin{array}{l}\underset{\left\{n,n_0\right\}}{min}{N}^2\left(1-\frac{n}{N}\right)\frac{{\widehat{\sigma }}^2}{n}\frac{n_0\left({\mathit{Nn}}_0-1\right)}{\left(N-1\right)}\left[1-0.032\right]-\frac{{\mathit{Nn}}_0\left(n_0-1\right)}{\left(N-1\right)}\left[1-\left(1-\widehat{\tau }\right)\left(1-\frac{3}{2}\frac{\sqrt{2n_0}}{\widehat{a}}+\frac{{\left(\sqrt{2n_0}\right)}^3}{2{\widehat{a}}^3}\right)\right]\\ s.t.C_0+C_w{\mathit{nn}}_0+C_k\sqrt{\mathit{\text{An}}}\le C\end{array}.$$

To find the solution to this problem, we use the cost function coefficients in [21] ($C=120000\text{\$}$, $C_0=20500\text{\$}$, $C_w=50\text{\$}$, and $C_k=1\text{\$}/\mathit{Km}$) and the nloptr package [22]. The solution is in Table 12.

Variable	Optimal solution
	Lower bound for the sampling unit size
	1	2	3	4	5	10
Sample size: $n$	1661	891	611	465	375	192
Sampling unit size: $n_0$	1	2	3	4	5	10
Minimum Anticipated Variance	1.1 e+9	1.4 e+9	1.7e+9	2.0e+9	2.3e+9	3.8e+9

Table 12.

SRS: Optimum solution for several values of lower bound for the sampling unit size.

The optimum sampling unit size is the lower bound allowed. The minimum anticipated variance increases from 1.1 e+9 to 3.8e+9, when the sampling unit size increases from 1 to 10. In other words, using sampling units of size 10, the sampling variance is 1.86 times higher than using sampling units of size 1.

6. Concluding remarks

We have shown how to integrate CTM and ground data to produce agricultural and environmental statistics, along with measures of their accuracy that take into account any source of uncertainty, including classification errors. We have developed prototypes for this integration at the national, provincial, and municipal levels. For large areas (national and provincial level), these prototypes are based on the design of a probabilistic scheme to select the sample from which to collect the ground data, and we use a statistical model to integrate this sample with CTM data. The resulting statistics are design-consistent and robust; that is, they are based on the sampling design rather than on the model. The sample design is key, and we have shown how to optimize it using CTM. For small areas (municipal level), the prototypes are based solely on the model so that the resulting statistics are model-dependent.

To facilitate the transference of these prototypes, we have included examples of their application to data provided by two National Statistical Offices: Spain and Senegal. The former is an example of countries where ground data are collected at the parcel level, and we use linear models for their integration with CTM. The latter is an example of countries where ground data are observed at the pixel level, and we use multinomial logit models for their integration with CTM. We show how to estimate the model parameters and how to use CTM data to compute the required statistics, as well as their accuracy measures in these two countries, at both large and small area levels.

We have evaluated the efficiency of CTM data to improve the methods currently used to produce official agricultural statistics in these two countries, and we have found it to be high: using CTM data, the current sample size at the level of large areas could be reduced to less than one-fifth in the case of Spain and between one-third and one-fifth in the case of Senegal, without loss of neither robustness nor accuracy. At the small area level, where the sampling fraction is almost null, we have shown how to produce estimates accurate enough to be labeled as official statistics, that is, useful for most applications, without increasing costs.

Acknowledgments

The research results presented in this chapter were supported by the European Space Agency, in the Sen4Stat project framework [Contract No. 4000127181/19/I-NS]. The national statistical offices of Senegal and Spain provided the ground data used in the examples.

A. Appendix 1—pixels classification

We want to divide the $i=1,2,\cdots ,N$ pixels into a set of $k=1,2,\cdots ,K$ groups, each one corresponding to a discriminable class on the terrain identified by a crop type or land use. Let $f\left(i\in k,{\mathbf{x}}_i\right)$ be the joint probability function of the events “pixel $i$ is in the terrain class $k$ $\left(i\in k\right)$, and the spectral measure in $i$ is x_i”. The decision rule is to “classify $i$ in the class $k$ if and only if $f\left(i\in k,{\mathbf{x}}_i\right)=\underset{k=1,2,\cdots ,K}{max}f\left(i\in k,{\mathbf{x}}_i\right)$.” Maximum likelihood, maximum entropy, and multinomial regression are three approaches to implement this rule.

A.1 Maximum likelihood

This approach looks for minimizing the expected loss due to classification errors. Let ${\lambda }_i\left(k^{\prime }\left|k\right.\right)$ be a loss function measuring the cost incurred when pixel $i$ is classified in class $k^{\prime }\ne k$ but is actually of class $k$. The expected loss is $E{\lambda }_i\left(k^{\prime }\left|k\right.\right)={\sum }_{k^{\prime }=1}^K{\lambda }_i\left(k^{\prime }\left|k\right.\right)f\left(i\in k^{\prime }\left|{\mathbf{x}}_i\right.\right)$, where $f\left(i\in k^{\prime }\left|{\mathbf{x}}_i\right.\right)=f\left(i\in k^{\prime }\right)f\left({\mathbf{x}}_i\left|i\in k^{\prime }\right.\right)/f\left({\mathbf{x}}_i\right)$. Using this loss function, we consider the set of $K$ discriminant functions, $\left\{g_{k^{\prime }}\left({\mathbf{x}}_i\right)=-{\sum }_{k^{\prime }=1}^K{\lambda }_i\left(k^{\prime }\left|k\right.\right)f\left(i\in k^{\prime }\right)f\left({\mathbf{x}}_i\left|i\in k^{\prime }\right.\right)/f\left({\mathbf{x}}_i\right);i=1,2,\cdots ,N;k^{\prime }=1,2,\cdots ,K\right\}$ and we decide to classify the pixel $i$ in the class $k^{\prime }$ if $g_{k^{\prime }}\left({\mathbf{x}}_i\right)=\underset{k=1,2,\cdots ,K}{max}{g}_{k^{\prime }}\left({\mathbf{x}}_i\right)=\underset{k=1,2,\cdots ,K}{min}-g_{k^{\prime }}\left({\mathbf{x}}_i\right)$. Note that this classifier minimizes the expected loss, as desired. A simple loss function is ${\lambda }_i\left(k^{\prime }\left|k\right.\right)=0\mathit{\text{if}}k=k^{\prime }$ and ${\lambda }_i\left(k^{\prime }\left|k\right.\right)=1$ otherwise. Using this loss function, the discriminant function is $g_{k^{\prime }}\left({\mathbf{x}}_i\right)=-{\sum }_{k^{\prime }\left(\ne k\right)=1}^{K}f\left({\mathbf{x}}_i\left|i\in k^{\prime }\right.\right)f\left(i\in k^{\prime }\right)/f\left({\mathbf{x}}_i\right)$, and since $f\left({\mathbf{x}}_i\right)$ is a constant for a given ${\mathbf{x}}_i$, an equivalent discriminant function is $g_{k^{\prime }}^{\prime }\left({\mathbf{x}}_i\right)=g_{k^{\prime }}\left({\mathbf{x}}_i\right)f\left({\mathbf{x}}_i\right)=-{\sum }_{k^{\prime }\left(\ne k\right)=1}^{K}f\left({\mathbf{x}}_i\left|i\in k^{\prime }\right.\right)f\left(i\in k^{\prime }\right)$. Note that $f\left({\mathbf{x}}_i\right)={\sum }_{k^{\prime }=1}^{K}f\left(i\in k^{\prime }\right)f\left({\mathbf{x}}_i\left|i\in k^{\prime }\right.\right)={\sum }_{k^{\prime }\left(\ne k\right)=1}^{K}f\left({\mathbf{x}}_i\left|i\in k^{\prime }\right.\right)f\left(i\in k^{\prime }\right)+f\left({\mathbf{x}}_i\left|i\in k\right.\right)f\left(i\in k\right)$, and as a result, $g_{k^{\prime }}^{\prime }\left({\mathbf{x}}_i\right)=-f\left({\mathbf{x}}_i\right)+f\left({\mathbf{x}}_i\left|i\in k\right.\right)f\left(i\in k\right)$. Since $f\left({\mathbf{x}}_i\right)$ is a constant for a given ${\mathbf{x}}_i$, an equivalent discriminant function is $g_{k^{\prime }}^{″}\left({\mathbf{x}}_i\right)=g_{k^{\prime }}^{\prime }\left({\mathbf{x}}_i\right)+f\left({\mathbf{x}}_i\right)=f\left({\mathbf{x}}_i\left|i\in k\right.\right)f\left(i\in k\right)$.

Thus, we consider the set of $K$ discriminant functions, $\left\{g_{k^{\prime }}^{″}\left({\mathbf{x}}_i\right)=f\left({\mathbf{x}}_i\left|i\in k\right.\right)f\left(i\in k\right);i=1,2,\cdots ,N;k=1,2,\cdots ,K\right\}$ and we decide to classify the pixel $i$ in the class $k^{\prime }$ if $g_{k^{\prime }}^{″}\left({\mathbf{x}}_i\right)=\underset{k=1,2,\cdots ,K}{max}{g}_{k^{\prime }}^{″}\left({\mathbf{x}}_i\right)=\underset{k=1,2,\cdots ,K}{max}f\left({\mathbf{x}}_i\left|i\in k\right.\right)f\left(i\in k\right)$.

Note that $\underset{k=1,2,\cdots ,K}{max}f\left({\mathbf{x}}_i\left|i\in k\right.\right)f\left(i\in k\right)=\underset{k=1,2,\cdots ,K}{max}logf\left({\mathbf{x}}_i\left|i\in k\right.\right)f\left(i\in k\right)$, where $logf\left({\mathbf{x}}_i\left|i\in k\right.\right)f\left(i\in k\right)=logf\left(i\in k\right)+logf\left({\mathbf{x}}_i\left|i\in k\right.\right)$, where $logf\left({\mathbf{x}}_i\left|i\in k\right.\right)=-\frac{L}{2}log\left(2\pi \right)-\frac{1}{2}log\left|det{\mathbf{V}}_k\right|-\left[\frac{1}{2}\left({\mathbf{x}}_i-{\mathbf{\mu }}_k\right){\mathbf{V}}_k^{-1}{\left({\mathbf{x}}_i-{\mathbf{\mu }}_k\right)}^T\right]$. Thus, an equivalent set of $K$ discriminant functions is $\left\{g_k^{‴}\left({\mathbf{x}}_i\right)=logf\left(i\in k\right)-\frac{1}{2}log\left|det{\mathbf{V}}_k\right|-\left[\frac{1}{2}\left({\mathbf{x}}_i-{\mathbf{\mu }}_k\right){\mathbf{V}}_k^{-1}{\left({\mathbf{x}}_i-{\mathbf{\mu }}_k\right)}^T\right],;,,i=1,,,,2,,,,\cdots ,N,,,,k=1,,,,2,,,,\cdots ,K\right\}$ and we decide to classify the pixel $i$ in the class $k$ if $g_k^{‴}\left({\mathbf{x}}_i\right)=\underset{k=1,2,\cdots ,K}{max}{g}_k^{‴}\left({\mathbf{x}}_i\right)$.

A.2 Maximum entropy

1. Multinomial distribution. The entropy $H\left(\mathbf{\mu }\right)$ can be derived from the multinomial distribution. Assume that the survey vector, ${\mathbf{y}}_i=\underset{1\le k\le K}{\mathit{col}}\left(y_{\mathit{ik}}\right)$, follows a multinomial distribution, $\mathit{MN}\left(1;{\mathbf{\mu }}_i\right)$. Then, ${\sum }_{i=1}^n{\mathbf{y}}_i=\underset{1\le k\le K}{\mathit{col}}\left(y_k={\sum }_{i=1}^n{y}_{\mathit{ik}}\right)$ follows a multinomial distribution $\mathit{MN}\left(n={\sum }_{k=1}^K{y}_k;\mathbf{\mu }=\underset{1\le k\le K}{\mathit{col}}\left({\mu }_k\right)\right)$, with ${\sum }_{k=1}^K{\mu }_k=1$. The number of ${\sum }_{i=1}^n{\mathbf{y}}_i$ realizations favorable to the observed sample ${\mathbf{y}}_n=\underset{1\le k\le K}{\mathit{col}}\left(y_k\right)$ is $W=\frac{n!}{{\prod }_{k=1}^K{y}_k!}$. Thus, $lnW=lnn!-{\sum }_{k=1}^{K}ln{y}_k!$ and, using the Stirling approximation, $lnW=nlnn-n-{\sum }_{k=1}^K{y}_{k}ln{y}_k+n=nlnn-{\sum }_{k=1}^K{y}_{k}ln{y}_k$. For large samples, and according to the Bernoulli theorem, $\frac{y_k}{n}\to {\mu }_k$ when $n\to \infty$, and as a result, $lnW\approx nlnn-{\sum }_{k=1}^{K}n{\mu }_{k}lnn{\mu }_k=nlnn-n\left({\sum }_{k=1}^K{\mu }_k\left(lnn+ln{\mu }_k\right)\right)=-n{\sum }_{k=1}^K{\mu }_{k}ln{\mu }_k$. Finally, $n^{-1}lnW\approx -{\sum }_{k=1}^K{\mu }_{k}ln{\mu }_k=H\left(\mathbf{\mu }\right)$. The entropy and multinomial logit approaches provide coincident class probability estimates because they are maximizing the same function.

2. Information gain measures. In the simple case where the available information reduces to sample ground data, $\mathbf{y}=\underset{1\le i\le n}{\mathit{col}}\underset{1\le k\le K}{\mathit{col}}\left(y_{\mathit{ik}}\right)$, $\mathbf{\mu }$ is considered compatible with the observed data ${\mathbf{y}}_n={\sum }_{i=1}^n{\mathbf{y}}_i=\underset{1\le k\le K}{\mathit{col}}\left(y_k={\sum }_{i=1}^n{y}_{\mathit{ik}}\right)$, when ${\mathbf{y}}_n=\mathbf{\mu }+{\mathbf{\epsilon }}_n$, where ${\mathbf{\epsilon }}_n=\underset{1\le k\le K}{\mathit{col}}\left({\sum }_{i=1}^n{\epsilon }_{\mathit{ik}}\right)$ and ${\epsilon }_{\mathit{ik}}$ is a random noise. The set of $\mathbf{\mu }=\underset{1\le k\le K}{\mathit{col}}\left({\mu }_k\right)$ values compatible with ${\mathbf{y}}_n$ is countless, and according to the maximum entropy principle, we choose as $\mathbf{\mu }$ value the one making $H\left(\mathbf{\mu }\right)$ maximum: $\widehat{\mathbf{\mu }}=\underset{\mathbf{\mu }}{max}H\left(\mathbf{\mu }\right)=\underset{\mathbf{\mu }}{max}-{\sum }_{k=1}^K{\mu }_{k}ln{\mu }_k=\underset{\mathbf{\mu }}{min}{\sum }_{k=1}^K{\mu }_{k}ln{\mu }_k$, subject to ${\sum }_{k=1}^K{\mu }_k=1$. To find $\widehat{\mathbf{\mu }}$, we solve the Lagrange function $\mathrm{\Psi }=-{\sum }_{k=1}^K{\mu }_{k}ln{\mu }_k+\lambda \left(1-{\sum }_{k=1}^K{\mu }_k\right)$. The $\mathrm{\Psi }$ derivative with respect to ${\mu }_k$ is $\frac{\mathrm{\partial \Psi }}{\partial {\mu }_k}=\frac{\partial }{\partial {\mu }_k}\left[-{\sum }_{k=1}^K{\mu }_{k}ln{\mu }_k+\lambda \left(1-{\sum }_{k=1}^K{\mu }_k\right)\right]=-\left({\mu }_k\frac{1}{{\mu }_k}+ln{\mu }_k\right)-\lambda =-ln{\mu }_k-1-\lambda$, for $k=1,2,\cdots ,K$ and with respect to $\lambda$ is $\frac{\mathrm{\partial \Psi }}{\partial \lambda }=\left(1-{\sum }_{k=1}^K{\mu }_k\right)$. The solution of this system of $K+1$ equations is $-ln{\widehat{\mu }}_k-1-\widehat{\lambda }=0\iff {\widehat{\mu }}_k=e^{-1-\widehat{\lambda }};k=1,2,\cdots ,K$, with ${\sum }_{k=1}^K{\widehat{\mu }}_k={\mathit{Ke}}^{-1-\widehat{\lambda }}=1$, so that $ln{\mathit{Ke}}^{-1-\widehat{\lambda }}=ln1=0\iff lnK-1-\widehat{\lambda }=0\iff \widehat{\lambda }=lnK-1$. The entropy is $H\left(\widehat{\mathbf{\mu }}\right)=-{\sum }_{k=1}^K{e}^{-1-\widehat{\lambda }}\left(-1-\widehat{\lambda }\right)={\mathit{Ke}}^{-1-\widehat{\lambda }}\left(1+\widehat{\lambda }\right)=\left(1+\widehat{\lambda }\right)=lnK$, and it is achieved when ${\widehat{\mu }}_k=e^{-1-\widehat{\lambda }}=e^{-lnK};k=1,2,\cdots ,K$, that is, when, $ln{\widehat{\mu }}_k=-lnK=ln\frac{1}{K}\iff {\widehat{\mu }}_k=\frac{1}{K};k=1,2,\cdots ,K$ so that the assigned probabilities $\widehat{\mathbf{\mu }}=\underset{1\le k\le K}{\mathit{col}}\left({\widehat{\mu }}_k\right)$ correspond to a uniform distribution.

When additional information, $\mathbf{X}=\underset{1\le i\le n}{\mathit{col}}\underset{1\le l\le L}{\mathit{row}}\left(x_{\mathit{il}}\right)$, is introduced, then a reduction of uncertainty is achieved, from $lnK$ to $-{\sum }_{k=1}^K{\widehat{\mu }}_{k}ln{\widehat{\mu }}_k$. This reduction, $lnK+{\sum }_{k=1}^K{\widehat{\mu }}_{k}ln{\widehat{\mu }}_k$, is called “information gain,” and the ratio between this reduction and the maximum uncertainty, $I\left(\widehat{\mathbf{\mu }}\right)=\frac{lnK+{\sum }_{k=1}^K{\widehat{\mu }}_{k}ln{\widehat{\mu }}_k}{lnK}=1-S\left(\widehat{\mathbf{\mu }}\right)$, is called “information index”, where $S\left(\widehat{\mathbf{\mu }}\right)=-\frac{{\sum }_{k=1}^K{\widehat{\mu }}_{k}ln{\widehat{\mu }}_k}{lnK}$. A value $S\left(\widehat{\mathbf{\mu }}\right)=0$, that is, $I\left(\widehat{\mathbf{\mu }}\right)=1$, corresponds to no uncertainty (${\widehat{\mu }}_k=1$ for some $k$ and ${\widehat{\mu }}_{k^{\prime }}=0$ for every $k^{\prime }\ne k$), and a value $S\left(\widehat{\mathbf{\mu }}\right)=1$, that is, $I\left(\widehat{\mathbf{\mu }}\right)=0$, corresponds to perfect uncertainty (${\widehat{\mu }}_k=\frac{1}{K};k=1,2,\cdots ,K$). The measure $S\left(\widehat{\mathbf{\mu }}\right)$ can be used to evaluate the effect of each piece of information. For instance, the effect of reducing the sample size in a unit can be measured by comparing $S\left(\left.\widehat{\mathbf{\mu }}\right|n\right)$ and $S\left(\left.\widehat{\mathbf{\mu }}\right|n-1\right)$: if $S\left(\left.\widehat{\mathbf{\mu }}\right|n\right)<S\left(\left.\widehat{\mathbf{\mu }}\right|n-1\right)$, then the unit in question reduced the uncertainty. In the same way, the effect of excluding an auxiliary variable can be measured using $S\left(\left.\widehat{\mathbf{\mu }}\right|L\right)$ and $S\left(\left.\widehat{\mathbf{\mu }}\right|L-1\right)$: if $S\left(\left.\widehat{\mathbf{\mu }}\right|L\right)<S\left(\left.\widehat{\mathbf{\mu }}\right|L-1\right)$, then the variable in question reduces the uncertainty.

B. Appendix 2 - Small area estimation based on multinomial logit mixed models

B.1 Parameter estimation

To estimate ${\mu }_{\mathit{dik}}$, it is sufficient to acquire estimates of ${\mathbf{\beta }}_k$ for $k=1,2,\cdots ,K-1$. The probability that crop $K$ covers pixel $i$ of $d$ is ${\mu }_{\mathit{diK}}=1-{\sum }_{k=1}^{K-1}{\mu }_{\mathit{dik}}$ and can be estimated using the estimates of ${\mathbf{\beta }}_k$ for $k=1,2,\cdots ,K-1$. We specify the vector $\mathbf{\eta }=\underset{1\le d\le D}{\mathit{col}}\underset{1\le i\le N_d}{\mathit{col}}\underset{1\le k\le K-1}{\mathit{col}}\left({\eta }_{\mathit{dik}}\right)$ where ${\eta }_{\mathit{dik}}={\mathbf{x}}_{\mathit{di}}{\mathbf{\beta }}_k+v_{\mathit{dk}}$ as a linear mixed model $\mathbf{\eta }=\mathbf{X\beta }+\mathbf{Zv}$, where $\mathbf{X}=\underset{1\le d\le D}{\mathit{col}}\underset{1\le i\le N_d}{\mathit{col}}\underset{1\le k\le K-1}{\mathit{diag}}\underset{1\le l\le L}{\mathit{row}}\left(x_{\mathit{dil}}\right)$ are known CTM values, $\mathbf{\beta }=\underset{1\le k\le K-1}{\mathit{col}}\underset{1\le l\le L}{\mathit{col}}\left({\beta }_{\mathit{kl}}\right)$ is a vector of parameters, and $\mathbf{Z}=\underset{1\le d\le D}{\mathit{diag}}\underset{1\le i\le N_d}{\mathit{col}}\underset{1\le k\le K-1}{\mathit{row}}\left(z_{\mathit{dik}}\right)$ is an indicator matrix where $z_{\mathit{dik}}=1$ if $k$ covers $i$ of $d$ and $z_{\mathit{dik}}=0$ otherwise. It is assumed that the vector of small-area random effect $\mathbf{v}=\underset{1\le d\le D}{\mathit{col}}\underset{1\le k\le K-1}{\mathit{col}}\left(v_{\mathit{dk}}\right)$ has mean $E\mathbf{v}=\mathbf{0}$ and covariance matrix $\mathbf{\Sigma }=\mathbf{\theta }\otimes {\mathbf{I}}_D$, where $\mathbf{\theta }=\underset{1\le k\le K-1}{\mathit{row}}\underset{1\le k^{\prime }\le K-1}{\mathit{col}}\left({\sigma }_{k{k}^{\prime }}^2\right)$ and ${\sigma }_{k{k}^{\prime }}^2={\sigma }_k^2$ if $k=k^{\prime }$.

The parameters in $\mathbf{\beta }$ and $\mathbf{\theta }$ are often estimated following a penalized quasi-likelihood approach [22, 23, 24]. As discussed by McCulloch and Searle [13], these methods are not completely satisfactory in practice, and these authors recommended to linearize the nonlinear mixed model instead, thus obtaining a working linear mixed model fitted by maximum likelihood. Here, we follow this recommendation and consider the link function applied to the ground data $g_k\left({\mathbf{y}}_{\mathit{di}}\right)=log\left(\frac{y_{\mathit{dik}}}{1-{\sum }_{k=1}^{K-1}{y}_{\mathit{dik}}}\right)$ for $k=1,2,\cdots ,K-1$, and its Taylor expansion about $E\left(\left.{\mathbf{y}}_{\mathit{di}}\right|\mathbf{v}\right)={\mathbf{\mu }}_{\mathit{di}}$, $g_k\left({\mathbf{y}}_{\mathit{di}}\right)\simeq g_k\left({\mathbf{\mu }}_{\mathit{di}}\right)+{\left.\frac{\partial g_k\left({\mathbf{y}}_{\mathit{di}}\right)}{\partial y_{\mathit{di}1}}\right|}_{{\mathbf{y}}_{\mathit{di}}={\mathbf{\mu }}_{\mathit{di}}}\left(y_{\mathit{di}1}-{\mu }_{\mathit{di}1}\right)+\cdots +{\left.\frac{\partial g_k\left({\mathbf{y}}_{\mathit{di}}\right)}{\partial y_{\mathit{diK}-1}}\right|}_{{\mathbf{y}}_{\mathit{di}}={\mathbf{\mu }}_{\mathit{di}}}\left(y_{\mathit{diK}-1}-{\mu }_{\mathit{diK}-1}\right)$, to obtain the working linear model $\mathbf{\xi }=\mathbf{X\beta }+\mathbf{Zv}+\mathbf{e}$ [24]. Here, $\mathbf{\xi }=\underset{1\le d\le D}{\mathit{col}}\underset{1\le i\le N_d}{\mathit{col}}\left({\mathbf{\xi }}_{\mathit{di}}\right)$ and ${\mathbf{\xi }}_{\mathit{di}}=\underset{1\le k\le K-1}{\mathit{col}}\left({\xi }_{\mathit{dik}}\right)$ is a working vector, $\mathbf{e}=\underset{1\le d\le D}{\mathit{col}}\underset{1\le i\le N_d}{\mathit{col}}\left({\mathbf{e}}_{\mathit{di}}\right)$ with ${\mathbf{e}}_{\mathit{di}}={\mathbf{H}}_{\mathit{di}}^{-1}\left({\mathbf{y}}_{\mathit{di}}-{\mathbf{\mu }}_{\mathit{di}}\right)$, and ${\mathbf{H}}_{\mathit{di}}=\frac{\partial h\left({\mathbf{\eta }}_{\mathit{di}}\right)}{\partial {\mathbf{\eta }}_{\mathit{di}}}=\underset{1\le k\le K-1}{\mathit{col}}\underset{1\le k^{\prime }\le K-1}{\mathit{row}}\left({\delta }_{k{k}^{\prime }}{\mu }_{\mathit{dik}}-{\mu }_{\mathit{dik}}{\mu }_{\mathit{di}{k}^{\prime }}\right)$ with ${\delta }_{k{k}^{\prime }}=1$ if $k=k^{\prime }$ and ${\delta }_{k{k}^{\prime }}=0$ otherwise, ${\mathbf{y}}_{\mathit{di}}=\underset{1\le k\le K-1}{\mathit{col}}\left(y_{\mathit{dik}}\right)$, and ${\mathbf{\mu }}_{\mathit{di}}=\underset{1\le k\le K-1}{\mathit{col}}\left({\mu }_{\mathit{dik}}\right)$. Here, $h\left({\mathbf{\eta }}_{\mathit{di}}\right)$ is the inverse of the link function $g\left({\mathbf{\mu }}_{\mathit{di}}\right)={\mathbf{\eta }}_{\mathit{di}}$, that is, ${\mathbf{\mu }}_{\mathit{di}}=g^{-1}\left({\mathbf{\eta }}_{\mathit{di}}\right)=h\left({\mathbf{\eta }}_{\mathit{di}}\right)$.

The working vector covariance matrix is $\mathbf{V}=\mathbf{V\xi }={\mathbf{Z\Sigma Z}}^T+\mathbf{Ve}$, where $\mathbf{Ve}=\underset{1\le d\le D}{\mathit{diag}}\underset{1\le i\le N_d}{\mathit{diag}}\left({\mathbf{Ve}}_{\mathit{di}}\right)$ and ${\mathbf{Ve}}_{\mathit{di}}={\mathbf{H}}_{\mathit{di}}^{-1}\mathbf{V}\left({\mathbf{y}}_{\mathit{di}}-{\mathbf{\mu }}_{\mathit{di}}\right){\mathbf{H}}_{\mathit{di}}^{-1}$, where $\mathbf{V}\left({\mathbf{y}}_{\mathit{di}}-{\mathbf{\mu }}_{\mathit{di}}\right)$ is the covariance matrix of${\mathbf{y}}_{\mathit{di}}\to \mathit{MN}\left(1,{\mathbf{\mu }}_{\mathit{di}}\right)$.

We assume that the $\mathbf{\xi }$ asymptotic distribution is normal, $\mathbf{\xi }\to \mathit{\text{AN}}\left(\mathbf{X\beta },\mathbf{V}\right)$, and we see $\mathbf{\xi }=\mathbf{X\beta }+\mathbf{Zv}+\mathbf{e}$ as a linear mixed model. ${\mathbf{Ve}}_{\mathit{di}}$ is complex because the derivative of the link function ${\mathbf{H}}_{\mathit{di}}$ depends on $\mathbf{V}$ through $\mathbf{\mu }$, and ${\mathbf{\mu }}_{\mathit{di}}$ also depends on $\mathbf{V}$. We consider a simplification of this problem; it consists of setting $\mathbf{v}=\mathbf{0}$ in ${\mathbf{\mu }}_{\mathit{di}}$ that reduces ${\mathbf{H}}_{\mathit{di}}$ to ${\mathbf{H}}_{\mathit{di}}^{\ast }=\underset{1\le k\le K-1}{\mathit{col}}\underset{1\le k^{\prime }\le K-1}{\mathit{row}}\left({\delta }_{k{k}^{\prime }}{\mu }_{\mathit{dik}}^{\ast }-{\mu }_{\mathit{dik}}^{\ast }{\mu }_{\mathit{di}{k}^{\prime }}^{\ast }\right)$, where ${\mu }_{\mathit{dik}}^{\ast }=\frac{exp\left({\mathbf{x}}_{\mathit{di}}{\mathbf{\beta }}_k\right)}{1+{\sum }_{k=1}^{K-1}exp\left({\mathbf{x}}_{\mathit{di}}{\mathbf{\beta }}_k\right)}$ no longer depends on the random effects $\mathbf{v}$ [13]. As a result, ${\mathbf{Ve}}_{\mathit{di}}^{\ast }={\mathbf{H}}_{\mathit{di}}^{\ast -1}\mathbf{V}\left({\mathbf{y}}_{\mathit{di}}-{\mathbf{\mu }}_{\mathit{di}}\right){\mathbf{H}}_{\mathit{di}}^{\ast -1}$, where ${\mathbf{\mu }}_{\mathit{di}}$ does depend on $\mathbf{v}$.

To determine the $\mathbf{\beta }$ and $\mathbf{\theta }$ maximum likelihood estimates, we follow an iterative procedure:

1. Let ${\mathbf{\theta }}^{\left(0\right)}$ and ${\mathbf{\beta }}^{\left(0\right)}$ be starting values of $\mathbf{\theta }$ and $\mathbf{\beta }$, respectively

2. Let ${\mathbf{v}}^{\left(\mathbf{0}\right)}=0$ and compute ${\mathbf{\xi }}^{\left(0\right)}={\mathbf{X\beta }}^{\left(0\right)}+{\mathbf{H}}^{\ast -1\left(0\right)}\left(\mathbf{y}-{\mathbf{\mu }}^{\left(0\right)}\right)$

3. Compute ${\mathbf{V\xi }}^{\left(0\right)}=\underset{1\le d\le D}{\mathit{diag}}\left({\mathbf{Ve}}_{\mathit{di}}^{\ast \left(0\right)}\right)$, where ${\mathbf{Ve}}_{\mathit{di}}^{\ast \left(0\right)}={\mathbf{H}}_{\mathit{di}}^{\ast -1\left(0\right)}\mathbf{V}\left({\mathbf{y}}_{\mathit{di}}-{\mathbf{\mu }}_{\mathit{di}}^{\left(0\right)}\right){\mathbf{H}}_{\mathit{di}}^{\ast -1\left(0\right)}$ with ${\mathbf{\mu }}_{\mathit{di}}^{\left(0\right)}=\underset{1\le k\le K-1}{\mathit{col}}\left({\mu }_{\mathit{dik}}^{\left(0\right)}\right)$ and ${\mu }_{\mathit{dik}}^{\left(0\right)}=\frac{exp\left({\mathbf{x}}_{\mathit{di}}{\mathbf{\beta }}_k^{\left(0\right)}\right)}{1+{\sum }_{k=1}^{K-1}exp\left({\mathbf{x}}_{\mathit{di}}{\mathbf{\beta }}_k^{\left(0\right)}\right)}$

4. Update ${\mathbf{\beta }}^{\left(0\right)}$using ${\widehat{\mathbf{\beta }}}^{\left(1\right)}={\left({\mathbf{X}}^T{\mathbf{V\xi }}^{{\left(0\right)}^{-1}}\mathbf{X}\right)}^{-1}{\mathbf{X}}^T{\mathbf{V\xi }}^{{\left(0\right)}^{-1}}{\mathbf{\xi }}^{\left(0\right)}$

5. Update ${\mathbf{\theta }}^{\left(0\right)}$ using ${\widehat{\mathbf{\theta }}}^{\left(1\right)}={\mathbf{\theta }}^{\left(0\right)}-{\left[\frac{{\partial }^{2}l\left(\mathbf{\beta }\right)}{\partial \mathbf{\theta }\mathbf{\partial }{\mathbf{\theta }}^T}\right]}_{\mathbf{\beta }={\widehat{\mathbf{\beta }}}^{\left(1\right)}}^{-1}{\left.\frac{\partial l\left(\mathbf{\beta }\right)}{\partial \mathbf{\theta }}\right|}_{\mathbf{\beta }={\widehat{\mathbf{\beta }}}^{\left(1\right)}}$, where $l\left(\mathbf{\beta }\right)$ is the profile log likelihood of the $\mathbf{\xi }$ distribution.

6. Update ${\widehat{\mathbf{V}}}^{\left(0\right)}$ using ${\widehat{\mathbf{v}}}^{\left(1\right)}={\widehat{\mathbf{\Sigma }}}^{\left(1\right)}{\mathbf{Z}}^T{\mathbf{V\xi }}^{{\left(1\right)}^{-1}}\left({\mathbf{\xi }}^{\left(0\right)}-\mathbf{X}{\widehat{\mathbf{\beta }}}^{\left(1\right)}\right)$ and compute ${\mathbf{\xi }}^{\left(1\right)}=\mathbf{X}{\widehat{\mathbf{\beta }}}^{\left(1\right)}+\mathbf{Z}{\widehat{\mathbf{v}}}^{\left(1\right)}+{\widehat{\mathbf{H}}}^{\ast -1\left(1\right)}\left(\mathbf{y}-{\widehat{\mathbf{\mu }}}^{\left(1\right)}\right)$

7. Update ${\widehat{\mathbf{\beta }}}^{\left(m\right)}$ using ${\widehat{\mathbf{\beta }}}^{\left(m+1\right)}={\left({\mathbf{X}}^T{\mathbf{V\xi }}^{{\left(m\right)}^{-1}}\mathbf{X}\right)}^{-1}{\mathbf{X}}^T{\mathbf{V\xi }}^{{\left(m\right)}^{-1}}{\mathbf{\xi }}^{\left(m\right)}$

8. Update ${\widehat{\mathbf{\theta }}}^{\left(m\right)}$ using ${\widehat{\mathbf{\theta }}}^{\left(m+1\right)}={\widehat{\mathbf{\theta }}}^{\left(m\right)}-{\left[\frac{{\partial }^{2}l\left(\mathbf{\beta }\right)}{\mathrm{\partial \theta \partial }{\mathbf{\theta }}^T}\right]}_{\mathbf{\beta }={\widehat{\mathbf{\beta }}}^{\left(m\right)}}^{-1}{\left.\frac{\partial l\left(\mathbf{\beta }\right)}{\mathrm{\partial \theta }}\right|}_{\mathbf{\beta }={\widehat{\mathbf{\beta }}}^{\left(m\right)}}$

   The gradient vector for $\mathbf{\theta }=\underset{1\le k\le K-1}{\mathit{row}}\underset{1\le k^{\prime }\le K-1}{\mathit{col}}\left({\sigma }_{k{k}^{\prime }}^2\right)$ is $\mathbf{b}=\frac{\partial l\left(\mathbf{\beta }\right)}{\mathrm{\partial \theta }}=\underset{1\le l\le L}{\mathit{col}}\left(\frac{\partial l}{\partial {\theta }_l}\right)$, where ${\theta }_l={\sigma }_k$ for $k=k^{\prime }$ and ${\theta }_l={\sigma }_{k{k}^{\prime }}^2$ for $k\ne k^{\prime }$ and $\frac{\partial l}{\partial {\theta }_l}=-\frac{1}{2}\mathit{tr}\left({\mathbf{V}}^{-1}{\mathbf{V}}_l\right)+\frac{1}{2}{\mathrm{r}}^T\left({\mathbf{V}}^{-1}{\mathbf{V}}_l{\mathbf{V}}^{-1}\right)\mathrm{r}$, where ${\mathbf{r}}^{\left(m+1\right)}={\mathbf{\xi }}^{\left(m\right)}-\mathit{X}{\mathbf{\beta }}^{\left(m+1\right)}$ and ${\mathbf{V}}_l=\frac{\partial \mathbf{V}}{\partial {\theta }_l}$, which for ${\theta }_l={\sigma }_k$ is equal to ${\mathbf{V}}_{{\sigma }_k}=\frac{\partial \mathbf{V}}{\partial {\sigma }_k}=\mathbf{Z}\left({\mathbf{I}}_D\otimes \frac{\partial \mathbf{\theta }}{\partial {\sigma }_k}\right){\mathbf{Z}}^T$ with $\frac{\partial \mathbf{\theta }}{\partial {\sigma }_k}=\underset{1\le k\le K-1}{\mathit{diag}}\left(2{\sigma }_k\right)$ and for ${\theta }_l={\sigma }_{k{k}^{\prime }}^2$ when $k\ne k^{\prime }$ it is ${\mathbf{V}}_{{\sigma }_{k{k}^{\prime }}^2}=\frac{\partial \mathbf{V}}{\partial {\sigma }_{k{k}^{\prime }}^2}=\mathbf{Z}\left({\mathbf{I}}_D\otimes \frac{\partial \mathbf{\theta }}{\partial {\sigma }_{k{k}^{\prime }}^2}\right){\mathbf{Z}}^T$ with $\frac{\partial \mathbf{\theta }}{\partial {\sigma }_{k{k}^{\prime }}^2}=\frac{\partial \underset{1\le k\le K-1}{\mathit{row}}\underset{1\le k^{\prime }\le K-1}{\mathit{col}}\left({\sigma }_{k{k}^{\prime }}^2\right)}{\partial {\sigma }_{k{k}^{\prime }}^2}$.

   The Hessian matrix for $\mathbf{\theta }$ is $\frac{{\partial }^{2}l\left(\mathbf{\beta }\right)}{\partial \mathbf{\theta }\partial {\mathbf{\theta }}^T}=\underset{1\le l\le L}{\mathit{row}}\underset{1\le l^{\prime }\le L}{\mathit{col}}\left(\frac{{\partial }^{2}l}{\partial {\theta }_l\partial {\theta }_{l^{\prime }}}\right)$, where we use $l$ instead of $l\left(\mathbf{\beta }\right)$. Here, $\frac{{\partial }^{2}l}{\partial {\theta }_l\partial {\theta }_{l^{\prime }}}=-\frac{1}{2}\mathit{tr}\left({\mathbf{V}}^{-1}{\mathbf{V}}_{l{l}^{\prime }}-{\mathbf{V}}^{-1}{\mathbf{V}}_l{\mathbf{V}}^{-1}{\mathbf{V}}_{l^{\prime }}\right)-\frac{1}{2}{\mathbf{e}}^T{\mathbf{V}}^{l{l}^{\prime }}\mathbf{e}$ with ${\mathbf{V}}_{l{l}^{\prime }}=\frac{{\partial }^2\mathbf{V}}{\partial {\theta }_l\partial {\theta }_{l^{\prime }}}$ and ${\mathbf{V}}^{l{l}^{\prime }}=\frac{{\partial }^2{\mathbf{V}}^{-1}}{\partial {\theta }_l\partial {\theta }_{l^{\prime }}}={\mathbf{V}}^{-1}\left({\mathbf{V}}_{l{l}^{\prime }}-{\mathbf{V}}_l{\mathbf{V}}^{-1}{\mathbf{V}}_{l^{\prime }}-{\mathbf{V}}_{l^{\prime }}{\mathbf{V}}^{-1}{\mathbf{V}}_l\right){\mathbf{V}}^{-1}$.

   For $K=3$ and $l,l^{\prime }={\sigma }_1,{\sigma }_2,{\sigma }_{12}^2$, the Hessian matrix is $\frac{{\partial }^{2}l\left(\mathbf{\beta }\right)}{\partial \mathbf{\theta }\partial {\mathbf{\theta }}^T}=\left[\begin{array}{l}\frac{{\partial }^{2}l}{\partial {\sigma }_1\partial {\sigma }_1}\frac{{\partial }^{2}l}{\partial {\sigma }_1\partial {\sigma }_2}\frac{{\partial }^{2}l}{\partial {\sigma }_1\partial {\sigma }_{12}^2}\\ \frac{{\partial }^{2}l}{\partial {\sigma }_2\partial {\sigma }_1}\frac{{\partial }^{2}l}{\partial {\sigma }_2\partial {\sigma }_2}\frac{{\partial }^{2}l}{\partial {\sigma }_2\partial {\sigma }_{12}^2}\\ \frac{{\partial }^{2}l}{\partial {\sigma }_{12}^2\partial {\sigma }_1}\frac{{\partial }^{2}l}{\partial {\sigma }_{12}^2\partial {\sigma }_2}\frac{{\partial }^{2}l}{\partial {\sigma }_{12}^2\partial {\sigma }_{12}^2}\end{array}\right]$, and ${\mathbf{V}}_{{\sigma }_1{\sigma }_1}=\frac{{\partial }^2\mathbf{V}}{\partial {\sigma }_1\partial {\sigma }_1}=\mathbf{Z}\left({\mathbf{I}}_D\otimes \frac{{\partial }^2\mathbf{\theta }}{\partial {\sigma }_1\partial {\sigma }_1}\right){\mathbf{Z}}^T=\mathbf{Z}\left({\mathbf{I}}_D\otimes \left[\begin{array}{l}20\\ 00\end{array}\right]\right){\mathbf{Z}}^T$, ${\mathbf{V}}_{{\sigma }_1{\sigma }_2}=\frac{{\partial }^2\mathbf{V}}{\partial {\sigma }_1\partial {\sigma }_2}=\frac{{\partial }^2\mathbf{V}}{\partial {\sigma }_1\partial {\sigma }_2}=\mathbf{Z}\left({\mathbf{I}}_D\otimes \frac{{\partial }^2\mathbf{\theta }}{\partial {\sigma }_1\partial {\sigma }_2}\right){\mathbf{Z}}^T=\mathbf{Z}\left({\mathbf{I}}_D\otimes \left[\begin{array}{l}00\\ 00\end{array}\right]\right){\mathbf{Z}}^T={\mathbf{0}}_{n\left(K-1\right)\times n\left(K-1\right)}$, ${\mathbf{V}}_{{\sigma }_1{\sigma }_{12}^2}=\frac{{\partial }^2\mathbf{V}}{\partial {\sigma }_1\partial {\sigma }_{12}^2}=\mathbf{Z}\left({\mathbf{I}}_D\otimes \frac{{\partial }^2\mathbf{\theta }}{\partial {\sigma }_1\partial {\sigma }_{12}^2}\right){\mathbf{Z}}^T=0_{n\left(K-1\right)\times n\left(K-1\right)}$, ${\mathbf{V}}_{{\sigma }_2{\sigma }_1}={\mathbf{V}}_{{\sigma }_1{\sigma }_2}$ ${\mathbf{V}}_{{\sigma }_2{\sigma }_2}=\frac{{\partial }^2\mathbf{V}}{\partial {\sigma }_2\partial {\sigma }_2}=\mathbf{Z}\left({\mathbf{I}}_D\otimes \frac{{\partial }^2\mathbf{\theta }}{\partial {\sigma }_2\partial {\sigma }_2}\right){\mathbf{Z}}^T=\mathbf{Z}\left({\mathbf{I}}_D\otimes \left[\begin{array}{l}00\\ 02\end{array}\right]\right){\mathbf{Z}}^T$, ${\mathbf{V}}_{{\sigma }_2{\sigma }_{12}^2}=\frac{{\partial }^2\mathbf{V}}{\partial {\sigma }_2\partial {\sigma }_{12}^2}=\mathbf{Z}\left({\mathbf{I}}_D\otimes \frac{{\partial }^2\mathbf{\theta }}{\partial {\sigma }_2\partial {\sigma }_{12}^2}\right){\mathbf{Z}}^T={\mathbf{0}}_{n\left(K-1\right)\times n\left(K-1\right)}$, ${\mathbf{V}}_{{\sigma }_{12}^2{\sigma }_1}={\mathbf{V}}_{{\sigma }_1{\sigma }_{12}^2}$, ${\mathbf{V}}_{{\sigma }_{12}^2{\sigma }_2}={\mathbf{V}}_{{\sigma }_2{\sigma }_{12}^2}$,${\mathbf{V}}_{{\sigma }_{12}^2{\sigma }_{12}^2}=\frac{{\partial }^2\mathbf{V}}{\partial {\sigma }_{12}^2\partial {\sigma }_{12}^2}=\mathbf{Z}\left({\mathbf{I}}_D\otimes \frac{{\partial }^2\mathbf{\theta }}{\partial {\sigma }_{12}^2\partial {\sigma }_{12}^2}\right){\mathbf{Z}}^T={\mathbf{0}}_{n\left(K-1\right)\times n\left(K-1\right)}$

9. Continue until convergence.

B.2 Prediction

We partition the linear predictor vector ${\mathbf{\eta }}_d=\underset{1\le i\le N_d}{\mathit{col}}\underset{1\le k\le K-1}{\mathit{col}}\left({\eta }_{\mathit{dik}}\right)$ into a part ${\mathbf{\eta }}_{\mathit{ds}}=\underset{1\le i_s\le n_d}{\mathit{col}}\underset{1\le k\le K-1}{\mathit{col}}\left({\eta }_{{\mathit{di}}_{s}k}\right)$ containing the $n_d$ known sample values and another part ${\mathbf{\eta }}_{\mathit{dr}}=\underset{1\le i_r\le N_d-n_d}{\mathit{col}}\underset{1\le k\le K-1}{\mathit{col}}\left({\eta }_{{\mathit{di}}_{r}k}\right)$ containing $N_d-n_d$ unknown values we want to predict. We partition the model ${\mathbf{\eta }}_d={\mathbf{X}}_d\mathbf{\beta }+{\mathbf{Z}}_d\mathbf{V}$ in the same way: $\left[\begin{array}{l}{\mathbf{\eta }}_{\mathit{ds}}\\ {\mathbf{\eta }}_{\mathit{dr}}\end{array}\right]=\left[\begin{array}{l}{\mathbf{X}}_{\mathit{ds}}\\ {\mathbf{X}}_{\mathit{dr}}\end{array}\right]\mathbf{\beta }+\left[\begin{array}{l}{\mathbf{Z}}_{\mathit{ds}}\\ {\mathbf{Z}}_{\mathit{dr}}\end{array}\right]\mathbf{V}=\left[\begin{array}{l}{\mathbf{X}}_{\mathit{ds}}\mathbf{\beta }+{\mathbf{Z}}_{\mathit{ds}}\mathbf{V}\\ {\mathbf{X}}_{\mathit{dr}}\mathbf{\beta }+{\mathbf{Z}}_{\mathit{dr}}\mathbf{V}\end{array}\right]$. We use ${\widehat{\mathbf{y}}}_{N_d}={\sum }_{i_s=1}^{n_d}{\mathbf{y}}_{{\mathit{di}}_s}+{\sum }_{i_r=1}^{N_d-n_d}{\widehat{\mathbf{\mu }}}_{{\mathit{di}}_r}$ as a predictor of ${\mathbf{y}}_{N_d}$, where ${\widehat{\mathbf{\mu }}}_{{\mathit{di}}_r}=\underset{1\le k\le K-1}{\mathit{col}}\left({\widehat{\mu }}_{{\mathit{di}}_{r}k}\right)$ with ${\widehat{\mu }}_{{\mathit{di}}_{r}k}=\frac{exp\left({\mathbf{x}}_{{\mathit{di}}_r}{\widehat{\mathbf{\beta }}}_k+{\widehat{v}}_{\mathit{dk}}\right)}{1+{\sum }_{k=1}^{K-1}exp\left({\mathbf{x}}_{{\mathit{di}}_r}{\widehat{\mathbf{\beta }}}_k+{\widehat{v}}_{\mathit{dk}}\right)}$ for crop $k=1,2,\cdots ,K-1$ and ${\widehat{\mu }}_{{\mathit{di}}_{r}K}=\frac{1}{1+{\sum }_{k=1}^{K-1}exp\left({\mathbf{x}}_{{\mathit{di}}_r}{\widehat{\mathbf{\beta }}}_k+{\widehat{v}}_{\mathit{dk}}\right)}$ for crop $K$. Here, ${\widehat{\mathbf{\beta }}}_k$ is the ${\mathbf{\beta }}_k$ estimate at the convergence iteration ${\widehat{\mathbf{\beta }}}^{\left(m\right)}=\underset{1\le k\le K-1}{\mathit{col}}\left({\widehat{\mathbf{\beta }}}_k^{\left(m\right)}\right)$ and ${\widehat{v}}_{\mathit{dk}}$ is the predictor of the small-area random effect $v_{\mathit{dk}}$ in $\mathbf{v}=\underset{1\le d\le D}{\mathit{col}}\underset{1\le k\le K-1}{\mathit{col}}\left(v_{\mathit{dk}}\right)$ defined by $\widehat{\mathbf{v}}=\underset{1\le d\le D}{\mathit{col}}\underset{1\le k\le K-1}{\mathit{col}}\left({\widehat{v}}_{\mathit{dk}}\right)=\widehat{\mathbf{\Sigma }}{\mathbf{Z}}^T\widehat{\mathbf{V}}{\mathbf{\xi }}^{-1}\left(\mathbf{\xi }-\mathbf{X}\widehat{\mathbf{\beta }}\right)$, where $\widehat{\mathbf{\Sigma }}=\widehat{\mathbf{\theta }}\otimes {\mathbf{I}}_{\mathbf{D}}$ and $\widehat{\mathbf{\theta }}$ is the $\mathbf{\theta }$ estimate at the convergence iteration ${\widehat{\mathbf{\theta }}}^{\left(m\right)}=\underset{1\le k\le K-1}{\mathit{row}}\underset{1\le k^{\prime }\le K-1}{\mathit{col}}\left({\widehat{\sigma }}_{k{k}^{\prime }}^{\left(m\right)2}\right)$.

B.3 Mean squared error

The mean squared error of predictor ${\widehat{\mathbf{y}}}_{N_d}$ is $\mathit{MSE}\left({\widehat{\mathbf{y}}}_{N_d}\right)=\mathit{MSE}\left({\widehat{\mathbf{\tau }}}_{\mathit{dr}}\right)+{\sum }_{i_r=1}^{N_d-n_d}{\mathbf{Vy}}_{{\mathit{di}}_r}$. The term $\mathit{MSE}\left({\widehat{\mathbf{\tau }}}_{\mathit{dr}}\right)$ can be approximated by a Taylor series expansion of ${\widehat{\mathbf{\mu }}}_{\mathit{di}}=h\left({\widehat{\mathbf{\eta }}}_{\mathit{di}}\right)$ around ${\mathbf{\mu }}_{\mathit{di}}=h\left({\mathbf{\eta }}_{\mathit{di}}\right)$, where ${\mathbf{\eta }}_{\mathit{di}}={\mathbf{X}}_{\mathit{di}}\mathbf{\beta }+{\mathbf{Z}}_{\mathit{di}}\mathbf{V}$.

As a result, we have ${\widehat{\mathbf{\mu }}}_{\mathit{di}}-{\mathbf{\mu }}_{\mathit{di}}\simeq {\mathbf{H}}_{\mathit{di}}\left({\widehat{\mathbf{\eta }}}_{\mathit{di}}-{\mathbf{\eta }}_{\mathit{di}}\right)$. Let ${\mathbf{\tau }}_{\mathit{dr}}^T={\sum }_{i_r=1}^{N_d-n_d}{\mathbf{H}}_{{\mathit{di}}_r}{\mathbf{\eta }}_{{\mathit{di}}_r}={\sum }_{i_r=1}^{N_d-n_d}\left({\mathbf{H}}_{{\mathit{di}}_r}{\mathbf{X}}_{{\mathit{di}}_r}\mathbf{\beta }+{\mathbf{H}}_{{\mathit{di}}_r}{\mathbf{Z}}_{{\mathit{di}}_r}\mathbf{V}\right)={\mathrm{{\rm M}}}_{\mathit{dr}}\mathbf{\beta }+{\mathrm{K}}_{\mathit{dr}}\mathbf{V}$, where ${\mathrm{{\rm M}}}_{\mathit{dr}}={\sum }_{i_r=1}^{N_d-n_d}{\mathbf{H}}_{{\mathit{di}}_r}{\mathbf{X}}_{{\mathit{di}}_r}$ and ${\mathrm{K}}_{\mathit{dr}}={\sum }_{i_r=1}^{N_d-n_d}{\mathbf{H}}_{{\mathit{di}}_r}{\mathbf{Z}}_{{\mathit{di}}_r}$. The approximation of $\mathit{MSE}\left({\widehat{\mathbf{\tau }}}_{\mathit{dr}}\right)=\mathit{MSE}\left({\widehat{\mathbf{\tau }}}_{\mathit{dr}}^T\right)$ is $\mathit{MSE}\left({\widehat{\mathbf{\tau }}}_{\mathit{dr}}\right)={\mathbf{G}}_{1d}\left(\mathbf{\theta }\right)+{\mathbf{G}}_{2d}\left(\mathbf{\theta }\right)+{\mathbf{G}}_{3d}\left(\mathbf{\theta }\right)+{\mathbf{G}}_{4d}\left(\mathbf{\theta }\right)$, where ${\mathbf{G}}_{1d}\left(\mathbf{\theta }\right)={\mathrm{K}}_d{\mathrm{TK}}_d^T$ ${\mathbf{G}}_{2d}\left(\mathbf{\theta }\right)={\mathrm{\Lambda }}_d{\mathrm{P\Lambda }}_d^T$, ${\mathbf{G}}_{3d}\left(\mathbf{\theta }\right)={\left(\frac{\partial {\mathrm{\Gamma }}_d}{\mathrm{\partial \theta }}\right)}^T\mathbf{V}\left(\frac{\partial {\mathrm{\Gamma }}_d}{\mathrm{\partial \theta }}\right){\mathbf{V}}_{\mathbf{\theta }}^{-1}$ and ${\mathbf{G}}_{4d}\left(\mathbf{\theta }\right)={\mathbf{Vy}}_{\mathit{dr}}={\sum }_{i_r=1}^{N_d-n_d}{\mathbf{Vy}}_{{\mathit{di}}_r}$, with $\mathrm{T}={\left({\mathbf{Z}}^T\mathrm{HZ}+{\mathbf{\Sigma }}^{-1}\right)}^{-1}$, $\mathbf{H}=\underset{1\le d\le D}{\mathit{diag}}\underset{1\le i\le N_d}{\mathit{diag}}\left({\mathbf{H}}_{\mathit{di}}\right)$, $\mathrm{P}=\mathbf{V}\widehat{\mathbf{\beta }}={\left({\mathbf{X}}^T{\mathbf{V}}^{-1}\mathbf{X}\right)}^{-1}$, $\mathbf{V}={\mathrm{Z\Sigma Z}}^T+\mathbf{Ve}$, ${\mathrm{\Lambda }}_d={\mathbf{m}}_d-{\mathrm{K}}_d{\mathrm{TZ}}^T\mathbf{VeX}$, and $\frac{\partial {\mathrm{\Gamma }}_d}{\mathrm{\partial \theta }}={\mathbf{V}}^{-1}{\mathrm{ZK}}_d^T$. Here, ${\mathbf{V}}_{\mathbf{\theta }}$ denotes the Fisher information of the variance component parameters $\mathbf{\theta }$.

As estimator of the mean squared error of predictor ${\widehat{\mathbf{y}}}_{N_d}$, we use the asymptotically unbiased estimator proposed in [23], $M\widehat{S}E\left({\widehat{\mathbf{y}}}_{N_d}\right)={\mathbf{G}}_{1d}\left(\widehat{\mathbf{\theta }}\right)+{\mathbf{G}}_{2d}\left(\widehat{\mathbf{\theta }}\right)+2{\mathbf{G}}_{3d}\left(\widehat{\mathbf{\theta }}\right)+{\mathbf{G}}_{4d}\left(\widehat{\mathbf{\theta }}\right)$, where $\widehat{\mathbf{\theta }}$ is the $\mathbf{\theta }$ estimate at the convergence iteration.