#### MARS --- We now really only use the 'earth' package

help(trees) ## trees data {standard R 'datasets'} :
## Girth, Height and Volume for Black Cherry Trees
##---
## Description:
##---
##      This data set provides measurements of the girth, height and
##      volume of timber in 31 felled black cherry trees.  Note that girth
##      is the diameter of the tree (in inches) measured at 4 ft 6 in
##      above the ground.
##---
##      A data frame with 31 observations on 3 variables.
##        'Girth'   numeric  Tree diameter in inches
##        'Height'  numeric  Height in ft
##        'Volume'  numeric  Volume of timber in cubic ft

library(earth)
Mfit <- earth(Volume ~ ., data = trees)
summary(Mfit)
## or as in lecture:
summary(Mfit, style = "pmax", details=TRUE)
## Call: .....

## Volume =
## 29.05995
## -  3.419806 * pmax(0,   14.2 -  Girth)
## +  6.229514 * pmax(0,  Girth -   14.2)
## + 0.5813644 * pmax(0, Height -     75)

## Number of cases: 31
## Selected 4 of 5 terms, and 2 of 2 predictors
## Termination condition: RSq changed by less than 0.001 at 5 terms
## Importance: Girth, Height
## Number of terms at each degree of interaction: 1 3 (additive model)
## GCV 11.25439    RSS 209.1139    GRSq 0.959692    RSq 0.9742029

## Now, read
?earth
## --> default 'degree' := maximal degree of interaction is 1  <==>  additive model
Mf2 <- earth(Volume ~ ., data = trees,
             degree = 2, # <- interaction up to 2
             trace = 3) # <- show algorithm progress
## and the result is the same *after* pruning.

predict(Mfit,
        data.frame(Girth= 5:15, Height= seq(60,80, length=11)))
## Larger but more interesting:
Girth. <- 5:15
Height. <- 60:80
newDat <- expand.grid(Girth = Girth.,
                      Height= Height.)
dim(newDat)#- 231 x 2 ; 231 = 11 x 21
predVol <- predict(Mfit, newDat)

persp(Girth., Height.,
      matrix(predVol, length(Girth.),length(Height.)),
      zlab = "predicted  Volume", col = "lightblue", shade = 0.85,
      ticktype = "detailed", theta = -40)

## or with dynamic 3D-Graphics
library(rgl)
persp3d(Girth., Height.,
        matrix(predVol, length(Girth.),length(Height.)),
        col = "lightblue", zlab = "predicted  Volume")
rgl.close()


if(FALSE) { ## --- old ---------------------

### MARS --- exists in package  "mda" --- as said in manuscript,
help(mars, package = "mda", offline=TRUE)

### but also in  Kooperberg's  "polspline" :
help(polymars, package = "polspline", offline=TRUE)

} ## old



### MARS (earth) for the ozone data :
### ---- MM: this is just ./mars-ozone.R
###                         ~~~~~~~~~~~~

require(earth)

### The example from 'plot.earth' :
data(ozone1)
earth.mod <- earth(O3 ~ ., data = ozone1, degree = 2)
plot(earth.mod)

summary(earth.mod)

## Ok, try alternatives:
summary( earth.mod <- earth(O3 ~ ., data = ozone1, degree = 2, nk = 40) )
## doing crossvalidation -- costs more!
summary( earth.modC <- earth(O3 ~ ., data = ozone1, degree = 2, nk = 50,
                             nfold=16, ncross = 4) )
plot(earth.modC)
plotmo(earth.mod)

if(FALSE) ## NB: update() does *not* work with  earth() fits  [FIXME ?]
    loz.mars <- update(earth.mod, log(O3) ~ .) #

### "playing around" -- trying versions of the above

summary(loz.earth <- earth(log(O3) ~ ., data = ozone1, degree = 2))
summary(loz.earth <- earth(log(O3) ~ ., data = ozone1, degree = 2, nk = 30))
summary(loz.earth) # still reached boundary 30
summary(loz.earth <- earth(log(O3) ~ ., data = ozone1, degree = 2, nk = 50))

summary( earth.L.oz <- earth(log(O3) ~ ., data = ozone1, degree = 2, nk = 50, nfold=16, ncross = 4) )
plotmo(earth.L.oz)
plotmo(earth.mod)

## Finally omitting "the outlier" (see gam() etc)
summary( earth.L.oz <- earth(log(O3) ~ ., data = ozone1, subset=-92, degree = 2, nk = 50,
                             nfold=16, ncross = 4) )