3
$\begingroup$

I am working to understand what is represented by the standard errors and confidence intervals of smooths of ordered factors. In this toy example, I wish to determine (1) the overall form of the trend (reference smooth), and (2) whether the trend shape depends on the factor level (difference smooth).

An example is the following (R code below). The model summary suggests that the by-smooth may capture differences with respect to the reference smooth.

Link function: identity 

Formula:
y ~ s(x) + s(x, by = fac) + fac

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.01639    0.01867   0.878    0.381
fac.L        0.01648    0.02641   0.624    0.533

Approximate significance of smooth terms:
            edf Ref.df       F p-value    
s(x)      8.335  8.866 192.290 < 2e-16 ***
s(x):fac2 2.631  3.347   4.773 0.00244 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.935   Deviance explained = 93.9%
GCV = 0.074581  Scale est. = 0.069746  n = 200

Following from Gavin Simpson's discussion of differences of smooths, I construct the confidence intervals:

Where the confidence interval excludes zero, we might infer significant differences between a pair of estimated smooths.

However, the confidence intervals of the difference smooth include zero everywhere. What, then, can we infer? Although some differences in the form of the smooths exists, we cannot identify them?

On a related note, although I think I understand the meaning of the confidence interval, I am not sure how to interpret the standard error of the smooth. What information do we obtain by seeing this plotted out? Does it matter where the standard error encompasses zero?

example smooth diff

require(mgcv)

# set up some toy data
set.seed(3)
n <- 100

x <- sort(rnorm(n, 10,2))

y1 <- scale( cos(x) ) + rnorm(n, 0,0.25)
y2 <- scale( cos(x/0.98) ) + rnorm(n, 0,0.25)

dat <- data.frame(x=rep(x,2), y=c(y1,y2), fac=factor(rep(c(1,2),each=n),ordered=TRUE))

# fit gam with parametric part of factor, reference smooth, and difference smooth  
mod <- gam(y~s(x)+s(x,by=fac)+fac, data=dat)                  
summary(mod)

# plot data, reference smooth, diff. smooth with se, and diff. smooth with CI
par(mfrow=c(1,4))
plot(y~x, col=fac, pch=16, data=dat,
     main='raw data')

modfits <- plot(mod, seWithMean=TRUE)
abline(h=0)
fits <- lapply( modfits, function(x) 
  {
    fit <- x$fit
se <- x$se
    upr <- fit+se*1.96
    lwr <- fit-se*1.96

    data.frame(x=x$x,fit,se,upr,lwr)

  } )
title(main='diff. smooth w/ se')
with(fits[[2]],
     {
       plot( fit~x, type='l', ylim=c(-3,3) )
       lines( upr~x, type='l', lty='dashed' )
       lines( lwr~x, type='l', lty='dashed' )
       abline(h=0)
     } )
title(main='diff. smooth w/ CI')
$\endgroup$
3
$\begingroup$

What you think is in se is not actually the standard error. It's an unfortunate choice of name for this component, but it actually contains the values se.mult * se.fit, with the latter being computed inside the plot method. se.mult is 2 by default as per the arguments of plot.gam.

Given that, the rest of your question is moot as you actually created an ~99.99557% confidence interval, which tells you something different from a 95% one and from the test in the output of summary().

plot.gam is a bit of a complicated function to read through as it does a lot of different things and handles many different types of smoothers. But you can confirm that se contains the requested ~95% confidence interval and not the standard errors of the smooth by following this example, which replicates the plot produced by plot.gam() for a single smoother in an ordered by factor smooth model.

library('mgcv')
set.seed(1)
dat <- gamSim(4)
dat <- transform(dat, fac = ordered(fac))

m <- gam(y ~ fac + s(x2) + s(x2, by = fac) + s(x0), data = dat, method = "REML")
plt_dat <- plot(m, pages = 1)

df <- plt_dat[[3]]
df$upr <- df$fit + df$se
df$lwr <- df$fit - df$se

## plot.gam does this by default, scale = -1;
## it looks at the range of fitted values +/- CI
## at sets limits as range over all smooths
foo <- function(x) {
  upr <- x$fit + x$se
  lwr <- x$fit - x$se
  range(upr, lwr)
}
ylim <- range(unlist(lapply(plt_dat, foo)))

layout(matrix(1:2, ncol = 2))
plot(m, select = 3)
plot(df$x, df$fit, type = "l", ylim = ylim, xlim = df$xlim)
lines(df$x, df$upr, lty = "dashed")
lines(df$x, df$lwr, lty = "dashed")
rug(df$raw, side = 1)
layout(1)

which produces

The plot on the left is the one produced by <code>plot.gam</code> and the one on the right is the one I created from the data returned by <code>plot.gam</code>.

The plot on the left is the one produced by plot.gam and the one on the right is the one I created from the data returned by plot.gam.

$\endgroup$
3
  • 1
    $\begingroup$ This simultaneously explains why I was reading references to the 'CI' returned by plot.gam(). Brilliant, thank you very much. $\endgroup$
    – noname
    May 22 '19 at 16:46
  • $\begingroup$ Follow-up: for plot.gam(mod,all.terms=TRUE), termplot() is called to plot the partial effects of the parametric terms. termplot() does plot the actual standard errors when se=TRUE. So the smooth nonparametric components will be plotted with the CIs, but the parametric components with the SE? $\endgroup$
    – noname
    May 22 '19 at 16:58
  • 1
    $\begingroup$ Yes; ?termplot indicates that termplot only knows how to plot point-wise standard errors. This is all a bit of a fudge; I get the impression this was just thrown in as an after-thought as there is no consistency here. $\endgroup$ May 22 '19 at 17:57

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.