Solved – Multi-adaptive Regression Splines (MARS) equation interpretation

marsnonlinear regressionrregressionsplines

I was curious about certain aspects of multi-adaptive regression splines (MARS or earth in R). Taking the equation of MARS from the wiki page:

  ozone = 5.2 +
          + 0.93max(0,temp - 58)
          - 0.64max(0,temp - 68)
          - 0.046max(0,234-ibt)
          - 0.016max(0,wind-7)max(0,200-vis)

My questions are:

1) How to interpret the above equation? If I understand it correctly, for the following value of temp = 52, ibt = 230, wind = 8 and vis = 190, the above equation will become:

      ozone = 5.2 +
          + 0.93max(0,52 - 58)
          - 0.64max(0,52- 68)
          - 0.046max(0,234- 230)
          - 0.016max(0,8-7)max(0,200-190)

    ozone = 5.2 +
          + 0.93max(0,-6)
          - 0.64max(0,-16)
          - 0.046max(0,4)
          - 0.016max(0,1)max(0,10)

    ozone = 5.2 +
          + 0
          - 0
          - 0.184
          - 0.16

Is this correct?

2) Before running earth (MARS), do the predictor variables have to uncorrelated? What is the way to go about it: remove the correlated predictor variables and then run MARS or run MARS using all predictors and somehow MARS deal with the correlation

3) In R, for earth equation, how do I take out the final equation in step 1 for my model

4) Lets say I run two MARS model: first model predictor x1, x2 and x3 are retained in the final model and in second model x1 and x2 are retained in the final model. Is there any way I can combine the two MARS equation to have a single equation presumably where parameters of x1 and x2 are average of the two models and x3 on its own.

EDIT

To elaborate on (4), let's say I collect some data and and develop a MARS model

 ozone1 = 5.2 + 0.93max(0,temp - 58) - 0.64max(0,temp - 68)- 0.046max(0,234-ibt)- 0.016max(0,wind-7)max(0,200-vis)

I go back and collect another set of data and develop a second model:

  ozone2 = 4.6 + 0.89max(0,temp - 58) - 0.033max(0,234-ibt)- 0.016max(0,wind-7)max(0,200-vis)

Can I combine the two equation in a single equation something like this:

   ozone.fin = (5.2 + 4.6)/2 + ((0.93 + 0.89)/2)max(0,temp - 58) - 0.64max(0,temp - 68) - (0.046+0.033/2)max(0,234-ibt) - 0.016max(0,wind-7)max(0,200-vis)

Thaks

Best Answer

1) This is correct. It would be best to interpret the functions locally, as different areas in the domain have different slopes.

2) MARS exhibits problems in choosing among predictor variables when multicollinearity is present. To improve the ability of MARS to deal with multicollinearity, principal components can be used to reduce the dimensionality of the input variables before invoking MARS. However, slightly correlated predictors should not give problems.

3) Use the summary function to see the results

4) I do not see exactly what you mean there? Can you explain it any further?

Related Solutions

Solved – Testing slopes in multivariate adaptive regression splines (MARS/earth)

The short answer: you can't quantify that precisely.

You can have confidence in the model AS A WHOLE (because of the GCV, see for example the FAQ in the earth package vignette "What is a GCV in simple terms?") but can't easily quantify the significance of each part of the model.

The problem is actually worse than quantifying confidence in the effect of two combined hinge functions, because we can't even quantify the effect of each hinge function in isolation. The standard error values you quote are unlikely to be valid. (I assume you got those SEs by doing a linear regression on the MARS basis matrix i.e. lm(y ~ earth.model$bx - 1) ).

If we were certain of the position of a knot, we can do a standard linear regression on the slope at that knot, and get standard linear regression SEs and p-values for that slope. But of course with a different sample, MARS might choose a somewhat different position of the knot. Since we aren't certain of the position of the knot, those p-values overstate certainty, often dramatically.

Another way of looking at is to remember that the MARS algorithm does a lot of work to select only good terms. Terms that would have a bad p-value are tossed out by the MARS forward pass during the knot selection process. So lm(y ~ earth.model$bx - 1) will always give good p-values, because all terms in bx have been pre-selected to be "good". But all the uncertainty in the selection process is ignored by lm. To quote the earth vignette Section 3: variable selection overstates the significance of the selected variables. This is true in general and also in the specific case of MARS.

To gain confidence in the model as a whole, one possibility is to do as Gavin Simpson says --- bootstrapping. The earth function has built in cross-validation facilities, so cross-validation would be an easier way to go than bootstrapping. Each bootstrap or cross-validation iteration will give slightly different hinge functions, so it doesn't help in formally quantifying confidence in the position of the hinges in the full model, although may give some informal confidence. Earth will soon have facilities for estimating prediction intervals that may help in terms of the big picture of what you are trying to do.

Solved – Adding a quadratic term flips the signs of the coefficients

This is a Q&A site for questions about statistics, not for questions about how R works or how to use it. Your question #1 is off topic here, and I don't see why you want us to say if what they told you on Stack Overflow was correct—that is the appropriate site for R coding questions. At any rate, you can get the canonical answer in the documentation for formulae for statistical models in the Intro to R manual:

y ~ poly(x,2)
...
Polynomial regression of y on x of degree 2. ...

...
y ~ A*B
y ~ A + B + A:B
...
Two factor non-additive model of y on A and B. ...

From these you can see that you can get an interaction between two polynomials like so:

poly(x,2)*poly(y,2)

Regarding your question about the sign changing, there are several ways to think about this. First of all, when you add another variable, the signs of existing variables can change (see: here, and here, e.g.), and from the model's perspective, a squared term is just another variable (see: Why is polynomial regression considered a special case of multiple linear regression?). On the other hand, I have argued that we should consider all polynomial terms together when interpreting a model (cf., Does it make sense to add a quadratic term but not the linear term to a model?). Still, it may seem weird to think that the squared term is correlated with the original variable, and it is legitimate to wonder what that means.

For what it's worth, if you use poly(x,2), you will get orthogonal polynomials and there won't be an issue; likewise if you center your x variable before squaring, cf., here. However, the issue remains if you simply square your x variable in the most intuitive way (as you did), or use poly(x, 2, raw=TRUE). So, what does it mean for the coefficient on x to become, say, negative where it was positive before? As @whuber explains in his answer here, when you have (only) an $x$ and an $x^2$ term, you will have a parabola, and the x-value of the location of the apex of the parabola will be related to those coefficients. More specifically, if we call the coefficients on $x$ and $x^2$ $\beta_1$ and $\beta_2$ (eliding the hats for simplicity), the x position of the apex is $\beta_1/-2\beta_2$. (This is complicated, because it is also contingent on the value of $\beta_2$. Moreover, $\beta_2$ will be positive when the parabola is concave up—and vice versa; without loss of generality, let's assume that is your case.) So here a negative value of $\beta_1$ makes the x location of the apex to the right of $0$. You can see this demonstrated in the example below (coded in R):

set.seed(981)                 # this makes the example exactly reproducible
x = runif(100, min=0, max=5)  # here I generate some data
y = 0 -1*x + .5*x**2 + rnorm(100, mean=0, sd=.1)

m.l = lm(y~x)                 # here I fit linear & quadratic models
m.p = lm(y~poly(x, 2, raw=TRUE))
summary(m.l)$coefficient      # these are the fitted coefficients:
#              Estimate Std. Error   t value     Pr(>|t|)
# (Intercept) -2.103753 0.18623242 -11.29639 1.911436e-19
# x            1.445807 0.06430754  22.48270 1.822846e-40
summary(m.p)$coefficient
#                             Estimate Std. Error     t value     Pr(>|t|)
# (Intercept)              0.009625853 0.02979478   0.3230718 7.473363e-01
# poly(x, 2, raw = TRUE)1 -1.008168121 0.02698804 -37.3561142 2.216150e-59
# poly(x, 2, raw = TRUE)2  0.501424686 0.00533965  93.9059119 4.859026e-97

xs   = seq(0,5, by=.1)        # this code makes the plot
apex = coef(m.p)[2] / (-2*coef(m.p)[3])
windows()
  plot(x, y)
  abline(coef(m.l), col="blue", lty=2)
  lines(xs, predict(m.p, data.frame(x=xs)), col="red")
  points(apex, predict(m.p, data.frame(x=apex)), col="red", pch="*", cex=5)
  legend("topleft", legend=c("linear fit", "quadratic fit", "apex"),
         lty=c(2,1,NA), pch=c(NA, NA, "*"), col=c("blue","red","red"))

Best Answer

Related Solutions

Solved – Testing slopes in multivariate adaptive regression splines (MARS/earth)

Solved – Adding a quadratic term flips the signs of the coefficients

Related Question