0

How could I plot the fitted model with the confident intervals for each predicted value (in x= 5, 10, 15, 20, 25, 30, 35) ?

For the next dataset

df<-data.frame(
  x = rep(c( 5, 10, 15,  20,  25, 30, 35), each=4 ), 
  y = c(0.2, 1.1, 1.5, 0.9, 
        2.1, 1.9, 2.75, 3.4, 
        5.15, 4.6, 4.75, 4.15, 
        7, 6.7, 6.7, 6.95,  
        7, 5.45, 6.15, 6.4, 
        0.001, 0.001, 0.5, 
        0.001, 0.001, 0.001, 0.001, 0.001)
)
head(df)

And the follow fitted model:

fun <-with(df,  
 y ~ Yopt*((x-Tmin)/(Topt-Tmin))^(b1*(Topt-Tmin)/(Tmax-Topt))*((Tmax-x)/(Tmax-Topt))^b1
    )

 starters <- expand.grid(Yopt = seq(4, 8, len = 4),
                         Tmin = seq(0, 5, len = 4), 
                         Topt = seq(15, 25, len = 4),
                         Tmax= seq(28, 38, len = 4),
                         b1 = seq(0, 4, len = 4))

fit <- nls2(fun, start = starters, algorithm = "brute-force")

summary(fit)

with(df, c(plot(y~x))); points(fitted(fit)~I(df$x), pch=19)

with(as.list(coef(fit)),
     curve(
 Yopt*((x-Tmin)/(Topt-Tmin))^(b1*(Topt-Tmin)/(Tmax-Topt))*((Tmax-x) / (Tmax-Topt)) ^ b1,
           add=TRUE, col="red"))

enter image description here

Juanchi
  • 1,147
  • 2
  • 18
  • 36
  • You should define a proper function instead of using `with` like that. Anyway, I was going to suggest bootstrapping of residuals, but apparently the model fit is not stable. When passing the coefficients found by `nls2` to `nls` it complains. – Roland Dec 13 '14 at 22:30
  • 1
    @Roland - Yes it is unstable, but he insists on using `nls2(...)`: http://stackoverflow.com/questions/26681385/nls-convergence-error. – jlhoward Dec 14 '14 at 01:06
  • I was surprised to see that predict.nls(...)` does *not* support the `se.fit` argument!! – jlhoward Dec 14 '14 at 01:07
  • @jlhoward Confidence intervals for a non-linear fit can only be calculated approximately. For some models extremely approximately. R sometimes provides you with the rope to hang yourself, but in this case it doesn't. – Roland Dec 14 '14 at 12:05

1 Answers1

1

You should take @Roland's comment on your use of with seriously.

There are much better ways of doing this, but building on your code, I would have done something like this.

library(ggplot2)
library(reshape2)
df$fitted <- fitted(fit)
df$upper <- df$fitted + 1 #I didn't bother to produce actual confidence band
df$lower <- df$fitted - 1
df <- melt(data = df, id.vars = c("x","upper","lower"))
coef <- coef(fit)
fit.fun <-function(x) coef[1]*((x-coef[2])/(coef[3]-coef[2]))^(coef[5]*(coef[3]-coef[2])/(coef[4]-coef[3]))*((coef[4]-x) / (coef[4]-coef[3])) ^ coef[5]

ggplot(df, aes(x=x)) + geom_point( aes(y=value, color=variable)) +
  geom_segment(aes(x=x, xend=x, y=upper, yend=lower)) +
  stat_function(fun=fit.fun, color="blue")

enter image description here

Khashaa
  • 7,293
  • 2
  • 21
  • 37
  • I think the hole point of OP's question is how to determine the CL of the fitted values. How would you do that?? – jlhoward Dec 14 '14 at 01:09
  • Double question mark suggested impending attack from expert on the subject, and made me think double. I am not sure what you mean by "CL of the fitted values"; guessing it is "CI of the fitted value" my quick answer is: there must be `predict` of some kind designed for that purpose, if not built-in. `as.lm.nls` gives wrong answer, so I would use the delta method. – Khashaa Dec 14 '14 at 03:26
  • I thought the OP's question was about putting vertical segment on a graph. – Khashaa Dec 14 '14 at 03:27
  • I just used the dataset and the fitted model to play with them and try to find out how to add the vertical CI for every fitted value of x. As @Khashaa explain in his plot, but with the right CL. By this way I would find the Tmin value when the dataset is not complete to find it. I would consider that CI that includes "0" as the Tmin. – Juanchi Dec 14 '14 at 23:21