> # test.varmmod.mgcv.R
> # mgcv has to be tested separately because of clashes between library(gam) and library(mgcv)
> # Stephen Milborrow Apr 2015 Berea
> 
> source("test.prolog.R")
> library(earth)
Loading required package: Formula
Loading required package: plotmo
Loading required package: plotrix
Loading required package: TeachingDemos
> options(warn=2)
> 
> printh <- function(caption)
+     cat("===", caption, "\n", sep="")
> 
> CAPTION <- NULL
> 
> multifigure <- function(caption, nrow=3, ncol=3)
+ {
+     CAPTION <<- caption
+     printh(caption)
+     par(mfrow=c(nrow, ncol))
+     par(cex = 0.8)
+     par(mar = c(3, 3, 5, 0.5)) # small margins but space for right hand axis
+     par(mgp = c(1.6, 0.6, 0))  # flatten axis elements
+     oma <- par("oma") # make space for caption
+     oma[3] <- 2
+     par(oma=oma)
+ }
> do.caption <- function() # must be called _after_ first plot on new page
+     mtext(CAPTION, outer=TRUE, font=2, line=1, cex=1)
> old.par <- par(no.readonly=TRUE)
> 
> library(mgcv)
Loading required package: nlme
This is mgcv 1.8-28. For overview type 'help("mgcv-package")'.
> 
> for(varmod.method in c("gam", "x.gam")) {
+ 
+     multifigure(sprint("varmod.method=\"%s\"", varmod.method), 2, 3)
+     par(mar = c(3, 3, 2, 3)) # space for right margin axis
+ 
+     set.seed(6)
+     earth.mod <- earth(Volume~Girth, data=trees, nfold=3, ncross=3,
+                        varmod.method=varmod.method,
+                        trace=if(varmod.method %in% c("const", "lm", "power")) .3 else 0)
+     printh(sprint("varmod.method %s: summary(earth.mod)", varmod.method))
+     printh("summary(earth.mod)")
+     print(summary(earth.mod))
+ 
+     # summary(mgcv) prints environment as hex address which messes up the diffs
+     printh("skipping summary(mgcv::gam) etc.\n")
+ 
+     printh(sprint("varmod.method %s: predict(earth.mod, interval=\"pint\")", varmod.method))
+     pints <- predict(earth.mod, interval="pint")
+     print(pints)
+ 
+     plotmo(earth.mod$varmod, do.par=FALSE, col.response=2, clip=FALSE,
+            main="plotmo residual model",
+            xlab="x", ylab="varmod residuals")
+ 
+     plotmo(earth.mod, level=.90, do.par=FALSE, col.response=1, clip=FALSE,
+            main="main model plotmo Girth")
+     do.caption()
+ 
+     plot(earth.mod, which=3, do.par=FALSE, level=.95)
+ 
+     # plot.varmod
+     plot(earth.mod$varmod, do.par=FALSE, which=1:3, info=(varmod.method=="earth"))
+ }
===varmod.method="gam"
===varmod.method gam: summary(earth.mod)
===summary(earth.mod)
Call: earth(formula=Volume~Girth, data=trees,
            trace=if(varmod.method%in%c("const","lm","power"))0.3els...),
            nfold=3, ncross=3, varmod.method=varmod.method)

              coefficients
(Intercept)      28.766764
h(13.8-Girth)    -3.427802
h(Girth-13.8)     6.570747

Selected 3 of 4 terms, and 1 of 1 predictors
Termination condition: RSq changed by less than 0.001 at 4 terms
Importance: Girth
Number of terms at each degree of interaction: 1 2 (additive model)
GCV 14.20145  RSS 309.6832  GRSq 0.949137  RSq 0.9617962  CVRSq 0.9133168

Note: the cross-validation sd's below are standard deviations across folds

Cross validation:   nterms 3.00 sd 0.00    nvars 1.00 sd 0.00

     CVRSq    sd     MaxErr   sd
     0.913 0.045       14.3 8.71

varmod: method "gam" (mgcv package)    min.sd 0.408    iter.rsq 0.319

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)     4.078983    0.505543           12
s(Volume).1     0.000000          NA           NA
s(Volume).2     0.000000          NA           NA
s(Volume).3     0.000000          NA           NA
s(Volume).4     0.000000          NA           NA
s(Volume).5     0.000000          NA           NA
s(Volume).6     0.000000          NA           NA
s(Volume).7     0.000000          NA           NA
s(Volume).8     0.000000          NA           NA
s(Volume).9     2.061017          NA           NA

                              mean   smallest    largest      ratio
95% prediction interval   15.98932    5.66953   38.03631   6.708899

                                         68%    80%    90%    95% 
response values in prediction interval   77     90     97     97  
===skipping summary(mgcv::gam) etc.

===varmod.method gam: predict(earth.mod, interval="pint")
         fit       lwr      upr
1   9.913855  7.079089 12.74862
2  10.942195  7.845491 14.03890
3  11.627755  8.356425 14.89909
4  17.455018 12.699366 22.21067
5  18.140578 13.210301 23.07086
6  18.483359 13.465768 23.50095
7  19.168919 13.976702 24.36114
8  19.168919 13.976702 24.36114
9  19.511699 14.232169 24.79123
10 19.854479 14.487637 25.22132
11 20.197259 14.743104 25.65142
12 20.540040 14.998571 26.08151
13 20.540040 14.998571 26.08151
14 21.568380 15.764972 27.37179
15 22.596721 16.531374 28.66207
16 25.681742 18.830578 32.53291
17 25.681742 18.830578 32.53291
18 27.052863 19.852447 34.25328
19 28.423983 20.874315 35.97365
20 28.766764 21.129782 36.40374
21 30.080913 22.109191 38.05263
22 31.395063 23.088601 39.70152
23 33.366287 24.557714 42.17486
24 43.222408 31.903282 54.54153
25 45.193632 33.372396 57.01487
26 51.764379 38.269442 65.25932
27 53.078529 39.248851 66.90821
28 55.706828 41.207669 70.20599
29 56.363903 41.697374 71.03043
30 56.363903 41.697374 71.03043
31 73.447846 54.429692 92.46600
===varmod.method="x.gam"
===varmod.method x.gam: summary(earth.mod)
===summary(earth.mod)
Call: earth(formula=Volume~Girth, data=trees,
            trace=if(varmod.method%in%c("const","lm","power"))0.3els...),
            nfold=3, ncross=3, varmod.method=varmod.method)

              coefficients
(Intercept)      28.766764
h(13.8-Girth)    -3.427802
h(Girth-13.8)     6.570747

Selected 3 of 4 terms, and 1 of 1 predictors
Termination condition: RSq changed by less than 0.001 at 4 terms
Importance: Girth
Number of terms at each degree of interaction: 1 2 (additive model)
GCV 14.20145  RSS 309.6832  GRSq 0.949137  RSq 0.9617962  CVRSq 0.9133168

Note: the cross-validation sd's below are standard deviations across folds

Cross validation:   nterms 3.00 sd 0.00    nvars 1.00 sd 0.00

     CVRSq    sd     MaxErr   sd
     0.913 0.045       14.3 8.71

varmod: method "x.gam" (mgcv package)    min.sd 0.398    iter.rsq 0.391

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)     3.981114    0.452568           11
s(Volume).1     0.000000          NA           NA
s(Volume).2     0.000000          NA           NA
s(Volume).3     0.000000          NA           NA
s(Volume).4     0.000000          NA           NA
s(Volume).5     0.000000          NA           NA
s(Volume).6     0.000000          NA           NA
s(Volume).7     0.000000          NA           NA
s(Volume).8     0.000000          NA           NA
s(Volume).9     1.712690          NA           NA

                              mean   smallest    largest      ratio
95% prediction interval   15.60568   4.844286   31.59344   6.521795

                                         68%    80%    90%    95% 
response values in prediction interval   81     90     97     97  
===skipping summary(mgcv::gam) etc.

===varmod.method x.gam: predict(earth.mod, interval="pint")
         fit       lwr      upr
1   9.913855  7.491712 12.33600
2  10.942195  8.193843 13.69055
3  11.627755  8.661930 14.59358
4  17.455018 12.640674 22.26936
5  18.140578 13.108762 23.17240
6  18.483359 13.342806 23.62391
7  19.168919 13.810893 24.52694
8  19.168919 13.810893 24.52694
9  19.511699 14.044937 24.97846
10 19.854479 14.278981 25.42998
11 20.197259 14.513024 25.88149
12 20.540040 14.747068 26.33301
13 20.540040 14.747068 26.33301
14 21.568380 15.449199 27.68756
15 22.596721 16.151331 29.04211
16 25.681742 18.257725 33.10576
17 25.681742 18.257725 33.10576
18 27.052863 19.193900 34.91183
19 28.423983 20.130075 36.71789
20 28.766764 20.364118 37.16941
21 30.080913 21.460795 38.70103
22 31.395063 22.557472 40.23265
23 33.366287 24.202487 42.53009
24 43.222408 32.427562 54.01725
25 45.193632 34.072577 56.31469
26 51.764379 39.555960 63.97280
27 53.078529 40.652637 65.50442
28 55.706828 42.845990 68.56767
29 56.363903 43.394329 69.33348
30 56.363903 43.394329 69.33348
31 73.447846 57.651125 89.24457
> par(old.par)
> 
> source("test.epilog.R")
