14
$\begingroup$

Reading mgcv::gam's help page:

confidence/credible intervals are readily available for any quantity predicted using a fitted model

However I can't figure a way to actually get one. I thought predict.gam would have a type=confidence and a level parameter but it doesn't. Can you help me on how to create it ?

$\endgroup$
33
$\begingroup$

In the usual way:

p <- predict(mod, newdata, type = "link", se.fit = TRUE)

Then note that p contains a component $se.fit with standard errors of the predictions for observations in newdata. You can then form CI by multipliying the SE by a value appropriate to your desired level. E.g. an approximate 95% confidence interval is formed as:

upr <- p$fit + (2 * p$se.fit)
lwr <- p$fit - (2 * p$se.fit)

You substitute in an appropriate value from a $t$ or Gaussian distribution for the interval you need.

Note that I use type = "link" as you don't say if you have a GAM or just an AM. In the GAM, you need to form the confidence interval on the scale of the linear predictor and then transform that to the scale of the response by applying the inverse of the link function:

upr <- mod$family$linkinv(upr)
lwr <- mod$family$linkinv(lwr)

Now note that these are very approximate intervals. In addition these intervals are point-wise on the predicted values and they don't take into account the fact that the smoothness selection was performed.

A simultaneous confidence interval can be computed via simulation from the posterior distribution of the parameters. I have an example of that on my blog.

If you want a confidence interval that is not conditional upon the smoothing parameters (i.e. one that takes into account that we do not know, but instead estimate, the values of the smoothness parameters), then add unconditional = TRUE to the predict() call.

Also, if you don't want to do this yourself, note that newer versions of mgcv have a plot.gam() function that returns an object with all data used to create the plots of the smooths and their confidence intervals. You can just save the output from plot.gam() in an obj

obj <- plot(model, ....)

and then inspect obj, which is a list with one component per smooth. Add seWithMean = TRUE to the plot() call to get confidence intervals that are not conditional upon smoothness parameter.

$\endgroup$
  • $\begingroup$ Doing the simultaneous CI and the parametric bootstrap get a bit more involved with the code so if you can get by with just the pointwise intervals great. If not I can provide further examples for each of those. $\endgroup$ – Reinstate Monica - G. Simpson Jul 30 '12 at 7:26
  • $\begingroup$ +1 for the answer. Impressive blog post indeed, I'm going to be studying it for a while to improve my graphics skills. $\endgroup$ – jbowman Jul 30 '12 at 17:11
  • $\begingroup$ Any way that I could get access to that impressive blog post (ucfagls.wordpress.com/2011/06/12/…)? Currently the blog requires a login. $\endgroup$ – geneorama Feb 3 '15 at 19:06
  • $\begingroup$ @geneorama I moved my blog away from Wordpress and for a year paid for redirects to the new one for all URLs but I let that lapse recently. Sorry about that. I've edited in the new link, and that doesn't require a login. (The login is to avoid two copies of same post and I've been too lazy to delete the pages from the Wordpress site as of yet.) $\endgroup$ – Reinstate Monica - G. Simpson Feb 3 '15 at 23:04
  • $\begingroup$ The original blog post (see the edit history of this Q&A) had a fundamental flaw in the way the simultaneous interval was created. The link in the current (as of Dec 2016) version of the Answer computes the simultaneous interval correctly. $\endgroup$ – Reinstate Monica - G. Simpson Dec 16 '16 at 18:00
5
$\begingroup$

If you just want to plot them the plot.gam function has shading that defaults to confidence intervals using the shade argument. Also see gam.vcomp for getting the intervals.

$\endgroup$
5
$\begingroup$

The package mgcv (newer than gam) readily plots credible intervals. This Bayesian approach is different from confidence intervals, but the results are almost the same, as numerical simulations have shown (see the paper by Marra and Wood linked in mgcv).

$\endgroup$
  • 2
    $\begingroup$ +1 A key result of Marra & Wood's paper is that they develop Nychka's understanding/explanation of why the empirical Bayes credible intervals also have pretty extraordinary frequentist interpretation/behaviour when viewed as "across-the-function" confidence intervals. You can treat the intervals in a Bayesian or frequentist manner and the coverage property implied by the $1-\alpha$ interval holds, approximately. $\endgroup$ – Reinstate Monica - G. Simpson Dec 16 '16 at 18:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.