0

I have a basic GAM model setup, with one predictor:

fit <- gam(response ~ s(predictor, k=12), data = data, )

I noticed I can easily pull out the fitted values with the fit$fitted.values accessor.

The documentation says the confidence intervals are computed in a bayesian manner, but I wasn't able to find them in the returned object. It seems to somehow be related to fit$hat.

I could feed the model back into predict but that just produces another regression model, does it not? I'd like to use the confidence bans from the mgcv package GAM model itself.

Learning stats by example
  • 211
  • 4
  • 17
  • ¿Does the function fitted(fit) produce anything? If so, try confint(fit). Next option is to try names(fit) and names(summary(fit)) to see what you might be able to access. – Gregg H Apr 21 '18 at 14:38
  • 3
    https://stats.stackexchange.com/questions/33327/confidence-interval-for-gam-model – D.G Apr 21 '18 at 14:39

1 Answers1

2

Don't use fit$fitted.values use fitted(fit) to access the fitted values.

$hat contains the leading diagonal of the hat (influence) matrix, i.e. the hat values. Again, access them using influence().

No; using the predict() method does not refit the model nor produce another model. It will produce a range of predictions depending on argument type.

As the post linked to in the comment shows, you need to generate predictions from the model over the range of predictor and compute a normal point-wise interval from the standard errors returned by predict() when se.fit = TRUE. If this is a GAM (i.e. using a non-Gaussian family) then this must be done with type = "link", the default.

You can do this quickly using the confint() method from the schoenberg package, currently on github:

## install.packages("devtools")
## devtools::install_github("gavinsimpson/schoenberg")
library("schoenberg")
library("mgcv")
set.seed(2)
dat <- gamSim(1, n = 400, dist = "normal", scale = 2)
mod <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = dat, method = "REML")
##'
## point-wise interval
ci <- confint(mod, parm = "x1", type = "confidence")
head(ci)

This produces:

>     head(ci)
  smooth           x1      est        se     crit    lower    upper
1  s(x1) 0.0006632213 5.888546 0.3134611 1.959964 5.274174 6.502919
2  s(x1) 0.0056813456 5.895688 0.3045623 1.959964 5.298757 6.492619
3  s(x1) 0.0106994698 5.902830 0.2958458 1.959964 5.322983 6.482677
4  s(x1) 0.0157175940 5.909974 0.2873297 1.959964 5.346818 6.473130
5  s(x1) 0.0207357183 5.917121 0.2790332 1.959964 5.370226 6.464016
6  s(x1) 0.0257538425 5.924272 0.2709761 1.959964 5.393169 6.455376

If, as here, you don't provide any newdata it will generate new data over the range of the indicated covariate (x1 in this case) to generate a fine enough set of values over the range of the spline to produce a nice smooth plot.

You can always provide it your data (via argument newdata) to generate prediction and confidence intervals for your observations.

Gavin Simpson
  • 37,567
  • 5
  • 110
  • 153