# \donttest{
require(ggplot2)
require(nlme)
data(Soybean)
SoyF <- subset(Soybean, Variety == "F" & Year == 1988)
fm1 <- nls(weight ~ SSlogis(Time, Asym, xmid, scal), data = SoyF)
## The SSlogis also supplies analytical derivatives
## therefore the predict function returns the gradient too
prd1 <- predict(fm1, newdata = SoyF)
## Gradient
head(attr(prd1, "gradient"))
## Prediction method using gradient
prds <- predict2_nls(fm1, interval = "conf")
SoyFA <- cbind(SoyF, prds)
ggplot(data = SoyFA, aes(x = Time, y = weight)) +
geom_point() +
geom_line(aes(y = Estimate)) +
geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "purple", alpha = 0.3) +
ggtitle("95% Confidence Bands")
## This is equivalent
fm2 <- nls(weight ~ Asym/(1 + exp((xmid - Time)/scal)), data = SoyF,
start = c(Asym = 20, xmid = 56, scal = 8))
## Prediction interval
prdi <- predict2_nls(fm1, interval = "pred")
SoyFA.PI <- cbind(SoyF, prdi)
## Make prediction interval plot
ggplot(data = SoyFA.PI, aes(x = Time, y = weight)) +
geom_point() +
geom_line(aes(y = Estimate)) +
geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "purple", alpha = 0.3) +
ggtitle("95% Prediction Band")
## For these data we should be using gnls instead with an increasing variance
fmg1 <- gnls(weight ~ SSlogis(Time, Asym, xmid, scal),
data = SoyF, weights = varPower())
IC_tab(fm1, fmg1)
prdg <- predict_gnls(fmg1, interval = "pred")
SoyFA.GPI <- cbind(SoyF, prdg)
## These prediction bands are not perfect, but they could be smoothed
## to eliminate the ragged appearance
ggplot(data = SoyFA.GPI, aes(x = Time, y = weight)) +
geom_point() +
geom_line(aes(y = Estimate)) +
geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), fill = "purple", alpha = 0.3) +
ggtitle("95% Prediction Band. NLS model which \n accomodates an increasing variance")
# }
Run the code above in your browser using DataLab