# Confidence Interval for non-smooth term in gam (mgcv)

I fitted a gam model in mgcv package and now want to get the confidence intervals for the non-smooth term.

require(datasets)
require(mgcv)
b = gam(Temp ~ s(Ozone) + Solar.R, data=airquality, family=gaussian)


So, when I wanted to print the confidence intervals, I ran:

summary(b)


and it prints only

Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 77.8580645  1.4091014  55.254   <2e-16 ***
Solar.R     -0.0003532  0.0069707  -0.051     0.96


So, I am interested in the confidence intervals for the non-smooth term Solar.R. How can I do this?

Thank you.

You need the standard error of the term, which is given in the output from summary(), but you can grab the standard errors from the diagonals of the variance-covariance matrix, which is extracted using vcov().

require(datasets)
require(mgcv)
b <- gam(Temp ~ s(Ozone) + Solar.R, data=airquality, family=gaussian)


Grab the coefficients:

beta <- coef(b)


...and the standard errors, which are the square roots of the diagonals of the variance-covariance matrix $$\mathrm{V}$$, noting that here I use the variance-covariance matrix that is adjusted to account for the selection of the smoothness parameter for the smooth in the model, $$\mathrm{V_b}$$

Vb <- vcov(b, unconditional = TRUE)
se <- sqrt(diag(Vb))


Next, identify the $$\hat{\beta}$$ etc for the term you want

i <- which(names(beta) == "Solar.R")


Then an approximate 95% confidence interval is

beta[i] + (c(-1,1) * (2 * se[i]))


This gives:

> beta[i] + (c(-1,1) * (2 * se[i]))
 -0.01429463  0.01358823


The 2 above comes from the 0.975th probability quantile of the Gaussian distribution.

> qnorm(0.975)
 1.959964


Often we take into account the additional uncertainty of having estimated the variance of the residuals by using the t distribution, for which we need the residual degrees of freedom

rdf <- df.residual(b)
qt(0.975, df = rdf)


which gives

> qt(0.975, df = rdf)
 1.982583


Which is close to 2 also. This distinction only really makes any difference when you have a small data set, so you'll generally see 95% intervals given as plus/minus 2 * standard error in mgcv.

• Thank you! Just one more question: as Pr(>|t|) is 0.96 , so above 0.05, I would have expected a positive value for the upper limit of the 95 % CI. Why isn't this the case? Thank you. – Nad Mar 25 '19 at 15:12
• Because I appear to have lost my sense. The Standard errors are the square roots of the diagonal elements of the variance-covariance matrix. Let me fix that. – Gavin Simpson Mar 25 '19 at 17:01
• Thank you! Worked perfectly! – Nad Mar 28 '19 at 10:27

Gavin's useful answer can be generalized by the following in the situation where you may have more than one coefficient of interest. Let's say these coefficients are called c1, c2, c3, c4 for example.

Replace

i <- which(names(beta) == "Solar.R")

with

i <- which(names(beta) %in% c("c1", "c2", "c3", "c4"))

Then replace

beta[i] + (c(-1,1) * (2 * se[i]))

with

my.ci <- data.frame(cbind(estimate = beta[i], lci = beta[i]-2se[i], uci = beta[i] + 2se[i]))

If the estimates were on a log scale, for example if the GAM was for a binomial outcome, then the estimates can be simply converted and rounded to 2 decimal places as follows

round(exp(my.ci), 2)