[R] Confidence intervals nls

[Walmes Zeviani]

The problem was because attr(predict(model, list(DOY=days))) didn't provide
any gradient, it returned NULL. To correct this we need to specify the
gradient by using deriv3() function. Then, when you use predict(model,
list(DOY=days)) will be provided the gradient at this new experimental
points.

# toy data ---------------------------------------------------------------
da <- data.frame(DOY=seq(1,120,l=30))
da$CET <- 9.5+(-6.5*sin(((2*pi)/365)*(da$DOY+65)))+rnorm(da$DOY,0,0.2)
plot(CET~DOY, data=da)

# to get the gradient ----------------------------------------------------
ff <- deriv3(~a+(b*sin(((2*pi)/365)*(DOY+t))),
             c("a","b","t"), function(a, b, t, DOY) NULL)

# model fit --------------------------------------------------------------
model <- nls(CET~ff(a,b,t,DOY), data=da, start=list(a=9.5, b=-6.5, t=65))

# prediction -------------------------------------------------------------
days <- seq(0,365,1)
str(predict(model, list(DOY=days)))
attr(predict(model, list(DOY=days)),"gradient")

se.fit <- sqrt(apply(attr(predict(model, list(DOY=days)),"gradient"), 1,
                     function(x) sum(vcov(model)*outer(x,x))))

# plot -------------------------------------------------------------------
matplot(days, predict(model,list(DOY = days))+
        outer(se.fit, qnorm(c(.5, .025,.975))),
        type="l", col=c(1,2,2), lty=c(1,2,2))

#-------------------------------------------------------------------------


Sincerely.
Walmes Zeviani.

-----
..oooO
..................................................................................................
..(....)... 0ooo...                              Walmes Zeviani
...\..(.....(.....)...     Master in Statistics and Agricultural
Experimentation
....\_)..... )../....       walmeszeviani at hotmail.com, Lavras - MG, Brasil
............
(_/............................................................................................
-- 
View this message in context: http://n4.nabble.com/Confidence-intervals-nls-tp1556487p1556702.html
Sent from the R help mailing list archive at Nabble.com.