Introduction
Fertilization is the single most expensive operation in oil palm cultivation.
The average proportion of fertilization cost to total maintenance cost per hectare
has in recent years exceeded 70% (Xaviar, 2000), which is 1620% higher than
that reported by Chan et al. (1993) and Ng (1977). On a whole farm basis,
fertilization represents 3050% of the total operating costs.
Oil palm fertilizer trials, from which response curves are derived, are routinely
performed on soils that exhibit significant fertility gradients. Spatial variations
in soil fertility can obscure treatment effects and hence lead to incorrect
fertilizer recommendations. Fertilizer trials are typically laid out as factorial
experiments with a randomized complete block design. The purpose of blocking
is to increase precision by ensuring that treatments are evaluated with respect
to similar environmental conditions within a block. Withinblock heterogeneity
distorts the estimation of treatment effects. According to Brownie et al.
(1993), the greater the heterogeneity within blocks, the greater the variation
in estimates of treatment effect and the poorer the precision of the study.
Two important aspects of experimental design that facilitate accurate interpretation
of trials conducted on spatially variable soils are replication and randomization
(Mulla et al., 1990). Replication of treatments enables estimation of
experimental error, reduces the standard deviation of the treatment mean and
increases the range of conditions under which treatment effects are evaluated
(Steel et al., 1997). Random allocation of treatments ensures that experimental
error estimates are unbiased (Steel et al., 1997) and to some extent
compensates for correlated errors (Mulla et al., 1990). Traditionally,
classical statistical techniques such as analysis of variance, analysis of covariance
and multiple/simple linear regression are used to interpret results from fertilizer
trials and evaluate treatment effects. However, these techniques do not allow
for complete removal of effects induced by soil spatial variability, especially
for designs that feature insufficient replication or improper size and shape
of trial plots/blocks (Mulla et al., 1990).
There are two approaches for dealing with the effects of spatial variability in field experiments. The first approach is to use special experimental designs such as Latin square (Myers, 1966) or nearest neighbor (Freeman, 1979). Both these designs, unfortunately, have elaborate requirements that translate into prohibitively large number of trial plots in situations where many treatments are employed. The second approach is to use alternative statistical methods that facilitate the removal of spatial effects so as to increase the precision of field experiments (Mulla et al., 1990; Samra et al., 1990; Bhatti et al., 1991; Grondona and Cressie, 1991; Hernandez and Mulla, 2002). Essentially, these methods allow the response variable, such as crop yield or plant growth, to be adjusted for spatial variability such that comparisons can be made between the true treatment means. A widely used method for improving the estimation of treatment effects is Nearest Neighbor Analysis (NNA).
This study was aimed at evaluating oil palm growth response following potassium (K) application. To achieve this objective, it was necessary to identify and remove spatial trends in growth response that could not be compensated using classical statistical methods. Specifically, the following issues were looked at: 1) effect of treatments on oil palm growth, 2) spatial variability of growth variables and 3) growth response to fertilizer treatments after removal of spatial effects.
Materials and Methods
Site Attributes
This study was carried out in Sungai Lilin, South Sumatra, Indonesia. Two
sites were selected based on uniformity of stand and palm age. The first site,
Sungai Pelepah Estate, is geographically located at 02° 28’
South and 104° 04’ East while
the second site, Sri Gunung Estate, is 02° 31’
South and 104° East. Sungai Pelepah has a slight rolling topography with
04% slope. In contrast, Sri Gunung is characterized by a stronger rolling topography
with 412% slope. Both sites had redyellow podzolic soils (Typic Paleudult).
The annual rainfall at both sites ranges between 2800 and 3200 mm. The cumulative
yield at Sungai Pelepah was 214 kg Fresh Fruit Bunches (FFB) while that at Sri
Gunung was 204 kg FFB.
Experimental
The trial was laid out using a classical Randomized Complete Block (RCB)
design with the following four treatments: 1) 2.5 kg palm^{1} (standard
plantation rate), 2) 5 kg palm^{1}, 3) 1.25 kg palm^{1} and
4) a control.

Notes: 
1) T1,…., T4  denote treatments, which were randomized
within each replicate 

2) Betweenrow distance = 8.7 m; betweenpalm distance = 9.1 m 
Fig. 1: 
Layout of experimental plot at both study sites 
Treatments were administered based on the plantation’s fertilization
schedule, i.e., two applications per year. Muriate of Potash (MOP) was used
as the nutrient carrier. MOP was broadcast manually within 2 m of the palm circle.
Treatment dates were spaced 23 weeks from the regular fertilization rounds.
Two strips of palms, amounting to 26 trees, were designated as a treatment plot
bearing a dimension of 8x119 m. Each treatment was replicated three times. The
experimental plot layout is given in Fig. 1. Ground maintenance,
such as weed and insect control and pruning of fronds, were carried out in conformity
to plantation policy.
The response variable in this study was plant growth, expressed as plant height and leaf length. Plant height was measured as the vertical distance between ground level and base of the first opened leaf. Leaf length was taken as distance between the rudimentary leaflet and the tip of rachis. Prior to imposing the treatments, growth measurements were made from each experimental plot to establish a baseline data set, which included leaf/soil nutrient levels. Posttreatment growth measurements were made at 6month intervals for 1 year. For each treatment unit, growth measurements were made on six palms. These measurements were averaged to give three observations per treatment unit. Measurement points were spaced 36.4 m between palms and 8.7 m between rows.
Data Analysis
Treatment effects on plant height and leaf length were investigated based
on a classical RCB using Analysis Of Variance (AOV). Essentially, the RCBAOV
featured the following model:
where: y_{ij} is the response variable in block i with treatment j,
μ_{..} is the overall mean, ρ_{i} is the effect of
the ith block, τ_{j} is the effect of the jth treatment, ε_{ij}
are the random errors (classical residuals) that are assumed to be independently
and identically distributed (iid) with mean zero and constant variance, I is
the number of blocks and J is the number of treatments.
The AOV was computed using the Generalized Linear Model (GLM) procedure (SAS Institute, Cary, NC). Growth response curves were constructed using mean values derived from the AOV based on an appropriate function.
Growth variables were assessed for spatial structure using semivariance analysis, computed using GS+ Version 5.1.1 (Gamma Software Design, Plainwell, MI). In constructing the semivariogram, the data were assumed to be stationary (trendfree). Stationarity of the data requires that the semivariance between any two locations in the study region depends only on the distance and direction of separation between the two locations and not their geographic location. The semivariogram was also assumed to be isotropic and omnidirectional, meaning that pairwise squared differences were averaged without regard to direction. The extent to which variability was attributable to spatial structure was defined using the nugget to sill ratio (Cambardella et al., 1994).
Growth variables were next subjected to Nearest Neighbor Analysis (NNA). The NNA involved a 3step procedure carried out in an iterative fashion. In the first step, nonclassical deviations (residuals) from treatment means were computed using the following formula:
where: r_{i} is the calculated residual at point i, g_{i} is the measured height/length at point i and τ_{j} is the treatment mean
In the second step, a unidirectional NNA covariate was computed based on neighboring residuals using the following convention:
where: c_{i} is the estimated covariate at point i; r_{i1} and r_{i+1} represent the residuals calculated from immediate neighbors
Covariates were computed along palm rows across treatments. For border points, which had only one immediate neighbor, the covariate was set equal to the calculated residual. Neighboring points that had unequal separation distances from the estimation point were assigned a weighting factor based on inverse distance using the following convention:
where: d_{i1} and d_{i+1} are separation distances (m) between
the estimation point and both immediate neighbors; d_{i1} + d_{i+1}
is the separation distance (m) between both neighbors.
In the third step, an Analysis Of Covariance (AOC) using the computed NNA covariate was performed to derive an adjusted response (height/length). The AOC had the following model:
where: y_{ij} is the response variable in block i with treatment j, μ_{j} is the estimated true treatment mean, δ_{i} represents the spatial trend component and ε_{ij} represents the effect of random errors
Using the adjusted height/length, residuals were recalculated and then assessed
for spatial structure. Covariates were regenerated and fed into the AOC to
derive a second set of adjusted height/length. The cycle of calculating residuals,
assessing the residuals for spatial structure, computing covariates and AOC
was repeated until the semivariogram exhibited no remaining spatial structure.
However, in situations where the structural variability showed a progressive
reduction but was not completely removed at the final iteration, a corrective
step was imposed. The corrective step, which essentially involved a trend removal
procedure, was performed using a multiple linear regression approach where the
final residual was taken as the response while the predictors comprised 1) distance
between palms (x direction), 2) distance between rows (y direction) and/or 3)
residuals of variable(s) that showed significant correlation with height/length.
The trend factor obtained from the regression model was subtracted from the
residuals generated at the final iteration and then assessed for spatial structure.
Treatment effects on the NNAadjusted height/length were examined using AOV
and compared with results obtained using the original height/length measurements.
Results and Discussion
Descriptive statistics leaf/soil nutrients are given in Table 1. Plant height increments at Sungai Pelepah and Sri Gunung were comparable at 26 and 32 cm, respectively, with CVs of 3031%. However, leaf length increment at Sungai Pelepah was about 6fold lower than that at Sri Gunung. Both sites registered a high CV for leaf length. Comparatively, the concentration of leaf and soil K were higher and less variable at Sungai Pelepah.
Effect of Treatments on Growth
Results of the RCB AOV are given in Table 2. The corresponding
growth response curves are shown in Fig. 2. Plant height and
leaf length differences among treatments were insignificant at both study sites.
Table 1: 
Descriptive statistics of crop yield, growth and leaf/soil
nutrient levels 

^{§}Unit increase (cm) between pretreatment
(Mar, 2002) and posttreatment (Sept, 2002), ^{Φ}Average concentration
reflecting a 1year period (Mar, 2002March, 2003); expressed as % for leaf
nutrients, mg kg^{1} for soil P and m.e. 100 g^{1} for
soil K, Mg, Ca 
Table 2: 
Treatment effects on growth at (a) Sungai Pelepah Estate
and (b) Sri Gunung Estate 

^{**} denotes significance at p=0.05 
(a) Sungai Pelepah Estate 

(b) Sri Gunung Estate 

Note: 
MSD = Minimum Significance Difference based on the WallerDuncan
Kratio t Test; means that are separated by values smaller or equal
to the MSD are not significantly different at p = 0.05 
Fig. 2: 
Growth response to treatments at (a) Sungai Pelepah Estate and (b) Sri
Gunung Estate 
However, plant height differences across replicates at Sri Gunung were significant
at the 95% probability level. At Sungai Pelepah, the linear models explained
58 and 33% of the variation in plant height and leaf length, respectively.
(a) Sungai Pelepah Estate


(b) Sri Gunung Estate


Fig. 3: 
Spatial structural attributes of growth parameters at (a)
Sungai Pelepah Estate and (b) Sri Gunung Estate 
Meanwhile, at Sri Gunung, the linear models explained 75 and 41% of the variation
in plant height and leaf length, respectively. Plant height and leaf length
responses to K at both sites were described using quadratic response functions
with R^{2} values approaching unity in most cases.
The CV for leaf length in all trial plots was relatively higher than that for plant height, indicating that plant height is possibly a more reliable growth indicator. Results clearly show that both plant height and leaf length were not responsive to treatments. It is possible that treatment effects on growth were masked by the spatial trends that have not been accounted for in the RCB AOV model.
Spatial Structure Assessment of Growth Variables
Semivariograms of measured differences in plant height and leaf length at
both study sites are given in Fig. 3. At Sungai Pelepah, leaf
length exhibited a moderate spatial dependence, while plant height showed a
random spatial structure. The spatial structure of leaf length was defined by
a spherical model with a short range at 38 m and a relatively low nugget variance.
At Sri Gunung, plant height exhibited a moderate spatial dependence, while leaf
length showed a random spatial distribution. The spatial structure of plant
height was defined by a spherical model with a short range at 57 m and a moderate
nugget variance.

^{§}Trend factor (^{3}r*)
was derived using the following regression model (R^{2} = 0.8):
^{3}r* = 8.10.003(x)0.11(y), where ^{3}r refers to leaf
length residuals at iteration 3, x is the distance between palms and y
is the distance between rows. The ^{3}r* values were subtracted
from the actual final residuals (^{3}r) and then assessed for
spatial structure 
Fig. 4: 
Semivariograms of leaf length residuals computed iteratively
from Sungai Pelepah 
Nearestneighbor Analysis on Growth Variables
Spatial structure assessment of growth (plant height and leaf length) residuals
was performed iteratively. Essentially, the semivariogram served as an indicator
of spatial effect removal/reduction. Only leaf length residuals at Sungai Pelepah
and plant height residuals at Sri Gunung exhibited a defined spatial structure.
Diagnostic model parameters such as R^{2}, Mean Square Error (MSE) and
Coefficient of Variation (CV) computed from the analysis of covariance (AOC)
for each iterative step are given below (Table 3). Semivariograms
of growth residuals for each iterative step are given by site in Fig.
4 and 5.
At Sungai Pelepah, the three iterations caused a 153% increase in R^{2}, a 9899% decrease in MSE and a 3542% decrease in CV. The structural variance of length residuals decreased by 70% while its range decreased by 36%. At Sri Gunung, which had three iterations, R^{2} increased by 52%, while MSE and CV decreased by 91 and 33%, respectively. The structural variance of height residuals decreased by 33% while the range decreased by 20%. For both sites, the remaining spatial structure on the final iteration was essentially removed by the corrective procedure.

^{§}Trend factor (^{3}r*)
was derived using the following regression model (R^{2} = 0.6):
^{3}r* = 17.90.05(x)0.17(y)+816.4(r_{lP}), where
^{3}r refers to plant height residuals at iteration 3, x is the
distance between palms, y is the distance between rows, and r_{lP}
is the residuals of leaf P. The ^{3}r* values were subtracted
from the actual final residuals (^{3}r) and then assessed for
spatial structure.

Fig. 5: 
Semivariograms of plant height residuals computed iteratively
from Sri Gunung 
Table 3: 
Diagnostic model parameters of growth variables derived from
analysis of variance 

^{§}The original growth data were used for iteration
1 to generate adjusted plant height/leaf length 
Clearly, the progression of iterative cycles produced an increase in R^{2} with an accompanying decrease in MSE and CV. These improvements corresponded with a substantial decrease in structural variance of growth residuals. Results also showed that two to four iteration cycles were sufficient to substantially remove the spatial effects.
Effect of K Treatments on Growth Variables after Adjustment for Spatial
Effects
Results of the RCB AOV performed on NNAadjusted growth variables are given
in Table 4. The corresponding growth response curves are shown
in Fig. 6.
Table 4: 
Treatment effects on NNAadjusted growth at a) Sungai Pelepah
Estate and b) Sri Gunung Estate 

^{**} denotes significance at p = 0.05 
(a) Sungai Pelepah Estate 

(b) Sri Gunung Estate 

Fig. 6: 
NNAadjusted growth response: (a) leaf length at Sungai Pelepah
Estate and (b) plant height at Sri Gunung Estate. 
Notes: 
1. MSD = Minimum Significance Difference based on the WallerDuncan Kratio
ttest; means that are separated by values smaller or equal to the MSD are
not significantly different at p = 0.05, 2. SE = Standard Error 
At both sites, treatment effects on the NNAadjusted growth variables were
significant. For the leaf length data at Sungai Pelepah, the R^{2} increased
from 0.33 to 0.98 while the MSE and CV were reduced by 96 and 43%, respectively.
A similar effect was observed with the plant height data at Sri Gunung, which
had a narrower increase in R^{2} (0.75 to 0.90) and a smaller reduction
in MSE (38%) as well as CV (4.5%).
The leaf length data from Sungai Pelepah exhibited a 2phase response with an initial decreasing trend followed by a gentle upward trend. At Sri Gunung, however, the plant height response clearly demonstrated that lower treatment rates gave higher height increments.
The NNA adjustment resulted in growth variables showing significant differences among treatments. At both study sites, the untreated control gave the highest increase in plant growth. At Sungai Pelepah, leaf K and leaf Mg were negatively correlated (r = 0.70). Based on Foster et al. (1988), leaf K values (mean of 0.96%) indicated mild deficiency while leaf Mg (mean of 0.45%) showed optimum levels. These observations suggest an antagonistic relationship between K and Mg, which could have affected plant growth. Both K and Mg are known to compete for entry into the plant. At Sri Gunung, positive correlations occurred between leaf K and leaf P (r = 0.67) and between plant height and leaf P (r = 0.65). Increased K rates could have imposed a negative impact on height increments by inducing higher demand for P.
Apart from separating the spatial trend from the true treatment effects, this approach also substantially increased the precision of the AOV model in terms of R^{2} and MSE.
Conclusions
Before removing spatial effects using the NNA, K treatments had no significant effect on leaf length and plant height. Both growth variables exhibited a significant spatial trend. A similar observation was found for the growth residuals. The NNA technique was found to substantially reduce structural variance present in the growth data sets, which enabled the assessment of true treatment effects. As a result of the NNA adjustment, leaf length and plant height differences among treatments became significant at Sungai Pelepah and Sri Gunung, respectively. At both study sites, the untreated control gave the highest increase in plant growth. The NNA adjustment also rendered improved precision to the linear model.
Acknowledgment
The authors are grateful to Cargill Inc. (USA) and its plantation subsidiary (PT Hindoli, South Sumatra, Indonesia) for financial support and technical collaboration.