> # test.varmod.R
> 
> 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)
> 
> multifigure("test predict.earth with pints", 2, 2)
===test predict.earth with pints
> 
> set.seed(2)
> earth.trees <- earth(Volume~Girth, data=trees, nfold=3, ncross=3, varmod.method="earth")
> 
> old.environment <- attr(earth.trees$varmod, ".Environment")
> stopifnot(is.environment(old.environment))
> # following necessary else print.default prints a different default environment hex address each time
> attr(earth.trees$varmod, ".Environment") <- NULL
> printh("print.default(earth.trees$varmod)")
===print.default(earth.trees$varmod)
> print.default(earth.trees$varmod)
$call
varmod(parent = parent, method = method, exponent = exponent, 
    conv = conv, clamp = clamp, minspan = minspan, trace = trace, 
    parent.x = x, parent.y = y, model.var = model.var, degree = ..1)

$parent
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.9451125

$method
[1] "earth"

$exponent
[1] 1

$lambda
[1] 1

$rmethod
[1] "hc12"

$converged
[1] TRUE

$iters
[1] 1

$residmod
Selected 2 of 4 terms, and 1 of 1 predictors
Termination condition: RSq changed by less than 0.001 at 4 terms
Importance: RHS
Number of terms at each degree of interaction: 1 1 (additive model)
GCV 3.4882    RSS 88.2177    GRSq 0.2106366    RSq 0.3123768

$min.sd
[1] 0.3800377

$model.var
              [,1]
 [1,] 0.4237080282
 [2,] 0.9796261257
 [3,] 0.3353434600
 [4,] 0.0014690631
 [5,] 0.1956131647
 [6,] 0.1656436269
 [7,] 0.4260785739
 [8,] 0.1660133056
 [9,] 0.1308601982
[10,] 0.4700957111
[11,] 0.0001419694
[12,] 0.1660851055
[13,] 0.5168044702
[14,] 0.2612776068
[15,] 0.2310163116
[16,] 0.5379486915
[17,] 1.6295066790
[18,] 0.4264025111
[19,] 0.0090402313
[20,] 0.0305532477
[21,] 1.2162723789
[22,] 1.7221628721
[23,] 1.5223938207
[24,] 0.1820720480
[25,] 0.2350058950
[26,] 0.3501834397
[27,] 2.7639280510
[28,] 2.5209256938
[29,] 3.8249680345
[30,] 0.2403475621
[31,] 9.1817721531

$abs.resids
         Volume
 [1,] 0.7975990
 [2,] 1.2429614
 [3,] 1.7509813
 [4,] 1.1484398
 [5,] 0.8405929
 [6,] 1.3783014
 [7,] 3.9078464
 [8,] 1.1225815
 [9,] 3.3502387
[10,] 0.6873876
[11,] 4.3101763
[12,] 0.6411992
[13,] 1.1719193
[14,] 0.5870355
[15,] 3.7951596
[16,] 3.8635479
[17,] 8.9363279
[18,] 0.7564978
[19,] 3.0376093
[20,] 4.3294863
[21,] 5.0370789
[22,] 1.3549887
[23,] 3.4526918
[24,] 5.3635770
[25,] 2.8631420
[26,] 4.0604341
[27,] 3.3523508
[28,] 3.3190359
[29,] 5.8240248
[30,] 6.0695913
[31,] 5.6473173

$parent.x
      Girth
 [1,]   8.3
 [2,]   8.6
 [3,]   8.8
 [4,]  10.5
 [5,]  10.7
 [6,]  10.8
 [7,]  11.0
 [8,]  11.0
 [9,]  11.1
[10,]  11.2
[11,]  11.3
[12,]  11.4
[13,]  11.4
[14,]  11.7
[15,]  12.0
[16,]  12.9
[17,]  12.9
[18,]  13.3
[19,]  13.7
[20,]  13.8
[21,]  14.0
[22,]  14.2
[23,]  14.5
[24,]  16.0
[25,]  16.3
[26,]  17.3
[27,]  17.5
[28,]  17.9
[29,]  18.0
[30,]  18.0
[31,]  20.6

$parent.y
      Volume
 [1,]   10.3
 [2,]   10.3
 [3,]   10.2
 [4,]   16.4
 [5,]   18.8
 [6,]   19.7
 [7,]   15.6
 [8,]   18.2
 [9,]   22.6
[10,]   19.9
[11,]   24.2
[12,]   21.0
[13,]   21.4
[14,]   21.3
[15,]   19.1
[16,]   22.2
[17,]   33.8
[18,]   27.4
[19,]   25.7
[20,]   24.9
[21,]   34.5
[22,]   31.7
[23,]   36.3
[24,]   38.3
[25,]   42.6
[26,]   55.4
[27,]   55.7
[28,]   58.3
[29,]   51.5
[30,]   51.0
[31,]   77.0

$iter.rsq
[1] 0.3123768

$iter.stderr
     (Intercept) bxh(25.6817-RHS) 
      0.40081837       0.06710405 

attr(,"class")
[1] "varmod"
> attr(earth.trees$varmod, ".Environment") <- old.environment
> 
> printh("summary(earth.trees)")
===summary(earth.trees)
> print(summary(earth.trees))
Call: earth(formula=Volume~Girth, data=trees, nfold=3, ncross=3,
            varmod.method="earth")

              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.9451125

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.945 0.017       9.83 7.52

varmod: method "earth"    min.sd 0.38    iter.rsq 0.312

stddev of predictions:
                  coefficients iter.stderr iter.stderr%
(Intercept)          4.9378880    0.502351           10
h(25.6817-Volume)   -0.3052613   0.0841025           28

                              mean   smallest    largest      ratio
95% prediction interval   14.92951    1.48972   19.35617   12.99315

                                         68%    80%    90%    95% 
response values in prediction interval   71     90     97     97  
> 
> # level arg not allowed with interval="se"
> expect.err(try(predict(earth.trees, interval="se", level=.8)))
Error : predict.varmod: the level argument is not allowed with type="se"
Got error as expected from try(predict(earth.trees, interval = "se", level = 0.8))
> 
> printh("predict(earth.trees, interval=\"se\")")
===predict(earth.trees, interval="se")
> stderrs <- predict(earth.trees, interval="se")
> print(stderrs)
 [1] 0.3800377 0.4384742 0.6477493 2.4265873 2.6358624 2.7404999 2.9497750
 [8] 2.9497750 3.0544125 3.1590500 3.2636876 3.3683251 3.3683251 3.6822377
[15] 3.9961503 4.9378880 4.9378880 4.9378880 4.9378880 4.9378880 4.9378880
[22] 4.9378880 4.9378880 4.9378880 4.9378880 4.9378880 4.9378880 4.9378880
[29] 4.9378880 4.9378880 4.9378880
> 
> # level arg not allowed with interval="abs.residual"
> expect.err(try(predict(earth.trees, interval="abs.res", level=.8)))
Error : predict.varmod: the level argument is not allowed with type="abs.residual"
Got error as expected from try(predict(earth.trees, interval = "abs.res", level = 0.8))
> 
> printh("predict(earth.trees, interval=\"abs.residual\")")
===predict(earth.trees, interval="abs.residual")
> stderrs <- predict(earth.trees, interval="abs.residual")
> print(stderrs)
 [1] 0.3032262 0.3498518 0.5168292 1.9361365 2.1031139 2.1866026 2.3535799
 [8] 2.3535799 2.4370686 2.5205572 2.6040459 2.6875346 2.6875346 2.9380006
[15] 3.1884666 3.9398646 3.9398646 3.9398646 3.9398646 3.9398646 3.9398646
[22] 3.9398646 3.9398646 3.9398646 3.9398646 3.9398646 3.9398646 3.9398646
[29] 3.9398646 3.9398646 3.9398646
> 
> expect.err(try(predict(earth.trees, newdata=trees, interval="cint")))
Error : predict.varmod: newdata is not allowed with interval="cint"
Got error as expected from try(predict(earth.trees, newdata = trees, interval = "cint"))
> 
> printh("predict(earth.trees, interval=\"cint\")")
===predict(earth.trees, interval="cint")
> cints <- predict(earth.trees, interval="cint")
> print(cints)
         fit       lwr      upr
1   9.913855  8.638058 11.18965
2  10.942195  9.002300 12.88209
3  11.627755 10.492763 12.76275
4  17.455018 17.379896 17.53014
5  18.140578 17.273722 19.00743
6  18.483359 17.685666 19.28105
7  19.168919 17.889558 20.44828
8  19.168919 18.370337 19.96750
9  19.511699 18.802690 20.22071
10 19.854479 18.510659 21.19830
11 20.197259 20.173906 20.22061
12 20.540040 19.741285 21.33879
13 20.540040 19.131039 21.94904
14 21.568380 20.566538 22.57022
15 22.596721 21.654680 23.53876
16 25.681742 24.244207 27.11928
17 25.681742 23.179806 28.18368
18 27.052863 25.773016 28.33271
19 28.423983 28.237630 28.61034
20 28.766764 28.424172 29.10936
21 30.080913 27.919372 32.24245
22 31.395063 28.822978 33.96715
23 33.366287 30.947979 35.78459
24 43.222408 42.386093 44.05872
25 45.193632 44.243492 46.14377
26 51.764379 50.604545 52.92421
27 53.078529 49.820076 56.33698
28 55.706828 52.594910 58.81875
29 56.363903 52.530698 60.19711
30 56.363903 55.403025 57.32478
31 73.447846 67.508873 79.38682
> 
> printh("predict(earth.trees, interval=\"pin\", level=.80)")
===predict(earth.trees, interval="pin", level=.80)
> news <- predict(earth.trees, interval="pin", level=.80)
> print(news)
         fit       lwr      upr
1   9.913855  9.426817 10.40089
2  10.942195 10.380268 11.50412
3  11.627755 10.797631 12.45788
4  17.455018 14.345221 20.56481
5  18.140578 14.762585 21.51857
6  18.483359 14.971267 21.99545
7  19.168919 15.388630 22.94921
8  19.168919 15.388630 22.94921
9  19.511699 15.597312 23.42609
10 19.854479 15.805994 23.90296
11 20.197259 16.014676 24.37984
12 20.540040 16.223357 24.85672
13 20.540040 16.223357 24.85672
14 21.568380 16.849403 26.28736
15 22.596721 17.475448 27.71799
16 25.681742 19.353584 32.00990
17 25.681742 19.353584 32.00990
18 27.052863 20.724705 33.38102
19 28.423983 22.095825 34.75214
20 28.766764 22.438605 35.09492
21 30.080913 23.752755 36.40907
22 31.395063 25.066904 37.72322
23 33.366287 27.038129 39.69444
24 43.222408 36.894250 49.55057
25 45.193632 38.865474 51.52179
26 51.764379 45.436221 58.09254
27 53.078529 46.750371 59.40669
28 55.706828 49.378670 62.03499
29 56.363903 50.035745 62.69206
30 56.363903 50.035745 62.69206
31 73.447846 67.119688 79.77600
> expect.err(try(predict(earth.trees, interval="none", level=.80)), "predict.earth: level=0.8 was specified but interval=\"none\"")
Error : predict.earth: level=0.8 was specified but interval="none"
Got error as expected from try(predict(earth.trees, interval = "none", level = 0.8))
> expect.err(try(predict(earth.trees, interval="pin", type="class")), "predict.earth: the interval argument is not allowed with type=\"class\"")
Error : predict.earth: the interval argument is not allowed with type="class"
Got error as expected from try(predict(earth.trees, interval = "pin", type = "class"))
> expect.err(try(predict(earth.trees, interval="pin", type="cl")), "predict.earth: the interval argument is not allowed with type=\"class\"")
Error : predict.earth: the interval argument is not allowed with type="class"
Got error as expected from try(predict(earth.trees, interval = "pin", type = "cl"))
> expect.err(try(predict(earth.trees, interval="pin", type="ter")), "predict.earth: the interval argument is not allowed with type=\"terms\"")
Error : predict.earth: the interval argument is not allowed with type="terms"
Got error as expected from try(predict(earth.trees, interval = "pin", type = "ter"))
> 
> printh("print.default(earth.trees$varmod$residmod)")
===print.default(earth.trees$varmod$residmod)
> # have to modify earth.trees because terms field stores the environment
> # as a hex address which messes up the diffs
> earth.trees$varmod$residmod$terms <- NULL
> print.default(earth.trees$varmod$residmod)
$rss
[1] 88.2177

$rsq
[1] 0.3123768

$gcv
[1] 3.4882

$grsq
[1] 0.2106366

$bx
      (Intercept) h(25.6817-RHS)
 [1,]           1      15.767888
 [2,]           1      14.739547
 [3,]           1      14.053987
 [4,]           1       8.226724
 [5,]           1       7.541164
 [6,]           1       7.198383
 [7,]           1       6.512823
 [8,]           1       6.512823
 [9,]           1       6.170043
[10,]           1       5.827263
[11,]           1       5.484483
[12,]           1       5.141702
[13,]           1       5.141702
[14,]           1       4.113362
[15,]           1       3.085021
[16,]           1       0.000000
[17,]           1       0.000000
[18,]           1       0.000000
[19,]           1       0.000000
[20,]           1       0.000000
[21,]           1       0.000000
[22,]           1       0.000000
[23,]           1       0.000000
[24,]           1       0.000000
[25,]           1       0.000000
[26,]           1       0.000000
[27,]           1       0.000000
[28,]           1       0.000000
[29,]           1       0.000000
[30,]           1       0.000000
[31,]           1       0.000000

$dirs
               RHS
(Intercept)      0
h(RHS-25.6817)   1
h(25.6817-RHS)  -1
h(RHS-43.2224)   1

$cuts
                    RHS
(Intercept)     0.00000
h(RHS-25.6817) 25.68174
h(25.6817-RHS) 25.68174
h(RHS-43.2224) 43.22241

$selected.terms
[1] 1 3

$prune.terms
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    1    3    0    0
[3,]    1    3    4    0
[4,]    1    2    3    4

$fitted.values
      abs.resids
 [1,] 0.09938581
 [2,] 0.34985182
 [3,] 0.51682916
 [4,] 1.93613655
 [5,] 2.10311389
 [6,] 2.18660256
 [7,] 2.35357990
 [8,] 2.35357990
 [9,] 2.43706857
[10,] 2.52055724
[11,] 2.60404591
[12,] 2.68753458
[13,] 2.68753458
[14,] 2.93800059
[15,] 3.18846660
[16,] 3.93986463
[17,] 3.93986463
[18,] 3.93986463
[19,] 3.93986463
[20,] 3.93986463
[21,] 3.93986463
[22,] 3.93986463
[23,] 3.93986463
[24,] 3.93986463
[25,] 3.93986463
[26,] 3.93986463
[27,] 3.93986463
[28,] 3.93986463
[29,] 3.93986463
[30,] 3.93986463
[31,] 3.93986463

$residuals
       abs.resids
 [1,]  0.69821318
 [2,]  0.89310955
 [3,]  1.23415211
 [4,] -0.78769675
 [5,] -1.26252101
 [6,] -0.80830120
 [7,]  1.55426655
 [8,] -1.23099838
 [9,]  0.91317016
[10,] -1.83316966
[11,]  1.70613040
[12,] -2.04633535
[13,] -1.51561531
[14,] -2.35096508
[15,]  0.60669298
[16,] -0.07631673
[17,]  4.99646323
[18,] -3.18336680
[19,] -0.90225532
[20,]  0.38962172
[21,]  1.09721426
[22,] -2.58487588
[23,] -0.48717284
[24,]  1.42371239
[25,] -1.07672265
[26,]  0.12056947
[27,] -0.58751387
[28,] -0.62082874
[29,]  1.88416021
[30,]  2.12972668
[31,]  1.70745270

$coefficients
               abs.resids
(Intercept)     3.9398646
h(25.6817-RHS) -0.2435633

$rss.per.response
[1] 88.2177

$rsq.per.response
[1] 0.3123768

$gcv.per.response
[1] 3.4882

$grsq.per.response
[1] 0.2106366

$rss.per.subset
[1] 128.29366  88.21770  81.31999  81.22837

$gcv.per.subset
[1] 4.419004 3.488200 3.729171 4.371666

$leverages
 [1] 0.24689453 0.21180022 0.19014329 0.06223841 0.05380006 0.05010267
 [7] 0.04375147 0.04375147 0.04109766 0.03879170 0.03683360 0.03522336
[13] 0.03522336 0.03247977 0.03286691 0.05281259 0.05281259 0.05281259
[19] 0.05281259 0.05281259 0.05281259 0.05281259 0.05281259 0.05281259
[25] 0.05281259 0.05281259 0.05281259 0.05281259 0.05281259 0.05281259
[31] 0.05281259

$pmethod
[1] "backward"

$nprune
NULL

$penalty
[1] 2

$nk
[1] 21

$thresh
[1] 0.001

$termcond
[1] 4

$weights
NULL

$call
earth(formula = abs.resids ~ ., data = data, weights = weights, 
    keepxy = TRUE, trace = trace, degree = ..1, minspan = minspan)

$namesx.org
[1] "RHS"

$namesx
[1] "RHS"

$xlevels
named list()

$data
   abs.resids       RHS
1   0.7975990  9.913855
2   1.2429614 10.942195
3   1.7509813 11.627755
4   1.1484398 17.455018
5   0.8405929 18.140578
6   1.3783014 18.483359
7   3.9078464 19.168919
8   1.1225815 19.168919
9   3.3502387 19.511699
10  0.6873876 19.854479
11  4.3101763 20.197259
12  0.6411992 20.540040
13  1.1719193 20.540040
14  0.5870355 21.568380
15  3.7951596 22.596721
16  3.8635479 25.681742
17  8.9363279 25.681742
18  0.7564978 27.052863
19  3.0376093 28.423983
20  4.3294863 28.766764
21  5.0370789 30.080913
22  1.3549887 31.395063
23  3.4526918 33.366287
24  5.3635770 43.222408
25  2.8631420 45.193632
26  4.0604341 51.764379
27  3.3523508 53.078529
28  3.3190359 55.706828
29  5.8240248 56.363903
30  6.0695913 56.363903
31  5.6473173 73.447846

$y
      abs.resids
 [1,]  0.7975990
 [2,]  1.2429614
 [3,]  1.7509813
 [4,]  1.1484398
 [5,]  0.8405929
 [6,]  1.3783014
 [7,]  3.9078464
 [8,]  1.1225815
 [9,]  3.3502387
[10,]  0.6873876
[11,]  4.3101763
[12,]  0.6411992
[13,]  1.1719193
[14,]  0.5870355
[15,]  3.7951596
[16,]  3.8635479
[17,]  8.9363279
[18,]  0.7564978
[19,]  3.0376093
[20,]  4.3294863
[21,]  5.0370789
[22,]  1.3549887
[23,]  3.4526918
[24,]  5.3635770
[25,]  2.8631420
[26,]  4.0604341
[27,]  3.3523508
[28,]  3.3190359
[29,]  5.8240248
[30,]  6.0695913
[31,]  5.6473173

attr(,"class")
[1] "earth"
> # prevent mistakes later where we try to use a modified earth.trees
> remove(earth.trees)
> 
> multifigure("test example for varmod help page", 2, 2)
===test example for varmod help page
> 
> data(ozone1)
> set.seed(1) # optional, for cross validation reproducibility
> 
> a <- earth(O3~temp, data=ozone1, nfold=10, ncross=3, varmod.method="earth")
> 
> print(summary(a)) # note additional info on the variance model
Call: earth(formula=O3~temp, data=ozone1, nfold=10, ncross=3,
            varmod.method="earth")

            coefficients
(Intercept)    7.5768304
h(58-temp)    -0.1572930
h(temp-58)     0.6090229

Selected 3 of 4 terms, and 1 of 1 predictors
Termination condition: RSq changed by less than 0.001 at 4 terms
Importance: temp
Number of terms at each degree of interaction: 1 2 (additive model)
GCV 22.43818  RSS 7181.917  GRSq 0.6514492  RSq 0.6598731  CVRSq 0.6168999

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

Cross validation:   nterms 3.37 sd 0.89    nvars 1.00 sd 0.00

     CVRSq    sd     MaxErr   sd
     0.617 0.107       17.2 11.1

varmod: method "earth"    min.sd 0.46    iter.rsq 0.154

stddev of predictions:
              coefficients iter.stderr iter.stderr%
(Intercept)      6.4704573    0.307898            5
h(16.1032-O3)   -0.3329966   0.0430647           13

                              mean   smallest    largest      ratio
95% prediction interval   18.01673   7.458628   25.36373   3.400589

                                         68%    80%    90%    95% 
response values in prediction interval   70     82     94     97  
> 
> old.mfrow <- par(mfrow=c(2,2))
> 
> # the variance model assumes residuals are symmetric, which is not
> # quite true in this example, so the lower band is a bit too big
> plotmo(a, do.par=FALSE, col.response=1, level=.95, main="earth model: O3~temp")
> 
> plot(a, which=1)            # model selection plot, same as ever
> plot(a, which=3, level=.95) # residual plot: note 95% pred and darker conf intervals
> plot(a, which=3, level=.95, standardize=TRUE) # standardize resids are approx homoscedastic
> 
> par(par=old.mfrow)
> 
> plot(a$varmod)              # plot the embedded variance model (this calls plot.varmod)
> 
> multifigure("test example for plot.varmod help page", 2, 2)
===test example for plot.varmod help page
> 
> # multivariate example (for univariate, see the example on the varmod help page)
> 
> data(ozone1)
> set.seed(1) # optional, for cross validation reproducibility
> 
> mod.temp.vh.doy <- earth(O3~temp+vh+vis+doy, data=ozone1, nfold=5, ncross=3, varmod.method="x.earth")
> 
> print(summary(mod.temp.vh.doy))  # note additional info on the variance model
Call: earth(formula=O3~temp+vh+vis+doy, data=ozone1, nfold=5, ncross=3,
            varmod.method="x.earth")

            coefficients
(Intercept)    9.8874670
h(temp-58)     0.4016529
h(vh-5530)     0.0210146
h(vis-17)     -0.0324305
h(vis-200)     0.0474193
h(106-doy)    -0.0972229
h(doy-151)    -0.0285859

Selected 7 of 15 terms, and 4 of 4 predictors
Termination condition: Reached nk 21
Importance: temp, vis, doy, vh
Number of terms at each degree of interaction: 1 6 (additive model)
GCV 17.87166  RSS 5442.137  GRSq 0.7223847  RSq 0.742267  CVRSq 0.6740982

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

Cross validation:   nterms 8.87 sd 1.77    nvars 4.00 sd 0.00

     CVRSq    sd     MaxErr   sd
     0.674 0.056       17.8 14.4

varmod: method "x.earth"    min.sd 0.418    iter.rsq 0.259

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)    8.4358492     1.31775           16
h(vh-5690)     0.0138052   0.0026229           19
h(vis-70)     -0.0222822    0.005732           26
h(vis-150)     0.0215700   0.0085761           40
h(doy-119)    -0.0218882   0.0058761           27
h(302-doy)    -0.0215648   0.0061368           28

                              mean   smallest    largest      ratio
95% prediction interval   16.38722   3.600676   31.47704   8.741982

                                         68%    80%    90%    95% 
response values in prediction interval   73     82     93     96  
> 
> plot(mod.temp.vh.doy, level=.95) # note 95% pred and darker conf intervals in resids plot
> 
> plot(mod.temp.vh.doy$varmod)     # plot the variance model (this calls plot.varmod)
> 
> plot(mod.temp.vh.doy, versus="", level=.9, caption="plot.earth versus=\"\"")
> plot(mod.temp.vh.doy, versus="v", level=.9, caption="plot.earth versus=\"v\" and versus=\"temp\"", do.par=2)
> plot(mod.temp.vh.doy, versus="temp", level=.9, caption="", main="temp on same page")
> 
> # plot.earth will silently not plots it cannot plot below, so 1:9 becomes c(3,5,6)
> plot(mod.temp.vh.doy, which=1:9, versus="v", info=T, caption='which=c(3,5) versus="v" info=T')
> 
> # versus="b:"
> 
> plot(mod.temp.vh.doy, versus="b:", level=.9,
+      caption="plot.earth versus=\"b:\"")
> 
> # versus="b:" and versus=1:4 with info
> plot(mod.temp.vh.doy, versus="b:", level=.8, info=TRUE,
+      caption="plot.earth versus=\"b:\" with info")
> 
> multifigure("versus=1:4", 3, 3)
===versus=1:4
> 
> plot(mod.temp.vh.doy, versus=1, caption="", do.par=FALSE, which=3)
> do.caption()
> plot(mod.temp.vh.doy, versus=2, caption="", do.par=FALSE)
> plot(mod.temp.vh.doy, versus=3, caption="", do.par=FALSE)
> 
> plot(mod.temp.vh.doy, versus=1, info=TRUE, caption="", do.par=FALSE, which=3)
> plot(mod.temp.vh.doy, versus=2, info=TRUE, caption="", do.par=FALSE)
> plot(mod.temp.vh.doy, versus=3, info=TRUE, caption="", do.par=FALSE)
> 
> plot(mod.temp.vh.doy, versus=1, info=TRUE, caption="", do.par=FALSE, level=.8, which=3)
> plot(mod.temp.vh.doy, versus=2, info=TRUE, caption="", do.par=FALSE, level=.8)
> plot(mod.temp.vh.doy, versus=3, info=TRUE, caption="", do.par=FALSE, level=.8)
> 
> expect.err(try(plot(mod.temp.vh.doy, versus=9)))
Error : versus must be an integer or a string:
  1    fitted (default)
  2    observation numbers
  3    response
  4    leverages
  ""   predictors
  "b:" basis functions
Got error as expected from try(plot(mod.temp.vh.doy, versus = 9))
> expect.err(try(plot(mod.temp.vh.doy, versus=1.2)))
Error : versus must be an integer or a string:
  1    fitted (default)
  2    observation numbers
  3    response
  4    leverages
  ""   predictors
  "b:" basis functions
Got error as expected from try(plot(mod.temp.vh.doy, versus = 1.2))
> expect.err(try(plot(mod.temp.vh.doy, versus=2:3)))
Error : illegal 'versus' argument (length of 'versus' must be 1 when 'versus' is numeric)
Got error as expected from try(plot(mod.temp.vh.doy, versus = 2:3))
> 
> # versus="b:doy"
> 
> plot(mod.temp.vh.doy, versus="b:doy", level=.9, caption="plot.earth versus=\"b:doy\"")
> 
> # test warnings from plotres about which
> options(warn=1) # print warnings as they occur
> plot(mod.temp.vh.doy, which=1, versus="b:doy")
Warning: which=1 is now empty because plots were removed because versus=b:doy
Warning: plotres: nothing to plot
> options(warn=2)
> 
> multifigure("test example in (very old) earth vignette", 2, 2)
===test example in (very old) earth vignette
> 
> data(ozone1)
> x <- ozone1$temp
> y <- ozone1$O3
> 
> set.seed(1) # optional, for cross validation reproducibility
> earth.mod <- earth(y~x, nfold=10, ncross=3, varmod.method="earth", trace=.1)
Full model GRSq 0.651 RSq 0.660, starting cross validation
CV fold  1.1  CVRSq  0.507   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  1.2  CVRSq  0.447   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  1.3  CVRSq  0.684   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  1.4  CVRSq  0.714   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  1.5  CVRSq  0.658   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  1.6  CVRSq  0.725   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  1.7  CVRSq  0.565   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  1.8  CVRSq  0.539   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  1.9  CVRSq  0.450   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  1.10 CVRSq  0.763   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  2.1  CVRSq  0.654   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  2.2  CVRSq  0.462   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  2.3  CVRSq  0.550   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  2.4  CVRSq  0.510   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  2.5  CVRSq  0.628   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  2.6  CVRSq  0.742   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  2.7  CVRSq  0.692   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  2.8  CVRSq  0.684   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  2.9  CVRSq  0.530   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  2.10 CVRSq  0.710   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  3.1  CVRSq  0.352   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  3.2  CVRSq  0.652   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  3.3  CVRSq  0.574   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  3.4  CVRSq  0.789   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  3.5  CVRSq  0.598   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  3.6  CVRSq  0.732   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  3.7  CVRSq  0.658   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  3.8  CVRSq  0.632   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  3.9  CVRSq  0.649   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV fold  3.10 CVRSq  0.657   n.oof 297 10%   n.infold.nz 297 100%   n.oof.nz  33 100%
CV all        CVRSq  0.617                   n.infold.nz 330 100%

> predict <- predict(earth.mod, interval="pint")
> 
> order <- order(x)
> x <- x[order]
> y <- y[order]
> predict <- predict[order,]
> 
> inconf <- y >= predict$lwr & y <= predict$upr
> 
> plot(x, y, pch=20, col=ifelse(inconf, 1, 2), main=sprint(
+     "Prediction intervals\n%.0f%% of the points are in the estimated band",
+     100 * sum(inconf) / length(y)))
> do.caption()
> 
> lines(x, predict$fit)
> lines(x, predict$lwr, lty=2)
> lines(x, predict$upr, lty=2)
> 
> # Plot the Residuals vs Fitted graph
> plot(earth.mod, which=3, level=.95)
> 
> # Plot the embedded residual model
> plot(earth.mod$varmod, do.par=F, which=1:2)
> 
> cat('head(residuals(earth.mod))\n')
head(residuals(earth.mod))
> print(head(residuals(earth.mod)))
           y
1 -1.7455563
2 -0.5320213
3 -1.9476584
4  2.0409087
5 -1.5320213
6 -3.1049514
> cat('head(residuals(earth.mod, type="standardize"))\n')
head(residuals(earth.mod, type="standardize"))
> print(head(residuals(earth.mod, type="standardize")))
           y
1 -0.6542917
2 -0.1810585
3 -0.5709062
4  0.8525427
5 -0.5213804
6 -0.8966962
> 
> multifigure("plot.earth varmod options", 2, 2)
===plot.earth varmod options
> 
> plot(earth.mod, which=3, level=.95, level.shade=0, main="plot.earth varmod options")
> do.caption()
> plot(earth.mod, which=3, level.shade="orange", level.shade2="darkgray", level=.99)
> plot(earth.mod, which=3, level=.95, level.shade=0, level.shade2="mistyrose4")
> 
> multifigure("plot.earth delever and standardize", 2, 2)
===plot.earth delever and standardize
> 
> set.seed(4)
> earth.mod1 <- earth(O3~temp, data=ozone1, nfold=5, ncross=3, varmod.method="lm", keepxy=T, trace=.1)
Full model GRSq 0.651 RSq 0.660, starting cross validation
CV fold  1.1  CVRSq  0.711   n.oof 264 20%   n.infold.nz 264 100%   n.oof.nz  66 100%
CV fold  1.2  CVRSq  0.642   n.oof 264 20%   n.infold.nz 264 100%   n.oof.nz  66 100%
CV fold  1.3  CVRSq  0.654   n.oof 264 20%   n.infold.nz 264 100%   n.oof.nz  66 100%
CV fold  1.4  CVRSq  0.571   n.oof 264 20%   n.infold.nz 264 100%   n.oof.nz  66 100%
CV fold  1.5  CVRSq  0.604   n.oof 264 20%   n.infold.nz 264 100%   n.oof.nz  66 100%
CV fold  2.1  CVRSq  0.646   n.oof 264 20%   n.infold.nz 264 100%   n.oof.nz  66 100%
CV fold  2.2  CVRSq  0.746   n.oof 264 20%   n.infold.nz 264 100%   n.oof.nz  66 100%
CV fold  2.3  CVRSq  0.579   n.oof 264 20%   n.infold.nz 264 100%   n.oof.nz  66 100%
CV fold  2.4  CVRSq  0.550   n.oof 264 20%   n.infold.nz 264 100%   n.oof.nz  66 100%
CV fold  2.5  CVRSq  0.602   n.oof 264 20%   n.infold.nz 264 100%   n.oof.nz  66 100%
CV fold  3.1  CVRSq  0.716   n.oof 264 20%   n.infold.nz 264 100%   n.oof.nz  66 100%
CV fold  3.2  CVRSq  0.717   n.oof 264 20%   n.infold.nz 264 100%   n.oof.nz  66 100%
CV fold  3.3  CVRSq  0.511   n.oof 264 20%   n.infold.nz 264 100%   n.oof.nz  66 100%
CV fold  3.4  CVRSq  0.550   n.oof 264 20%   n.infold.nz 264 100%   n.oof.nz  66 100%
CV fold  3.5  CVRSq  0.648   n.oof 264 20%   n.infold.nz 264 100%   n.oof.nz  66 100%
CV all        CVRSq  0.630                   n.infold.nz 330 100%

> plot(earth.mod1, which=3, ylim=c(-16,20), info=TRUE, level=.95)
> do.caption()
> plot(earth.mod1, which=3, ylim=c(-16,20), delever=TRUE, level=.95)
> plot(earth.mod1, which=3, standardize=TRUE, info=TRUE,    level=.95)
> # the standardize and delever arguments cannot both be set
> expect.err(try(plot(earth.mod1, which=3, standardize=TRUE, delever=TRUE, level=.95)))
Error : the standardize and delever arguments cannot both be set
Got error as expected from try(plot(earth.mod1, which = 3, standardize = TRUE, delever = TRUE,     level = 0.95))
> 
> multifigure("plot.earth which=5 and which=6", 2, 3)
===plot.earth which=5 and which=6
> plot(earth.mod1, which=5, info=T,            main="which=5, info=T")
> plot(earth.mod1, which=5, standardize=T, info=T, main="which=5, standardize=T, info=T")
> plot(earth.mod1, which=5, standardize=T,         main="which=5, standardize=T")
> do.caption()
> plot(earth.mod1, which=6, info=T,            main="which=6, info=T")
> plot(earth.mod1, which=6, standardize=T, info=T, main="which=6, standardize=T, info=T")
> plot(earth.mod1, which=6, standardize=T,         main="which=6, standardize=T")
> 
> multifigure("plot.earth which=7", 2, 3)
===plot.earth which=7
> plot(earth.mod1, which=7, info=T,            main="which=7, info=T")
> plot(earth.mod1, which=7, standardize=T, info=T, main="which=7, standardize=T, info=T")
> plot(earth.mod1, which=7, standardize=T,         main="which=7, standardize=T")
> do.caption()
> 
> multifigure("plot.earth which=8 and which=9", 2, 3)
===plot.earth which=8 and which=9
> plot(earth.mod1, which=8, info=T,            main="which=8, info=T")
> plot(earth.mod1, which=8, standardize=T, info=T, main="which=8, standardize=T, info=T")
> plot(earth.mod1, which=8, standardize=T,         main="which=8, standardize=T")
> do.caption()
> plot(earth.mod1, which=9, info=T,            main="which=9, info=T")
> plot(earth.mod1, which=9, standardize=T, info=T, main="which=9, standardize=T, info=T")
> plot(earth.mod1, which=9, standardize=T,         main="which=9, standardize=T")
> 
> multifigure("plot.earth versus=4, which=3 and which=5", 2, 3)
===plot.earth versus=4, which=3 and which=5
> plot(earth.mod1, versus=4, which=3,                    main="versus=4, which=3")
> plot(earth.mod1, versus=4, which=3, standardize=T, info=T, main="versus=4, which=3, standardize=T, info=T")
> plot(earth.mod1, versus=4, which=3, standardize=T,         main="versus=4, which=3, standardize=T")
> do.caption()
> plot(earth.mod1, versus=4, which=5,                    main="versus=4, which=5")
> plot(earth.mod1, versus=4, which=5, standardize=T, info=T, main="versus=4, which=5, standardize=T, info=T")
> plot(earth.mod1, versus=4, which=5, standardize=T,         main="versus=4, which=5, standardize=T")
> 
> cat("summary(earth.mod1, newdata=ozone1)\n")
summary(earth.mod1, newdata=ozone1)
> print(summary(earth.mod1, newdata=ozone1))
RSq 0.660 on newdata (330 cases)

                                 68%    80%    90%    95% 
newdata in prediction interval   71     83     93     97  
> 
> cat("summary(earth.mod1, newdata=ozone1[1:100,]:)\n")
summary(earth.mod1, newdata=ozone1[1:100,]:)
> print(summary(earth.mod1, newdata=ozone1[1:100,]))
RSq 0.199 on newdata (100 cases)

                                 68%    80%    90%    95% 
newdata in prediction interval   78     85     92     95  
> 
> # earth.default
> O3 <- ozone1$O3
> temp <- ozone1$temp
> set.seed(4)
> earth.default <- earth(temp, O3, nfold=5, ncross=3, varmod.method="lm")
> cat("summary(earth.default)\n")
summary(earth.default)
> print(summary(earth.default))
Call: earth(x=temp, y=O3, nfold=5, ncross=3, varmod.method="lm")

            coefficients
(Intercept)    7.5768304
h(58-temp)    -0.1572930
h(temp-58)     0.6090229

Selected 3 of 4 terms, and 1 of 1 predictors
Termination condition: RSq changed by less than 0.001 at 4 terms
Importance: temp
Number of terms at each degree of interaction: 1 2 (additive model)
GCV 22.43818  RSS 7181.917  GRSq 0.6514492  RSq 0.6598731  CVRSq 0.6297344

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

Cross validation:   nterms 3.47 sd 0.92    nvars 1.00 sd 0.00

     CVRSq    sd     MaxErr   sd
      0.63 0.071       17.7 13.8

varmod: method "lm"    min.sd 0.472    iter.rsq 0.147

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)     1.575154    0.336099           21
O3              0.267135   0.0354983           13

                              mean   smallest    largest      ratio
95% prediction interval   18.50547   8.673159   36.42941   4.200247

                                         68%    80%    90%    95% 
response values in prediction interval   71     83     93     97  
> cat("summary(earth.default, newdata=ozone1[1:100,]:)\n")
summary(earth.default, newdata=ozone1[1:100,]:)
> print(summary(earth.default, newdata=ozone1[1:100,]))
RSq 0.199 on newdata (100 cases)

                                 68%    80%    90%    95% 
newdata in prediction interval   78     85     92     95  
> plot(earth.default, level=.80, caption="earth.default")
> plotmo(earth.default, level=.80, col.response=3, caption="earth.default\nlevel = .80")
> 
> # TODO the following give err msg as expected, but do not give a try error
> # expect.err(try(summary(earth.mod1, newdata=c(1,2,3))))
> # expect.err(try(summary(earth.mod1, newdata=ozone1[1:100,1:3])))
> 
> multifigure("plot(earth.mod1)", 2, 2)
===plot(earth.mod1)
> set.seed(5)
> earth.mod2 <- earth(y~x, nfold=10, ncross=5, varmod.method="earth")
> plot(earth.mod2, caption="plot(earth.mod2)", level=.95)
> do.caption()
> 
> multifigure("plot(earth.mod2) with standardize=TRUE", 2, 2)
===plot(earth.mod2) with standardize=TRUE
> plot(earth.mod2, standardize=TRUE, level=.95,
+      caption="plot(earth.mod2, standardize=TRUE, level=.95)")
> do.caption()
> 
> multifigure("plot.varmod by calling plot(earth.mod2$varmod)", 2, 2)
===plot.varmod by calling plot(earth.mod2$varmod)
> plot(earth.mod2$varmod)
> 
> multifigure("embedded earth model by calling plot(earth.mod2$varmod$residmod)", 2, 2)
===embedded earth model by calling plot(earth.mod2$varmod$residmod)
> plot(earth.mod2$varmod$residmod, caption="embedded earth model")
> do.caption()
> 
> # test varmod.* args like varmod.conv
> 
> # cat("test varmod.exponent=.5\n")
> # set.seed(1)
> # (earth(Volume~Girth, data=trees, nfold=3, ncross=3, varmod.method="lm", trace=.3, varmod.exponent=.5))
> 
> # cat("test varmod.lambda=2/3\n")
> # set.seed(1)
> # (earth(Volume~Girth, data=trees, nfold=3, ncross=3, varmod.method="lm", trace=.3, varmod.lambda=2/3))
> 
> cat("test varmod.conv=50%\n")
test varmod.conv=50%
> set.seed(1)
> (earth(Volume~Girth, data=trees, nfold=3, ncross=3, varmod.method="lm", trace=.3, varmod.conv=50))
Full model GRSq 0.949 RSq 0.962, starting cross validation
CV fold  1.1  CVRSq  0.965   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  1.2  CVRSq  0.908   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  1.3  CVRSq  0.969   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  2.1  CVRSq  0.950   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  2.2  CVRSq  0.935   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  2.3  CVRSq  0.943   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  3.1  CVRSq  0.918   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  3.2  CVRSq  0.961   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  3.3  CVRSq  0.934   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV all        CVRSq  0.943                  n.infold.nz 31 100%


varmod method="lm" rmethod="hc12" lambda=1 exponent=1 conv=50 clamp=0.1 minspan=-3:
     iter weight.ratio coefchange% (Intercept) Volume
        1           11           0        1.07  0.063
        2           22          36        0.59  0.080

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.9425888
> 
> cat("test varmod.conv=-5\n")
test varmod.conv=-5
> set.seed(1)
> (earth(Volume~Girth, data=trees, nfold=3, ncross=3, varmod.method="lm", trace=.3, varmod.conv=-5))
Full model GRSq 0.949 RSq 0.962, starting cross validation
CV fold  1.1  CVRSq  0.965   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  1.2  CVRSq  0.908   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  1.3  CVRSq  0.969   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  2.1  CVRSq  0.950   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  2.2  CVRSq  0.935   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  2.3  CVRSq  0.943   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  3.1  CVRSq  0.918   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  3.2  CVRSq  0.961   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  3.3  CVRSq  0.934   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV all        CVRSq  0.943                  n.infold.nz 31 100%


varmod method="lm" rmethod="hc12" lambda=1 exponent=1 conv=-5 clamp=0.1 minspan=-3:
     iter weight.ratio coefchange% (Intercept) Volume
        1           11         0.0        1.07  0.063
        2           22        36.3        0.59  0.080
        3           27        14.9        0.45  0.086
        4           28         5.2        0.42  0.087
        5           29         1.6        0.40  0.088

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.9425888
> 
> cat("test varmod.clamp\n")
test varmod.clamp
> set.seed(1)
> a.noclamp <- earth(Volume~Girth, data=trees, nfold=3, ncross=3, varmod.method="lm")
> plot(a.noclamp$varmod, which=1:2, caption="a.noclamp and a.clamp", do.par=FALSE)
> set.seed(1)
> a.clamp <- earth(Volume~Girth, data=trees, nfold=3, ncross=3, varmod.method="lm", varmod.clamp=.6)
> plot(a.clamp$varmod, which=1:2, caption="", do.par=FALSE)
> 
> cat("test varmod.minspan=-5\n")
test varmod.minspan=-5
> set.seed(1)
> a.varmod.minspan.minus5 <- earth(Volume~Girth, data=trees, nfold=3, ncross=3, varmod.method="earth", trace=.3, varmod.minspan=-5)
Full model GRSq 0.949 RSq 0.962, starting cross validation
CV fold  1.1  CVRSq  0.965   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  1.2  CVRSq  0.908   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  1.3  CVRSq  0.969   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  2.1  CVRSq  0.950   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  2.2  CVRSq  0.935   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  2.3  CVRSq  0.943   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  3.1  CVRSq  0.918   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  3.2  CVRSq  0.961   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  3.3  CVRSq  0.934   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV all        CVRSq  0.943                  n.infold.nz 31 100%


varmod method="earth" rmethod="hc12" lambda=1 exponent=1 conv=1 clamp=0.1 minspan=-5:
     iter weight.ratio coefchange% (Intercept) h(28.7668-Volume)
        1           34           0           4             -0.18

> print(coef(a.varmod.minspan.minus5$varmod))
      (Intercept) h(28.7668-Volume) 
         4.975862         -0.227696 
> cat("test varmod.minspan=1\n")
test varmod.minspan=1
> set.seed(1)
> a.varmod.minspan1 <- earth(Volume~Girth, data=trees, nfold=3, ncross=3, varmod.method="earth", trace=.3, varmod.minspan=1)
Full model GRSq 0.949 RSq 0.962, starting cross validation
CV fold  1.1  CVRSq  0.965   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  1.2  CVRSq  0.908   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  1.3  CVRSq  0.969   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  2.1  CVRSq  0.950   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  2.2  CVRSq  0.935   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  2.3  CVRSq  0.943   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  3.1  CVRSq  0.918   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  3.2  CVRSq  0.961   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  3.3  CVRSq  0.934   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV all        CVRSq  0.943                  n.infold.nz 31 100%


varmod method="earth" rmethod="hc12" lambda=1 exponent=1 conv=1 clamp=0.1 minspan=1:
     iter weight.ratio coefchange% (Intercept) h(25.6817-Volume)
        1           25           0         3.8             -0.23

> print(coef(a.varmod.minspan1$varmod))
      (Intercept) h(25.6817-Volume) 
        4.7942482        -0.2863875 
> 
> # gam and y.gam are repeated below and on the repeat we will use the mgcv not gam package
> use.mgcv.package <- FALSE
> 
> for(varmod.method in c(earth:::VARMOD.METHODS, "gam", "x.gam")) {
+ 
+     multifigure(sprint("varmod.method=\"%s\"", varmod.method), 2, 3)
+     par(mar = c(3, 3, 2, 3)) # space for right margin axis
+ 
+     if(varmod.method %in% c("gam", "x.gam")) {
+         if(use.mgcv.package) {
+             # TODO with R 3.2.1 unload(gam) no longer works
+             cat("skipping mgcv tests\n")
+             next                                # NOTE next
+             cat("library(mgcv)\n")
+             library(mgcv)
+         } else
+             library(gam)
+     }
+     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))
+ 
+     if(use.mgcv.package && (varmod.method == "x.gam" || varmod.method == "gam")) {
+         # summary(mgcv) prints environment as hex address which messes up the diffs
+         printh("skipping summary(mgcv::gam) etc.\n")
+     } else {
+         printh("earth.mod$varmod")
+         print(earth.mod$varmod, style="unit")
+ 
+         printh("summary(earth.mod$varmod)")
+         print(summary(earth.mod$varmod))
+ 
+         printh("summary(earth.mod$varmod$residmod)")
+         print(summary(earth.mod$varmod$residmod))
+     }
+     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"))
+ 
+     # on second use of gam and y.gam we want to use the mgcv package
+     if(varmod.method == "x.gam" && !use.mgcv.package) {
+         use.mgcv.package <- TRUE
+         cat("detach(\"package:gam\", unload=TRUE)\n")
+         detach("package:gam", unload=TRUE)
+     }
+ }
===varmod.method="const"
Full model GRSq 0.949 RSq 0.962, starting cross validation
CV fold  1.1  CVRSq  0.864   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  1.2  CVRSq  0.932   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  1.3  CVRSq  0.888   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  2.1  CVRSq  0.950   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  2.2  CVRSq  0.966   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  2.3  CVRSq  0.842   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  3.1  CVRSq  0.945   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  3.2  CVRSq  0.884   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  3.3  CVRSq  0.950   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV all        CVRSq  0.913                  n.infold.nz 31 100%


varmod method="const" rmethod="hc12" lambda=1 exponent=1 conv=1 clamp=0.1 minspan=-3:
     iter weight.ratio coefchange% (Intercept)
        1            1           0         3.2

===varmod.method const: 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 "const"    min.sd 0.4

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)     3.995368     0.48554           12

                              mean   smallest    largest   ratio
95% prediction interval   15.66155   15.66155   15.66155       1

                                         68%    80%    90%    95% 
response values in prediction interval   81     94     97     97  
===earth.mod$varmod
method "const"    min.sd 0.4

stddev of predictions (scaled by unit):
            coefficients iter.stderr iter.stderr%
(Intercept)            1       0.122           12
unit                   4          NA           NA

                          mean   smallest   largest   ratio
95% prediction interval     16         16        16       1

                                         68%    80%    90%    95% 
response values in prediction interval   81     94     97     97  
===summary(earth.mod$varmod)
Parent model: 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)

method "const"    min.sd 0.4

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)            4        0.49           12

                          mean   smallest   largest   ratio
95% prediction interval     16         16        16       1

                                         68%    80%    90%    95% 
response values in prediction interval   81     94     97     97  

Regression submodel (Abs Residuals):

Call:
lm(formula = abs.resids ~ 1, data = data, weights = weights, 
    y = TRUE)

Coefficients:
(Intercept)  
        3.2  

===summary(earth.mod$varmod$residmod)

Call:
lm(formula = abs.resids ~ 1, data = data, weights = weights, 
    y = TRUE)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.7649 -2.0649  0.1714  1.2442  5.6925 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   3.1878     0.3874   8.229 3.48e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.157 on 30 degrees of freedom

===varmod.method const: predict(earth.mod, interval="pint")
         fit       lwr      upr
1   9.913855  2.083078 17.74463
2  10.942195  3.111418 18.77297
3  11.627755  3.796979 19.45853
4  17.455018  9.624241 25.28579
5  18.140578 10.309802 25.97136
6  18.483359 10.652582 26.31414
7  19.168919 11.338142 26.99970
8  19.168919 11.338142 26.99970
9  19.511699 11.680922 27.34248
10 19.854479 12.023703 27.68526
11 20.197259 12.366483 28.02804
12 20.540040 12.709263 28.37082
13 20.540040 12.709263 28.37082
14 21.568380 13.737603 29.39916
15 22.596721 14.765944 30.42750
16 25.681742 17.850965 33.51252
17 25.681742 17.850965 33.51252
18 27.052863 19.222086 34.88364
19 28.423983 20.593207 36.25476
20 28.766764 20.935987 36.59754
21 30.080913 22.250136 37.91169
22 31.395063 23.564286 39.22584
23 33.366287 25.535510 41.19706
24 43.222408 35.391631 51.05318
25 45.193632 37.362855 53.02441
26 51.764379 43.933603 59.59516
27 53.078529 45.247752 60.90931
28 55.706828 47.876051 63.53760
29 56.363903 48.533126 64.19468
30 56.363903 48.533126 64.19468
31 73.447846 65.617069 81.27862
===varmod.method="power"
Full model GRSq 0.949 RSq 0.962, starting cross validation
CV fold  1.1  CVRSq  0.864   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  1.2  CVRSq  0.932   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  1.3  CVRSq  0.888   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  2.1  CVRSq  0.950   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  2.2  CVRSq  0.966   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  2.3  CVRSq  0.842   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  3.1  CVRSq  0.945   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  3.2  CVRSq  0.884   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  3.3  CVRSq  0.950   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV all        CVRSq  0.913                  n.infold.nz 31 100%


varmod method="power" rmethod="hc12" lambda=1 exponent=1 conv=1 clamp=0.1 minspan=-3:

     start: (Intercept)=0.66    coef=0.084    exponent=1

     iter weight.ratio coefchange% (Intercept) coef exponent
        1           67         0.0        -9.5 6.31     0.21
        2           41       109.7        -1.7 0.79     0.55
        3           44         8.4        -2.0 0.87     0.53
        4           44         2.2        -1.9 0.84     0.54
        5           44         0.2        -1.9 0.85     0.54

===varmod.method power: 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 "power"    min.sd 0.4    iter.rsq 0.373

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)   -2.3895955     7.82659          328
coef           1.0600809     4.01852          379
exponent       0.5366486    0.758238          141

                              mean   smallest    largest      ratio
95% prediction interval   15.68284   4.864412   32.31915   6.643998

                                         68%    80%    90%    95% 
response values in prediction interval   81     90     97     97  
===earth.mod$varmod
method "power"    min.sd 0.4    iter.rsq 0.373

stddev of predictions (scaled by unit):
            coefficients iter.stderr iter.stderr%
(Intercept)        -2.25         7.4          328
coef                1.00         3.8          379
exponent            0.54         0.8          141
unit                1.06          NA           NA

                          mean   smallest   largest   ratio
95% prediction interval     16        4.9        32     6.6

                                         68%    80%    90%    95% 
response values in prediction interval   81     90     97     97  
===summary(earth.mod$varmod)
Parent model: 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)

method "power"    min.sd 0.4    iter.rsq 0.373

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)        -2.39         7.8          328
coef                1.06           4          379
exponent            0.54         0.8          141

                          mean   smallest   largest   ratio
95% prediction interval     16        4.9        32     6.6

                                         68%    80%    90%    95% 
response values in prediction interval   81     90     97     97  

Regression submodel (Abs Residuals):

Nonlinear regression model
  model: abs.resids ~ `(Intercept)` + coef * RHS^exponent
   data: data
(Intercept)        coef    exponent 
      -1.91        0.85        0.54 
 weighted residual sum-of-squares: 57

Number of iterations to convergence: 3 
Achieved convergence tolerance: 2.6e-10
===summary(earth.mod$varmod$residmod)

Formula: abs.resids ~ `(Intercept)` + coef * RHS^exponent

Parameters:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -1.9066     6.2447  -0.305    0.762
coef          0.8458     3.2063   0.264    0.794
exponent      0.5366     0.7582   0.708    0.485

Residual standard error: 1.427 on 28 degrees of freedom

Number of iterations to convergence: 3 
Achieved convergence tolerance: 2.588e-10

===varmod.method power: predict(earth.mod, interval="pint")
         fit       lwr      upr
1   9.913855  7.481648 12.34606
2  10.942195  8.122956 13.76143
3  11.627755  8.559808 14.69570
4  17.455018 12.498875 22.41116
5  18.140578 12.983071 23.29809
6  18.483359 13.226492 23.74023
7  19.168919 13.715863 24.62198
8  19.168919 13.715863 24.62198
9  19.511699 13.961768 25.06163
10 19.854479 14.208458 25.50050
11 20.197259 14.455914 25.93860
12 20.540040 14.704117 26.37596
13 20.540040 14.704117 26.37596
14 21.568380 15.453028 27.68373
15 22.596721 16.208047 28.98539
16 25.681742 18.505934 32.85755
17 25.681742 18.505934 32.85755
18 27.052863 19.541368 34.56436
19 28.423983 20.584599 36.26337
20 28.766764 20.846559 36.68697
21 30.080913 21.854919 38.30691
22 31.395063 22.869409 39.92072
23 33.366287 24.401834 42.33074
24 43.222408 32.224413 54.22040
25 45.193632 33.815804 56.57146
26 51.764379 39.172844 64.35592
27 53.078529 40.253006 65.90405
28 55.706828 42.421244 68.99241
29 56.363903 42.964885 69.76292
30 56.363903 42.964885 69.76292
31 73.447846 57.288273 89.60742
===varmod.method="power0"
===varmod.method power0: 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 "power0"    min.sd 0.403    iter.rsq 0.329

stddev of predictions:
         coefficients iter.stderr iter.stderr%
coef        0.1819598    0.136285           75
exponent    0.9124346    0.225676           25

                              mean   smallest    largest      ratio
95% prediction interval   15.80737    5.78442   35.96148   6.216955

                                         68%    80%    90%    95% 
response values in prediction interval   81     90     97     97  
===earth.mod$varmod
method "power0"    min.sd 0.403    iter.rsq 0.329

stddev of predictions (scaled by unit):
         coefficients iter.stderr iter.stderr%
coef             0.20       0.149           75
exponent         0.91       0.226           25
unit             0.91          NA           NA

                          mean   smallest   largest   ratio
95% prediction interval     16        5.8        36     6.2

                                         68%    80%    90%    95% 
response values in prediction interval   81     90     97     97  
===summary(earth.mod$varmod)
Parent model: 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)

method "power0"    min.sd 0.403    iter.rsq 0.329

stddev of predictions:
         coefficients iter.stderr iter.stderr%
coef             0.18       0.136           75
exponent         0.91       0.226           25

                          mean   smallest   largest   ratio
95% prediction interval     16        5.8        36     6.2

                                         68%    80%    90%    95% 
response values in prediction interval   81     90     97     97  

Regression submodel (Abs Residuals):

Nonlinear regression model
  model: abs.resids ~ coef * RHS^exponent
   data: data
    coef exponent 
    0.15     0.91 
 weighted residual sum-of-squares: 66

Number of iterations to convergence: 4 
Achieved convergence tolerance: 1.2e-10
===summary(earth.mod$varmod$residmod)

Formula: abs.resids ~ coef * RHS^exponent

Parameters:
         Estimate Std. Error t value Pr(>|t|)    
coef       0.1452     0.1087   1.335 0.192217    
exponent   0.9124     0.2257   4.043 0.000356 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.506 on 29 degrees of freedom

Number of iterations to convergence: 4 
Achieved convergence tolerance: 1.212e-10

===varmod.method power0: predict(earth.mod, interval="pint")
         fit       lwr      upr
1   9.913855  7.021645 12.80606
2  10.942195  7.777452 14.10694
3  11.627755  8.282580 14.97293
4  17.455018 12.608892 22.30114
5  18.140578 13.121078 23.16008
6  18.483359 13.377388 23.58933
7  19.168919 13.890424 24.44741
8  19.168919 13.890424 24.44741
9  19.511699 14.147146 24.87625
10 19.854479 14.404001 25.30496
11 20.197259 14.660985 25.73353
12 20.540040 14.918096 26.16198
13 20.540040 14.918096 26.16198
14 21.568380 15.690173 27.44659
15 22.596721 16.463317 28.73012
16 25.681742 18.788655 32.57483
17 25.681742 18.788655 32.57483
18 27.052863 19.824756 34.28097
19 28.423983 20.862342 35.98562
20 28.766764 21.121961 36.41157
21 30.080913 22.118082 38.04374
22 31.395063 23.115417 39.67471
23 33.366287 24.613579 42.11899
24 43.222408 32.138291 54.30652
25 45.193632 33.649178 56.73809
26 51.764379 38.697712 64.83105
27 53.078529 39.709517 66.44754
28 55.706828 41.735075 69.67858
29 56.363903 42.241858 70.48595
30 56.363903 42.241858 70.48595
31 73.447846 55.467108 91.42858
===varmod.method="lm"
Full model GRSq 0.949 RSq 0.962, starting cross validation
CV fold  1.1  CVRSq  0.864   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  1.2  CVRSq  0.932   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  1.3  CVRSq  0.888   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  2.1  CVRSq  0.950   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  2.2  CVRSq  0.966   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  2.3  CVRSq  0.842   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  3.1  CVRSq  0.945   n.oof 20 35%   n.infold.nz 20 100%   n.oof.nz 11 100%
CV fold  3.2  CVRSq  0.884   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV fold  3.3  CVRSq  0.950   n.oof 21 32%   n.infold.nz 21 100%   n.oof.nz 10 100%
CV all        CVRSq  0.913                  n.infold.nz 31 100%


varmod method="lm" rmethod="hc12" lambda=1 exponent=1 conv=1 clamp=0.1 minspan=-3:
     iter weight.ratio coefchange% (Intercept) Volume
        1           21        0.00        0.66  0.084
        2           40       44.82        0.20  0.101
        3           45       20.39        0.13  0.104
        4           46        4.64        0.11  0.104
        5           46        0.72        0.11  0.104

===varmod.method lm: 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 "lm"    min.sd 0.408    iter.rsq 0.321

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)    0.1419992    0.759816          535
Volume         0.1306286   0.0352932           27

                              mean   smallest    largest      ratio
95% prediction interval   16.00581   5.633063   38.16594   6.775345

                                         68%    80%    90%    95% 
response values in prediction interval   77     90     97     97  
===earth.mod$varmod
method "lm"    min.sd 0.408    iter.rsq 0.321

stddev of predictions (scaled by unit):
            coefficients iter.stderr iter.stderr%
(Intercept)         1.09         5.8          535
Volume              1.00         0.3           27
unit                0.13          NA           NA

                          mean   smallest   largest   ratio
95% prediction interval     16        5.6        38     6.8

                                         68%    80%    90%    95% 
response values in prediction interval   77     90     97     97  
===summary(earth.mod$varmod)
Parent model: 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)

method "lm"    min.sd 0.408    iter.rsq 0.321

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)         0.14        0.76          535
Volume              0.13        0.04           27

                          mean   smallest   largest   ratio
95% prediction interval     16        5.6        38     6.8

                                         68%    80%    90%    95% 
response values in prediction interval   77     90     97     97  

Regression submodel (Abs Residuals):

Call:
lm(formula = abs.resids ~ ., data = data, weights = weights, 
    y = TRUE)

Coefficients:
(Intercept)          RHS  
       0.11         0.10  

===summary(earth.mod$varmod$residmod)

Call:
lm(formula = abs.resids ~ ., data = data, weights = weights, 
    y = TRUE)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-1.8911 -0.9998 -0.2203  0.7893  5.1212 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.11330    0.60625   0.187 0.853051    
RHS          0.10423    0.02816   3.701 0.000895 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.497 on 29 degrees of freedom
Multiple R-squared:  0.3208,	Adjusted R-squared:  0.2974 
F-statistic:  13.7 on 1 and 29 DF,  p-value: 0.0008945

===varmod.method lm: predict(earth.mod, interval="pint")
         fit       lwr      upr
1   9.913855  7.097323 12.73039
2  10.942195  7.862380 14.02201
3  11.627755  8.372418 14.88309
4  17.455018 12.707742 22.20229
5  18.140578 13.217781 23.06338
6  18.483359 13.472800 23.49392
7  19.168919 13.982838 24.35500
8  19.168919 13.982838 24.35500
9  19.511699 14.237857 24.78554
10 19.854479 14.492876 25.21608
11 20.197259 14.747895 25.64662
12 20.540040 15.002914 26.07717
13 20.540040 15.002914 26.07717
14 21.568380 15.767971 27.36879
15 22.596721 16.533028 28.66041
16 25.681742 18.828200 32.53528
17 25.681742 18.828200 32.53528
18 27.052863 19.848276 34.25745
19 28.423983 20.868352 35.97961
20 28.766764 21.123371 36.41016
21 30.080913 22.101063 38.06076
22 31.395063 23.078754 39.71137
23 33.366287 24.545291 42.18728
24 43.222408 31.877975 54.56684
25 45.193632 33.344512 57.04275
26 51.764379 38.232968 65.29579
27 53.078529 39.210659 66.94640
28 55.706828 41.166042 70.24761
29 56.363903 41.654888 71.07292
30 56.363903 41.654888 71.07292
31 73.447846 54.364874 92.53082
===varmod.method="rlm"
===varmod.method rlm: 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 "rlm"    min.sd 0.374    iter.rsq 0.308

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)    0.1702205    0.666088          391
Volume         0.1183674   0.0307541           26

                              mean   smallest    largest     ratio
95% prediction interval   14.66633   5.267199   34.74645   6.59676

                                         68%    80%    90%    95% 
response values in prediction interval   74     90     97     97  
===earth.mod$varmod
method "rlm"    min.sd 0.374    iter.rsq 0.308

stddev of predictions (scaled by unit):
            coefficients iter.stderr iter.stderr%
(Intercept)         1.44         5.6          391
Volume              1.00         0.3           26
unit                0.12          NA           NA

                          mean   smallest   largest   ratio
95% prediction interval     15        5.3        35     6.6

                                         68%    80%    90%    95% 
response values in prediction interval   74     90     97     97  
===summary(earth.mod$varmod)
Parent model: 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)

method "rlm"    min.sd 0.374    iter.rsq 0.308

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)         0.17        0.67          391
Volume              0.12        0.03           26

                          mean   smallest   largest   ratio
95% prediction interval     15        5.3        35     6.6

                                         68%    80%    90%    95% 
response values in prediction interval   74     90     97     97  

Regression submodel (Abs Residuals):

Call:
rlm(formula = abs.resids ~ ., data = data, weights = weights, 
    method = "MM")
Converged in 6 iterations

Coefficients:
(Intercept)         RHS 
      0.136       0.094 

Degrees of freedom: 31 total; 29 residual
Scale estimate: 1.32 
===summary(earth.mod$varmod$residmod)

Call: rlm(formula = abs.resids ~ ., data = data, weights = weights, 
    method = "MM")
Residuals:
     Min       1Q   Median       3Q      Max 
-1.71205 -0.81715 -0.06093  0.96673  5.34892 

Coefficients:
            Value  Std. Error t value
(Intercept) 0.1358 0.5315     0.2556 
RHS         0.0944 0.0245     3.8488 

Residual standard error: 1.324 on 29 degrees of freedom
===varmod.method rlm: predict(earth.mod, interval="pint")
         fit       lwr      upr
1   9.913855  7.280255 12.54745
2  10.942195  8.070025 13.81437
3  11.627755  8.596538 14.65897
4  17.455018 13.071900 21.83814
5  18.140578 13.598413 22.68274
6  18.483359 13.861669 23.10505
7  19.168919 14.388182 23.94966
8  19.168919 14.388182 23.94966
9  19.511699 14.651439 24.37196
10 19.854479 14.914696 24.79426
11 20.197259 15.177952 25.21657
12 20.540040 15.441209 25.63887
13 20.540040 15.441209 25.63887
14 21.568380 16.230978 26.90578
15 22.596721 17.020748 28.17269
16 25.681742 19.390057 31.97343
17 25.681742 19.390057 31.97343
18 27.052863 20.443084 33.66264
19 28.423983 21.496110 35.35186
20 28.766764 21.759366 35.77416
21 30.080913 22.768639 37.39319
22 31.395063 23.777911 39.01221
23 33.366287 25.291819 41.44075
24 43.222408 32.861360 53.58346
25 45.193632 34.375269 56.01200
26 51.764379 39.421630 64.10713
27 53.078529 40.430902 65.72616
28 55.706828 42.449446 68.96421
29 56.363903 42.954082 69.77372
30 56.363903 42.954082 69.77372
31 73.447846 56.074621 90.82107
===varmod.method="earth"
===varmod.method earth: 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 "earth"    min.sd 0.4    iter.rsq 0.441

stddev of predictions:
                  coefficients iter.stderr iter.stderr%
(Intercept)          4.6397622    0.560983           12
h(25.6817-Volume)   -0.2734766   0.0859684           31
h(Volume-43.2224)    0.1299760   0.0612411           47

                              mean   smallest    largest     ratio
95% prediction interval   15.67065   1.566155   33.58729   21.4457

                                         68%    80%    90%     95% 
response values in prediction interval   77     94     94     100  
===earth.mod$varmod
method "earth"    min.sd 0.4    iter.rsq 0.441

stddev of predictions (scaled by unit):
                  coefficients iter.stderr iter.stderr%
(Intercept)              16.97        2.05           12
h(25.6817-Volume)        -1.00        0.31           31
h(Volume-43.2224)         0.48        0.22           47

                          mean   smallest   largest   ratio
95% prediction interval     16        1.6        34      21

                                         68%    80%    90%     95% 
response values in prediction interval   77     94     94     100  
===summary(earth.mod$varmod)
Parent model: 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)

method "earth"    min.sd 0.4    iter.rsq 0.441

stddev of predictions:
                  coefficients iter.stderr iter.stderr%
(Intercept)               4.64        0.56           12
h(25.6817-Volume)        -0.27        0.09           31
h(Volume-43.2224)         0.13        0.06           47

                          mean   smallest   largest   ratio
95% prediction interval     16        1.6        34      21

                                         68%    80%    90%     95% 
response values in prediction interval   77     94     94     100  

Regression submodel (Abs Residuals):

Selected 3 of 4 terms, and 1 of 1 predictors
Termination condition: RSq changed by less than 0.001 at 4 terms
Importance: RHS
Number of terms at each degree of interaction: 1 2 (additive model)
GCV 3.6    RSS 78    GRSq 0.26    RSq 0.44
===summary(earth.mod$varmod$residmod)
Call: earth(formula=abs.resids~., data=data, weights=weights, keepxy=TRUE,
            trace=trace, degree=..1, minspan=minspan)

               coefficients
(Intercept)       3.7019947
h(25.6817-RHS)   -0.2182027
h(RHS-43.2224)    0.1037059

Selected 3 of 4 terms, and 1 of 1 predictors
Termination condition: RSq changed by less than 0.001 at 4 terms
Importance: RHS
Number of terms at each degree of interaction: 1 2 (additive model)
GCV 3.578019    RSS 78.0239    GRSq 0.2557634    RSq 0.4409956
===varmod.method earth: predict(earth.mod, interval="pint")
         fit       lwr      upr
1   9.913855  9.130777 10.69693
2  10.942195  9.748888 12.13550
3  11.627755 10.066985 13.18853
4  17.455018 12.770810 22.13923
5  18.140578 13.088907 23.19225
6  18.483359 13.247956 23.71876
7  19.168919 13.566053 24.77179
8  19.168919 13.566053 24.77179
9  19.511699 13.725101 25.29830
10 19.854479 13.884150 25.82481
11 20.197259 14.043198 26.35132
12 20.540040 14.202247 26.87783
13 20.540040 14.202247 26.87783
14 21.568380 14.679393 28.45737
15 22.596721 15.156538 30.03690
16 25.681742 16.587975 34.77551
17 25.681742 16.587975 34.77551
18 27.052863 17.959096 36.14663
19 28.423983 19.330217 37.51775
20 28.766764 19.672997 37.86053
21 30.080913 20.987146 39.17468
22 31.395063 22.301296 40.48883
23 33.366287 24.272520 42.46005
24 43.222408 34.128641 52.31617
25 45.193632 35.597699 54.78956
26 51.764379 40.494560 63.03420
27 53.078529 41.473932 64.68313
28 55.706828 43.432676 67.98098
29 56.363903 43.922362 68.80544
30 56.363903 43.922362 68.80544
31 73.447846 56.654200 90.24149
===varmod.method="gam"
Loading required package: splines
Loading required package: foreach
Loaded gam 1.16

===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" (gam package)    min.sd 0.395    iter.rsq 0.403

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)    0.2023532    0.670707          331
s(Volume)      0.1227670   0.0312472           25

                              mean   smallest    largest      ratio
95% prediction interval   15.46709   5.446277   31.90195   5.857571

                                         68%    80%    90%    95% 
response values in prediction interval   84     90     97     97  
===earth.mod$varmod
method "gam" (gam package)    min.sd 0.395    iter.rsq 0.403

stddev of predictions (scaled by unit):
            coefficients iter.stderr iter.stderr%
(Intercept)         1.65         5.5          331
s(Volume)           1.00         0.3           25
unit                0.12          NA           NA

                          mean   smallest   largest   ratio
95% prediction interval     15        5.4        32     5.9

                                         68%    80%    90%    95% 
response values in prediction interval   84     90     97     97  
===summary(earth.mod$varmod)
Parent model: 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)

method "gam" (gam package)    min.sd 0.395    iter.rsq 0.403

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)         0.20        0.67          331
s(Volume)           0.12        0.03           25

                          mean   smallest   largest   ratio
95% prediction interval     15        5.4        32     5.9

                                         68%    80%    90%    95% 
response values in prediction interval   84     90     97     97  

Regression submodel (Abs Residuals):

Call:
gam::gam(formula = form, data = data, weights = weights)

Degrees of Freedom: 30 total; 26 Residual
Residual Deviance: 50.02 
===summary(earth.mod$varmod$residmod)

Call: gam::gam(formula = form, data = data, weights = weights)
Deviance Residuals:
    Min      1Q  Median      3Q     Max 
-1.7965 -0.8764 -0.1362  0.5448  3.8774 

(Dispersion Parameter for gaussian family taken to be 1.9238)

    Null Deviance: 83.84 on 30 degrees of freedom
Residual Deviance: 50.0183 on 25.9999 degrees of freedom
AIC: 128.7914 

Number of Local Scoring Iterations: 2 

Anova for Parametric Effects
          Df Sum Sq Mean Sq F value    Pr(>F)    
s(RHS)     1 29.696 29.6960  15.436 0.0005626 ***
Residuals 26 50.018  1.9238                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova for Nonparametric Effects
            Npar Df  Npar F Pr(F)
(Intercept)                      
s(RHS)            3 0.71485 0.552
===varmod.method gam: predict(earth.mod, interval="pint")
         fit       lwr      upr
1   9.913855  7.190716 12.63699
2  10.942195  8.050124 13.83427
3  11.627755  8.625781 14.62973
4  17.455018 13.300703 21.60933
5  18.140578 13.752619 22.52854
6  18.483359 13.966702 23.00002
7  19.168919 14.374580 23.96326
8  19.168919 14.374580 23.96326
9  19.511699 14.570431 24.45297
10 19.854479 14.761459 24.94750
11 20.197259 14.947249 25.44727
12 20.540040 15.127820 25.95226
13 20.540040 15.127820 25.95226
14 21.568380 15.635922 27.50084
15 22.596721 16.109219 29.08422
16 25.681742 17.653837 33.70965
17 25.681742 17.653837 33.70965
18 27.052863 18.530731 35.57499
19 28.423983 19.522256 37.32571
20 28.766764 19.785799 37.74773
21 30.080913 20.847286 39.31454
22 31.395063 21.969618 40.82051
23 33.366287 23.700148 43.03243
24 43.222408 32.576322 53.86849
25 45.193632 34.346557 56.04071
26 51.764379 39.964352 63.56441
27 53.078529 41.042569 65.11449
28 55.706828 43.166481 68.24717
29 56.363903 43.693855 69.03395
30 56.363903 43.693855 69.03395
31 73.447846 57.496870 89.39882
===varmod.method="x.lm"
===varmod.method x.lm: 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.lm"    min.sd 0.398    iter.rsq 0.391

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)   -3.3689284     1.40128           42
Girth          0.5547878     0.12848           23

                              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  
===earth.mod$varmod
method "x.lm"    min.sd 0.398    iter.rsq 0.391

stddev of predictions (scaled by unit):
            coefficients iter.stderr iter.stderr%
(Intercept)        -6.07        2.53           42
Girth               1.00        0.23           23
unit                0.55          NA           NA

                          mean   smallest   largest   ratio
95% prediction interval     16        4.8        32     6.5

                                         68%    80%    90%    95% 
response values in prediction interval   81     90     97     97  
===summary(earth.mod$varmod)
Parent model: 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)

method "x.lm"    min.sd 0.398    iter.rsq 0.391

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)        -3.37         1.4           42
Girth               0.55        0.13           23

                          mean   smallest   largest   ratio
95% prediction interval     16        4.8        32     6.5

                                         68%    80%    90%    95% 
response values in prediction interval   81     90     97     97  

Regression submodel (Abs Residuals):

Call:
lm(formula = abs.resids ~ ., data = data, weights = weights, 
    y = TRUE)

Coefficients:
(Intercept)        Girth  
      -2.69         0.44  

===summary(earth.mod$varmod$residmod)

Call:
lm(formula = abs.resids ~ ., data = data, weights = weights, 
    y = TRUE)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-1.8562 -1.0072  0.0099  0.6783  4.4187 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -2.6880     1.1181  -2.404 0.022822 *  
Girth         0.4427     0.1025   4.318 0.000168 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.355 on 29 degrees of freedom
Multiple R-squared:  0.3913,	Adjusted R-squared:  0.3704 
F-statistic: 18.65 on 1 and 29 DF,  p-value: 0.0001677

===varmod.method x.lm: 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
===varmod.method="x.rlm"
===varmod.method x.rlm: 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.rlm"    min.sd 0.376    iter.rsq 0.385

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)   -3.0621287     1.29393           42
Girth          0.5146341    0.117998           23

                              mean   smallest    largest      ratio
95% prediction interval   14.72302   4.740503   29.55364   6.234285

                                         68%    80%    90%    95% 
response values in prediction interval   81     90     97     97  
===earth.mod$varmod
method "x.rlm"    min.sd 0.376    iter.rsq 0.385

stddev of predictions (scaled by unit):
            coefficients iter.stderr iter.stderr%
(Intercept)        -5.95        2.51           42
Girth               1.00        0.23           23
unit                0.51          NA           NA

                          mean   smallest   largest   ratio
95% prediction interval     15        4.7        30     6.2

                                         68%    80%    90%    95% 
response values in prediction interval   81     90     97     97  
===summary(earth.mod$varmod)
Parent model: 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)

method "x.rlm"    min.sd 0.376    iter.rsq 0.385

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)        -3.06        1.29           42
Girth               0.52        0.12           23

                          mean   smallest   largest   ratio
95% prediction interval     15        4.7        30     6.2

                                         68%    80%    90%    95% 
response values in prediction interval   81     90     97     97  

Regression submodel (Abs Residuals):

Call:
rlm(formula = abs.resids ~ ., data = data, weights = weights, 
    method = "MM", y.ret = TRUE)
Converged in 6 iterations

Coefficients:
(Intercept)       Girth 
      -2.44        0.41 

Degrees of freedom: 31 total; 29 residual
Scale estimate: 1.31 
===summary(earth.mod$varmod$residmod)

Call: rlm(formula = abs.resids ~ ., data = data, weights = weights, 
    method = "MM", y.ret = TRUE)
Residuals:
    Min      1Q  Median      3Q     Max 
-1.7566 -0.9024  0.1548  0.8198  4.6098 

Coefficients:
            Value   Std. Error t value
(Intercept) -2.4432  1.0324    -2.3665
Girth        0.4106  0.0941     4.3614

Residual standard error: 1.31 on 29 degrees of freedom
===varmod.method x.rlm: predict(earth.mod, interval="pint")
         fit       lwr      upr
1   9.913855  7.543603 12.28411
2  10.942195  8.269344 13.61505
3  11.627755  8.753172 14.50234
4  17.455018 12.865705 22.04433
5  18.140578 13.349533 22.93162
6  18.483359 13.591447 23.37527
7  19.168919 14.075274 24.26256
8  19.168919 14.075274 24.26256
9  19.511699 14.317188 24.70621
10 19.854479 14.559102 25.14986
11 20.197259 14.801015 25.59350
12 20.540040 15.042929 26.03715
13 20.540040 15.042929 26.03715
14 21.568380 15.768670 27.36809
15 22.596721 16.494412 28.69903
16 25.681742 18.671635 32.69185
17 25.681742 18.671635 32.69185
18 27.052863 19.639290 34.46644
19 28.423983 20.606945 36.24102
20 28.766764 20.848859 36.68467
21 30.080913 21.961276 38.20055
22 31.395063 23.073692 39.71643
23 33.366287 24.742317 41.99026
24 43.222408 33.085442 53.35937
25 45.193632 34.754067 55.63320
26 51.764379 40.316150 63.21261
27 53.078529 41.428567 64.72849
28 55.706828 43.653400 67.76026
29 56.363903 44.209608 68.51820
30 56.363903 44.209608 68.51820
31 73.447846 58.671024 88.22467
===varmod.method="x.earth"
===varmod.method x.earth: 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.earth"    min.sd 0.4    iter.rsq 0.401

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)    2.4649863    0.516488           21
h(Girth-11)    0.6082285    0.138056           23

                              mean   smallest   largest      ratio
95% prediction interval   15.66155   9.662569    32.551   3.368773

                                         68%    80%    90%    95% 
response values in prediction interval   81     90     97     97  
===earth.mod$varmod
method "x.earth"    min.sd 0.4    iter.rsq 0.401

stddev of predictions (scaled by unit):
            coefficients iter.stderr iter.stderr%
(Intercept)            4        0.85           21
h(Girth-11)            1        0.23           23

                          mean   smallest   largest   ratio
95% prediction interval     16        9.7        33     3.4

                                         68%    80%    90%    95% 
response values in prediction interval   81     90     97     97  
===summary(earth.mod$varmod)
Parent model: 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)

method "x.earth"    min.sd 0.4    iter.rsq 0.401

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)         2.46        0.52           21
h(Girth-11)         0.61        0.14           23

                          mean   smallest   largest   ratio
95% prediction interval     16        9.7        33     3.4

                                         68%    80%    90%    95% 
response values in prediction interval   81     90     97     97  

Regression submodel (Abs Residuals):

Selected 2 of 5 terms, and 1 of 1 predictors
Termination condition: RSq changed by less than 0.001 at 5 terms
Importance: Girth
Number of terms at each degree of interaction: 1 1 (additive model)
GCV 3.3    RSS 84    GRSq 0.31    RSq 0.4
===summary(earth.mod$varmod$residmod)
Call: earth(formula=abs.resids~., data=data, weights=weights, keepxy=TRUE,
            trace=trace, degree=..1, minspan=minspan)

            coefficients
(Intercept)    1.9667745
h(Girth-11)    0.4852961

Selected 2 of 5 terms, and 1 of 1 predictors
Termination condition: RSq changed by less than 0.001 at 5 terms
Importance: Girth
Number of terms at each degree of interaction: 1 1 (additive model)
GCV 3.306152    RSS 83.61366    GRSq 0.3123124    RSq 0.4009477
===varmod.method x.earth: predict(earth.mod, interval="pint")
         fit       lwr      upr
1   9.913855  5.082570 14.74514
2  10.942195  6.110911 15.77348
3  11.627755  6.796471 16.45904
4  17.455018 12.623734 22.28630
5  18.140578 13.309294 22.97186
6  18.483359 13.652074 23.31464
7  19.168919 14.337635 24.00020
8  19.168919 14.337635 24.00020
9  19.511699 14.561204 24.46219
10 19.854479 14.784774 24.92418
11 20.197259 15.008343 25.38618
12 20.540040 15.231913 25.84817
13 20.540040 15.231913 25.84817
14 21.568380 15.902622 27.23414
15 22.596721 16.573330 28.62011
16 25.681742 18.585456 32.77803
17 25.681742 18.585456 32.77803
18 27.052863 19.479735 34.62599
19 28.423983 20.374013 36.47395
20 28.766764 20.597583 36.93594
21 30.080913 21.673311 38.48852
22 31.395063 22.749039 40.04109
23 33.366287 24.362632 42.36994
24 43.222408 32.430594 54.01422
25 45.193632 34.044186 56.34308
26 51.764379 39.422828 64.10593
27 53.078529 40.498556 65.65850
28 55.706828 42.650013 68.76364
29 56.363903 43.187877 69.53993
30 56.363903 43.187877 69.53993
31 73.447846 57.172345 89.72335
===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" (gam package)    min.sd 0.395    iter.rsq 0.410

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)    -3.345908     1.49309           45
s(Volume)       0.543721    0.135621           25

                              mean   smallest    largest      ratio
95% prediction interval   15.48763   5.590087   31.47631   5.630736

                                         68%    80%    90%     95% 
response values in prediction interval   84     90     97     100  
===earth.mod$varmod
method "x.gam" (gam package)    min.sd 0.395    iter.rsq 0.410

stddev of predictions (scaled by unit):
            coefficients iter.stderr iter.stderr%
(Intercept)        -6.15        2.75           45
s(Volume)           1.00        0.25           25
unit                0.54          NA           NA

                          mean   smallest   largest   ratio
95% prediction interval     15        5.6        31     5.6

                                         68%    80%    90%     95% 
response values in prediction interval   84     90     97     100  
===summary(earth.mod$varmod)
Parent model: 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)

method "x.gam" (gam package)    min.sd 0.395    iter.rsq 0.410

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)        -3.35        1.49           45
s(Volume)           0.54        0.14           25

                          mean   smallest   largest   ratio
95% prediction interval     15        5.6        31     5.6

                                         68%    80%    90%     95% 
response values in prediction interval   84     90     97     100  

Regression submodel (Abs Residuals):

Call:
gam::gam(formula = form, data = data, weights = weights)

Degrees of Freedom: 30 total; 26 Residual
Residual Deviance: 47.92 
===summary(earth.mod$varmod$residmod)

Call: gam::gam(formula = form, data = data, weights = weights)
Deviance Residuals:
    Min      1Q  Median      3Q     Max 
-1.7707 -0.8448 -0.1895  0.4960  3.5435 

(Dispersion Parameter for gaussian family taken to be 1.843)

    Null Deviance: 81.1946 on 30 degrees of freedom
Residual Deviance: 47.9171 on 26 degrees of freedom
AIC: 127.7746 

Number of Local Scoring Iterations: 2 

Anova for Parametric Effects
          Df Sum Sq Mean Sq F value    Pr(>F)    
s(RHS)     1 29.622  29.622  16.073 0.0004565 ***
Residuals 26 47.917   1.843                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova for Nonparametric Effects
            Npar Df  Npar F  Pr(F)
(Intercept)                       
s(RHS)            3 0.66113 0.5834
===varmod.method x.gam: predict(earth.mod, interval="pint")
         fit       lwr      upr
1   9.913855  7.118811 12.70890
2  10.942195  7.996429 13.88796
3  11.627755  8.591126 14.66438
4  17.455018 13.557720 21.35232
5  18.140578 13.992212 22.28895
6  18.483359 14.190289 22.77643
7  19.168919 14.558207 23.77963
8  19.168919 14.558207 23.77963
9  19.511699 14.733724 24.28967
10 19.854479 14.904647 24.80431
11 20.197259 15.069130 25.32539
12 20.540040 15.226289 25.85379
13 20.540040 15.226289 25.85379
14 21.568380 15.638894 27.49787
15 22.596721 15.982827 29.21061
16 25.681742 17.233368 34.13012
17 25.681742 17.233368 34.13012
18 27.052863 18.124009 35.98172
19 28.423983 19.173304 37.67466
20 28.766764 19.452897 38.08063
21 30.080913 20.656244 39.50558
22 31.395063 21.872191 40.91793
23 33.366287 23.695357 43.03722
24 43.222408 32.616128 53.82869
25 45.193632 34.353740 56.03352
26 51.764379 39.963839 63.56492
27 53.078529 41.057784 65.09927
28 55.706828 43.225535 68.18812
29 56.363903 43.764749 68.96306
30 56.363903 43.764749 68.96306
31 73.447846 57.709693 89.18600
detach("package:gam", unload=TRUE)
===varmod.method="gam"
skipping mgcv tests
===varmod.method="x.gam"
skipping mgcv tests
> # test varmod.exponent
> set.seed(6)
> earth.exponent <- earth(Volume~Girth, data=trees, nfold=3, ncross=3,
+                         varmod.method="lm", varmod.exponent=.5)
> printh("summary(earth.exponent)")
===summary(earth.exponent)
> print(summary(earth.exponent))
Call: earth(formula=Volume~Girth, data=trees, nfold=3, ncross=3,
            varmod.method="lm", varmod.exponent=0.5)

              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 "lm"  exponent 0.500  min.sd 0.399  iter.rsq 0.377

stddev of predictions:
             coefficients iter.stderr iter.stderr%
(Intercept)     -2.780989     1.31593           47
sqrt(Volume)     1.273693    0.304316           24

                              mean   smallest    largest      ratio
95% prediction interval   15.65879    4.81914   31.88774   6.616895

                                         68%    80%    90%    95% 
response values in prediction interval   81     94     97     97  
> 
> par(old.par)
> 
> source("test.epilog.R")
