Answer to: Obtain prediction confidence intervals for GLS model predictions

After some more digging I found another solution using the marginaleffects package: library(marginaleffects) GLSout <- gls(follicles ~ sin(2*pi*Time) + cos(2*pi*Time), weights = varIdent(form = ~ 1|Mare), data = Ovary) predictOut3 <- marginaleffects::predictions(GLSout, vcov = TRUE, by = FALSE, type = "response") The marginaleffects CIs on predictions seem more conservative -- I've plotted them (in red) along side with the ones calculated with boostrapping (in blue)