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
|
# test.pre.R: test the "pre" package with plotmo and plotres
source("test.prolog.R")
library(pre)
library(plotmo)
library(earth) # for ozone1
options(warn=1) # print warnings as they occur
data(airquality)
airq <- airquality[complete.cases(airquality), (c("Ozone", "Wind", "Temp"))]
# prevent confusion caused by integer rownames which don't match row numbers
rownames(airq) <- NULL
airq <- airq[1:50, ] # small set of data for quicker test
coef.glmnet <- glmnet:::coef.glmnet # TODO workaround required for glmnet 3.0
predict.cv.glmnet <- glmnet:::predict.cv.glmnet
set.seed(2018)
pre.mod <- pre(Ozone~., data=airq, ntrees=10) # ntrees=10 for faster test
plotres(pre.mod) # variable importance and residual plots
plotres(pre.mod, which=3, main="pre.mod residuals") # which=3 for just the residual vs fitted plot
plotmo(pre.mod) # plot model surface with background variables held at their medians
# sanity check: compare model surface to to randomForest
# (commented out to save test time)
#
# library(randomForest)
# set.seed(2018)
# rf.mod <- randomForest(Ozone~., data=airq)
# plotmo(rf.mod)
# compare singleplot and plotmo
par(mfrow=c(2,2)) # 4 plots per page
singleplot(pre.mod, varname="Temp", main="Temp\n(singleplot)")
plotmo(pre.mod,
pmethod="partdep", # plot partial dependence plot,
degree1="Temp", degree2=0, # plot only Temp, no degree2 plots
do.par=FALSE, # don't automatically set par(), use above par(mfrow)
main="Temp\n(plotmo partdep)")
# test penalty.par.val="lambda.min"
singleplot(pre.mod, varname="Temp",
main="penalty.par.val=lambda.min\n(singleplot)",
penalty.par.val="lambda.min")
plotmo(pre.mod,
pmethod="partdep",
degree1="Temp", degree2=0,
do.par=FALSE,
main="penalty.par.val=lambda.min\n(plotmo partdep)",
predict.penalty.par.val="lambda.min") # use "predict." to pass it on to predict.pre
par(org.par)
# compare pairplot and plotmo
par(mfrow=c(2,3)) # 6 plots per page
pairplot(pre.mod, c("Temp", "Wind"), main="pairplot")
plotmo(pre.mod, main="plotmo partdep",
pmethod="partdep",
degree1=0, degree2="Temp",
do.par=FALSE)
# Compare to pmethod="apartdep". An approximate partdep plot is
# faster than a full partdep plot (plotmo vignette Section 9.2).
plotmo(pre.mod, main="plotmo apartdep",
pmethod="apartdep",
degree1=0, degree2="Temp",
do.par=FALSE)
# plot contour and image plots with plotmo
plotmo(pre.mod, type2="contour",
degree1=0, degree2="Temp", do.par=FALSE)
plotmo(pre.mod, type2="image",
degree1=0, degree2="Temp", do.par=FALSE)
par(org.par)
# test gpe models
set.seed(2018)
gpe.mod <- gpe(Ozone~., data=airq,
base_learners=list(gpe_linear(), gpe_trees(), gpe_earth()))
plotmo(gpe.mod) # by default no degree2 plots because importance(gpe) not available
plotmo(gpe.mod, all2=TRUE, # force degree2 plot(s) by specifying all2=TRUE
persp.ticktype="detailed", persp.nticks=2) # optional (these get passed on to persp)
plotmo(gpe.mod, degree1=0, degree2=c("Wind", "Temp"), SHOWCALL=TRUE) # explictly specify degree2 plot
# which=3 below for only the residuals-vs-fitted plot
# optional info=TRUE to plot some extra information (RSq etc.)
plotres(gpe.mod, which=3, info=TRUE, main="gpe.mod residuals")
# multinomial response
set.seed(2018)
pre.iris <- pre(Species~., data=iris, ntrees=10) # ntrees=10 for faster testoptions(warn=2) # treat warnings as errors
options(warn=2) # treat warnings as errors
expect.err(try(plotmo(pre.iris)), "Defaulting to nresponse=1, see above messages")
options(warn=1) # print warnings as they occur
plotmo(pre.iris, all2=TRUE, nresponse="virginica", trace=1)
source("test.epilog.R")
|
