inst/slowtests/test.pre.R

In plotmo: Plot a Model's Residuals, Response, and Partial Dependence Plots

# 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")