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
|
> # test.pre.R: test the "pre" package with plotmo and plotres
>
> source("test.prolog.R")
> library(pre)
> library(plotmo)
Loading required package: Formula
Loading required package: plotrix
> 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
plotmo grid: Wind Temp
10.3 75
>
> # 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)")
calculating partdep for Temp
>
> # 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
calculating partdep for Temp
>
> par(org.par)
>
> # compare pairplot and plotmo
>
> par(mfrow=c(2,3)) # 6 plots per page
>
> pairplot(pre.mod, c("Temp", "Wind"), main="pairplot")
Loading required namespace: interp
> plotmo(pre.mod, main="plotmo partdep",
+ pmethod="partdep",
+ degree1=0, degree2="Temp",
+ do.par=FALSE)
calculating partdep for Wind:Temp 01234567890
>
> # 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)
calculating apartdep for Wind:Temp 01234567890
>
> # 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 grid: Wind Temp
10.3 75
> 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 grid: Wind Temp
10.3 75
> 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")
predict.pre[3,3]:
setosa versicolor virginica
1 0.9746686 0.01299582 0.01233561
2 0.9746686 0.01299582 0.01233561
3 0.9750720 0.01300120 0.01192680
predict.pre returned multiple columns (see above) but nresponse is not specified
Use the nresponse argument to specify a column.
Example: nresponse=2
Example: nresponse="versicolor"
Error : (converted from warning) Defaulting to nresponse=1, see above messages
Got expected error from try(plotmo(pre.iris))
> options(warn=1) # print warnings as they occur
> plotmo(pre.iris, all2=TRUE, nresponse="virginica", trace=1)
importance: Petal.Length Petal.Width
stats::predict(pre.object, data.frame[3,4], type="response")
stats::fitted(object=pre.object)
fitted() was unsuccessful, will use predict() instead
assuming "Species" in the model.frame is the response, because terms(object) did not return the terms
nresponse=3 but for plotmo_y using nresponse=1 because ncol(y) == 1
assuming "Species" in the model.frame is the response, because terms(object) did not return the terms
got model response from model.frame(Species ~ .,
data=object$data, na.action="na.fail")
plotmo grid: Sepal.Length Sepal.Width Petal.Length Petal.Width
5.8 3 4.35 1.3
>
> source("test.epilog.R")
|
