Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
code for plotting bootstrap predictions from a nonlinear model
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.