The result of the forward pass is the set of terms defined by $dirs and $cuts in earth’s return value.
Note that GCVs (via GRSq) are used during the forward pass only as one of the stopping conditions and in trace prints. Changing the penalty argument does not change the knot positions.
The various stopping conditions mean that the actual number of terms created by the forward pass may be less than nk. There are some other reasons why the actual number of terms may be less than nk
-
The forward pass discards one side of each term pair if it adds nothing to the model—but the forward pass counts terms as if they were actually created in pairs
-
As a final step, the forward pass deletes linearly dependent terms, if any, so all terms in $dirs and $cuts are independent. And remember that the pruning pass will further discard terms.
The Pruning Pass
The pruning pass is handed the sets of terms created by the forward pass. Its job is to find the subset of these terms that gives the lowest GCV. The pruning pass works like this: it determines the subset of terms (using pmethod) with the lowest RSS for each model size in 1:nprune (see the Force.xtx.prune argument above for some details). It saves the RSS and term numbers for each such subset in rss.per.subset and prune.terms. It then applies the Get. crit function with ppenalty to rss.per.subset to yield gcv.per.subset. It chooses the model with the lowest value in gcv.per.subset, and puts its term numbers into selected.terms. Finally, it runs lm to determine the fitted.values, residuals, and coefficients, by regressing the response y on the selected.terms of bx.
Set trace>=3 to trace the pruning pass.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%. In practice a dataset larger than the one below should be used for splitting. Also, remember that the test set should not be used for parameter tuning— use GCVs or separate validation sets for that.
train.subset <- sample(1:nrow(trees), .8 * nrow(trees))
test.subset <- (1:nrow(trees))[-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
Establishing Variable Importance
Establishing predictor importance is in general a tricky and even controversial problem.
Running plotmo with ylim=NULL (the default) gives an idea of which predictors make the largest changes to the predicted value.
You can also use drop1 (assuming you are using the formula interface to earth). Calling drop1(my.earth.model) will delete each predictor in turn from your model, rebuild the model from scratch each time, and calculate the GCV each time. You will get warnings that the earth library function extractAIC.earth is returning GCVs instead of AICs — but that is what you want so you can ignore the warnings. The column labeled AIC in the printed response from drop1 will actually be a column of GCVs not AICs. The Df column is not much use in this context.
You will get lots of output from drop1 if you built your original earth model with trace>0. You can set trace=0 by updating your model before calling drop1. Do it like this:
my.model <- update.earth(my.model, trace=0).
Remember that these techniques only tell you how important a variable is with the other variables already in the model. There are alternative ways of measuring variable importance (using resampling) but they are not yet implemented.
Which Predictors Were Added To the model first?
You can see the forward pass adding terms with trace=2 or higher. But remember, pruning will remove some of the terms. Another approach is to use
summary(my.model, decomp=”none”)
which will list the basis functions remaining after pruning, in the order they were added by the forward pass.
Which Predictors Are Actually Used In The Model?
The following function will give you a a list of the predictors in the model:
get.used.pred.names <- function(obj) # obj is an earth object
names(which(apply(obj$dirs[obj$selected.terms,,drop=FALSE],2,any)))
Why Do Are There Fewer Terms Than nk, Even With prune=”none”?
See the section above on the forward pass.
Multiple Response Models
If y has K columns then earth builds K simultaneous models. Each model has the same set of basis functions (i.e. same bx and selected.terms) but different coefficients (the returned coefficients will have K columns). The models are built and pruned as usual but with the GCVs and RSSs averaged across all K responses.
Since earth attempts to optimize for all models simultaneously, the results will not be as “good” as building the models independently. i.e. the GCV of the combined model will not be as good as the GCVS for independent models, on the whole. However, the combined model may be a better model in other senses, depending on what you are trying to achieve.
For more details on using GCVs averaged over multiple responses see section 4.1 of Hastie, Tibshirani, and Buja Flexible Discriminant Analysis by Optimal Scoring, JASA, December 1994
http://www-stat.stanford.edu/~hastie/Papers/fda.pdf.
Using Eearth with FDA and MDA
Earth can be used with fda and mda in the mda package. Earth will generate a multiple response model, as described above. Use keep.fitted=TRUE if you want to call plot.earth later (actually only necessary for large datasets, see the description of keep.fitted in fda). Use keepxy=TRUE if you want to call update or plotmo later. Use trace>=5 to see the call to earth generated by fda or mda. Example: library(mda)
(a <- fda(Species ~ ., data=iris, keep.fitted=TRUE, method=earth, keepxy=TRUE))
plot(a)
summary(a$fit) # examine earth model embedded in fda model
plot(a$fit)
plotmo(a$fit, ycolumn=1, ylim=c(-1.5,1.5), clip=FALSE)
plotmo(a$fit, ycolumn=2, ylim=c(-1.5,1.5), clip=FALSE)
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.
Author(s)
Stephen Milborrow, derived from mda::mars by Trevor Hastie and Robert Tibshirani.
References
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 the leaps package which is based on techniques in Miller.
Faraway Extending the Linear Model with R http://www.maths.bath.ac.uk/~jjf23
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 http://www-stat.stanford.edu/research/index.html
Hastie, Tibshirani, and Friedman (2001) The Elements of Statistical Learning http://www-stat.stanford.edu/~hastie/pub.htm
Miller, Alan (1990, 2nd ed. 2002) Subset Selection in Regression
