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

Answer 1

Try plot.drc in the drc package.

library(drc)

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

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