Answer to: Obtain prediction confidence intervals for GLS model predictions

I think I might I've found a solution using the boot package, but would love if someone who's savvier in GLS and bootstrapping could check if I'm not messing something up. Here is what I came up with: library(boot) library(nlme) gls.fun <- function(data, idx, x.pred) { dt <- data[idx,] mod <- gls(follicles ~ sin(2*pi*Time) + cos(2*pi*Time), weights = varIdent(form = ~ 1|Mare), data = dt) c(predict(mod, newdata = x.pred)) } set.seed(123) bootstrap <- boot(Ovary, gls.fun, R = 1000, x.pred = Ovary) # Display the result of boot function CIs <- sapply(1:length(bootstrap$t0), FUN = function(i) { CI <- boot.ci(boot.out = bootstrap, index = i, type = "norm")$normal colnames(CI) <- c("conf", "lower", "upper") CI }, simplify = FALSE) CIs <- Reduce(rbind, CIs) DF <- data.frame(Time = Ovary$Time, follicles = Ovary$follicles, fitted = bootstrap$t0, lowerCI = CIs[, "lower"], upperCI = CIs[, "upper"]) library(ggplot2) ggplot(DF, aes(x = Time)) + geom_point(aes(y = follicles)) + geom_line(aes(y = fitted)) + geom_ribbon(aes(ymin = lowerCI, ymax = upperCI), alpha = 0.5)