graphing confidence intervals nls r

Question

I'm in the process of putting some incidence data together for a proposal. I know that the data takes on a sigmoid shape overall so I fit it using NLS in R. I was trying to get some confidence intervals to plot as well so I used bootstrapping for the parameters, made three lines and here's where I'm having my problem. The bootstrapped CIs give me three sets of values, but because of equation the lines they are crossing.

Picture of Current Plot with "Ideal" Lines in Black

NLS is not my strong suit so perhaps I'm not going about this the right way. I've used mainly a self start function to this point just to get something down on the plot. The second NLS equation will give the same output, but I've put it down now so that I can alter later if needed.

Here is my code thus far:

data <- readRDS(file = "Incidence.RDS")

inc <- nls(y ~ SSlogis(x, beta1, beta2, beta3), 
           data = data, 
           control = list(maxiter = 100))

b1 <- summary(inc)$coefficients[1,1]
b2 <- summary(inc)$coefficients[2,1]
b3 <- summary(inc)$coefficients[3,1]

inc2 <- nls(y ~ phi1 / (1 + exp(-(x - phi2) / phi3)), 
            data = data, 
            start = list(phi1 = b1, phi2 = b2, phi3 = b3), 
            control = list(maxiter = 100))

inc2.boot <- nlsBoot(inc2, niter = 1000)    

phi1 <- summary(inc2)$coefficients[1,1]
phi2 <- summary(inc2)$coefficients[2,1]
phi3 <- summary(inc2)$coefficients[3,1]

phi1_L <- inc2.boot$bootCI[1,2]
phi2_L <- inc2.boot$bootCI[2,2]
phi3_L <- inc2.boot$bootCI[3,2]

phi1_U <- inc2.boot$bootCI[1,3]
phi2_U <- inc2.boot$bootCI[2,3]
phi3_U <- inc2.boot$bootCI[3,3]

#plot lines

age <- c(20:95)
mean_incidence <- phi1 / (1 + exp(-(age - phi2) / phi3))
lower_incidence <- phi1_L / (1 + exp(-(age - phi2_L) / phi3_L))
upper_incidence <- phi1_U / (1 + exp(-(age - phi2_U) / phi3_U))

inc_line <- data.frame(age, mean_incidence, lower_incidence, upper_incidence)


p <- ggplot()
p <- (p
      + geom_point(data = data, aes(x = x, y = y), color = "darkgreen")
      + geom_line(data = inc_line, 
                  aes(x = age, y = mean_incidence), 
                  color = "blue", 
                  linetype = "solid")

      + geom_line(data = inc_line, 
                  aes(x = age, y = lower_incidence), 
                  color = "blue", 
                  linetype = "dashed")

      + geom_line(data = inc_line, 
                  aes(x = age, y = upper_incidence), 
                  color = "blue", 
                  linetype = "dashed")

      + geom_ribbon(data = inc_line, 
                    aes(x = age, ymin = lower_incidence, ymax = upper_incidence), 
                    fill = "blue", alpha = 0.20)

      + labs(x = "\nAge", y = "Incidence (per 1,000 person years)\n")
      )

print(p)

Here's a link to the data.

Any help on what to do next or if this is even possible given my current set up would be appreciated.

Thanks

There is a difference between confidence intervals of parameters and confidence intervals of predictions. You can't calculate the latter from the former like you tried to do, e.g., because parameters are usually correlated. I suggest that you bootstrap the predictions. — Roland, Apr 07 '15 at 16:54
At one point the authors of `predict.nls` apparently planned to deliver SE's around predictions, but the `se.fit` argument is currently ignored. http://stackoverflow.com/questions/12425920/predict-cannot-display-the-standard-errors-of-the-predictions-with-se-fit-true — IRTFM, Apr 07 '15 at 17:25
Code for Roland's suggestion: http://www.r-bloggers.com/predictnls-part-1-monte-carlo-simulation-confidence-intervals-for-nls-models/ — IRTFM, Apr 07 '15 at 17:31

score 1 · Accepted Answer · answered Apr 08 '15 at 13:03

1

Try plot.drc in the drc package.

library(drc)

fm <- drm(y ~ x, data = data, fct = LL.3())
plot(fm, type = "bars")

screenshot

P.S. Please include the library calls in your questions so that the code is self contained and complete. In the case of the question here: library(ggplot2); library(nlstools) .

answered Apr 08 '15 at 13:03

G. Grothendieck

254,981
17
203
341

Thank for your answer. Have not heard of the drc package but it looks as though this will also be useful going forward. – MJH Apr 08 '15 at 17:09

graphing confidence intervals nls r

1 Answers1