Published online November 20, 2006

Sensitivity Analysis of the Nonparametric Nearest Neighbor Technique to Estimate Soil
Water Retention
A. Nemes,* W. J. Rawls, Ya. A. Pachepsky, and M. Th. van Genuchten

A. Nemes, Univ. of California, Dep. of Environmental Sciences, Riverside, CA 92521; A. Nemes and W.J. Rawls, USDA-ARS Hydrology
and Remote Sensing Lab., 10300 Baltimore Ave., Bldg. 007, BARCWest, Beltsville, MD 20705; Ya.A. Pachepsky, USDA-ARS Environmental Microbial Safety Lab., Powder Mill Road, Bldg. 173,
BARC-East, Beltsville, MD 20705; and M.Th. Van Genuchten,
USDA-ARS George E. Brown Jr. Salinity Lab., 450 West Big Springs
Rd., Riverside, CA 92501. Received 25 Jan. 2006. *Corresponding
author (anemes@hydrolab.arsusda.gov).

the case of small sample size in the development data
set. The equations need to be redeveloped and republished, should new data become available, and users
are not able to simply include any additional data sets
to improve performance for their site-specific range ofsoil properties.
An alternative approach for such estimations is theuse of nonparametric techniques. Such techniques are
based on pattern recognition rather than on fitting equations to data. One of these techniques is the (k-)nearest neighbor (k-NN), which has been widely applied inpattern recognition and statistical classification tasks(Dasarathy, 1991). This technique belongs to the groupof “lazy learning algorithms”. It is “lazy” in that it passively stores the data until the time of application; allcalculations are performed only when estimations need
to be generated. Applications of this technique can befound in the literature of many fields; some recent applications in the fields of meteorology and hydrology are
Yakowitz (1993), Lall and Sharma (1996), Tarboton
et al. (1998), Sharma and O’Neill (2002), Harrold et al.
(2003a, 2003b), and Mehrotra and Sharma (2006).
In a soil physical and hydrophysical context, asimilarity-based k-NN type technique has been applied
successfully by Nemes et al. (1999) to interpolate soil
particle-size distributions. They found this technique toperform well while estimating the missing 50-mm particle fraction for many European soils, which later served
as input to soil hydraulic PTFs. Jagtap et al. (2004)
introduced a dynamic k-NN technique to estimate thedrained upper limit and lower limit of plant water availability from field-measured soil water retention information. They compared their model to existing soil
hydraulic PTFs to make estimations for their data set
and concluded that the k-NN technique performs better
than three published regression-type parametric PTFs.
Most recently, Nemes et al. (2006) estimated soil watercontent at 233 and 21500 kPa matric potentials using ak-NN variant and a hierarchical set of input data thatoriginated from the USA. In this study, the k-NN technique was found to be competitive with neural network(NNet) models that are cited as probably the most advanced and accurate PTF technique of the day.
Application of the k-NN technique means identifying
and retrieving the nearest (most similar) stored objects
to the target object. The quality of such estimations
depends on, among others, which objects are considered
to be the nearest to the target object. An understanding of how this technique works suggests, however, that

Published in Vadose Zone Journal 5:1222–1235 (2006).
Original Researchdoi:10.2136/vzj2006.0017ª Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

Abbreviations: OM, organic matter; SSC, sand, silt, and clay contents(soil texture); NNet, neural network; k-NN, k-nearest neighbor technique; PTF, pedotransfer function; RMSR, root mean squared residual; u33, water retention at 233kPa matric potential; u1500, water
retention at 21500 kPa matric potential.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

ABSTRACTA k-nearest neighbor (k-NN) nonparametric algorithm variant was
earlier applied successfully to estimate soil water retention. In thisstudy, we tested the sensitivity of that k-NN variant to different data
and algorithm options, such as: (i) estimations made to soils with differing distribution of properties; (ii) the use of different sample weighting methods; (iii) the number of ensembles we developed; (iv) data
density in the reference data set; (v) the presence of outliers in thereference data set; (vi) unequal weighting of input attributes; and (vii)the addition of locally specific data to the reference data set. We used ahierarchical set of input attributes and data set sizes to develop
ensembles of predictions using multiple randomized subset selections.
The k-NN technique performed comparably well as neural networkmodels developed on the same data. Using .50 ensemble members
did not improve the results any further. The k-NN technique showed
little sensitivity to the choice of sample weighting methods and to
suboptimal weighting of input attributes. Differences in data density inparts of the reference data set did not substantially impact estimation
errors. Estimations substantially improved for locally specific datawhen some local samples were included in the reference data set, while
estimations for other samples remained almost unaffected. The k-NN
technique shows a large degree of stability and insensitivity to different
settings and options, can easily adopt new data without the need toredevelop equations, and is an effective alternative to other techniques
to estimate soil water retention.

M

ODELING WATER and solute transport has become
an important part of simulating agricultural productivity as well as environmental quality. The use of models, however, is often hindered by the lack of information
on soil hydraulic properties. For many applications, theestimation of those properties using pedotransfer functions (PTFs) is a feasible alternative to costly and timeconsuming measurements.
One common feature of today’s PTFs is that they are
all based on some parametric approach, i.e., they areequations with parameters found from fitting thoseequations to data. Identifying the right equation andensuring that the associated probability distributionsof errors will be similar across the data space is not always easy. Estimation results can be heavily biased in

1222

1223

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

there are several factors that may potentially have agreat influence on which objects are ruled to be nearest,and how their influence on the final estimate is accounted for. A standard k-NN does not perform attribute selection; it allows irrelevant or interacting inputs
to have as much effect on the distance calculation as anyother useful inputs. Some inputs may also have a (considerably) wider numerical range of data than others. A
unit change in one input variable may have a muchlarger influence on the distance measure than the samechange in the other. Such concerns lead to the introduction of data normalization systems and different attribute weighting systems in more recent k-NN variants(e.g., Wettschereck et al., 1997; Mehrotra and Sharma,2006). Different recommendations exist for the weighting of the retrieved nearest objects while formulating theoutput of the k-NN technique. Simple averaging of theoutput attributes of the retrieved neighbors is probably
the simplest solution, but it does not consider any difference in the resemblance of the retrieved neighbors to
the target object. Lall and Sharma (1996) recommended
a method that assigns weights based on the rank of eachof the retrieved k neighbors in being the nearest to thetarget object. Yates et al. (2003) derived realizations of
the output variable by the selection of a single neighborfrom the retrieved k neighbors using a random number
generator. Nemes et al. (2006) used a weighting methodthat assigns weights based on the distances of the retrieved k neighbors from the target object. Such methods have not previously been tested against each otherwhile being applied to estimate unsaturated soil hydraulic properties.
It is unknown to what extent estimations are sensitive
to local data density in the underlying development(reference) data set. In other words, it has not been
quantified how reliable the estimates are for a type ofsoil that can be found in the reference data set but
is relatively poorly represented. The above weighting
methods may perform differently, and the estimations
may be biased or generally be less accurate for cases thatare poorly represented in the reference data set. Pedotransfer functions are expected to give reliable estimatesfor soils that are from the same population as the data
set used to develop the PTF; however, PTFs are frequently sought and tested that work well for soils thatoriginate from a population different than that used todevelop the PTF. Schaap and Leij (1998) showed thegreat extent to which PTFs could be dependent on theorigin (and distribution) of the soils in the development
and application data sets. In the study of Nemes et al.
(2006), the k-NN technique was tested using soils thatwere not in the reference data set, but that originatedfrom the same data source and distribution.
The approach of using an ensemble of estimation orforecast techniques is widely used in meteorology (e.g.,
Molteni et al., 1996; Houtemaker et al., 1996; Palmer
et al., 2004). The approach essentially means averaging
the predictions of a number of models that are applied
to the data simultaneously. The rationale behind suchapproach is that the use of a particular single model isoften not justifiable. Using multiple models, featuresthat are consistent among ensemble members will bepreserved through averaging, while those that are inconsistent will be reduced in amplitude. In the meantime,
the output of each member can be viewed as a potential
sample from the outcome of the estimations and be used
to calculate estimation uncertainty. In subsurface hydrology, Ye et al. (2004) suggested averaging of thespatial variability models in unsaturated fractured tuff.
Guber et al. (2006) tested the performance of an ensemble of 22 PTFs to estimate soil water retention
against measured data and subsequently used such datain soil water flow simulations with success.
An alternative form of ensemble predictions is when
one uses multiple realizations of the PTF development
data set, obtained by, for example, randomized subset selection or bootstrapping. In this way, multiple
subestimates are obtained that can be interpreted as
samples from the statistical distribution of estimation
outcomes and can subsequently be used to characterize
estimation uncertainty. Examples of such estimations
can be found in, e.g., Schaap and Leij (1998), Schaap
et al. (1999), Nemes et al. (2003, 2006), and Baker(2005). It is not known, however, how many ensembles
are minimally needed to obtain estimates that are notsignificantly changed by the addition of an additional
ensemble member. Generating an excessive number ofensembles may not be beneficial in terms of improving
the estimations, while it increases computation time.
The objective of this study was to test the sensitivityof the k-NN variant introduced by Nemes et al. (2006)
to different algorithm options and to differences in theproperties of the underlying data. We also compared theperformance of the k-NN technique to the performanceof NNet models developed using the same data.
MATERIALS AND METHODS
Soil Data
This work comprises seven case studies. In most case studies, we used data from two data sets described below. In Case
Study 7 we used an additional third data set. The first dataset encompasses 2125 soil horizons that were selected fromthe NRCS Soil Characterization Database (Soil Survey Staff,
1997), according to the following criteria: (i) mineral soil horizons were selected from the contiguous USA having horizon notation A, A1, and Ap (and their derivatives), with thecondition that the top of the horizon was at the soil surface;(ii) organic matter (OM) content of the selected soils waslimited to 1 to 15%, and (iii) bulk density (Db) was limited to
0.5 to 2.0 g cm23. Selected soil properties were the following:
sand (50–2000 mm), silt (2–50 mm), and clay content (,2 mm)
according to the USDA classification system (Soil Survey Staff,
1951), Db, OM content, and retained (volumetric) water at
233 and 21500 kPa matric potentials, (u33 and u1500, respectively), with no missing data allowed in any of the fields.
Such matric potentials were chosen as those are often used
to approximate field capacity and the wilting point when calculating plant-available water, and thus are often preferred
points in water retention curve (WRC) determinations in thelaboratory. Measured WRC data at those matric potentials canbe found frequently in many soil hydraulic databases worldwide. No entries were allowed that showed obvious inconsistency in physical or hydraulic data (sand 1 silt 1 clay 5/ 1; u33 ,Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1224

VADOSE ZONE J., VOL. 5, NOVEMBER 2006

u1500; [(1 2 [Db]/2.65) 2 u33] , 0). This data set is referred to
below as NRCS, and has been used as reference data as well as
to provide data for testing the estimations. Samples have been
randomly drawn to be either a member of the reference data
set or a test data set. We elected to use 435 samples, i.e., |20%of all data, as test data in each case. We used different sizesof reference data sets to evaluate the effect of the size of thereference data set on each examined factor. Samples weredrawn to be members of reference data sets of 100, 200, 400,
800, and 1600 samples. All random data selections were repeated 200 times to allow the development of an ensembleof PTF estimations. By using a sufficiently large number ofensembles, the impact on the final estimation results of any
single ensemble (i.e., any particular data set division) can be
minimized. Optimization of the number of ensembles was partof this study.
A second data set from the European HYPRES database(Wösten et al., 1999) served as an alternative test data set inmost of the case studies. These data originated from a different
geographical area and show a different data distribution thanthe NRCS data set. The HYPRES data were donated by researchers of some 20 institutions of different countries in
Europe. The selection process and criteria were the same as for
the NRCS data set, with one exception: there are differences in
the water retention data reported by the different sources. We
allowed water retention data reported in the ranges of 230 to
234 kPa and 21500 to 21600 kPa matric potentials to represent u33 and u1500, respectively. Altogether, we selected
435 samples from HYPRES that were used solely as an independent second test data set, referred to as HYPRES.
An additional third data set was used in Case Study 7 below.
This data collection originated from Brazil and has been used
before by, e.g., Tomasella et al. (2000, 2003). To avoid thepresence of outliers, the same limitations as for the NRCS data
set were imposed on the source data collection. This yielded a
set of 428 samples, which are referred to as BRAZ.
Table 1 shows the summary statistics of selected soil attributes of the selected data sets. The data sets contain data on awide range of soils, in terms of the shown soil attributes. The
average sand content of the HYPRES soils is about 10% larger
than that of the NRCS soils, whereas its silt content is less by
Table 1. Summary statistics of selected soil attributes† in the data
sets.
Sand

Silt
kg kg

21

Clay

Db

g cm
NRCS

23

OM
%

u33

u1500
3

23

m m

Min.
Max.
Mean
SD
Median

0.004
0.955
0.280
0.231
0.211

0.034
0.922
0.492
0.194
0.491

0.002
0.520
0.811
1.890
0.228
1.362
0.133
0.186
0.205
1.380
HYPRES

1.000
14.861
3.082
2.063
2.500

0.051
0.724
0.316
0.083
0.325

0.022
0.725
0.171
0.094
0.153

Min.
Max.
Mean
SD
Median

0.007
0.931
0.383
0.289
0.291

0.040
0.791
0.407
0.212
0.397

0.023
0.899
0.670
1.700
0.210
1.402
0.138
0.170
0.188
1.430
BRAZ

1.000
13.700
2.688
1.944
2.090

0.047
0.583
0.292
0.105
0.298

0.039
0.422
0.170
0.090
0.139

Min.
Max.
Mean
SD
Median

0.000
0.910
0.271
0.225
0.180

0.034
0.710
0.230
0.165
0.190

0.040
0.810
0.499
0.214
0.528

1.000
11.016
2.730
1.532
2.300

0.080
0.548
0.321
0.093
0.333

0.023
0.414
0.229
0.080
0.238

0.720
1.760
1.173
0.204
1.176

† Db 5 bulk density; OM 5 organic matter content; u33 5 water retention
at 233 kPa matric potential; u1500 5 water retention at 21500 kPa
matric potential.

close to the same amount. The NRCS data set represents a
substantially wider range of soils in terms of silt and clay
content than the HYPRES data set. Differences in other listed
properties are less noticeable. We converted gravimetric watercontents stored in the NRCS database to volumetric watercontents to remain compatible with most existing PTFs andbetween data sets that we used. Different Db values are stored
in the NRCS database—measured at different states of wetness—that had to be used to convert 233 kPa and 21500 kPa
gravimetric water contents to their respective volumetric water content values. While the NRCS and HYPRES data setsconsist of soils that originated almost exclusively from areas
with temperate (continental or maritime) and Mediterranean
climates, most soils in the BRAZ data set originated fromareas with a tropical climate. A known physical characteristicof many soils in the tropics is the bimodality of their particlesize distribution (MacLean and Yager, 1972; Tomasella et al.,2003). This results in a composition of the solid phase that isrich in sand and clay particles while, in many soils, the silt
fraction is almost completely missing. This characteristic canbe identified in the BRAZ data set in Table 1; this set has, on
average, a substantially larger proportion of clay and lower
proportion of silt contents than the other two data sets. Thisdata set also has a considerably lower average Db than theother two sets.

The k-Nearest Neighbor Technique
Unlike classic PTFs, the k-NN technique does not use any
predefined mathematical functions to estimate a certain attribute. A “reference” data set—analogous to the development or training data sets used to develop classic PTFs—issearched for samples that are most similar to the target sample, based on selected input attributes. In most classic k-NN
variants, the “distance” measure is calculated as the classical Euclidean distance between the target and the known
instances (Wettschereck et al., 1997). In a simple case with only
two input attributes, e.g., sand and clay content, selection of
the nearest (or most similar) soil(s) can be represented geometrically using Pythagoras’ theorem, as demonstrated by,e.g., Jagtap et al. (2004). The “distance” of each soil from thetarget soil can be calculated as the square root of the sum ofsquared differences in sand and clay content between thetarget soil and each of the soils of the reference data set. Soilsof the reference data set will then be sorted in ascending orderof their distance from the target soil. The estimated value ofthe output attribute is calculated as the weighted average ofthe output attribute of a preselected number of the nearest soils.
Most PTFs, however, use information on more than two
input attributes. For such cases, the generalized form of

di 5

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
J2
j¼1 DaijO[1]provides a sufficient solution, where di is the “distance” of theith soil from the target soil, Daij represents the differenceof the ith soil from the target soil in the jth soil attribute, andJ is the total number of soil attributes considered as inputs.
The term distance does not refer to actual (physical) distance, but to a measure of similarity; the distance will besmaller for soils that are more similar to the target soil in their
input attributes.
A unit difference in one attribute may, however, not be as
influential as the same unit difference in another attribute. For
example, sand content, if given as a percentage, can take values anywhere between 0 and 100, whereas Db content ranges
theoretically from 0 to a maximum of 2.65 g cm23 in soils. In

1225

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

real data sets, the range of such values is usually narrower. A
unit difference in Db is expected to be more influential than thesame unit difference in sand content. To avoid bias toward one
attribute or the other, the data need to be normalized before
they are used to calculate “distance” using Eq. [1].
Nemes et al. (2006) suggested the following normalization.
First, all input attributes are transformed to obtain temporary
variables with distribution having zero mean and a standard
deviation of 1:

the distribution of distances of the k nearest neighbors fromthe target object as follows:aij(temp) 5 [(aij ) 2 aj ] / s(aj )where k is the number of neighbors considered, di is the
distance of the ith selected neighbor calculated using Eq. [1],and p is a power term that was optimized to provide the best
estimation results. The p term accounts for the weight/distance
relationship and was suggested to be set at p 5 0.767N 0.049,
where N is the number of samples in the reference data set.
Nemes et al. (2006) used different (hierarchical) combinations of the following input attributes: USDA sand, silt, andclay content (SSC), Db, and OM content. They assumed thatthese attributes are all equally relevant and important in theestimation of the output attributes. Four different sets of input
attributes were used to estimate u33 and u1500 from data of
the NRCS data set. The simplest model used only SSC as
predictors. In the following two models, either Db, or OM content was added to SSC as a predictor (SSCBD and SSCOM,
respectively). In the fourth model, all of these inputs were usedas predictors (SSCBDOM). In this study, we implemented thisapproach to avoid a possible bias while applying one particular
set of input attributes, and to account for different levels ofdata availability for potential future users.

[2]

where aij(temp) represents the temporary value of the jth
attribute of the ith soil, and aj and s(aj) represent the meanand standard deviation of the observed values of the jth
attribute in the reference data set. The difference between theminimum and maximum of those temporary variables is thenexamined, and the attribute that shows the widest range oftransformed (temporary) values is identified. The ratio of thewidest range of transformed values and the range of transformed values of each attribute were used as a scaling factor to
obtain zero mean and the same minimum–maximum range inthe data of all attributes:

aij(trans) 5 aij(temp) (Max{range[aj51(temp) ],…,
range[aj5x(temp) ]}) / range[aj(temp) ][3]

where aj(temp) represents the temporary data of the jth soil
attributes normalized using Eq. [2], and aij(trans) represents the
final transformed values of the jth attribute of the ith soil thatare to be used as input.
Nemes et al. (2006) experimented with the number of soils
used to obtain the estimate of the output attribute of the target
soil (k). They suggested using k 5 0.655N 0.493 to calculate theoptimal value of k based on the size of the reference data (N)set. Alternatives exist to this method. Lall and Sharma (1996),
for instance, suggested k 5 N1/2 for N . 100, based on theirexperience under certain conditions, where N is the length ofthe observed sample record, i.e., the number of known
instances in the reference data set. Both studies note, however,
the relatively low sensitivity of the technique to the choice of k.
In this study, we applied the formula suggested by Nemes
et al. (2006).
The user of the k-NN technique also has to decide how toweight each selected soil while calculating the estimate of theoutput attribute. As one solution, the simple average of their
output attribute can be calculated. In such case, the weight[w(i)] of each selected soil will equal 1/k, where k is as defined
above. The calculated “distance” of each soil from the target
object (see Eq. [1]) will be different, however, and it can beargued that a soil closer to the target object should have moreweight in calculating the estimated value than a soil that isfurther from it. Weighting methods that allow either rank- or
distance-dependent weighting of soils offer alternative solutions. One solution mentioned in the literature is that of Lalland Sharma (1996), who calculated weights for each selectedneighbor as

w(i) 5

1/i

O 1/i
k
i51

[4]

where w(i) is the weight associated with the ith nearest neighbor and k is the number of neighbors considered. This method
considers the rank of each sample in being the nearest neighbor to the target object, and does not consider the relative
distances of the selected k neighbors from the target object.
Nemes et al. (2006) used a weighting method that accounts forwi 5 di(rel) /O dki51 i(rel)[5]where k is the number of nearest neighbors retrieved, wi is theassigned weight, and di(rel) is the relative distance of the ithnearest neighbor, calculated asdi(rel) 5O d /d
k
i51 i

p

i

[6]

The Artificial Neural Network Technique
Recently, artificial NNet models have been used successfully in PTF development (e.g., Pachepsky et al., 1996; Tamari
et al., 1996; Schaap et al., 1998; Koekkoek and Booltink, 1999;
Minasny et al., 1999; Schaap and Leij, 2000; Minasny and
McBratney, 2002; Nemes et al., 2003). Most studies found thatthe predictive capabilities of NNet PTFs were equivalent or
superior to different regression-type PTFs. For this reason, we
chose the NNet technique to serve as the basis of comparison
for the k-NN technique in Case Study 1.
A NNet model consists of many simple computing elements(termed neurons or nodes), that are organized into subgroups(layers) and are interconnected as a network by weights. A
model typically consists of an input layer, an output layer, andone (or more) “hidden” layer(s) that connect(s) the input andoutput layers. The number of nodes in the input and output
layers correspond to the number of input and output variablesof the model; the number of hidden nodes can be varied freely.
Data flow goes from the input layer through the hidden
layer(s) to the output layer. A node in the hidden and output
layers receives multiple inputs—typically from all nodes of theprevious layer. Within the node, each input is weighted andcombined to produce a single value as the output of that node,
which is then directed to all the nodes of the next layer, or
outputted if it was a node of the output layer. The weightmatrices are obtained through a calibration (training) procedure, which can then be used to make estimations for independent data. For a more thorough description on NNets, see
Hecht-Nielsen (1990) or Haykin (1994).
Following Nemes et al. (2006), we used a three-layer backpropagation NNet model. There are different approaches to
set the number of nodes in the hidden layer. We elected to
calculate it as half of the total number of input and output
variables, rounded up to the nearest integer. Four different
models, each using a different set of input attributes, weredeveloped to estimate u33 and u1500 separately from data ofReproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1226

VADOSE ZONE J., VOL. 5, NOVEMBER 2006

the NRCS data set. To allow direct comparison with theperformance of the k-NN models, inputs to the four NNetmodels were the same as inputs to the four k-NN models outlined above. We transformed all data, before being presented
to the NNets, to take up the interval [0,1].
The NNets were combined with the data selection procedure of the bootstrap method (Efron and Tibshirani, 1993)
to generate internal calibration–validation data set pairs for an
early stopping procedure. We generated 10 bootstrap replica
data sets, each of which was used to calibrate the NNet models.
This procedure provided 10 subestimates that could be slightly
different from each other. The estimate of a PTF from one
particular ensemble data set—for each single value—was thencalculated by averaging the 10 subestimates of the value. Application of the bootstrap method took place internally in the
NNet program to derive the best estimates from each ensemble’s development data set, and was performed independently
within each of the 200 PTF ensembles described above. All
NNet modeling was performed using the Neural Network
Toolbox in MATLAB (Demuth and Beale, 1992).

Evaluation Criteria
In Case Study 1, k-NN estimations were compared with
NNet estimations. Other case studies provided comparisons
made only between k-NN models using different settings. In allcases, the goal was to report on or to minimize the estimation
errors for the test data set(s) at the population level, which wascharacterized by two measures. Root mean squared residual(RMSR) of the estimations is defined assﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

O (u 2 û )N

RMSR 5

(1/N)

i

i51

2

[7]

i

and the mean residual value (MR) is calculated as

O (u 2 û )N

MR 5 (1/N)

i51

i

[8]

i

In Eq. [7] and [8], N is the number of samples in the test dataset, u and û are measured and estimated water contents, respectively. The MR can quantify systematic errors between
measurements and estimations and the RMSR can give theaccuracy of the estimations in terms of standard deviations.
For Case Study 4, however, we do not directly report on MR
or RMSR values, as the actual values have no particular importance. Rather, we examined the correlation between theestimation error (residual) for the individual test soils and the
distance of the retrieved neighbors from the target object,using R2, the coefficient of determination, calculated as

R2 5

{N

O

O
O O
Ox) ]}{N O[ y 2 (Oy) ]}[N( xy) 2 ( x)( y)]2
[x2 2 (

2

2

2

[9]

where x and y represent the independent and dependent variables in the equation, respectively, and N is the number of samples. The R2 shows what proportion of the dependent variable
can be attributed to the effect(s) of independent variable(s).

set as the reference data set. The rationale behind thiscase study was that PTFs usually produce worse estimations for data that originate from a different geographical area than for the data used to develop or train the
PTF. We compared the degree of such loss in accuracy
by the k-NN technique and an alternative technique.
The k-NN technique was applied using each of fivereference data set sizes (N 5 100, 200, 400, 800, and1600) and four input attribute sets (SSC, SSCBD, SSCOM,
and SSCBDOM) to estimate two output attributes (u33
and u1500). To demonstrate the capabilities of this technique, NNet models were also applied to estimate thesame output attributes, using the same sets and the same
input attributes as for the k-NN technique. This way weshow the pure difference between the two techniques, without being affected by differences originating from the underlying data.

Comparison of Methods to Weight the Retrieved
k Neighbors
The user of the k-NN technique has to make a choice
how to weight the retained neighbors while forming theestimate of the output attribute. We examined threealternative methods: (i) simple averaging of the output
attributes of the retrieved neighbors; (ii) “rank”-based
weighting according to Lall and Sharma (1996) (Eq. [4]);and (iii) distance-based weighting according to Nemes
et al. (2006) (Eq. [5] and [6]). These weighting methods
were implemented in the algorithm separately, eachcombined with each of the reference data set sizes, input
attribute sets and output attributes, as outlined for Case
Study 1. The NRCS data set was used as the reference
data set, and estimations were made for soils of both the
NRCS and HYPRES test data sets.

Sensitivity of the Estimation Accuracy to the
Number of Model Ensembles
When multiple realizations are developed from thesame master data set, usually the number of ensembles
is set arbitrarily. It is rarely examined—and thus remains
a question—whether the outcome of the estimations
would be significantly changed if the developer chose
to use more or fewer ensembles. This may result in suboptimal estimation result—when the number chosen is
too small—or in unnecessarily long computations—when the number chosen is too large. We examined theeffect of the number of ensembles using each of the fivereference data set sizes, four input attribute sets, andtwo output attributes as outlined for Case Study 1. The
NRCS data set was used as the reference data set, andestimations were made for soils in both the NRCS and
HYPRES test data sets. Two hundred ensembles wererun and running RMSRs were recorded after the completion of estimations using each ensemble.

CASE STUDIES
Application of the k-Nearest Neighbor Technique
to Data from an Alternative Data Source

Sensitivity of the k-Nearest Neighbor Technique toData Density

The k-NN technique was applied to estimate u33 andu1500 for the HYPRES data set, using the NRCS data

One of the reported advantages of the k-NN technique is that, unlike parametric PTFs, it uses information

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

that is specific to the target object. It does so as the“nearest neighbors” are selected in terms of their properties, meaning that they are similar to the target object.
Parametric PTFs do not work this way; one (set of) parametric equation(s) describes the entire data space. It canbe argued, however, that estimation errors may be muchlarger when the selected “nearest” neighbors are still
distant from the target object, i.e., when the data density
in the reference data set is small at some locations in thedata domain, meaning that some of the target objects
are not well represented. We hypothesized that larger
distances calculated according to Eq. [1] will lead tolarger estimation errors. To test this, we correlated estimation errors with the calculated distances between the
retrieved nearest neighbors and the target objects.

Sensitivity of the k-Nearest Neighbor Technique tothe Potential Presence of Outliers
Input data are normalized before being presented tothe calculation algorithm. This is done to assure that thedifferent input attributes will receive equal weight in the
distance calculations (Eq. [1]). It may happen, however,that the k-NN technique there are just a few outliers in
the unique data set, which may expand the data range of
one or more of the input attributes. When an outlier ispresent, it may not add much to the characterization of
the whole data set, but may mask the “true” effective
data range in the data set. The significance of this phenomenon in the proposed k-NN variant is that while
performing the data normalization (i.e., Eq. [2] and [3]),an input attribute with outlier(s) may eventually get
lower weight in the distance calculations than it would
based on its effective data range. We examined the effectof the presence of outliers in the reference and test data
sets by simulating the presence of outliers. In this casestudy, we used only the NRCS data set that was also
used in the preceding case studies. Specifically for thiscase study, we imposed further limitations on this data
set. We left out the samples with their sand content ,0.1
or .0.8 kg kg21, and samples with OM content ,2% or.14%. As a result, we obtained a set of 650 samples. We
then put only two of those discarded samples back in thedata set in a controlled manner. We recorded whether
those samples appeared (i) only in the reference dataset, (ii) randomly in the reference or test data sets, or(iii) only in the test data set. We also ran a fourth option,in which neither appeared in any of the data sets. To stayconsistent with the established reference data set size
options, we elected to use 400 samples as a reference
data set in each of the ensembles. All other samples
were used as test data.

Unequal Input Attribute Weighting
It can be argued that equal weights assigned to each
input attribute—as introduced by applying Eq. [2] and[3]—may not provide the best possible results, and thatparticular attribute(s) should receive more weight. As asimple example, we refer to the known dominance ofclay content over that of sand content in the characterization of the dry end of the water retention curve, e.g.,

1227

u1500 in our study. This is because clay particles have asignificantly larger role in determining the presence anddistribution of finer pores in the soil, which control soil
hydraulic properties in the dry range. To test the performance and sensitivity of the k-NN technique to unequal attribute weighting, we introduced additional
scaling to each of the input attributes after normalizing them first according to Eq. [2] and [3]. This concept
has the same logic as the influence weight concept ofMehrotra and Sharma (2006). By applying all combinations of weight factors of 1, 5, 10, 30, 100, and 200, we
introduced scaling to cover ratios of 1, 2, 3, 3.33, 5, 6,
6.66, 10, 20, 30, 40, 100, and 200 to 1. All combinations of
all the above ratios were applied to all input attributecombinations in the algorithm that used all listed input
attributes (SSCBDOM) and the reference data set sizeof 1600 samples. We used the NRCS test data set to test
the performance of each of the weight combinations.
Root mean squared residual values were ordered andcombinations of attribute weights that yielded thesmallest RMSRs were logged. Models were comparedwith each other and with the base model with 1:1 attribute weights.

Estimations using Locally Specific Data
An advantage of this nonparametric technique is that
it is capable of easily adopting new, locally specific data
into a general reference data set that had been previously established. Should new data become available,
the user is able to include those in the reference data set
without the need to redevelop or republish any equations or calculation matrices. A user will presumably beable to improve estimations for specific local samples by
incorporating existing local information in the reference
data set, without affecting estimations for other sectionsof the data space. We tested this potential feature of thek-NN technique by introducing to the study an additional data set that originated from Brazil. This datacollection originated primarily from areas with a tropicalclimate, as opposed to our previously used two data sets,
NRCS and HYPRES.
It was hypothesized that the addition of “locallyspecific” samples would only have noticeable impact on
estimations made to other samples from this location
without causing significant alteration or degradation in
the performance for other samples. This was assumed
because the addition of locally specific data changes or
improves data density locally in that specific part of thedata space. Using the k-NN technique means that weselected samples from the close neighborhood of thetarget object—in terms of its properties—from the reference data set. For soils with other textures, the originalreference data set (without locally specific data) would
supposedly provide the same estimations as before, as
the locally specific samples would not be selected as thenearest neighbors because of the larger differences intheir properties. We used the default NRCS data set toprovide samples for the reference data set, similar to theother case studies. We used 1600 samples. As an option,
additional data from the BRAZ data set was also used asReproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1228

VADOSE ZONE J., VOL. 5, NOVEMBER 2006

input. Figure 1 shows the particle-size distribution of thesoils in all three data sets. We divided the BRAZ data set
into two parts. The section above the divider lines(labeled [OUT]) contains soils the texture of which arepractically not found in the NRCS and HYPRES data
sets. This is despite the fact that simple summary statistics did not reveal that (c.f., Table 1). The divider lines
are defined by the CLAY(kg kg21) . 0.7–0.75 3
SAND(kg kg21) and CLAY(kg kg21) . 0.1 equations.
To account for the presence of locally specific data in thereference data set, the BRAZ[OUT] data were randomly split, and part of the data (108 in each case) were
optionally added to the original reference data sets. The
rest of the BRAZ[OUT] data set (N 5 160) was used fortesting, along with the NRCS test set, and HYPRES andBRAZ[IN] (N 5 160) data sets. Ensemble estimationsof both u33 and u1500 water contents were made using areference data set of 1600 NRCS soils with the 108
BRAZ[OUT] soils optionally added. We used two different sets of input attributes: (i) textural properties only
(SSC); and (ii) textural properties plus Db and OM
content (SSCBDOM).

RESULTS AND DISCUSSION
Application of the k-Nearest Neighbor Technique
to Data from an Alternative Data Source
Using the design parameters k and p, calculated after
Nemes et al. (2006), we applied the k-NN technique,using all five data set sizes and four input attribute sets,
and made estimations for both output attributes for the
HYPRES data set. We also performed the same estimations on identical data using NNet models. Results, interms of RMSR, are summarized in Table 2. Root meansquared residual values are shown separately for each
output attribute, estimation technique, input attributeset, and each development or reference data set size.
Trends that can be observed in Table 2 correspond to
those reported by Nemes et al. (2006): more accurate
estimations of u1500 than of u33; a slight improvementin the estimations when more input attributes were used;
and the mostly insignificant loss of estimation accuracy
with a smaller number of samples used in the reference
data set. These apply to both techniques. Root meansquared residuals for the HYPRES data set are, in general, 0.01 to 0.016 m3 m23 worse than those reported by
Nemes et al. (2006) for the NRCS data set. This is ob-

served using both techniques, and is consistent with thefindings of, e.g., Schaap and Leij (1998). Estimations
made for a data set having different data characteristics
than the development data set is expected (and wasshown) to be worse than for a data set having the samecharacteristics. This was also the case in this study, independent of which estimation technique was used. Whenthe two techniques are compared pairwise, the NNetmodel resulted in 0.001 to 0.005 m3m23 smaller average
RMSRs; such differences are only 0.001 m3 m23 larger
than those reported by Nemes et al. (2006). Lesser estimation accuracy for the HYPRES data set does notseem to be due to lesser capabilities and suboptimal
settings of the k-NN technique because the NNet model
lost accuracy comparably.

Comparison of Methods to Weight the Retrieved
k Neighbors
We compared three weighting methods that are different in their fundamentals. Simple averaging naturallymeans that for a given estimation, any of the selected kneighbors will carry equal weights. The method of Lalland Sharma (1996) assigns weight to each of the selected kneighbors based on their rank in similarity to the target
object in their input attributes. The method of Nemes et al.
(2006) accounted for the distribution of their actual
degree of similarity to the target object. Results have been
averaged by the input and output attributes, and arepresented for each weighting method and reference data
set size in Fig. 2. The differences between RMSRs at the
population level using the three weighting methods are
statistically insignificant within each test data set andreference data set size. This is somewhat surprising, because it suggests that no differentiation is really necessary
among the k number of selected samples, regardless of
their actual resemblance to the target object. This may still
be understandable, given their closeness; the selected ksamples all have significant relevance to the target object.
There are special cases, however, when the target sample falls close to the edge of the data domain in one ormore properties, or the combination of its properties may
simply not be well represented in the reference data set. In
these cases, some of the selected k neighbors may fall a lot
farther from the target sample in their properties thanother selected samples. Neighbors that are significantly
farther than others should theoretically get much lessFig. 1. Distribution of samples in the NRCS, HYPRES, and BRAZ data sets according to the NRCS textural triangle.

1229

www.vadosezonejournal.org

Table 2. Root mean squared residuals for the HYPRES data set using the k-nearest neighbor technique with settings according to Nemes
et al. (2006), and the neural network models.
Sample size of the pedotransfer function development data set (N)Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

N 5 1600

N 5 800

N 5 400

N 5 200

N 5 100

Estimated attribute

Estimation method

Input attributes†

Mean

SD

Mean

SD

Mean

SD

Mean

SD

Mean

SD

Water retention at 233 kPa

Nearest neighbor

SSC
SSC, Db
SSC, OM
SSC, Db, OM
SSC
SSC, Db
SSC, OM
SSC, Db, OM
SSC
SSC, Db
SSC, OM
SSC, Db, OM
SSC
SSC, Db
SSC, OM
SSC, Db, OM

0.069
0.067
0.066
0.065
0.067
0.065
0.063
0.063
0.051
0.051
0.051
0.051
0.050
0.050
0.048
0.049

,0.001,0.001,0.001,0.001,0.001
0.001
0.001
0.001
,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001

0.069
0.067
0.066
0.065
0.067
0.066
0.064
0.063
0.051
0.051
0.051
0.051
0.049
0.050
0.048
0.048

0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001

0.070
0.068
0.066
0.066
0.067
0.065
0.065
0.064
0.052
0.052
0.052
0.052
0.049
0.050
0.048
0.048

0.001
0.001
0.001
0.001
0.001
0.001
0.002
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001

0.071
0.068
0.067
0.067
0.067
0.067
0.065
0.065
0.052
0.053
0.052
0.052
0.049
0.050
0.048
0.049

0.002
0.002
0.002
0.002
0.001
0.002
0.002
0.002
0.001
0.001
0.001
0.001
0.001
0.001
0.002
0.002

0.072
0.070
0.069
0.068
0.068
0.067
0.067
0.066
0.053
0.054
0.054
0.054
0.050
0.050
0.048
0.049

0.003
0.003
0.002
0.002
0.002
0.003
0.004
0.003
0.002
0.002
0.002
0.002
0.001
0.002
0.003
0.003

Neural network

Water retention at 215000 kPA

Nearest neighbor

Neural network† SSC 5 sand, silt, and clay content; Db 5 bulk density; OM 5 organic matter content.

weight when calculating the final output. Differences inthe estimations for such special cases may be the reason
why the weighting method of Nemes et al. (2006) is, by anarrow and insignificant margin, still the most accurate
method of the three.

It was still desirable, however, to establish a minimum
number of ensembles to be used to obtain stable RMSR,
i.e., quasi-flat lines in Fig. 3. The case shown in Fig. 3

Sensitivity of the Estimation Accuracy to the
Number of Model Ensembles
We plotted the running RMSR values against the total
number of ensembles after each replication data set had
been applied to make estimations. We show a representative example of such plots in Fig. 3, for the SSCBDOM
model and u1500 as output. In Fig. 3, the average RMSR
obtained after using M ensembles is shown, meaningthat the last record on each line (at M 5 200) in Fig. 3bmatches the mean values shown in Table 2, line 12. For
any M and M 1 1 ensemble number pairs, differencesamong their performance were too small for the t-test
to show statistically significant differences at p 5 0.95.
When M is small, e.g., M , 20, changes in average
RMSR are visible in Fig. 3; however, due to large standard deviations, changes are not statistically significant.

Fig. 2. Estimation accuracy using three different weighting methods tomake estimations for the NRCS (below) and HYPRES (above) test
data sets.

Fig. 3. Running root mean squared residuals (RMSR) for the (a)NRCS and (b) HYPRES test data sets for up to 200 ensembles
using sand, silt, clay, bulk density, and organic matter content as
input and water retention at 21500 kPa matric potential as output.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1230

VADOSE ZONE J., VOL. 5, NOVEMBER 2006shows, for both test data sets, that the amplitude ofchanges is relatively larger using the first 30 (NRCS) to50 (HYPRES) ensembles, and that there is practicallyno change after those cutoff points in any of the curves.
For larger reference data set sizes, this flat tail of thegraphs emerges after a smaller number of ensembles.
Such behavior could be expected, since with larger subdata sets picked from the same master data set, theoverlap between subdata sets will be more expressed,
meaning that a smaller difference in estimations is expected, leading to less amplitude caused by an individual
ensemble member. We have applied various criteria todetermine the cutoff point at which we consider the estimations unchanged by the addition of further ensemble members. Table 3 lists cutoff ensemble numbers thatwere found after setting an absolute and a relative criterion to consider the estimations unchanged. In relative
terms, we were seeking the ensemble number after
which the change in RMSR did not exceed 0.01% of its
value. Depending on the particular test data set andinput–output combination used, it took the use of 21 to
55 ensembles to reach stability in this respect. In absolute terms, we were searching for the ensemble number
from which the RMSR would remain within 0.001 m3 m23of the RMSR obtained using 200 ensemble members.
This was achieved by using as few as two to seven ensemble members in many cases, and the maximum numberof ensembles to meet this criterion was 26. Each valuein Table 3 is the largest out of five values, representing
the five different reference data set sizes. Typically, thepresented values come from the models that used 100 or
200 samples in the reference data set because less overlap between data sets leads to larger variability in theestimations. Overall, it seems to be safe to state that, forthis application, approximately 50 ensemble members areenough to minimize estimation errors.

Sensitivity of the k-Nearest Neighbor Technique toData Density
We hypothesized that larger distances calculated according to Eq. [1] would lead to larger estimation errors.
We correlated the absolute value of the estimation er-

rors with the calculated distances of the nearest soils
neighbors found in the reference data set to each individual target object. We used the absolute value of
estimation errors in the analysis to make sure positive
and negative errors did not cancel each other. Figure 4
shows the summary of our findings, where R2 values
have been averaged across reference data set sizes andinput attribute sets used.
Correlation between the absolute value of errors andthe distance of the kth neighbors is generally weak; themaximum value of R2 was just above 0.1 when u1500 was
estimated for the NRCS test data set. Such correlations were consistently smaller for the HYPRES dataset, regardless of which output attribute was estimated(maximum R2 5 0.056). Another general observation isthat it was the first nearest neighbor’s distance thatmatters least for the magnitude of the estimation errors,
regardless of which data set or output was used. This issuggested by the smallest R2 found for k 5 1 on eachcurve in Fig. 4. The importance of the distance of theother selected soils (k 5 2 | 25) did not vary much,except for u1500 of the NRCS set. There are substantial
differences between R2 values for the different kth
neighbors, and the clustering of the points is noticeable.
The reason for such clustering is demonstrated inFig. 5. As was mentioned above, points obtained forFig. 4 are a result of averaging the outcomes of using
different reference data set sizes. Figure 5 shows theresults separately for the five reference data set sizes
that yielded the “NRCS 1500” points in Fig. 4. As wefollowed the recommendation of Nemes et al. (2006) tooptimize the number of selected soils (k), k is a functionof the size of the reference data set that was used. Using the suggested equation—k 5 0.655N 0.493, as above—
we obtained values of 6, 9, 13, 18, and 25 for k for reference
data set sizes of 100, 200, 400, 800, and 1600, respectively,
which is reflected in Fig. 5. For this test data set (NRCS)
and output variable (u1500), we found larger variation in
R2 values obtained using different reference data set sizes
than for other cases in Fig. 4. The R2 values are more
evenly distributed among neighbors for larger reference

Table 3. Number of ensembles needed for each input–output combination to reach a relative stability of ,0.01% of the actualroot mean squared residuals value if one more ensemble isadded (Columns a) or an absolute stability of ,0.001 m3m23
compared with the accuracy obtained using 200 ensembles(Columns b) for two different data sets. Each entry is the worst
performing of five input data set sizes.
NRCS
u33‡

HYPRES
u1500‡

u33

u1500

Input attributes†

a

b

a

b

a

b

a

b

SSC
SSC, Db
SSC, OM
SSC, Db,OM

36
52
21
30

18
7
2
16

54
36
30
36

14
7
14
7

55
52
39
39

21
17
26
18

25
29
24
27

4
6
6
6

† SSC 5 sand, silt, and clay content; Db 5 bulk density; OM 5 organic
matter content.
‡ u33 5 water retention at 233 kPa matric potential; u1500 5 water
retention at 21500 kPa matric potential.

Fig. 4. Correlation between the absolute values of the estimation errors and the actual distance values (di) of each of the k neighbors.
Different data set sizes are averaged.

1231

www.vadosezonejournal.org

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Sensitivity of the k-Nearest Neighbor Technique tothe Potential Presence of OutliersFig. 5. Correlation between the absolute values of the estimation errors and the actual distance values (di) of each of the k neighbors.
Different data set sizes are expanded for the NRCS data set for
estimating water retention at 21500 kPa matric potential.

data set sizes, where more soils were selected. The largest
correlation was found for the last (sixth) selected neighborof the smallest reference data set (N 5 100), R2 5 0.125.
Figure 6 shows the scatter plot and the regressed line thatyielded this value in Fig. 5; all other data set input–output
combinations yielded weaker correlations. Overall, the kNN technique is not very sensitive to data density in termsof the estimation errors. Distant, but probably evenly distributed (surrounding), neighbors do not result in a muchbiased estimate for the individual soil. One selected neighbor may have a substantially larger value for the output
attribute, but the other may have a smaller value by thesame margin. The two values, when averaged, will yield areasonable estimate and estimation error for the target
soil, despite the large actual distances calculated for theneighbors, and the potentially largely different individual
values in terms of the output attribute.

Fig. 6. Correlation between the absolute values (ABS) of the estimation errors and the distance values (di) of the sixth neighbor.
Estimations are shown for estimating water retention at 21500 kPa
matric potential for the NRCS test data set using 100 samples in thereference data set and all described inputs (sand, silt, clay content,bulk density, and organic matter content).

Figure 7 shows the sand and OM content of the soils
used in this case study. The two outliers were added oromitted in each case according to the four different outlier settings. Figure 8 shows the RMSR values that wereobtained using different input–output combinations.
Outlier Setting 4 can be considered the point of reference in each case, as there were no outliers in any of thedata sets in that setting. For u33, there is practically no
difference in the estimation results between any of theother outlier settings when cases using the same set of
input attributes are compared. For u1500, there is aslight (,0.002 m3 m23) increase in the estimation errors
with the outliers always being in the reference data set
or being allowed to randomly be included in any of thetwo data sets. This difference is insignificant, however,compared with the variation within the ensemble
RMSRs. In this case study, results in Outlier Setting 3—
i.e., when the outliers are in the test data set—were notexpected to differ much from the reference case. This isbecause the outliers represent only two out of 250 soils,
and are thus practically a negligible part of the RMSR
calculations. When they are present in the reference
data set, however—in all cases or only randomly (i.e.,
Settings 1 and 2)—they have an impact on all calculations for all test soils because they have an impact on theratio between different input attributes in the data
standardization (see Eq. [2] and [3]). Our study, however, only tested the potential presence of outliers
whose properties differed from the properties of thereference data set to a reasonable and realistic extent.

Unequal Input Attribute Weighting
Figure 9 shows the best combination of weights tominimize RMSR for the NRCS test data set. For both
output variables, we show the running average weights
to each input attribute, obtained considering the best
1 to 50 weight combinations. We examine the runningaverages of multiple models rather than individual mod-Fig. 7. Sand and organic matter (OM) contents in the limited NRCS
data subset used in Case Study 5. Data points symbolized by open
circles (N 5 2) are considered hypothetical outliers.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1232

VADOSE ZONE J., VOL. 5, NOVEMBER 2006

Fig. 8. Mean root mean squared residual (RMSR) values and theirstandard deviations for estimating water retention at (above)
233 kPa and (below) 21500 kPa matric potential obtained using
different scenarios to include outliers in the data. Outlier settings:
1 5 outliers always in the reference data; 2 5 outliers mixed in
randomly; 3 5 outliers always in the test data; 4 5 outliers not in thedata. Vertical bars represent 61 standard deviation, based on
200 ensembles. SSC 5 sand, silt, and clay content; SSCBD 5 sand,silt, clay, and bulk density; SSCOM 5 sand, silt, clay, and organic
matter contents; SSCBDOM 5 sand, silt, clay, and organic matter
contents and bulk density.

els to counter possible outlying cases and to monitor
trends in the weighting of input attributes. It may also
happen that the best combination of weights falls somewhere between our picks of scaling factors. In Fig. 9,
weights are presented in a standardized way, so that theminimum weight assigned to any attribute will be equal
to 1. That attribute will carry the least weight in the best
performing models. In Fig. 9a, in the estimation of u33,
sand content was the least important input attribute,having its value at 1, regardless of how many of the first
50 best models we averaged. The weight of silt and OM
content were also very close to 1. Clay content and Db
were two input attributes that had significantly larger
weights assigned in the best models. They both started
around a weight factor of 4, with Db receiving graduallyless weight but clay content assuming gradually more
average weight when more of the best ranked models
were considered. The weights of 3 to 5 that these two
attributes assumed significantly differ from the weight of
1. In the estimation of u1500, sand content was again theleast important input attribute, having its value fixed at
1 throughout the examined spectrum of models. Silt
content started with a weight close to 3, which quickly
decreased to close to 1, as in the case for u33. Bulk
density and OM content carried practically the same
weight throughout the examined spectrum of models.
The weight of these attributes decreased quickly from
3.5 to a steady 2.7 to 2.8. Such quick decrease could also
be observed in the case of clay content, the weight ofwhich became steady at around 7.
In the estimation of both properties, the superior roleof clay content became evident. As expected, it wasmore emphasized in the estimation of u1500 than in theestimation of u33. Bulk density is known to play an

Fig. 9. Running input attribute weights to minimize estimation root
mean squared residuals (RMSR) of (a) water retention at 233 kPa
and (b) 21500 kPa matric potential. The model used sand, silt, andclay content, bulk density (Db), and organic matter (OM) content
as input and 1600 samples in the reference data set.

important role in determining the water-holding capacity of the soil closer to the wet end of the water retention
curve, where soil structure has more influence throughthe formation of macro- and mesopores. Our findings
correspond to this, with Db having a larger impact on theestimation of u33 than on the estimation of u1500. Siltand sand content played a less important role; however,these properties are strongly correlated to clay content,
the input attribute with the most weight.
Table 4 summarizes the results of estimations using unequal weighting of input attributes quantitatively. Oncethe best performing combination of attribute weights
had been found for the NRCS test data set, we appliedthat combination to make estimations for the HYPRES
test data set as well. The gain in estimation accuracy
using unequal attribute weighting was marginal compared with the originally used equal weighting case,regardless of which test data set was used and which
output attribute was estimated. The improvement in1233

www.vadosezonejournal.org

Table 4. Comparison of root mean squared residuals (RMSR)
obtained by the best combination of unequal input attributeweights and by equal input attribute weights for two different
data sets.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

NRCS

RMSR minimum
RMSR with 1:1 weights

addition of BRAZ[OUT] soils to the reference data set.
A very small improvement (0.002–0.003 m3 m23) was
found when estimations were made for samples in the
BRAZ[IN] data subset. The values in absolute terms arecomparable—only slightly worse—to the values for the
HYPRES test data set. This suggests that the k-NN
technique is applicable for data sets that originate fromother geographical areas, but whose data ranges are
represented in the reference data set. The difference
between the accuracy for the HYPRES and BRAZ[IN]data sets is expectedly due to the major difference inclimatic conditions under which the BRAZ soils andsoils of the NRCS and HYPRES data sets developed.
For the BRAZ[OUT] test data set, estimations weresignificantly worse than for any other examined test dataset, using only the NRCS data set as reference data. This
was the case in terms of RMSR as well as MR. An MR
of 0.1 m3 m23 is considered very large. With the inclusion of some BRAZ[OUT] data in the reference dataset, estimation results changed dramatically for this test
data set, as estimation accuracy in terms of RMSR became almost as good as for the NRCS data set (not significantly different at p 5 0.95), and in terms of MR the
improvement is 75 to 80% compared with the original values, depending on how many input attributes were used.
The hypothesis that using local data will improve estimations was not rejected by our experiment with these
data. We showed that, while large improvements can beachieved for “locally specific” data, estimations for test
data originating from other sources could remain practically unchanged. The fact that changes in the estimations
were very limited for the BRAZ[IN] data set confirmsthat the change is not data origin specific, but specific forthe actual soil properties. Test data set samples that fall tothe edge of properties of nonspecific and locally specificsoils or to the edge of the entire data domain still denotespecial cases. Presumably, many of the samples that fallto the edge of the data domain are the ones that show themost bias in the estimations. It is because these test
samples are not uniformly surrounded by neighbors oneach side within the data domain.
When parametric PTFs are used and (additional) locally specific data becomes available, a user either needs
to develop his or her own independent PTF—for which
there may not be enough data—or needs to add the data

HYPRES

u33†

u1500†

u33

u1500

0.04790.04910.03290.03380.0641
0.06490.04970.0504† u33 5 water retention at 233 kPa matric potential; u1500 5 water
retention at 21500 kPa matric potential.

model performance always remained ,0.0012 m3 m23
and did not yield a statistically significant improvement.
This is an indicator of the relative insensitivity of thistechnique to suboptimal attribute weighting. Nemes
et al. (2006) noted similar insensitivity when they examined other perspectives of finding the optimal settings
for the k-NN model. Finding the technique insensitive tounequal weighting to a large extent also helps explain
our findings in Case Study 5. Including outliers in thereference data set in practical terms means that weslightly change the weight of particular input attribute(s)compared with other input attributes. Finding relative
insensitivity to unequal weighting points in the same
direction as our findings in Case Study 5 about thesimulated presence of outliers.

Estimations using Locally Specific Data
Results for all four test data sets are shown in Table 5in terms of RMSR and MR for estimating u33. Results
for the estimation of u1500 were similar and thus are notshown. Underlined data are data that were publishedearlier by Nemes et al. (2006) or in this study (c.f.,Table 2). Estimation accuracy for the NRCS test data set
did not change significantly with the addition of the
BRAZ[OUT] data to the reference data set. Regardless of the number of input attributes used, changes to
RMSR and MR remained well within the variation
among ensemble members. Trends in RMSR and MR
for the HYPRES data set were similar, except thatthe absolute values were somewhat larger than for the
NRCS data set. This was, however, expected based on
the findings of Case Study 1. Estimation accuracy for the
HYPRES data set was not affected significantly by the

Table 5. Summary of results, in terms of root-mean-squared residuals and mean residuals, of the estimation of water retention at 233 kPa
matric potential for the NRCS, HYPRES, and BRAZ data sets. Underlined are numbers of comparison that have been published by
Nemes et al. (2006) (for the NRCS set) or are included elsewhere in this study (for the HYPRES set, c.f., Table 2).
NRCS
Reference data set

Input attributes†

Mean

BRAZ(OUT)
SD

Mean
3

BRAZ(IN)

HYPRES

SD

Mean

SD

Mean

SD

23

RMSR, m m
NRCS
NRCS 1 BRAZ(OUT)
NRCS
NRCS 1 BRAZ(OUT)

SSC
SSC
SSC, Db, OM
SSC, Db, OM

0.054
0.054
0.050
0.051

0.003
0.113
0.003
0.055
0.002
0.112
0.002
0.053
3 23
MR, m m

0.005
0.003
0.004
0.003

0.072
0.070
0.072
0.069

0.001
0.001
0.001
0.001

0.069
0.069
0.065
0.065

,0.001,0.001,0.001,0.001NRCS
NRCS 1 BRAZ(OUT)
NRCS
NRCS 1 BRAZ(OUT)

SSC
SSC
SSC, Db, OM
SSC, Db, OM

0.000
0.001
0.002
0.003

0.003
0.003
0.002
0.002

20.10020.02020.09720.025

0.005
0.005
0.004
0.005

20.01220.00920.00920.005

0.001
0.001
0.002
0.002

0.004
0.004
0.008
0.009

0.001
0.001
0.001
0.001

† SSC 5 sand, silt, and clay content; Db 5 bulk density; OM 5 organic matter content.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

1234

VADOSE ZONE J., VOL. 5, NOVEMBER 2006

and redevelop or readjust the original PTF. In the latter case, however, addition of data with differing datacharacteristics will change the final form of the relationships between inputs and the output. Such relationships—i.e., the equations of parametric PTFs—are valid
globally for the entire range of soils on which they weredeveloped. This way, estimations made for the entire
range of soils are affected by the addition of somespecific soils. This is not the case when the k-NN technique is used, given its design to use samples of thereference data set only from the neighborhood of eachtarget object.

CONCLUSIONS
Nemes et al. (2006) introduced a k-NN techniqueto make estimations of water contents at 233 and21500 kPa matric potentials. They found the performance of the technique to be comparable to that of
NNets. In this study, we further tested the performanceof the k-NN technique under different settings. Estimations were made for different data sets; the technique
was exposed to different weighting methods, and its
performance was monitored for different specific partsof the data space. We also examined the performance of
this technique using different numbers of ensembles.
We found that the k-NN technique is, in general,
insensitive to potentially suboptimal settings in many
aspects. The k-NN technique performed almost equallywell as NNet models developed on the same data to
make estimations for data sets of different origin. The
use of approximately 50 ensemble members resulted inestimation results that were not significantly affected by
the addition of new ensemble members. The k-NN technique showed little sensitivity to the choice of applied
sample weighting methods and to potential suboptimal
settings in terms of input attribute weighting. Differences in data density in parts of the reference data set
did not seem to substantially impact the estimation
error. Substantial improvement was achieved for locallyspecific data when some local samples were included in
the reference data set, while estimations for other samples remained almost unaffected.
The k-NN technique appears to be a competitive
alternative to other techniques to develop PTFs. The
technique shows a large degree of stability and insensitivity to nonoptimal algorithm settings and the use ofdifferent options. Differences introduced by such suboptimal settings and options in the algorithm or in data
weighting caused only marginal changes in estimation
accuracy. Such small differences are unlikely to havenoticeable impact on simulation results that use one
water retention estimate or the other. This technique
can easily adopt new data without the need to redevelop
equations, and can be developed into an “umbrella
PTF” tool to make geographically or climatically specific estimations of soil hydraulic properties.
ACKNOWLEDGMENTS
We thank Dr. Javier Tomasella for making the Brazilian
data collection available for this study.

REFERENCES
Baker, L. 2005. Optimisation of pedotransfer function models forsoil hydraulic properties using an artificial neural network ensemble method. Ph.D. diss. Univ. of Abertay, Dundee, UK.
Dasarathy, B.V. (ed.). 1991. Nearest neighbor (NN) norms: NN pattern classification techniques. IEEE Comput. Soc. Press, Los
Alamitos, CA.
Demuth, H., and M. Beale. 1992. Neural network toolbox manual.
MathWorks, Natick, MA.
Efron, B., and R.J. Tibshirani. 1993. An introduction to the bootstrap.
Monographs on statistics and applied probability. Chapman andHall, New York.
Guber, A.K., Ya.A. Pachepsky, M.Th. van Genuchten, W.J. Rawls, J.
Simunek, D. Jacques, T.J. Nicholson, and R.E. Cady. 2006. Fieldscale water flow simulations using ensembles of pedotransfer functions for soil water retention. Vadose Zone J. 5:234–247.
Harrold, T.I., A. Sharma and S.J. Sheather. 2003a. A nonparametricmodel for stochastic generation of daily rainfall occurrence. Water
Resour. Res. 39:1300. doi:10.1029/2003WR002182.
Harrold, T.I., A. Sharma and S.J. Sheather. 2003b. A nonparametricmodel for stochastic generation of daily rainfall amounts. Water
Resour. Res. 39:1343. doi:10.1029/2003WR002570.
Haykin, S. 1994. Neural networks, a comprehensive foundation. Macmillan College Publ. Co., New York.
Hecht-Nielsen, R. 1990. Neurocomputing. Addison-Wesley, Reading, MA.
Houtemaker, P.L., L. Lefaivre, J. Derome, H. Ritchie, and H.L.
Mitchell. 1996. A system simulation approach to ensemble
prediction. Monthly Weather Rev. 124:1225–1242.
Jagtap, S.S., U. Lall, J.W. Jones, A.J. Gijsman, and J.T. Ritchie. 2004.
Dynamic nearest-neighbor method for estimating soil water parameters. Trans. ASAE 47:1437–1444.
Koekkoek, E.J.W., and H. Booltink. 1999. Neural network models to
predict soil water retention. Eur. J. Soil Sci. 50:489–495.
Lall, U., and A. Sharma. 1996. A nearest-neighbor bootstrap forresampling hydrologic time series. Water Resour. Res. 32:679–693.
MacLean, A.H., and T.U. Yager. 1972. Available water capacities ofZambian soils in relation to pressure plate measurements andparticle size analysis. Soil Sci. 113:23–29.
Mehrotra, R., and A. Sharma. 2006. Conditional resampling of hydrologic time series using multiple predictor variables: A K-nearest
neighbour approach. Adv. Water Resour. 29:987–999.
Minasny, B., and A.B. McBratney. 2002. The neuro-m method forfitting neural network parametric pedotransfer functions. Soil Sci.
Soc. Am. J. 66:352–361.
Minasny, B., A.B. McBratney, and K.L. Bristow. 1999. Comparison ofdifferent approaches to the development of pedotransfer functions
for water-retention curves. Geoderma 93:225–253.
Molteni, F., R. Buizza, T.N. Palmer, and T. Petroliagis. 1996. The
ECMWF ensemble prediction system: Methodology and validation.
Q. J. R. Meteorol. Soc. 122:73–119.
Nemes, A., W.J. Rawls, and Ya.A. Pachepsky. 2006. Use of a nonparametric nearest-neighbor technique to estimate soil water retention. Soil Sci. Soc. Am. J. 70:327–336.
Nemes, A., M.G. Schaap, and J.H.M. Wösten. 2003. Functional
evaluation of pedotransfer functions derived from different scalesof data collection. Soil Sci. Soc. Am. J. 67:1093–1102.
Nemes, A., J.H.M. Wösten, A. Lilly, and J.H. Oude Voshaar. 1999.
Evaluation of different procedures to interpolate the cumulative
particle-size distribution to achieve compatibility within a soil database. Geoderma 90:187–202.
Pachepsky, Ya.A., D. Timlin, and G. Várallyay. 1996. Artificial neuralnetworks to estimate soil water retention from easily measurable
data. Soil Sci. Soc. Am. J. 60:727–773.
Palmer, T.N., A. Alessandri, U. Andersen, P. Cantelaube, M. Davey, P.
Délécluse, M. Déqué, E. Dı́ez, F.J. Doblas-Reyes, H. Feddersen, R.
Graham, S. Gualdi et al. 2004. Development of a European multimodel ensemble system for seasonal to inter-annual prediction(DEMETER). Bull. Am. Meteorol. Soc. 85:853–872.
Schaap, M.G., and F.J. Leij. 1998. Database-related accuracy anduncertainty of pedotransfer functions. Soil Sci. 163:765–779.
Schaap, M.G., and F.J. Leij. 2000. Improved prediction of unsaturated
hydraulic conductivity with the Mualem–van Genuchten model.
Soil Sci. Soc. Am. J. 64:843–851.
Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 1998. NeuralReproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

network analysis for hierarchical prediction of soil hydraulic properties. Soil Sci. Soc. Am. J. 62:847–855.
Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 1999. A bootstrap–
neural network approach to predict soil hydraulic parameters.
p. 1237–1250. In M.Th. van Genuchten et al. (ed.) Proc. Int.
Worksh., Characterization and Measurements of the Hydraulic
Properties of Unsaturated Porous Media, Riverside, CA. 22–24
Oct. 1997. Univ. of California, Riverside.
Sharma, A., and R. O’Neill. 2002. A nonparametric approach forrepresenting interannual dependence in monthly streamflow sequences. Water Resour. Res. 38:5. doi:10.1029/2001WR000953.
Soil Survey Staff. 1951. Soil survey manual. U.S. Dep. Agric. Handbook no. 18. U.S. Gov. Print. Office, Washington, DC.
Soil Survey Staff. 1997. National characterization data. Soil Survey
Lab., Natl. Soil Survey Center, Lincoln, NE.
Tamari, S., J.H.M. Wösten, and J.C. Ruiz-Suárez. 1996. Testing an
artificial neural network for predicting soil hydraulic conductivity.
Soil Sci. Soc. Am. J. 60:771–774.
Tarboton, D.G., A. Sharma, and U. Lall. 1998. Disaggregation procedures for stochastic hydrology based on nonparametric density
estimation. Water Resour. Res. 34:107–119.

1235

Tomasella, J., M.G. Hodnett, and L. Rossato. 2000. Pedotransfer
functions for the estimation of soil water retention in Brazilian soils.
Soil Sci. Soc. Am. J. 64:327–338.
Tomasella, J., Ya.A. Pachepsky, S. Crestana, and W.J. Rawls. 2003.
Comparison of two techniques to develop pedotransfer functions
for water retention. Soil Sci. Soc. Am. J. 67:1085–1092.
Wettschereck, D., D.W. Aha, and T. Mohri. 1997. A review andempirical evaluation of feature weighting methods for a class of lazylearning algorithms. Artif. Intell. Rev. 11(1–5):273–314.
Wösten, J.H.M., A. Lilly, A. Nemes, and C. Le Bas. 1999. Development
and use of a database of hydraulic properties of European soils.
Geoderma 90:169–185.
Yakowitz, S. 1993. Nearest-neighbor estimation for null-recurrent
Markov time series. Stoch. Proc. Appl. 48:311–318.
Yates, D., S. Gangopadhyay, B. Rajagopalan, and K. Strzepek. 2003. A
technique for generating regional climate scenarios using a nearestneighbor algorithm. Water Resour. Res. 39(7):1199. doi:10.1029/
2002WR001769.
Ye, M., S.P. Neuman, and P.D. Meyer. 2004. Maximum likelihood Bayesian averaging of spatial variability models in unsaturated fractured
tuff. Water Resour. Res. 40:W05113. doi:10.1029/2003WR002557.