davharris · June 12, 2014 20:17 · Jun 12, 2014
diff --git a/boot.R b/boot.R
@@ -0,0 +1,38 @@
+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)
No results found