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
....\_)..... )../.... walmeszevi...@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.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
