I am using earth packageearth: Multivariate Adaptive Regression Spline Models regression to get a constant piecewise approximation of my data. I want to plot a band of confidence around it. Does this make sense to estimate a confidence interval of the smoothed function? If yes how can I do this?

I know that confidence intervals cannot be calculated directly, since it is a non parametric regression but I want to have a coherent result like (plot a band of confidence around my smoothed function) what I can get using lm or loess smoothing.


I add some code to clarify my idea, Here what I would do if I am using linear regression or loess.

x <- rnorm(15)
dat <- data.frame(x = x, y = c(x + rnorm(15)))

Using lm

mod <- lm(y ~ x,data=dat)
predfit <- predict(mod,se=TRUE,interval="confidence")$fit

Using loess

mod <- loess(y~x,data=dat)
pred <- predict(mod, se = TRUE)
y = pred$fit
    ci <- pred$se.fit * qt(level / 2 + .5, pred$df)
data.frame(ymin = y - ci,
           ymax = y + ci)

My question how to get ymin and ymax if I use earth :

mod = earth(y~x,data=dat)
  • $\begingroup$ please edit your source data.frame "dat". When I try to create it a error occurs: (longer object length is not a multiple of shorter object length).Thx $\endgroup$ Jan 20 '14 at 15:24
  • $\begingroup$ Yes, it makes sense to estimate a confidence band of the smoothed function. A reference is the scb function for simultaneously confidence band in package locfit. $\endgroup$
    – Randel
    Jan 23 '14 at 4:15
  • $\begingroup$ Wasserman (2006) stated that, "For more complicated methods like trees, MARS and projection pursuit regression, I am not aware of rigorous results that lead to valid bands." Maybe we need to wait for the breakthrough in theory. $\endgroup$
    – Randel
    Jan 23 '14 at 19:07
  • $\begingroup$ @Randel thanks for your interest. I read this also, but I think we can get some practical results using bootstrapping methods... $\endgroup$
    – agstudy
    Jan 23 '14 at 19:17
  • $\begingroup$ I consulted a professor about this, and he told me that it is possible to construct the confidence interval for a fixed point, but impossible to construct the confidence band. But a non-rigorous variant of confidence band, variability band, can be constructed. $\endgroup$
    – Randel
    Feb 6 '14 at 20:17

Earth will estimate prediction intervals for you. Do it like this:

mod = earth(y~x,data=dat, ncross=30, nfold=3, varmod.method="lm")
plotmo(mod, pt.col=1, level=.95) # show prediction intervals

summary(mod) will print this:

varmod: method "lm"    min.sd 0.124    iter.rsq 0.016

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)       1.2265      0.2059           17
y                 0.0883      0.1918          217

                          mean   smallest   largest   ratio
95% prediction interval   4.86       3.91       5.1     1.3

The above code illustrates the principle, but a problem with your small sample size is that there is not enough data to estimate prediction intervals reliably. (We know this because the iter.stderr% are very big, and the iter.rsq is very small. There is a vignette that comes with the earth package that explains in detail how to get prediction intervals, and some of the potential pitfalls.)


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.