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
At one point the authors of predict.nls apparently planned to deliver SE's around predictions, but the se.fit argument is currently ignored. stackoverflow.com/questions/12425920/… — IRTFM

G. Grothendieck G. Grothendieck · Accepted Answer · 2015-04-08T13:03:20

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) .

graphing confidence intervals nls r

1 Answers