Hi,
I'm having trouble with the confidence interval of the nls function.
I did my home work, and searched acros the support list until I came up with
following solution of Peter Dalgaard:
example(predict.nls)
se.fit <- sqrt(apply(attr(predict(fm,list(Time =
tt)),"gradient"),1,
function(x) sum(vcov(fm)*outer(x,x))))
matplot(tt, predict(fm,list(Time = tt))+
outer(se.fit,qnorm(c(.5, .025,.975))),type="l")
points(demand ~ Time, data = BOD)
One slight issue is that it doesn't work if "newdata"
is omitted, but then you can easily get the gradient from fm$m$gradient()
I tried this with my own data:
Data <- data.frame(Temp=rep(c(25,40),each=3), Mnd = c(1:3),Degr =
c(0.057,0.077,0.108,0.148,0.198,0.223))
model <- nls(Degr~exp((A/(Temp)+log(Mnd))*B),Data,start=list(A=-10,B=1))
Months <- c(1,2,3,6,9,12,24,48)
se.fit <- sqrt(apply(attr(predict(model,list(Temp =
25,Mnd=Months)),"gradient"),1, function(x) sum(vcov(fm)*outer(x,x))))
But unfortunately I get an error ( Error in apply(attr(predict(model, list(Temp
= 25, Mnd = Months)), "gradient"), :
dim(X) must have a positive length)
When I try using the gradient of the model instead of using the new data then
there is no problem:
se.fit <- sqrt(apply(model$m$gradient(),1, function(x)
sum(vcov(model)*outer(x,x))))
matplot(Data$Mnd, predict(model,list(Temp =
Data$Temp,Mnd=Data$Mnd))+outer(se.fit,qnorm(c(.5,
025,.975))),type="l")
But how about calculating confidence intervals of new data? How do I get an
gradient for these values?
I'm using Windows XP, R 2.4.1.
Thanks
Bart
[[alternative HTML version deleted]]
