if (FALSE) {
library(agridat)
data(miguez.biomass)
dat <- miguez.biomass
dat <- subset(dat, doy > 141)
libs(lattice)
xyplot(yield ~ doy | crop*input, data = dat,
main="miguez.biomass",
groups = crop,
type=c('p','smooth'),
auto.key=TRUE)
# ----------
# Archontoulis et al fit some nonlinear models.
# Here is a simple example which does NOT account for crop/input
# Slow, so dont run
if(0){
dat2 <- transform(dat, eu = paste(block, input, crop))
dat2 <- groupedData(yield ~ doy | eu, data = dat2)
fit.lis <- nlsList(yield ~ SSfpl(doy, A, B, xmid, scal),
data = dat2,
control=nls.control(maxiter=100))
print(plot(intervals(fit.lis)))
libs(nlme)
# use all data to get initial values
inits <- getInitial(yield ~ SSfpl(doy, A, B, xmid, scal), data = dat2)
inits
xvals <- 150:325
y1 <- with(as.list(inits), SSfpl(xvals, A, B, xmid, scal))
plot(yield ~ doy, dat2)
lines(xvals,y1)
# must have groupedData object to use augPred
dat2 <- groupedData(yield ~ doy|eu, data=dat2)
plot(dat2)
# without 'random', all effects are included in 'random'
m1 <- nlme(yield ~ SSfpl(doy, A, B, xmid,scale),
data= dat2,
fixed= A + B + xmid + scale ~ 1,
# random = B ~ 1|eu, # to make only B random
random = A + B + xmid + scale ~ 1|eu,
start=inits)
fixef(m1)
summary(m1)
plot(augPred(m1, level=0:1),
main="miguez.biomass - observed/predicted data") # only works with groupedData object
}
}
Run the code above in your browser using DataLab