How can I get confidence intervals for an nls(broken stick) class object

Question

I am trying to figure out how to obtain/plot confidence bounds for nls objects in R.

For example here is a nls model

bstick.lm.mean <- nls(TCTmean ~ cbind("intercept" = 1, 
                                      "l2Flow" = l2Flow, 
                                      "l2FlowBr" = ifelse(l2Flow > Br, 
                                                          l2Flow - Br, 0)),
                      start = list(Br = 6), 
                      algorithm = "plinear", 
                      data = flow.new.sum.dat)

So bstick.lm.mean is a nls class object.

new.seq4 = seq(min(flow.new.sum.dat$l2Flow), max(flow.new.sum.dat$l2Flow), length = 200)
new.seq4 = data.frame(new.seq4)
names(new.seq4) = 'l2Flow'

pz = predict(bstick.lm.mean, newdata = new.seq4, 
             interval = 'confidence', se.fit = TRUE, level = 0.95)
test.frame2 = data.frame(new.seq4,pz)

ggplot(data = test.frame2)+
  geom_point(mapping = aes(x = l2Flow, y = pz),
             shape = 1, col = 'red') +
  geom_point(data = test.frame,
             aes(x = l2Flow, y = TCTmean),
             shape = 0) +
  theme(panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black")) +
  xlab("Log2 Flow (KL)") +
  ylab("Mean Transformed Ct")

Produces the following plot:

My question is, how can I obtain confidence bands for this non linear plot? I am referring to the bands around the regression line.

Answer 1

I extracted data from the scatterplot for analysis, and performed an equation search. I found that several different sigmoidal equations were fitting the data best, and for individual sigmoidal equations it should be simpler to determine the confidence intervals with standard statistical software. As an example, here are my results for a Hyperbolic Logistic sigmoidal equation, "y = (a * pow(x, b)) / (c + pow(x, b))", with fitted parameters a = 1.6177820755100655E+01, b = -1.5270446610701983E+01, and c = 4.2601082365916449E-12 yielding RMSE = 2.58 and R-squared = 0.85. Again, there were several equally "good" sigmoidal equations to choose from.

Answer 2

Estimates of uncertainty around change points (point of break) are notoriously hard to do analytically. Try the R package mcp which takes a computational (Bayesian) approach:

library(mcp)
model = list(
  y ~ 1 + x,  # Slope
  ~ 0  # joined plateau
)

fit = mcp(model, df)
plot(fit, q_fit = TRUE)  # Plot with quantiles

The red lines are the highest-density interval, the grey lines random posterior draws, and the blue curve is the posterior distribution of the change point location. Use plot_pars(fit) and summary(fit) to get parameter-wise summaries and plots, including uncertainty.