1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
|
> # test.varmmod.mgcv.R
> # mgcv has to be tested separately because of clashes between library(gam) and library(mgcv)
> # Stephen Milborrow Apr 2015 Berea
>
> source("test.prolog.R")
> library(earth)
Loading required package: plotmo
Loading required package: plotrix
Loading required package: TeachingDemos
> options(warn=2)
>
> printh <- function(caption)
+ cat("===", caption, "\n", sep="")
>
> CAPTION <- NULL
>
> multifigure <- function(caption, nrow=3, ncol=3)
+ {
+ CAPTION <<- caption
+ printh(caption)
+ par(mfrow=c(nrow, ncol))
+ par(cex = 0.8)
+ par(mar = c(3, 3, 5, 0.5)) # small margins but space for right hand axis
+ par(mgp = c(1.6, 0.6, 0)) # flatten axis elements
+ oma <- par("oma") # make space for caption
+ oma[3] <- 2
+ par(oma=oma)
+ }
> do.caption <- function() # must be called _after_ first plot on new page
+ mtext(CAPTION, outer=TRUE, font=2, line=1, cex=1)
> old.par <- par(no.readonly=TRUE)
>
> library(mgcv)
Loading required package: nlme
This is mgcv 1.8-26. For overview type 'help("mgcv-package")'.
>
> for(varmod.method in c("gam", "x.gam")) {
+
+ multifigure(sprint("varmod.method=\"%s\"", varmod.method), 2, 3)
+ par(mar = c(3, 3, 2, 3)) # space for right margin axis
+
+ set.seed(6)
+ earth.mod <- earth(Volume~Girth, data=trees, nfold=3, ncross=3,
+ varmod.method=varmod.method,
+ trace=if(varmod.method %in% c("const", "lm", "power")) .3 else 0)
+ printh(sprint("varmod.method %s: summary(earth.mod)", varmod.method))
+ printh("summary(earth.mod)")
+ print(summary(earth.mod))
+
+ # summary(mgcv) prints environment as hex address which messes up the diffs
+ printh("skipping summary(mgcv::gam) etc.\n")
+
+ printh(sprint("varmod.method %s: predict(earth.mod, interval=\"pint\")", varmod.method))
+ pints <- predict(earth.mod, interval="pint")
+ print(pints)
+
+ plotmo(earth.mod$varmod, do.par=FALSE, col.response=2, clip=FALSE,
+ main="plotmo residual model",
+ xlab="x", ylab="varmod residuals")
+
+ plotmo(earth.mod, level=.90, do.par=FALSE, col.response=1, clip=FALSE,
+ main="main model plotmo Girth")
+ do.caption()
+
+ plot(earth.mod, which=3, do.par=FALSE, level=.95)
+
+ # plot.varmod
+ plot(earth.mod$varmod, do.par=FALSE, which=1:3, info=(varmod.method=="earth"))
+ }
===varmod.method="gam"
===varmod.method gam: summary(earth.mod)
===summary(earth.mod)
Call: earth(formula=Volume~Girth, data=trees,
trace=if(varmod.method%in%c("const","lm","power"))0.3els...),
nfold=3, ncross=3, varmod.method=varmod.method)
coefficients
(Intercept) 28.766764
h(13.8-Girth) -3.427802
h(Girth-13.8) 6.570747
Selected 3 of 4 terms, and 1 of 1 predictors
Termination condition: RSq changed by less than 0.001 at 4 terms
Importance: Girth
Number of terms at each degree of interaction: 1 2 (additive model)
GCV 14.20145 RSS 309.6832 GRSq 0.949137 RSq 0.9617962 CVRSq 0.9133168
Note: the cross-validation sd's below are standard deviations across folds
Cross validation: nterms 3.00 sd 0.00 nvars 1.00 sd 0.00
CVRSq sd MaxErr sd
0.913 0.045 14.3 8.71
varmod: method "gam" (mgcv package) min.sd 0.408 iter.rsq 0.319
stddev of predictions:
coefficients iter.stderr iter.stderr%
(Intercept) 4.078983 0.505543 12
s(Volume).1 0.000000 NA NA
s(Volume).2 0.000000 NA NA
s(Volume).3 0.000000 NA NA
s(Volume).4 0.000000 NA NA
s(Volume).5 0.000000 NA NA
s(Volume).6 0.000000 NA NA
s(Volume).7 0.000000 NA NA
s(Volume).8 0.000000 NA NA
s(Volume).9 2.061017 NA NA
mean smallest largest ratio
95% prediction interval 15.98932 5.66953 38.03631 6.708899
68% 80% 90% 95%
response values in prediction interval 77 90 97 97
===skipping summary(mgcv::gam) etc.
===varmod.method gam: predict(earth.mod, interval="pint")
fit lwr upr
1 9.913855 7.079089 12.74862
2 10.942195 7.845491 14.03890
3 11.627755 8.356425 14.89909
4 17.455018 12.699366 22.21067
5 18.140578 13.210301 23.07086
6 18.483359 13.465768 23.50095
7 19.168919 13.976702 24.36114
8 19.168919 13.976702 24.36114
9 19.511699 14.232169 24.79123
10 19.854479 14.487637 25.22132
11 20.197259 14.743104 25.65142
12 20.540040 14.998571 26.08151
13 20.540040 14.998571 26.08151
14 21.568380 15.764972 27.37179
15 22.596721 16.531374 28.66207
16 25.681742 18.830578 32.53291
17 25.681742 18.830578 32.53291
18 27.052863 19.852447 34.25328
19 28.423983 20.874315 35.97365
20 28.766764 21.129782 36.40374
21 30.080913 22.109191 38.05263
22 31.395063 23.088601 39.70152
23 33.366287 24.557714 42.17486
24 43.222408 31.903282 54.54153
25 45.193632 33.372396 57.01487
26 51.764379 38.269442 65.25932
27 53.078529 39.248851 66.90821
28 55.706828 41.207669 70.20599
29 56.363903 41.697374 71.03043
30 56.363903 41.697374 71.03043
31 73.447846 54.429692 92.46600
===varmod.method="x.gam"
===varmod.method x.gam: summary(earth.mod)
===summary(earth.mod)
Call: earth(formula=Volume~Girth, data=trees,
trace=if(varmod.method%in%c("const","lm","power"))0.3els...),
nfold=3, ncross=3, varmod.method=varmod.method)
coefficients
(Intercept) 28.766764
h(13.8-Girth) -3.427802
h(Girth-13.8) 6.570747
Selected 3 of 4 terms, and 1 of 1 predictors
Termination condition: RSq changed by less than 0.001 at 4 terms
Importance: Girth
Number of terms at each degree of interaction: 1 2 (additive model)
GCV 14.20145 RSS 309.6832 GRSq 0.949137 RSq 0.9617962 CVRSq 0.9133168
Note: the cross-validation sd's below are standard deviations across folds
Cross validation: nterms 3.00 sd 0.00 nvars 1.00 sd 0.00
CVRSq sd MaxErr sd
0.913 0.045 14.3 8.71
varmod: method "x.gam" (mgcv package) min.sd 0.398 iter.rsq 0.391
stddev of predictions:
coefficients iter.stderr iter.stderr%
(Intercept) 3.981114 0.452568 11
s(Volume).1 0.000000 NA NA
s(Volume).2 0.000000 NA NA
s(Volume).3 0.000000 NA NA
s(Volume).4 0.000000 NA NA
s(Volume).5 0.000000 NA NA
s(Volume).6 0.000000 NA NA
s(Volume).7 0.000000 NA NA
s(Volume).8 0.000000 NA NA
s(Volume).9 1.712690 NA NA
mean smallest largest ratio
95% prediction interval 15.60568 4.844286 31.59344 6.521795
68% 80% 90% 95%
response values in prediction interval 81 90 97 97
===skipping summary(mgcv::gam) etc.
===varmod.method x.gam: predict(earth.mod, interval="pint")
fit lwr upr
1 9.913855 7.491712 12.33600
2 10.942195 8.193843 13.69055
3 11.627755 8.661930 14.59358
4 17.455018 12.640674 22.26936
5 18.140578 13.108762 23.17240
6 18.483359 13.342806 23.62391
7 19.168919 13.810893 24.52694
8 19.168919 13.810893 24.52694
9 19.511699 14.044937 24.97846
10 19.854479 14.278981 25.42998
11 20.197259 14.513024 25.88149
12 20.540040 14.747068 26.33301
13 20.540040 14.747068 26.33301
14 21.568380 15.449199 27.68756
15 22.596721 16.151331 29.04211
16 25.681742 18.257725 33.10576
17 25.681742 18.257725 33.10576
18 27.052863 19.193900 34.91183
19 28.423983 20.130075 36.71789
20 28.766764 20.364118 37.16941
21 30.080913 21.460795 38.70103
22 31.395063 22.557472 40.23265
23 33.366287 24.202487 42.53009
24 43.222408 32.427562 54.01725
25 45.193632 34.072577 56.31469
26 51.764379 39.555960 63.97280
27 53.078529 40.652637 65.50442
28 55.706828 42.845990 68.56767
29 56.363903 43.394329 69.33348
30 56.363903 43.394329 69.33348
31 73.447846 57.651125 89.24457
> par(old.par)
>
> source("test.epilog.R")
|
