Crop input response functions with stochastic plateaus.

where equation (19) describes the locus of points at which the crop receives none of the K inputs in excess amounts. Note that, although equation (18) seems to indicate that [x.sub.j] has to be known in advance, imposing constraint equation (19) on the optimization process makes this requirement unnecessary. Furthermore, since equation (19) is an equality constraint, this problem can be converted to an unconstrained equivalent by substituting the constraint (equation 19) into equation (18). Thus, the objective function becomes

(20)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The first-order condition for expected profit maximization can then be obtained by differentiating equation (20) with respect to [x.sub.j]. Using the chain rule and the Liebnitz integral rule (McDonald and Moffitt 1980), the first-order condition is

(21)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Notice that equation (8) is a special case of equation (21) and arises when [h.sub.k](*) is linear and k = {1}.

In general, if [h.sub.k](*) is linear in [x.sub.k], [for all] k, equation (21) will also be linear and could be solved analytically for [x.sub.j, as long as [F.sup.-1](*) can be computed, analytically or numerically. To illustrate this, suppose

(22) [h.sub.k]([x.sub.k]) = [[beta].sub.0k] + [[beta].sub.1k][x.sub.k], [for all] k

where [[beta].sub.0k] and [[beta].sub.1k] are the intercept and slope for the kth factor. It can be shown, by substituting equation (22) into equation (21) and rearranging terms, that

(23) [x.sub.j] = 1/[[beta].sub.1j]{[F.sup.-1](1 - 1/p [summation over (k)] [r.sub.k]/[[beta].sub.1k]) - [[beta].sub.0j]}.

Solution values for each of the other inputs can be obtained by substituting the optimal [x.sub.j] into equation (19):

(24) [x.sub.k] = 1/[[beta].sub.1k]([[beta].sub.0j] + [[beta].sub.1j][x.sub.j] - [[beta].sub.0k]), [for all] k [not equal to] j.

If [h.sub.k](*) is nonlinear in [x.sub.k], equation (21) will be nonlinear and its solution will typically require nonlinear optimization techniques. As an example, suppose [h.sub.k](*)= h(*) is a univariate quadratic function:

(25) h(x) = [[beta].sub.0] + [[beta].sub.1]x + [[beta].sub.2][x.sup.2].

Then equation (21) becomes

(26)

[partial derivative]E([pi] | x)/[partial derivative]x = p([[beta].sub.1] + 2[[beta].sub.2]x)[1 F([[beta].sub.0] + [[beta].sub.1]x + [[beta].sub.2][x.sup.2])] - r = 0

which does not have an explicit analytical solution and can only be solved numerically.

Maximum Likelihood Estimation

This section gives an overview of the procedures used to estimate equation (1). Year random effects associated with the intercept [u.sub.t] and the random error term [[epsilon].sub.it] enter equation (1) linearly but year random effects associated with the plateau, [v.sub.t], enter nonlinearly. We propose maximizing the marginal log-likelihood function directly using the theory of nonlinear mixed-effects models (Wolfinger 1993; Wolfinger 1999; SAS Institute Inc 2004).

Let [f.sub.y]([y.sub.it] | [x.sub.it], [[beta].sub.0], [[beta].sub.1], [[sigma].sup.2.sub.[epsilon]], [v.sub.t], [u.sub.t]), [f.sub.v]([v.sub.t] | [[sigma].sup.2.sub.v]), and [f.sub.u]([u.sub.t]|[[sigma].sup.2.sub.u]) denote the conditional probability density function (pdf) of [y.sub.it], the pdf of v, and pdf of u in equation (1). Let the symbol [theta] represent the vector of unknown parameters, defined as [theta] = [([[beta].sub.0], [[beta].sub.1], [[sigma].sup.2.sub.[epsilon]], [[sigma].sup.2.sub.u], [[sigma].sup.2.sub.v])].sup./]. Then the joint probability density function is

(27)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [v.sub.t] and [u.sub.t] are assumed to have mean zero and variance [[sigma].sup.2.sub.u] and [[sigma].sup.2.sub.u]. The marginal likelihood function of [y.sub.it] is obtained by integrating (27) with respect to [v.sub.t] and [u.sub.t] and taking the product over t and i

(28)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where a and b are the limits of integration for the distribution of [v.sub.t], c and d are the limits of the distribution of [u.sub.t], t is the number of years under consideration, and [I.sub.t] is the number of cross-sectional observations in year t, [f.sub.u](*) is the pdf of u, and [f.sub.y](*) is the pdf of [y.sub.it] conditional on both [v.sub.t] and [u.sub.t].

Theoretically, (28) is the function whose logarithmic transformation is supposed to be maximized in maximum likelihood estimation. However, because [v.sub.t] enters nonlinearly, equation (28) has no closed form and can only be approximated numerically. As with most nonlinear optimization, convergence of equation (28) is not assured. Making sure that convergence occurs at a global maximum rather than a local maximum poses an additional challenge. Several approaches are available for approximating such integrals, including Monte Carlo integration and Gaussian quadrature. Of all the numerical integration techniques, Gaussian quadrature is believed to offer the highest degree of accuracy (Liu 1997; Stiegert and Hertel 1997; Ghomi and Hashemin 1999).

We tried various combinations of starting values, optimization algorithms, and quadrature methods available in the SAS NLMIXED (SAS Institute 2004) procedure. Scaling was used so that the diagonal elements of the Hessian did not differ too much. The largest likelihood value was obtained using first-order approximation starting values, and then non-adaptive Gaussian quadrature. Convergence was slow with nonadaptive Gaussian quadrature, but with adaptive Gaussian quadrature the estimation failed.

Data and Empirical Procedures

Data are from a long-term experiment conducted at the North Central Oklahoma research station located near Lahoma, Oklahoma in the Southern Great Plains of the United States. The study was established in 1970 to investigate winter wheat grain yield response to fertilizer application, using a randomized complete block design (Westerman et al. 1994; Raun et al. 1998). The treatments include a control (no nitrogen) and five levels of nitrogen (20, 40, 60, 80, and 100 pounds per acre). Each treatment was replicated four times. Data from 36 years (1971-2006) were used for estimation. Parameter estimates of the stochastic plateau, nonstochastic plateau, and the switching regression models are estimated using SAS NLMIXED. Then, these estimates are used in MAPLE (2002) to determine the maximum expected profit for the stochastic and nonstochastic plateau models. The analytical solution was verified using a grid search. The expected profits for the switching regression model were estimated using Monte Carlo integration and grid search in SIMETAR (2006) and verified with MAPLE.

The data from this experiment are from plots that are relatively close together and selected for soil uniformity. With other data sets, it might be appropriate to further decompose the random effects to include spatial random effects within a field. Adding such a spatial random effect to the entire equation would not affect the derivation of the optimal level of fertilizer in (14). But the standard deviation of a spatial random effect in the plateau error would need to be added to [[sigma].sub.v] in (14). The SAS NLMIXED procedure used here cannot currently estimate a nonlinear model with multilevel random effects, so a different estimation procedure would need to be used. The possibility of spatial random effects is important in developing nitrogen recommendations from precision sensing. Raun et al. (2002) have developed algorithms that make recommendations for applying the same nitrogen level to the whole field as well as applying a different level to each square meter. If there is spatial variability in the plateau, then the whole field precision-sensing recommendation would have some remaining variability in the plateau that should be considered in determining the optimal amount of nitrogen to apply.

Results

The estimation results for the linear response plateau function, linear response stochastic plateau function, and Maddala and Nelson switching regression models are reported in table 1. All parameters and variance components are significant at the 1% level, except the variance correlation p in the switching regression model. The hypothesis that the plateau is non- stochastic ([H.sub.0] : [[sigma].sup.2.sub.v] = 0) is rejected at the 1% level of significance based on a likelihood ratio test. The likelihood dominance criterion for testing competing nonnested models (Pollack and Wales 1991) indicated that the stochastic plateau model is about 668 times more likely than the switching regression model.

The expected plateau wheat grain yield is about 42, 41.8, and 39.7 bushels per acre for the nonstochastic plateau, stochastic plateau, and switching regression models. The threshold level of nitrogen is 70.6, 57.7, and 38.7 pounds per acre for the switching regression, nonstochastic plateau, and stochastic plateau models. The key difference is that the estimated marginal productivity of nitrogen is higher with the stochastic plateau model and so less nitrogen is needed. Nitrogen productivity is the lowest with the switching regression model and so using it would suggest more nitrogen is needed. There may be some attenuation bias in the nonstochastic plateau and switching regression models that causes the low estimates of nitrogen productivity.

[FIGURE 2 OMITTED]