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

I am using the earth package in R to estimate the number of breakpoints in a curve. There is only a single predictor.

I was hoping there was a sensible way to compare a slope for some segment against 0. Normally I would use the estimate and standard error, but it's not obvious (to me) what to do for segments that are combinations of multiple slopes:

        estimate    SE
(intercept) 2.78        .015
h(age-4)    -.01        .0029
h(4-age)    -.35        .015
h(age-11)   -.076       .0055


So the slope after age 11 is -.086. How do I quantify my confidence that this slope really is negative?

• MARS is a useful algorithm, but terrible for this application. {segmented} should be better. Or if you know the number of breakpoints segmented regression Aug 11 '14 at 18:29

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.

To some extent, MARS has already "tested" the slopes in so far as predictive ability via the CV procedure used. Of course, this assumes the CV has been done in a way that matches the experimental design or dependencies in the data.

If you wanted to get better estimates of the standard errors, I would suggest a bootstrap approach, which is computationally demanding but will give you SEs or a distribution from which to compute a confidence interval

You'll need to bootstrap the whole process of course. Start with the boot package and write a wrapper to the MARS function accordingly to fit the model to a given bootstrap sample and return the slopes from the wrapper.

• Resampling is the approach which seems sensible here. I would only caution that you need to think through your application and determine whether or not your response+predictor sets are statistically dependent upon one another or not. If they are, then you want to use something like the Politis and Romano STATIONARY BOOTSTRAP rather than a vanilla one which will treat these as independent. The Politis and Romano algorithm is available in the R package TSERIES within the function TSBOOTSTRAP, using type == "stationary" as an option. There are other parameters you will need to check sensibly. Dec 6 '16 at 4:58