hernandez.nitrogen: Multi-environment trial of corn with nitrogen fertilizer at 5...

Description Format Details Source Examples

Description

Corn response to nitrogen fertilizer at 5 sites.

Format

A data frame with 136 observations on the following 5 variables.

site

site factor, 5 levels

loc

location name

rep

rep, 4 levels

nitro

nitrogen, kg/ha

yield

yield, Mg/ha

Details

Experiment was conducted in 2006 at 5 sites in Minnesota.

Source

Hernandez, J.A. and Mulla, D.J. 2008. Estimating uncertainty of economically optimum fertilizer rates, Agronomy Journal, 100, 1221-1229. https://doi.org/10.2134/agronj2007.0273

Electronic data kindly supplied by Jose Hernandez.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
data(hernandez.nitrogen)
dat <- hernandez.nitrogen
cprice <- 118.1 # $118.1/Mg or $3/bu
nprice <- 0.6615 # $0.66/kg N or $0.30/lb N

# Hernandez optimized yield with a constraint on the ratio of the prices.
# Simpler to just calculate the income and optimize that.
dat <- transform(dat, inc = yield * cprice - nitro * nprice)
require(lattice)
xyplot(inc ~ nitro|site, dat, groups=rep, auto.key=list(columns=4),
       xlab="nitrogen", ylab="income", main="hernandez.nitrogen")

# Site 5 only
dat1 <- subset(dat, site=='S5')

# When we optimize on income, a simple quadratic model works just fine,
# and matches the results of the nls model below.
# Note, 'poly(nitro)' gives weird coefs
lm1 <- lm(inc ~ 1 + nitro + I(nitro^2), data=dat1) 
c1 <- coef(lm1)
-c1[2] / (2*c1[3])
##    nitro
## 191.7198    # Optimum nitrogen is 192 for site 5


## Not run: 
# Use the delta method to get a conf int
require("car")
del1 <- deltaMethod(lm1, "-b1/(2*b2)", parameterNames= paste("b", 0:2, sep=""))
# Simple Wald-type conf int for optimum
del1$Est +  c(-1,1) * del1$SE * qt(1-.1/2, nrow(dat1)-length(coef(lm1)))
## 118.9329 264.5067


# Nonlinear regression
# Reparameterize b0 + b1*x + b2*x^2 using th2 = -b1/2b2 so that th2 is optimum
nls1 <- nls(inc ~ th11- (2*th2*th12)*nitro + th12*nitro^2,
          data = dat1, start = list(th11 = 5, th2 = 150, th12 =-0.1),)
summary(nls1)
# Wald conf int
wald <- function(object, alpha=0.1){
  nobs <- length(resid(object))
  npar <- length(coef(object))
  est <- coef(object)
  stderr <- summary(object)$parameters[,2]
  tval <- qt(1-alpha/2, nobs-npar)
  ci <- cbind(est - tval * stderr, est + tval * stderr)
  colnames(ci) <- paste(round(100*c(alpha/2, 1-alpha/2), 1), "pct", sep= "")
  return(ci)
}
round(wald(nls1),2)
##          5
## th11 936.44 1081.93
## th2  118.93  264.51   # th2 is the optimum
## th12  -0.03   -0.01


# Likelihood conf int
require(MASS)
round(confint(nls1, "th2", level = 0.9),2)
##       5
## 147.96 401.65


# Bootstrap conf int
require(boot)
dat1$fit <- fitted(nls1)
bootfun <- function(rs, i) { # bootstrap the residuals
  dat1$y <- dat1$fit + rs[i]
  coef(nls(y ~ th11- (2*th2*th12)*nitro + th12*nitro^2, dat1, start = coef(nls1)))
}
res1 <- scale(resid(nls1), scale = FALSE) # remove the mean.  Why? It is close to 0.
set.seed(1) # Sometime the bootstrap fails, but this seed works
boot1 <- boot(res1, bootfun, R = 500)
boot.ci(boot1, index = 2, type = c("perc"), conf = 0.9)
## Level     Percentile
## 90


## End(Not run)

agridat documentation built on May 2, 2019, 4:01 p.m.