Earth: Multivariate Adaptive Regression Splines

Build a regression model using the techniques in Friedman's papers Multivariate Adaptive Regression Splines and Fast MARS.

models, regression, smooth
## S3 method for class 'formula':
earth(formula, data, \dots)

## S3 method for class 'default': earth(x = stop("no 'x' arg"), y = stop("no 'y' arg"), subset = NULL, weights = NULL, na.action =, penalty = if(degree > 1) 3 else 2, trace = 0, degree = 1, nk = max(21, 2 * NCOL(x) + 1), thresh = 0.001, minspan = 0, newvar.penalty = 0, fast.k = 20, fast.beta = 1, fast.h = NULL, pmethod = "backward", ppenalty = penalty, nprune = NULL, Object = NULL, Get.crit = get.gcv, Eval.model.subsets = eval.model.subsets, Print.pruning.pass = print.pruning.pass, ...)

Model formula.
Data frame.
Matrix containing the independent variables.
Vector containing the response variable. If the y values are very big or very small, you may get better results if you scale y first.
Index vector specifying which rows in x and elements of y to use. Default is NULL, meaning all.
Weight vector (not yet supported).
NA action. Default is, and only is supported.
GCV penalty per knot. Default is if(degree>1) 3 else 2. A value of 0 penalises only terms, not knots. The value -1 is a special case, meaning no penalty, so GCV=RSS/n. Theory suggests values in the range of about 2 to 3.
Trace earth's execution. Default is 0. Values:

0 none 1 overview 2 forward 3 pruning 4 more pruning 5 ... The following arguments are for the forward pass

Maximum degree of interaction (Friedman's $mi$). Default is 1, meaning build an additive model.
Maximum number of model terms before pruning. Includes the intercept. Default is max(21,2*NCOL(x)+1). The number of terms created by the forward pass will be less than nk if there are linearly dependent terms
Forward stepping threshold. This is one of the arguments used to decide when forward stepping should terminate. See the section below on the forward pass. Default is 0.001.
Minimum distance between knots. Set trace>=2 to see the calculated value. Values: <0< code=""> add to the internally calculated min span (i.e. decrease span). 0 (default) use internally calculated min span as per Fri
Penalty for adding a new variable in the forward pass (Friedman's $gamma$, equation 74 in the MARS paper). This argument can mitigate the effects of collinearity or concurvity in the input data. Default is 0. Useful non-zero value
Maximum number of considered parent terms, as as described in Friedman's Fast MARS paper section 3.0. Default is 20. The special value -1 is equivalent to infinity, meaning no Fast MARS. Typical values, apart from -1, range from a
Fast MARS ageing coefficient, as described in the Fast MARS paper section 3.1. Default is 1. A value of 0 sometimes gives better results.
Fast MARS $h$, as described in the Fast MARS paper section 4.0. (not yet implemented). The following arguments are for the pruning pass
Pruning method. One of: backward none exhaustive forward seqrep. Default is "backward". Model subset evaluation for pruning uses the leaps package. Pruning can ta
Like penalty but for the pruning pass. Default is penalty.
Maximum number of terms (including intercept) in the pruned model. Default is NULL, meaning all terms. Use this to reduce exhaustive search time, or to enforce a maximum model size. Often used with
Earth object to be updated, for use by
Criterion function for model selection during pruning. By default a function that returns the GCV. See the section below on the pruning pass.
Function used to evaluate model subsets --- see notes in source code.
Function used to print pruning pass results. --- see notes in source code.
earth.formula: arguments passed to earth.default.

earth.default: unused, but provided for generic/method consistency.


  • An object of class earth which is a list with the components listed below. Term refers to a term created during the forward pass (each line of the output from is a term). Term number 1 is always the intercept.
  • fitted.valuesFitted values
  • residualsResiduals
  • coefficientsLeast squares coefficients for columns in bx. Each value corresponds to a selected term. coefficients[1] is the intercept.
  • rssResidual sum-of-squares of the model. Equal to rssVec[length(selected.terms)]. See also rssVec below.
  • rsq1-rss/rss.null. R-Squared of the model. A measure of how well the model fits the training data.
  • gcvGeneralised Cross Validation value (GCV) of the model. Equal to gcvVec[length(selected.terms)]. See also gcvVec below. For details of the GCV calculation, see equation 30 in Friedman's MARS paper and earth:::get.gcv.
  • grsq1-gcv/gcv.null. An estimate of the predictive power of the model.

    Unlike rsq, grsq can be negative. A negative grsq would indicate a severely over parameterised model --- a model that would not generalise well even though it may be a good fit to the training data. Example of a negative grsq:

    earth(mpg ~ ., data = mtcars, pmethod = "none", trace = 4)

  • bxMatrix of basis functions applied to x. Each column corresponds to a selected term. Each row corresponds to a row in in the input matrix x, after taking subset. See for an example of bx handling. Example:(Intercept) h(Girth-12.9) h(12.9-Girth) h(Girth-12.9)*h(... [1,] 1 0.0 4.6 0 [2,] 1 0.0 4.3 0 [3,] 1 0.0 4.1 0 ...
  • dirsMatrix with $ij$-th element equal to 1 if term $i$ has a factor of the form $x_j > c$, equal to $-1$ if term $i$ has a factor of the form $x_j \le c$, and to 0 if $x_j$ is not in term $i$. This matrix includes all terms generated by the forward.pass, including those not in selected.terms. Note that the terms may not be in pairs, because the forward pass deletes linearly dependent terms before handing control to the pruning pass.

    Example:Girth Height (Intercept) 0 0 #no factors in intercept h(Girth-12.9) 1 0 #2nd term uses Girth h(12.9-Girth) -1 0 #3rd term uses Girth h(Girth-12.9)*h(Height-76) 1 1 #4th term uses Girth and Height ...

  • cutsMatrix with $ij$-th element equal to the cut point for variable $j$ in term $i$. This matrix includes all terms generated by the forward.pass, including those not in selected.terms. Note that the terms may not be in pairs, because the forward pass deletes linearly dependent terms before handing control to the pruning pass.

    Example:Girth Height (Intercept) 0.0 0 #intercept, no cuts h(Girth-12.9) 12.9 0 #2nd term has cut at 12.9 h(12.9-Girth) 12.9 0 #3rd term has cut at 12.9 h(Girth-12.9)*h(Height-76) 12.9 76 #4th term has two cuts ...

  • selected.termsVector of term numbers in the best model. Can be used as a row index vector into cuts and dirs. The first element selected.terms[1] is always 1, the intercept.
  • rssVecResidual sum-of-squares for each model size considered by the pruning pass. The length of rssVec is nprune. The null RSS (i.e. the RSS of an intercept only-model) is rssVec[1]. The RSS of the selected model is rssVec[length(selected.terms)].
  • gcvVecGCV for each model in prune.terms. The length of gcvVec is nprune. The null GCV (i.e. the GCV of an intercept-only model) is gcvVec[1]. The GCV of the selected model is gcvVec[length(selected.terms)].
  • prune.termsThe row index of prune.terms is the model size (the model size is the number of terms in the model). Each row is a vector of term numbers for the best model of that size. An element is 0 if the term is not in the model, thus prune.terms is a lower triangular matrix, with dimensions nprune x nprune. The model selected by the pruning pass is at row length(selected.terms). Example:[1,] 1 0 0 0 0 0 0 #intercept-only model [2,] 1 2 0 0 0 0 0 #best 2 term model uses terms 1,2. [3,] 1 2 4 0 0 0 0 #best 3 term model uses terms 1,2,4 [4,] 1 2 9 8 0 0 0 ...
  • ppenaltyThe GCV penalty used during pruning. A copy of earth's ppenalty argument.
  • callThe call used to invoke earth.
  • termsModel frame terms. This component exists only if the model was built using earth.formula.


Standard Model Functions

Standard model functions such as case.names are provided for earth objects and are not explicitly documented.

Other Implementations

The results are similar to but not identical to other Multivariate Adaptive Regression Splines implementations. The differences stem from the forward pass where very small implementation differences (or perturbations of the input data) can cause rather different selection of terms and knots. The backward passes give identical or near identical results, given the same forward pass results.

The source code of earth is derived from mars in the mda package written by by Trevor Hastie and Robert Tibshirani. Unlike earth, mda::mars allows multiple responses. See also

The term MARS is trademarked and licensed exclusively to Salford Systems Their implementation uses an engine written by Friedman and offers more features than earth.


Multiple responses are not yet supported.

There is no special support for factors.

The following aspects of MARS are mentioned in Friedman's papers but not implemented in earth: i) Piecewise cubic models ii) Specifying which predictors must enter linearly iii) Specifying which predictors can interact iv) Model slicing (plotmo goes part way) v) Handling missing variables vi) Logistic regression and special handling of categorical predictors vii) Fast MARS h parameter.

The Forward Pass

The forward pass adds terms in pairs until the first of the following conditions is met: i) reach maximum number of terms (nterms>=nk). ii) reach DeltaRSq threshold (DeltaRSq. DeltaRSq is the difference in R-Squared caused by adding the current term pair. iii) reach max RSq (RSq>1-thresh). iv) reach min GRSq (GRSq< -10).

Thus you can change when forward stepping stops by adjusting nk, thresh, or penalty. Set trace>=2 to see the stopping condition.

The result of the forward pass is the set of terms defined by $dirs and $cuts. As a final step, the forward pass deletes linearly dependent terms, if any, so all terms in $dirs and $cuts are independent.

Note that GCVs (via GRSq) are used during the forward pass only as one of the stopping conditions and in trace prints.

The Pruning Pass

The pruning pass is handed the sets of terms created by the forward pass and works like this: it determines the subset of terms (using pmethod) with the lowest RSS for each model size in 1:nprune. It saves the RSS and term numbers for each such subset in rssVec and prune.terms. It then applies the Get.crit function with ppenalty to rssVec to yield gcvVec. Finally, it chooses the model with lowest value in gcvVec, and puts its term numbers into selected.terms.

By default Get.crit is earth:::get.gcv. Alternative Get.crit functions can be defined. See the source code of get.gcv for an example.

Testing on New Data

This example demonstrates one way to train on 80% of the data and test on the remaining 20%. Repeated runs of the code show the high variance of R-Squared associated with a model built from a small dataset from which many parameters have to be estimated. train.subset <- sample(1:nrow(ozone), .8 * nrow(ozone)) test.subset <- (1:nrow(ozone))[-train.subset] a <- earth(Volume~., data=trees[train.subset, ]) yhat <- predict(a, newdata=trees[test.subset, ]) y <- trees$Volume[test.subset] print(1 - sum((y - yhat)^2)/sum((y - mean(y))^2)) # print R-Squared Large Models and Execution Time

For a given set of input data, the following can increase the speed of the forward pass: i) increasing fast.k ii) decreasing nk iii) decreasing degree iv) increasing threshold v) increasing min.span. The backward pass is normally much faster than the forward pass, unless pmethod="exhaustive". Reducing npune reduces exhaustive search time. One strategy is to first do a forward pass with pmethod="none" and then use to adjust pruning parameters.

For big models, earth is much faster than mda::mars.

Using fast.k

With a low fast.k (say 5), plotmo is faster. In general, with a high fast.k, or with fast.k disabled (set to -1), plotmo builds a better model, but this is not always the case. The idea is to use high values of fast.k for exploratory work, and to disable fast.k for the best model possible. You will need to experiment using your data.

Warning and Error Messages

Earth prints most error and warning messages without printing the call. If you are mystified by a warning message, try setting options(warn=2) and using traceback.


  • regression
  • mars
  • Friedman


The primary references are the Friedman papers. Readers may find the MARS section in Hastie, Tibshirani, and Friedman a more accessible introduction. Faraway takes a hands-on approach, using the ozone data to compare mda::mars with other techniques. (If you use Faraway's examples with earth instead of mars, use $bx instead of $x). Earth's pruning pass uses leaps which is based on techniques in Miller.

Faraway Extending the Linear Model with R

Friedman (1991) Multivariate Adaptive Regression Splines (with discussion) Annals of Statistics 19/1, 1--141

Friedman (1993) Fast MARS Stanford University Department of Statistics, Technical Report 110

Hastie, Tibshirani, and Friedman (2001) The Elements of Statistical Learning

Miller, Alan (1990, 2nd ed. 2002) Subset Selection in Regression

See Also,, get.nused.preds.per.subset,,, ozone1,,, plotmo,,,,

  • earth
  • earth.default
  • earth.formula
a <- earth(Volume ~ ., data = trees)
summary(a, digits = 2)

# yields:
#    Call:
#    earth(formula = Volume ~ ., data = trees)
#    Expression:
#      23 
#      +  5.7 * pmax(0,  Girth -     13) 
#      -  2.9 * pmax(0,     13 -  Girth) 
#      + 0.72 * pmax(0, Height -     76) 
#    Number of cases: 31
#    Selected 4 of 5 terms, and 2 of 2 predictors
#    Number of terms at each degree of interaction: 1 3 (additive model)
#    GCV: 11     RSS: 213     GRSq: 0.96     RSq: 0.97
Documentation reproduced from package earth, version 0.1-2, License: GPL2

Community examples

Looks like there are no examples yet.