Agricultural Statistics & Data Science · Statistical computing, data science and practice · also Geostatistics, Spatial Data Analysis, Remote Sensing Statistics
Statistics for data with a location: spatial autocorrelation, variograms and kriging, spatial regression, point patterns, area-frame sampling, accuracy assessment of crop maps and satellite-based area and yield estimation. Remote-sensing physics and GIS operations are taught in the Agricultural Engineering track; this is the statistics.
In the exams: AICE-JRF Agricultural Physics 'spatial statistics: variogram and interpolation' (present); national crop-estimation programmes (applied); light in national exams
Almost every agricultural observation has a location: a soil sample taken at a grid point, the yield of a field, the count of diseased plants in a district, the reflectance of a satellite pixel. Observations close together tend to be more alike than observations far apart, a pattern that geographers summarise as Tobler's first law. Classical statistics assumes independent observations; with spatial data that assumption fails, and methods that ignore it give standard errors that are too small and maps that are wrong. Spatial statistics is the set of methods that describe, model and exploit this dependence.
This subject covers the three main kinds of spatial data and their methods: geostatistical data (variograms and kriging), areal data (spatial autocorrelation and spatial regression) and point patterns. It then turns to remote sensing as a statistical problem: area-frame sampling with satellite stratification, image classification, accuracy assessment of maps, and the estimation of crop area and yield by combining satellite data with ground surveys. Remote-sensing physics, image processing and GIS operations are taught in the Agricultural Engineering track (Remote Sensing, GIS and GNSS for Land and Water, and Precision Agriculture and System Management); this subject is the statistics behind them.
In examinations the subject appears as variograms and spatial interpolation in agricultural physics and soil science papers, as area-frame and crop-estimation methods in sample-survey papers, and in master's and PhD courses on spatial statistics. In practice it is used daily in soil mapping, precision agriculture, field-trial analysis, disease surveillance and national crop-area estimation. The related theory is found in this track's Multivariate Analysis (discriminant analysis and clustering), Mixed Models (spatial analysis of field trials), Crop Estimation Surveys (sampling) and Crop Forecasting and Growth Models (vegetation-index yield models).
Following Cressie's classification, spatial data fall into three types according to how locations arise. Geostatistical data are measurements of a variable that exists everywhere in a region, taken at a set of fixed sample locations: soil organic carbon at grid points, nitrate in groundwater wells, rainfall at stations. The aim is usually to predict the variable at unsampled places and to map it. Areal (lattice) data are values attached to fixed areal units that together cover the region: crop yield per district, disease incidence per village, the reflectance of each pixel. The aim is to describe and explain patterns across units, often with neighbouring units influencing each other. Point-pattern data record the locations of events themselves, which are random: the positions of weeds in a field, diseased trees in an orchard, outbreaks of a livestock disease. The aim is to say whether events are clustered, regular or random, and why. The same phenomenon can be recorded in more than one way, and the type determines the method.
The support of an observation is the size, shape and orientation of the area or volume it represents: a 5 cm soil core, a 10 m by 10 m plot, a 30 m pixel, a district. Values on different supports are not directly comparable: averages over larger supports vary less than point values, and correlations between variables can change with the level of aggregation, the modifiable areal unit problem. Combining point soil samples with 250 m satellite pixels, or district statistics with field data, raises this change-of-support problem, which block kriging and careful aggregation address. Locations must be recorded in a defined coordinate system; distances computed from latitude and longitude in degrees are distorted, so data are usually projected to a planar system (such as UTM, in metres) before computing variograms or neighbour distances.
Spatial models need assumptions that allow dependence to be estimated from a single realisation. A spatial process is second-order stationary if its mean is constant and the covariance between two locations depends only on the separation vector h between them, not on their absolute positions. A weaker assumption, intrinsic stationarity, requires only that differences Z(x + h) - Z(x) have zero mean and a variance that depends only on h; this is what the variogram needs. A process is isotropic if dependence depends only on the distance |h|, not on direction. Soil properties along a slope or along the direction of ploughing, or disease spread with the prevailing wind, are often anisotropic. A trend in the mean, for example soil pH rising steadily across a valley, violates stationarity and should be modelled explicitly (a trend surface, or universal kriging) before the spatial dependence in the residuals is studied.
For areal data, neighbourhood is defined by a spatial weights matrix W, whose element w(i,j) expresses how strongly unit j is a neighbour of unit i; w(i,i) = 0. Contiguity weights set w(i,j) = 1 if units share a boundary: rook contiguity counts shared edges only, queen contiguity counts shared edges or corners. Distance-based weights use a threshold distance, k nearest neighbours, or an inverse-distance function. Weights are often row-standardised so that each row sums to 1, which makes a weighted sum of neighbours a weighted average and makes coefficients comparable across units with different numbers of neighbours. Results can depend on the choice of W, so the choice should be justified (for example, districts that share a market or a river) and sensitivity to it checked.
Moran's I is the most widely used measure of global spatial autocorrelation. With n units, values x(i), mean x-bar and weights w(i,j) summing to W over all pairs, I = (n/W) x [sum over i and j of w(i,j)(x(i) - x-bar)(x(j) - x-bar)] / [sum over i of (x(i) - x-bar) squared]. It resembles a correlation coefficient between each value and its neighbours. Under no spatial autocorrelation its expected value is -1/(n - 1), slightly negative, not zero. Values above the expectation indicate positive autocorrelation (similar values cluster: high yields near high yields), values below indicate negative autocorrelation (a chequerboard pattern). Significance is assessed by a normal approximation under the assumption of normality or of randomisation, or more robustly by a permutation test that reassigns the observed values randomly to locations many times and compares the observed I with the resulting distribution.
Geary's C is based on squared differences between neighbours: C = (n - 1) x [sum over i and j of w(i,j)(x(i) - x(j)) squared] / [2W x sum over i of (x(i) - x-bar) squared]. Its expected value under no autocorrelation is 1; values below 1 indicate positive autocorrelation (neighbours are similar, so their differences are small), and values above 1 negative autocorrelation. C is more sensitive to differences between particular neighbouring pairs, I to the overall pattern. Global statistics summarise the whole region in one number and can miss local clusters. Local indicators of spatial association (LISA), introduced by Anselin, decompose Moran's I into a contribution for each unit, identifying hot spots (high values surrounded by high values), cold spots (low among low) and spatial outliers (high among low or the reverse). The Moran scatterplot, which plots each value against the weighted average of its neighbours, shows the same information graphically; its slope with row-standardised weights equals Moran's I. Because many local tests are made, significance levels must be adjusted for multiple testing.
Positive spatial autocorrelation means each observation carries less independent information than it appears to. A sample of 100 autocorrelated soil cores may contain the information of far fewer independent samples, so ordinary standard errors are too small and tests of differences between treatments or regions are too liberal. In field experiments, fertility gradients make neighbouring plots similar, which is why blocking, and the spatial analysis of field trials in the Mixed Models subject, improve precision. In regression with spatial data, autocorrelated residuals signal a missing spatially structured variable or the need for a spatial model (Chapter 5). Testing residuals with Moran's I is a routine diagnostic.
For geostatistical data the basic tool is the semivariogram (often simply called the variogram), gamma(h) = (1/2) Var[Z(x + h) - Z(x)], half the expected squared difference between values separated by h. If the process is second-order stationary with variance C(0) and covariance function C(h), then gamma(h) = C(0) - C(h): the semivariance rises as covariance falls. The empirical (Matheron) estimator groups pairs of sample points into distance classes (lags) and computes gamma-hat(h) = [1/(2N(h))] x sum over the N(h) pairs at lag h of (z(i) - z(j)) squared. Practical rules: use lags with at least 30 to 50 pairs, and compute the variogram only up to about half the maximum distance between samples, beyond which few pairs remain and estimates are unstable. Robust estimators (Cressie and Hawkins) reduce the influence of outliers, to which the squared differences are very sensitive.
A typical variogram rises from near zero and levels off. The nugget is the apparent semivariance at distance zero (the intercept): it combines measurement error and variation at distances shorter than the closest sample spacing. The sill is the plateau, equal to the total variance of a stationary process; the partial sill is the sill minus the nugget, the spatially structured part. The range is the distance at which the variogram reaches the sill: beyond it, observations are effectively uncorrelated. The range tells a practitioner how far apart samples can be while still informing one another, which guides sampling design for soil mapping and the cell size of management zones. The nugget-to-sill ratio is often used to describe the strength of spatial dependence; in soil science, ratios below about 25 per cent are commonly called strong, 25 to 75 per cent moderate and above 75 per cent weak dependence, a convention rather than a test. A variogram that keeps rising without a sill indicates a trend or a non-stationary process.
The empirical variogram is replaced by a fitted model, because kriging needs values at all distances and the model must be conditionally negative definite so that kriging variances are never negative; arbitrary curves do not guarantee this. Common valid models, each with nugget c0, partial sill c and a distance parameter a, are: the spherical model, which rises nearly linearly and reaches the sill exactly at the range a; the exponential model, gamma(h) = c0 + c[1 - exp(-h/a)], which approaches the sill asymptotically and reaches about 95 per cent of it at a practical range of 3a; the Gaussian model, gamma(h) = c0 + c[1 - exp(-h squared / a squared)], parabolic near the origin and suited to very smooth variables, with a practical range of about the square root of 3 times a; and unbounded models such as the linear and power models for non-stationary variables. Models are fitted by weighted least squares (weights proportional to the number of pairs, divided by the squared model value, following Cressie), by maximum likelihood or by REML. Several structures can be added (nested models) when variation occurs at more than one scale.
Directional variograms, computed for pairs in angular sectors (for example north-south and east-west), reveal anisotropy. In geometric anisotropy, the range differs with direction but the sill does not, and a coordinate transformation makes the process isotropic. In zonal anisotropy, the sill differs with direction, often because of a trend or layered structure. The fitted variogram is checked by leave-one-out cross-validation: each sample is removed in turn and predicted by kriging from the rest. A good model gives a mean error near zero (no bias), a root mean square error as small as possible, and a mean squared standardised error (error squared divided by kriging variance) near 1, which shows that the kriging variances correctly describe the actual errors.
Kriging, named after the South African mining engineer D. G. Krige and formalised by Matheron, predicts the value at an unsampled location x0 as a weighted average of nearby observations, Z-hat(x0) = sum of lambda(i) Z(x(i)), with weights chosen to make the prediction the best linear unbiased predictor (BLUP). In ordinary kriging the mean is assumed constant but unknown; unbiasedness requires the weights to sum to 1, and minimising the prediction-error variance subject to that constraint, using a Lagrange multiplier, gives a set of linear equations in the variogram values between sample points and between each sample point and x0. The solution gives the weights and the kriging variance, the expected squared prediction error at x0, which depends on the configuration of samples and the variogram but not on the observed values themselves. Kriging weights depend on distance and on clustering: two samples close together share information and receive less weight each than an isolated sample at the same distance would. Weights can be negative, and with no nugget effect kriging is an exact interpolator, reproducing observed values at sample points.
Simple kriging assumes the mean is known; it has no unbiasedness constraint, and the weights need not sum to 1. Universal kriging (kriging with a trend) models the mean as a function of coordinates, such as a linear or quadratic trend surface, and kriges the residuals simultaneously; the variogram used must be that of the residuals, not of the raw data, which contain the trend. Kriging with an external drift lets the mean depend on an auxiliary variable known everywhere, such as elevation or a satellite index. Regression kriging fits the trend by regression on covariates and kriges the residuals separately, a widely used approach in digital soil mapping. Block kriging predicts the average over an area (a field or management zone) rather than at a point; block kriging variances are smaller than point kriging variances because the averaging removes the nugget and short-range variation. Indicator kriging transforms the data to 0/1 indicators of exceeding a threshold, for example a critical nitrate level, and kriges the indicators to estimate the probability of exceeding it.
When a secondary variable correlated with the target is measured more densely than the target, co-kriging uses both. Soil organic carbon may be measured at a few hundred points, while apparent electrical conductivity from a proximal sensor or a satellite index is available almost everywhere. Co-kriging requires the variogram of each variable and their cross-variogram, gamma12(h) = (1/2) E[(Z1(x + h) - Z1(x))(Z2(x + h) - Z2(x))], fitted together under a linear model of coregionalisation that keeps the joint model valid. Co-kriging improves predictions most when the secondary variable is strongly correlated with the target and much more densely sampled; when both are measured at the same points, the gain is usually small. In practice regression kriging and kriging with external drift are often preferred because they are simpler to fit.
Inverse distance weighting (IDW) predicts with weights proportional to 1/d to the power p; it is simple and needs no model, but its weights ignore the spatial structure of the data and the clustering of samples, and it gives no measure of prediction error. Splines and other deterministic interpolators share the second limitation. Kriging's advantages are optimal weights under its model and a prediction-error variance, which can be mapped to show where predictions are reliable and where more samples are needed. Its prediction is only as good as the variogram; with few samples (fewer than about 100 is often cited as a minimum for a reliable variogram) the variogram is uncertain and kriging may do no better than simpler methods. Model-based geostatistics fits the covariance parameters by likelihood or Bayesian methods and carries parameter uncertainty into predictions. For precision agriculture, maps of soil properties and yield made by kriging feed the management-zone and variable-rate decisions described in Precision Agriculture and System Management in the Agricultural Engineering track.
Consider district maize yields regressed on fertiliser use and rainfall. If yields in neighbouring districts influence one another (through shared markets, spillover of technology or unmeasured soil and climate), or if omitted variables are spatially smooth, the residuals of ordinary least squares (OLS) are spatially autocorrelated. Moran's I applied to the residuals detects this. The consequences depend on the form. If the dependence is in the response itself, OLS estimates are biased and inconsistent. If it is only in the errors, OLS estimates remain unbiased but are inefficient, and their standard errors are wrong, usually too small. Spatial regression models make the dependence explicit.
Simultaneous autoregressive (SAR) models use the weights matrix W. The spatial lag model is y = rho W y + X beta + e: each unit's value depends directly on the weighted average of its neighbours' values, with rho measuring the strength of this spillover. Because W y is correlated with the error, OLS is inconsistent; the model is fitted by maximum likelihood or by instrumental variables (spatial two-stage least squares). In the lag model, a change in X in one unit affects its neighbours through the feedback, so the marginal effect is not beta alone; direct, indirect (spillover) and total effects are computed from the inverse of (I - rho W). The spatial error model is y = X beta + u, with u = lambda W u + e: the dependence is a nuisance in the errors, typically from omitted spatially structured variables. Anselin's Lagrange multiplier tests, computed from OLS residuals, test for lag and error dependence, and their robust forms help choose between the two. The spatial Durbin model adds W X, the neighbours' covariates, and nests both.
Conditional autoregressive (CAR) models, associated with Besag, specify the distribution of each unit's value conditional on its neighbours: the conditional mean of y(i) given the others is a weighted combination of its neighbours' values, and the conditional variance is inversely related to the number of neighbours. CAR models lead to a joint distribution with a sparse precision (inverse covariance) matrix, which makes them computationally convenient in large problems and natural in Bayesian hierarchical models fitted by Markov chain Monte Carlo or INLA. The intrinsic CAR prior, combined with an unstructured random effect in the Besag-York-Mollie (BYM) model, is the standard model for disease mapping: smoothing the incidence of a livestock disease or crop pest across villages so that small units with few animals or fields do not show extreme rates by chance. SAR and CAR are different models and are not interchangeable, although both describe areal dependence.
Global regression assumes the same relationship everywhere, but the response of yield to nitrogen or to rainfall may differ between regions with different soils and climates. Geographically weighted regression (GWR), developed by Fotheringham, Brunsdon and Charlton, fits a separate weighted regression at each location, with weights given by a kernel that declines with distance from that location (Gaussian or bi-square kernels), so that nearby observations count more. The bandwidth controls smoothness and is chosen by cross-validation or a corrected AIC. The result is a map of local coefficients showing where an effect is strong or weak. GWR is exploratory: local coefficients are often highly correlated with one another, multiple testing is extensive, and apparent variation in coefficients can come from mis-specification rather than real differences. Multiscale GWR allows each covariate its own bandwidth. In field trials, related ideas appear as nearest-neighbour adjustment (Papadakis's method) and the separable autoregressive (AR1 x AR1) residual models treated in Mixed Models.
A point pattern is a set of event locations in a study region: weed plants in a field, trees with a disease in a plantation, farms with an outbreak. The benchmark is complete spatial randomness (CSR), the homogeneous Poisson process with intensity lambda (expected events per unit area): the number of events in any region of area A is Poisson with mean lambda A, and counts in non-overlapping regions are independent. Departures from CSR are of two kinds. Clustering (aggregation) occurs when events attract each other or share a patchy environment: weeds spread from seed patches, a soil-borne disease spreads along rows. Regularity (inhibition) occurs when events repel each other: trees planted on a grid, territorial animals, plants competing for space. Distinguishing true interaction from response to an uneven environment (inhomogeneity) is a central difficulty.
Quadrat analysis divides the region into equal cells and counts events in each. Under CSR the counts are Poisson, so their variance equals their mean; the index of dispersion, variance/mean, is about 1 under CSR, greater than 1 for clustering and less than 1 for regularity, and (number of quadrats - 1) times the index follows approximately a chi-square distribution with (number of quadrats - 1) degrees of freedom. Results depend on quadrat size, because a pattern can be clustered at one scale and regular at another. Nearest-neighbour methods use distances from each event to its closest neighbour. The Clark-Evans ratio R compares the observed mean nearest-neighbour distance with its expected value under CSR, 1/(2 times the square root of lambda): R near 1 suggests randomness, R below 1 clustering, and R above 1 regularity, with a maximum of about 2.15 for a perfect hexagonal arrangement. Edge effects, events near the boundary whose true nearest neighbour lies outside the region, bias these statistics and are corrected by guard areas or weighting.
Second-order methods describe the pattern at all scales at once. Ripley's K function is defined so that lambda K(r) is the expected number of further events within distance r of a typical event. Under CSR, K(r) = pi r squared. Estimated K above pi r squared at distance r indicates clustering at that scale, and below it indicates regularity. The transformation L(r) = square root of (K(r)/pi) makes the CSR benchmark a straight line, L(r) = r, and stabilises the variance; plots of L(r) - r against r are common. Because the sampling distribution of K is complicated, significance is judged with simulation envelopes: many CSR patterns with the same number of events are simulated in the same region, and the observed curve is compared with the range of simulated curves. Pointwise envelopes invite multiple-testing errors; global envelope tests address this. The pair-correlation function, the derivative of K scaled appropriately, shows at which distances events are more or less frequent than random. Cross-K functions describe the relation between two types of events, for example diseased and healthy plants.
When CSR is rejected, models explain the pattern. The inhomogeneous Poisson process lets intensity vary with location or covariates (lambda(x) = exp(beta0 + beta1 times soil moisture), for example), fitted by maximum likelihood as a log-linear model; this separates environmental heterogeneity from interaction. Cluster processes, such as Neyman-Scott and Thomas processes, generate parent points at random and daughter events scattered around them, a natural model for weed patches or disease foci from initial infections. Inhibition processes, such as hard-core and Strauss processes, prevent events from occurring too close together. Kernel density estimation produces a smooth map of intensity from the events, with a bandwidth that controls smoothness. In epidemiology, comparing the intensity of cases with that of controls (for example all farms) gives a relative-risk surface and a test for clustering of disease beyond the clustering of the population at risk.
A sample survey needs a frame from which units are selected. A list frame is a list of farms or holdings; it is efficient when up to date, but in many countries farm lists are incomplete or out of date, and they miss new and informal holdings. An area frame divides the land itself into units that together cover the region without overlap: every piece of land, and therefore every crop field, belongs to exactly one unit. Area frames are complete by construction and do not age as quickly as lists, which makes them the basis of crop-area surveys in many countries. The United States national agricultural statistics service has long used an area frame for its June Agricultural Survey, and the European Union's LUCAS survey uses a systematic grid of points. Area frames are costlier per unit to enumerate, because enumerators must visit the land, and they are less efficient for rare items (such as specialised livestock farms), for which a combined (multiple-frame) design with a list frame is used.
Area-frame units may be segments or points. Segments are pieces of land, either with physical boundaries (roads, rivers, field edges) that enumerators can recognise on the ground, or square segments defined on a map grid. Within a selected segment, the enumerator records the crop on every field, giving area by crop. Points are small locations (for example a 3 m circle) at which the land cover is recorded; the proportion of points under a crop, multiplied by the total area, estimates the crop area. Points are cheaper and avoid boundary problems, but each carries less information. Segments can also be linked to farms (the open, closed and weighted segment estimators) to estimate items such as livestock or production that belong to holdings rather than to land. Sampling designs typically select segments with probability proportional to size within strata, or systematically, and may replicate selections to give simple variance estimates and to rotate units over years.
Stratifying the area frame greatly improves efficiency. Satellite images and maps of land use are used to divide the region into strata of different agricultural intensity: for example, more than 75 per cent cultivated, 50 to 75 per cent cultivated, 15 to 50 per cent, urban and non-agricultural land. Because the variable of interest (area under a crop per segment) varies much more between such strata than within them, stratified sampling with appropriate allocation reduces the variance compared with simple random sampling of segments. More samples are allocated to strata with high cultivation and high variance (Neyman allocation), and few to non-agricultural strata. Strata must be defined before sampling and remain fixed; reclassifying units after seeing data destroys the design. The frame and its stratification are typically rebuilt every 10 to 20 years as land use changes. Satellite data also serve after the survey, as an auxiliary variable for regression estimation, the subject of Chapter 10.
A multispectral satellite image records, for each pixel, the reflectance in several spectral bands (blue, green, red, near-infrared, short-wave infrared and others). Each pixel is therefore a vector observation in a feature space whose axes are the bands, and each land-cover class (wheat, mustard, fallow, water, built-up) occupies a region of that space described by its spectral signature. Classification assigns each pixel to a class, which is exactly the discriminant and classification problem of Multivariate Analysis. Multi-temporal images, which add the same bands on several dates, greatly improve separation of crops that look alike on one date but follow different growth calendars; time series of vegetation indices are especially useful. Texture, shape and context add further features.
Supervised classification uses training samples, pixels or polygons whose class is known from ground visits or very-high-resolution imagery, to estimate each class's distribution. The maximum likelihood classifier assumes each class is multivariate normal with its own mean vector and covariance matrix, computes the likelihood of each pixel under each class, and assigns the pixel to the class with the highest likelihood or, when prior probabilities (expected class proportions) are used, the highest posterior probability by Bayes' rule. It is the quadratic discriminant rule of multivariate statistics; with equal covariance matrices it reduces to the linear rule. Minimum-distance classifiers use only class means. Non-parametric classifiers, especially random forests and support-vector machines (treated in Statistical and Machine Learning for Agriculture), do not assume normality and handle many correlated features; they have largely replaced maximum likelihood in operational crop mapping. Training samples must be representative of each class's variability across the region, and enough must be held back, or collected separately, for accuracy assessment.
Unsupervised classification groups pixels by clustering in feature space, for example with k-means or the ISODATA algorithm, which splits and merges clusters according to their spread and size. The analyst then labels the clusters by comparing them with ground knowledge. It needs no training data at the start but may produce clusters that do not match the classes of interest. Object-based image analysis first segments the image into objects (groups of adjacent, similar pixels, ideally corresponding to fields) and then classifies objects using their mean spectra, texture and shape, which reduces salt-and-pepper noise and suits high-resolution images. In smallholder landscapes, fields may be smaller than a pixel: mixed pixels contain several cover types, and spectral unmixing estimates the fraction of each cover within a pixel by modelling its reflectance as a combination of pure endmember spectra.
Classification errors arise from spectral confusion between similar crops, mixed pixels at field edges, cloud and haze, differences in sowing date that put the same crop at different stages, unrepresentative training data, and geometric misregistration between images and ground data. Errors are rarely random: they cluster at boundaries and in particular landscapes, and they often favour one class over another. This is why a classified map cannot be treated as a census of land cover; its errors must be measured with an independent sample (Chapter 9) and accounted for when the map is used to estimate area (Chapter 10).
Accuracy is assessed by comparing the map with reference data, the best available determination of the true class, for a probability sample of locations independent of the training data. Results are summarised in an error (confusion) matrix. Using the common convention (Congalton), rows are map classes and columns are reference classes; the diagonal holds correctly classified samples. Overall accuracy is the sum of the diagonal divided by the total number of samples. User's accuracy for a class is the diagonal count divided by its row total: the probability that a pixel mapped as wheat is really wheat on the ground, what a map user cares about; one minus user's accuracy is the commission error. Producer's accuracy is the diagonal count divided by its column total: the probability that real wheat is mapped as wheat, what the map producer cares about; one minus producer's accuracy is the omission error. Because conventions for rows and columns vary, a report must state which is which.
Cohen's kappa coefficient was long reported with every error matrix: kappa = (p(o) - p(e)) / (1 - p(e)), where p(o) is the observed overall accuracy and p(e) is the chance agreement, computed as the sum over classes of the product of the row and column proportions. It was presented as accuracy corrected for chance. Pontius and Millones and others have argued that kappa is not useful for map accuracy: the chance agreement it removes is not a sensible baseline for a map, kappa is highly correlated with overall accuracy and adds little information, and it does not reveal whether errors come from wrong quantities of each class or from wrong locations. Current good practice (Olofsson and colleagues, 2014) recommends reporting the error matrix in terms of estimated area proportions, with overall, user's and producer's accuracies and their confidence intervals, rather than relying on kappa. Students should know how to compute kappa, because examinations ask it, and also why it is now discouraged.
Accuracy statistics are estimates and need a probability sampling design. Simple random and systematic samples are easy but may contain few samples of rare classes. Stratified random sampling, with the map classes as strata, guarantees enough samples in each class; a common recommendation is a minimum of about 50 to 100 samples per class, with more for large or important classes. When strata are sampled at different rates, the error matrix must be weighted by the stratum (map-class) area proportions before computing accuracies; an unweighted matrix gives biased producer's and overall accuracies. Reference data must be of higher quality than the map, collected at the right time, and interpreted by people who have not seen the map label, to avoid confirmation bias.
Counting the pixels of a class and multiplying by pixel area gives the mapped area, which is biased whenever omission and commission errors do not balance; for a crop that is under-mapped, it is an underestimate. The accuracy sample allows an unbiased estimate. With map classes as strata, the estimated proportion of area truly in class j is p-hat(j) = sum over map classes i of W(i) n(i,j)/n(i), where W(i) is the mapped area proportion of map class i, n(i,j) is the number of samples mapped as i and found to be j, and n(i) is the number of samples in map class i. The estimated area is p-hat(j) times the total area, and a standard error follows from stratified sampling formulae, giving a confidence interval. In this approach the map serves as a stratification and the reference sample provides the estimate: the map improves precision, the sample removes bias. This is the same logic as the regression estimators of Chapter 10.
The classical statistical design for crop-area estimation combines an area-frame ground survey with a classified satellite image. In each sampled segment, the enumerated crop area y is recorded, and the classified image gives the number (or area) of pixels classified as the crop, x. For the whole region, the total classified area X is known. The regression estimator of the mean crop area per segment is y-bar-reg = y-bar + b(X-bar - x-bar), where y-bar and x-bar are the sample means, X-bar is the population mean of classified area per segment and b is the least-squares slope of y on x in the sample. Its variance is approximately (1 - rho squared) times that of the direct ground estimate y-bar, where rho is the correlation between ground and classified areas across segments; the relative efficiency is about 1/(1 - rho squared). With rho = 0.9, the satellite data are worth about 5.3 times as many ground segments. The approach, developed at the US Department of Agriculture and used in the EU's MARS work, is unbiased in large samples even if the classification is imperfect, because the ground survey defines the target and the image only reduces variance.
Two related approaches use the error matrix to correct pixel counts. Direct calibration uses the proportion of each mapped class found in each true class (from the accuracy sample) to convert mapped areas into estimates of true areas; inverse calibration uses the proportions in the other direction. Inverse calibration is generally preferred when the sample is a random sample of the region, because it does not require inverting the matrix. Satellite data are also valuable for small-area estimation. National surveys give reliable estimates for large regions but too few samples for each district or block. Area-level models such as the Fay-Herriot model combine the direct survey estimate for each small area with a regression prediction from satellite covariates, weighting each according to its precision: areas with many samples lean on the survey, areas with few lean on the model. Unit-level models do the same at the level of segments or fields. These methods are increasingly used to produce district and sub-district crop statistics.
Satellite-based yield estimation relates yield measured on the ground, usually from crop-cutting experiments, to vegetation indices, weather and other covariates. Field-level models regress crop-cutting yields on indices such as peak or integrated NDVI computed for each sampled field, and apply the model to all fields of the crop; aggregated to a district, the result can be compared with or combined with the survey estimate. Smart sampling uses satellite-derived yield proxies to select crop-cutting experiments (stratifying fields by expected yield) or to reduce their number, and combines them with model predictions through regression or small-area estimators. The modelling side, including vegetation-index regressions, crop simulation and machine learning, is treated in Crop Forecasting and Growth Models and in Statistical and Machine Learning for Agriculture.
Every satellite-based estimate should come with a standard error or confidence interval derived from a probability sample, and the assumptions behind any model should be stated. Pitfalls include using the same ground data to train a classifier and to assess or correct it; ignoring spatial autocorrelation in model residuals, which makes model-based intervals too narrow; mismatches between the date of imagery and the date of ground observation; and treating the area of a map as the area of a crop. In crop insurance and policy, where money depends on estimates, the design must be documented and auditable. The guiding principle of this subject is that remote sensing improves estimates most when it is combined with sound sampling, not when it replaces it.
A student must be able to compute Moran's I or Geary's C for a small example and interpret it, sketch and read a variogram (nugget, sill, range), name valid variogram models, set up the ordinary kriging system, explain the difference between spatial lag and spatial error models, compute an index of dispersion or Clark-Evans ratio, and compute overall, user's and producer's accuracy and kappa from an error matrix. In a thesis the student fits and cross-validates a variogram for soil data, or assesses a crop map and estimates crop area with confidence intervals.
Five years into work at a soil survey, a remote-sensing centre or a precision-agriculture company, the statistician designs sampling grids from variogram ranges, produces kriged maps of soil nutrients with uncertainty maps, builds and validates crop classification models, runs accuracy assessments with stratified reference samples, and turns maps into error-adjusted area estimates. The practitioner explains to clients why a beautiful map can still give a biased area estimate.
A senior official or programme head with 25 years of experience designs national systems: the area frame and its stratification, the integration of satellite classification with ground surveys, district-level estimates for insurance and policy, and standards for reporting uncertainty. Decisions at this level concern budgets for ground data versus imagery, the frequency of frame renewal, auditability of estimates that trigger payments, and the governance of methods across agencies.
A professor or research leader works on model-based geostatistics, Bayesian spatial and spatio-temporal models for disease and yield, small-area estimation with satellite covariates, change-of-support problems in combining field and satellite data, and design-based inference for map accuracy. The research leader supervises theses, reviews remote-sensing papers whose accuracy claims rest on non-probability samples, and advises agencies on statistically defensible use of Earth observation.
Spatio-temporal models, Bayesian geostatistics fitted with INLA or MCMC, multivariate geostatistics beyond the linear model of coregionalisation, and change-of-support theory are introduced only in outline. Remote-sensing physics, atmospheric correction and GIS operations are left to the Agricultural Engineering track. Thresholds such as the nugget-to-sill classes and minimum samples per class are conventions from the literature, not standards, and current national area-frame designs should be read from the agencies' own documents.
Original practice questions written to the syllabus and pattern of the exams named in each question's examTags. They are not past papers and do not reproduce any board's questions. Every numerical answer was computed with a script; the values needed are stated in the question.
These notes are free to read here, with no account.
What the Laisha app adds:
Free for ten days, then a subscription. We would rather write that here than have you find it later.