Abstract: Spatial variations in soil fertility can obscure treatment effects and hence lead to incorrect fertilizer recommendations. This study was aimed at evaluating oil palm growth response to K application. The response variable in this study was plant growth, expressed as plant height and leaf length. Treatment effects on plant height and leaf length were investigated using Analysis Of Variance (AOV). Both growth variables were assessed for spatial structure using variography. This was followed by Nearest-Neighbor Analysis (NNA) to derive adjusted growth data. The NNA involved a 3-step procedure carried out in an iterative fashion. Treatment effects on the NNA-adjusted growth data were examined using AOV and compared with those obtained using the original growth measurements. Results showed that before removing spatial trends, the effect of treatments on plant growth were not significant. Growth variables exhibited a significant spatial trend. A corresponding observation was found for 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. Following the NNA adjustment, growth variables varied significantly among treatments with the untreated control giving the highest increase in plant growth. The NNA adjustment also rendered improved precision to the linear model, computed using AOV.
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 16-20% higher than that reported by Chan et al. (1993) and Ng (1977). On a whole farm basis, fertilization represents 30-50% 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. Within-block 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
0-4% slope. In contrast, Sri Gunung is characterized by a stronger rolling topography
with 4-12% slope. Both sites had red-yellow 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) Between-row distance = 8.7 m; between-palm 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 2-3 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. Post-treatment growth measurements were made at 6-month 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 RCB-AOV
featured the following model:
| (1) |
where: yij 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 (trend-free). 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 3-step procedure carried out in an iterative fashion. In the first step, non-classical deviations (residuals) from treatment means were computed using the following formula:
| (2) |
where: ri is the calculated residual at point i, gi 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:
| (3) |
where: ci is the estimated covariate at point i; ri-1 and ri+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:
| (4) |
where: di-1 and di+1 are separation distances (m) between the estimation point and both immediate neighbors; di-1 + di+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:
| (5) |
where: yij 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 re-calculated and then assessed for spatial structure. Covariates were re-generated 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 NNA-adjusted 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 30-31%. However, leaf length increment at Sungai Pelepah was about 6-fold 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 pre-treatment (Mar, 2002) and post-treatment (Sept, 2002), ΦAverage concentration reflecting a 1-year period (Mar, 2002-March, 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 Waller-Duncan K-ratio 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 R2 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 (3r*)
was derived using the following regression model (R2 = 0.8): 3r* = 8.1-0.003(x)-0.11(y), where 3r refers to leaf length residuals at iteration 3, x is the distance between palms and y is the distance between rows. The 3r* values were subtracted from the actual final residuals (3r) and then assessed for spatial structure |
|
| Fig. 4: | Semivariograms of leaf length residuals computed iteratively from Sungai Pelepah |
Nearest-neighbor 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 R2, 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 R2, a 98-99% decrease in MSE and a 35-42% 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, R2 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 (3r*)
was derived using the following regression model (R2 = 0.6): 3r* = -17.9-0.05(x)-0.17(y)+816.4(rlP), where 3r refers to plant height residuals at iteration 3, x is the distance between palms, y is the distance between rows, and rlP is the residuals of leaf P. The 3r* values were subtracted from the actual final residuals (3r) 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 R2 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 NNA-adjusted growth variables are given
in Table 4. The corresponding growth response curves are shown
in Fig. 6.
| Table 4: | Treatment effects on NNA-adjusted 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: | NNA-adjusted 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 Waller-Duncan K-ratio t-test; 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 NNA-adjusted growth variables were significant. For the leaf length data at Sungai Pelepah, the R2 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 R2 (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 2-phase 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 R2 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.