code for plotting bootstrap predictions from a nonlinear model

boot.R

set.seed(1)

# Create fake data
x = runif(100, 0, 5)
y = .25 * x^3 - x^2 + rnorm(length(x))
data = data.frame(x = x, y = y)

# Identify 500 points to include in the plots
x.sequence = seq(0, 5, length = 500)

# I'm fitting GAMs because I'm more familiar with them.
# Procedure should be identical for nls.
library(mgcv)

# Fit models to 100 bootstrap replicates of the data
predictions = replicate(
  100,{
    boot = data[sample.int(nrow(data), replace = TRUE), ]
    model = gam(y ~ s(x), data = boot)
    # Output predictions at each point that we'll want to plot later
    predict(model, data.frame(x = x.sequence))
  }
)

# "spaghetti plot" of all 100 predicted means at each of the 500 plotting points
matplot(x.sequence, predictions, type = "l")


# Plot raw data
plot(y ~ x)

# Add bootstrap mean at each of the x.sequence points
# Alternatively, you could plot the predictions of the full model.
lines(rowMeans(predictions) ~ x.sequence, lwd = 3, type = "l")

# Plot bootstrap CI at each of the x.sequence points
lines(apply(predictions, 1, quantile, .975) ~ x.sequence)
lines(apply(predictions, 1, quantile, .025) ~ x.sequence)