> # test.plotmo3.R: extra tests for plotmo version 3 and higher
> 
> source("test.prolog.R")
> library(earth)
Loading required package: Formula
Loading required package: plotmo
Loading required package: plotrix
> data(ozone1)
> data(etitanic)
> options(warn=1) # print warnings as they occur
> 
> # check check.numeric.scalar
> 
> xtest <- NA
> expect.err(try(plotmo:::check.numeric.scalar(xtest)), "'xtest' is NA")
Error : 'xtest' is NA
Got expected error from try(plotmo:::check.numeric.scalar(xtest))
> xtest <- NULL
> expect.err(try(plotmo:::check.numeric.scalar(xtest)), "'xtest' is NULL")
Error : 'xtest' is NULL
Got expected error from try(plotmo:::check.numeric.scalar(xtest))
> expect.err(try(plotmo:::check.numeric.scalar(NA)), "argument is NA")
Error : argument is NA
Got expected error from try(plotmo:::check.numeric.scalar(NA))
> expect.err(try(plotmo:::check.numeric.scalar(NULL)), "argument is NULL")
Error : argument is NULL
Got expected error from try(plotmo:::check.numeric.scalar(NULL))
> expect.err(try(plotmo:::check.numeric.scalar(try)), "'try' must be numeric (whereas its current class is \"function\")")
Error : 'try' must be numeric (whereas its current class is "function")
Got expected error from try(plotmo:::check.numeric.scalar(try))
> expect.err(try(plotmo:::check.numeric.scalar('try')), "\"try\" must be numeric (whereas its current class is \"character\")")
Error : "try" must be numeric (whereas its current class is "character")
Got expected error from try(plotmo:::check.numeric.scalar("try"))
> expect.err(try(plotmo:::check.numeric.scalar(NULL)), "argument is NULL")
Error : argument is NULL
Got expected error from try(plotmo:::check.numeric.scalar(NULL))
> expect.err(try(plotmo:::check.numeric.scalar(1234, min=2, max=3)), "argument=1234 but it should be between 2 and 3")
Error : argument=1234 but it should be between 2 and 3
Got expected error from try(plotmo:::check.numeric.scalar(1234, min = 2, max = 3))
> expect.err(try(plotmo:::check.numeric.scalar(0.1234, min=2, max=3)), "argument=0.1234 but it should be between 2 and 3")
Error : argument=0.1234 but it should be between 2 and 3
Got expected error from try(plotmo:::check.numeric.scalar(0.1234, min = 2, max = 3))
> 
> expect.err(try(plotmo:::check.numeric.scalar(.1234, min=2, max=3)), "argument=0.1234 but it should be between 2 and 3")
Error : argument=0.1234 but it should be between 2 and 3
Got expected error from try(plotmo:::check.numeric.scalar(0.1234, min = 2, max = 3))
> expect.err(try(plotmo:::check.numeric.scalar(+1234, min=2, max=3)), "argument=1234 but it should be between 2 and 3")
Error : argument=1234 but it should be between 2 and 3
Got expected error from try(plotmo:::check.numeric.scalar(+1234, min = 2, max = 3))
> expect.err(try(plotmo:::check.numeric.scalar(-1234, min=2, max=3)), "argument=-1234 but it should be between 2 and 3")
Error : argument=-1234 but it should be between 2 and 3
Got expected error from try(plotmo:::check.numeric.scalar(-1234, min = 2, max = 3))
> expect.err(try(plotmo:::check.numeric.scalar(+.1234, min=2, max=3)), "argument=0.1234 but it should be between 2 and 3")
Error : argument=0.1234 but it should be between 2 and 3
Got expected error from try(plotmo:::check.numeric.scalar(+0.1234, min = 2, max = 3))
> expect.err(try(plotmo:::check.numeric.scalar(-.1234, min=2, max=3)), "argument=-0.1234 but it should be between 2 and 3")
Error : argument=-0.1234 but it should be between 2 and 3
Got expected error from try(plotmo:::check.numeric.scalar(-0.1234, min = 2, max = 3))
> expect.err(try(plotmo:::check.numeric.scalar("", min=0, max=3)), "\"\" must be numeric (whereas its current class is \"character\"")
Error : "" must be numeric (whereas its current class is "character")
Got expected error from try(plotmo:::check.numeric.scalar("", min = 0, max = 3))
> 
> x.numeric.scalar <- 1234
> expect.err(try(plotmo:::check.numeric.scalar(x.numeric.scalar, min=0, max=3)), "x.numeric.scalar=1234 but it should be between 0 and 3")
Error : x.numeric.scalar=1234 but it should be between 0 and 3
Got expected error from try(plotmo:::check.numeric.scalar(x.numeric.scalar, min = 0,     max = 3))
> stopifnot(identical(plotmo:::check.numeric.scalar(x.numeric.scalar, min=2, max=1235), 1234))
> stopifnot(identical(plotmo:::check.numeric.scalar(1234, min=2, max=1235), 1234))
> 
> # check check.integer.scalar
> 
> xtest <- NA
> expect.err(try(plotmo:::check.integer.scalar(xtest)), "'xtest' is NA")
Error : 'xtest' is NA
Got expected error from try(plotmo:::check.integer.scalar(xtest))
> xtest <- NULL
> expect.err(try(plotmo:::check.integer.scalar(xtest)), "'xtest' is NULL")
Error : 'xtest' is NULL
Got expected error from try(plotmo:::check.integer.scalar(xtest))
> expect.err(try(plotmo:::check.integer.scalar(NA)), "argument is NA")
Error : argument is NA
Got expected error from try(plotmo:::check.integer.scalar(NA))
> expect.err(try(plotmo:::check.integer.scalar(NA, null.ok=TRUE)), "argument is NA")
Error : argument is NA
Got expected error from try(plotmo:::check.integer.scalar(NA, null.ok = TRUE))
> expect.err(try(plotmo:::check.integer.scalar(NULL)), "argument is NULL")
Error : argument is NULL
Got expected error from try(plotmo:::check.integer.scalar(NULL))
> expect.err(try(plotmo:::check.integer.scalar(xtest, na.ok=TRUE)), "'xtest' is NULL")
Error : 'xtest' is NULL
Got expected error from try(plotmo:::check.integer.scalar(xtest, na.ok = TRUE))
> expect.err(try(plotmo:::check.integer.scalar("xyz", na.ok=TRUE)), "\"xyz\" is a string but it should be an integer, or NA, or TRUE or FALSE")
Error : "xyz" is a string but it should be an integer, or NA, or TRUE or FALSE
Got expected error from try(plotmo:::check.integer.scalar("xyz", na.ok = TRUE))
> expect.err(try(plotmo:::check.integer.scalar("TRUE", na.ok=TRUE)), "\"TRUE\" is a string but it should be an integer, or NA, or TRUE or FALSE")
Error : "TRUE" is a string but it should be an integer, or NA, or TRUE or FALSE
Got expected error from try(plotmo:::check.integer.scalar("TRUE", na.ok = TRUE))
> stopifnot(identical(plotmo:::check.integer.scalar(TRUE), TRUE))
> stopifnot(identical(plotmo:::check.integer.scalar(NA, na.ok=TRUE), NA))
> x.integer.scalar <- 1234L
> expect.err(try(plotmo:::check.integer.scalar(x.integer.scalar, min=0, max=3)), "x.integer.scalar=1234 but it should be between 0 and 3")
Error : x.integer.scalar=1234 but it should be between 0 and 3
Got expected error from try(plotmo:::check.integer.scalar(x.integer.scalar, min = 0,     max = 3))
> stopifnot(identical(plotmo:::check.integer.scalar(x.integer.scalar, min=2, max=1235), 1234L))
> stopifnot(identical(plotmo:::check.integer.scalar(1234, min=2, max=1235), 1234))
> stopifnot(identical(plotmo:::check.integer.scalar(x.integer.scalar, min=2, max=1235), 1234L))
> stopifnot(identical(plotmo:::check.integer.scalar(1234, min=2, max=1235), 1234))
> xtest <- 1.234
> expect.err(try(plotmo:::check.integer.scalar(xtest, min=0, max=3)), "xtest=1.234 but it should be an integer, or TRUE or FALSE")
Error : xtest=1.234 but it should be an integer, or TRUE or FALSE
Got expected error from try(plotmo:::check.integer.scalar(xtest, min = 0, max = 3))
> 
> # check check.vec
> xtest <- "x"
> expect.err(try(plotmo:::check.vec(xtest, "xtest", na.ok=TRUE)), "'xtest' is not numeric")
Error : 'xtest' is not numeric
Got expected error from try(plotmo:::check.vec(xtest, "xtest", na.ok = TRUE))
> xtest <- as.double(NA)
> print(plotmo:::check.vec(xtest, "xtest", na.ok=TRUE))
NULL
> xtest <- as.double(1:3)
> print(plotmo:::check.vec(xtest, "xtest", na.ok=TRUE))
NULL
> xtest <- c(1,2,3,1/0,5,6,7)
> expect.err(try(plotmo:::check.vec(xtest, "xtest", na.ok=TRUE)), "non-finite value in xtest")
Error : non-finite value in xtest
       xtest[4] is Inf
Got expected error from try(plotmo:::check.vec(xtest, "xtest", na.ok = TRUE))
> xtest <- c(1,2,3,NA,5,6,7)
> expect.err(try(plotmo:::check.vec(xtest, "xtest")), "NA in xtest")
Error : NA in xtest
       xtest[4] is NA
Got expected error from try(plotmo:::check.vec(xtest, "xtest"))
> xtest <- c(1,2,3)
> expect.err(try(plotmo:::check.vec(xtest, "xtest", expected.len=2)), "'xtest' has the wrong length 3, expected 2")
Error : 'xtest' has the wrong length 3, expected 2
Got expected error from try(plotmo:::check.vec(xtest, "xtest", expected.len = 2))
> print(plotmo:::check.vec(c(TRUE, FALSE), "c(TRUE, FALSE)"))
NULL
> 
> plotmo1 <- function(object, ..., trace=0, SHOWCALL=TRUE, caption=NULL) {
+     if(is.null(caption))
+         caption <- paste(deparse(substitute(object)), collapse=" ")
+     call <- match.call(expand.dots=TRUE)
+     call <- strip.space(paste(deparse(substitute(call)), collapse=" "))
+     printf("%s\n", call)
+     plotmo(object, trace=trace, SHOWCALL=SHOWCALL, caption=caption, ...)
+ }
> plotres1 <- function(object, ..., trace=0, SHOWCALL=TRUE, caption=NULL) {
+     if(is.null(caption))
+         caption <- paste(deparse(substitute(object)), collapse=" ")
+     call <- match.call(expand.dots=TRUE)
+     call <- strip.space(paste(deparse(substitute(call)), collapse=" "))
+     printf("%s\n", call)
+     plotres(object, trace=trace, SHOWCALL=SHOWCALL, caption=caption, ...)
+ }
> # basic tests of plotmo on abbreviated titanic data
> 
> get.tita <- function()
+ {
+     tita <- etitanic
+     pclass <- as.character(tita$pclass)
+     # change the order of the factors so not alphabetical
+     pclass[pclass == "1st"] <- "first"
+     pclass[pclass == "2nd"] <- "class2"
+     pclass[pclass == "3rd"] <- "classthird"
+     tita$pclass <- factor(pclass, levels=c("class2", "classthird", "first"))
+     # log age is so we have a continuous predictor even when model is age~.
+     set.seed(2015)
+     tita$logage <- log(tita$age) + rnorm(nrow(tita))
+     tita$parch <- NULL
+     # by=12 gives us a small fast model with an additive and a interaction term
+     tita[seq(1, nrow(etitanic), by=12), ]
+ }
> tita <- get.tita()
> 
> mod.lm.age <- lm(age~., data=tita)
> plotmo1(mod.lm.age)
plotmo1(object=mod.lm.age)
 plotmo grid:    pclass survived  sex sibsp  logage
             classthird        0 male     0 3.06991
> plotmo1(mod.lm.age, level=.95)
plotmo1(object=mod.lm.age,level=0.95)
 plotmo grid:    pclass survived  sex sibsp  logage
             classthird        0 male     0 3.06991
> plotmo1(mod.lm.age, level=.95, col.resp=3)
plotmo1(object=mod.lm.age,level=0.95,col.resp=3)
 plotmo grid:    pclass survived  sex sibsp  logage
             classthird        0 male     0 3.06991
> 
> sexn <- as.numeric(tita$sex)
> mod.lm.sexn <- lm(sexn~.-sex, data=tita)
> plotmo1(mod.lm.sexn)
plotmo1(object=mod.lm.sexn)
 plotmo grid:    pclass survived  sex age sibsp  logage
             classthird        0 male  30     0 3.06991
> plotmo1(mod.lm.sexn, level=.95)
plotmo1(object=mod.lm.sexn,level=0.95)
 plotmo grid:    pclass survived  sex age sibsp  logage
             classthird        0 male  30     0 3.06991
> 
> set.seed(2020)
> mod.earth.age <- earth(age~., data=tita, degree=2, nfold=3, ncross=3, varmod.method="lm")
> plotmo1(mod.earth.age)
plotmo1(object=mod.earth.age)
 plotmo grid:    pclass survived  sex sibsp  logage
             classthird        0 male     0 3.06991
> plotmo1(mod.earth.age, level=.9, degree2=0)
plotmo1(object=mod.earth.age,level=0.9,degree2=0)
 plotmo grid:    pclass survived  sex sibsp  logage
             classthird        0 male     0 3.06991
> 
> # tita[,4] is age
> set.seed(2020)
> mod.earth.tita.age <- earth(tita[,-4], tita[,4], degree=2, nfold=3, ncross=3, trace=.5, varmod.method="lm")
Model with pmethod="backward": GRSq 0.335 RSq 0.512 nterms 6
CV fold  1.1  CVRSq -0.047   n.oof 58 34%   n.infold.nz 58 100%   n.oof.nz 30 100%
CV fold  1.2  CVRSq -0.022   n.oof 59 33%   n.infold.nz 59 100%   n.oof.nz 29 100%
CV fold  1.3  CVRSq -0.045   n.oof 59 33%   n.infold.nz 59 100%   n.oof.nz 29 100%
CV fold  2.1  CVRSq  0.133   n.oof 58 34%   n.infold.nz 58 100%   n.oof.nz 30 100%
CV fold  2.2  CVRSq  0.338   n.oof 59 33%   n.infold.nz 59 100%   n.oof.nz 29 100%
CV fold  2.3  CVRSq  0.149   n.oof 59 33%   n.infold.nz 59 100%   n.oof.nz 29 100%
CV fold  3.1  CVRSq  0.419   n.oof 58 34%   n.infold.nz 58 100%   n.oof.nz 30 100%
CV fold  3.2  CVRSq  0.107   n.oof 59 33%   n.infold.nz 59 100%   n.oof.nz 29 100%
CV fold  3.3  CVRSq  0.307   n.oof 59 33%   n.infold.nz 59 100%   n.oof.nz 29 100%
CV all        CVRSq  0.149                  n.infold.nz 88 100%


varmod method="lm" rmethod="hc12" lambda=1 exponent=1 conv=1 clamp=0.1 minspan=-3:
     iter weight.ratio coefchange% (Intercept) tita[, 4]
        1          1.4         0.0          13    -0.032
        2          1.2         7.1          12    -0.018
        3          1.3         3.0          13    -0.024
        4          1.3         1.2          13    -0.022
        5          1.3         0.5          13    -0.023

> cat("\nsummary(mod.earth.tita.age)\n")

summary(mod.earth.tita.age)
> print(summary(mod.earth.tita.age))
Call: earth(x=tita[,-4], y=tita[,4], trace=0.5, degree=2, nfold=3, ncross=3,
            varmod.method="lm")

                            coefficients
(Intercept)                    25.664968
pclassfirst                     9.028974
h(sibsp-1)                    -12.096706
h(1.68119-logage)              -7.502937
sexmale * h(logage-2.48137)     5.062358
sibsp * h(logage-1.68119)       3.280947

Selected 6 of 14 terms, and 4 of 6 predictors
Termination condition: Reached nk 21
Importance: logage, sexmale, pclassclassthird-unused, sibsp, pclassfirst, ...
Number of terms at each degree of interaction: 1 3 2
GCV 174.7603  RSS 11022.31  GRSq 0.335155  RSq 0.5124778  CVRSq 0.1487371

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

Cross validation:   nterms 3.89 sd 1.05    nvars 3.22 sd 0.97

     CVRSq    sd     MaxErr   sd
     0.149 0.174      -39.1 32.3

varmod: method "lm"    min.sd 1.49    iter.rsq 0.001

stddev of predictions:
            coefficients iter.stderr iter.stderr%
(Intercept)   15.7287403     2.77398           18
tita[, 4]     -0.0283536   0.0837154          295

                              mean   smallest    largest     ratio
95% prediction interval   58.24711   55.23254   62.56685   1.13279

                                         68%    80%    90%    95% 
response values in prediction interval   84     90     97     99  
> plotmo1(mod.earth.tita.age)
plotmo1(object=mod.earth.tita.age)
 plotmo grid:    pclass survived  sex sibsp  logage
             classthird        0 male     0 3.06991
> plotmo1(mod.earth.tita.age, level=.9, degree2=0)
plotmo1(object=mod.earth.tita.age,level=0.9,degree2=0)
 plotmo grid:    pclass survived  sex sibsp  logage
             classthird        0 male     0 3.06991
> 
> set.seed(2020)
> a.earth.sex <- earth(sex~., data=tita, degree=2, nfold=3, ncross=3, varmod.method="lm")
> plotmo1(a.earth.sex)
plotmo1(object=a.earth.sex)
 plotmo grid:    pclass survived age sibsp  logage
             classthird        0  30     0 3.06991
> plotmo1(a.earth.sex, level=.9)
plotmo1(object=a.earth.sex,level=0.9)
 plotmo grid:    pclass survived age sibsp  logage
             classthird        0  30     0 3.06991
> plotmo1(a.earth.sex, type="class")
plotmo1(object=a.earth.sex,type="class")
 plotmo grid:    pclass survived age sibsp  logage
             classthird        0  30     0 3.06991
> expect.err(try(plotmo1(a.earth.sex, level=.9, degree2=0, type="class")), "predicted values are strings")
plotmo1(object=a.earth.sex,level=0.9,degree2=0,type="class")
Error : the level argument is not allowed when the predicted values are strings
Got expected error from try(plotmo1(a.earth.sex, level = 0.9, degree2 = 0, type = "class"))
> 
> # tita[,3] is sex
> set.seed(2020)
> mod.earth.tita <- earth(tita[,-3], tita[,3], degree=2, nfold=3, ncross=3, varmod.method="lm")
> plotmo1(mod.earth.tita)
plotmo1(object=mod.earth.tita)
 plotmo grid:    pclass survived age sibsp  logage
             classthird        0  30     0 3.06991
> plotmo1(mod.earth.tita, level=.9, degree2=0)
plotmo1(object=mod.earth.tita,level=0.9,degree2=0)
 plotmo grid:    pclass survived age sibsp  logage
             classthird        0  30     0 3.06991
> plotmo1(mod.earth.tita, type="class")
plotmo1(object=mod.earth.tita,type="class")
 plotmo grid:    pclass survived age sibsp  logage
             classthird        0  30     0 3.06991
> expect.err(try(plotmo1(mod.earth.tita, level=.9, degree2=0, type="class")), "predicted values are strings")
plotmo1(object=mod.earth.tita,level=0.9,degree2=0,type="class")
Error : the level argument is not allowed when the predicted values are strings
Got expected error from try(plotmo1(mod.earth.tita, level = 0.9, degree2 = 0, type = "class"))
> 
> set.seed(2020)
> mod.earth.sex <- earth(sex~., data=tita, degree=2, nfold=3, ncross=3, varmod.method="earth", glm=list(family=binomial))
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
> plotmo1(mod.earth.sex)
plotmo1(object=mod.earth.sex)
 plotmo grid:    pclass survived age sibsp  logage
             classthird        0  30     0 3.06991
> plotmo1(mod.earth.sex, type="link")
plotmo1(object=mod.earth.sex,type="link")
 plotmo grid:    pclass survived age sibsp  logage
             classthird        0  30     0 3.06991
> plotmo1(mod.earth.sex, type="class")
plotmo1(object=mod.earth.sex,type="class")
 plotmo grid:    pclass survived age sibsp  logage
             classthird        0  30     0 3.06991
> plotmo1(mod.earth.sex, level=.9, type="earth")
plotmo1(object=mod.earth.sex,level=0.9,type="earth")
 plotmo grid:    pclass survived age sibsp  logage
             classthird        0  30     0 3.06991
> 
> # tita[,3] is sex
> set.seed(2020)
> mod.earth.tita <- earth(tita[,-3], tita[,3], degree=2, nfold=3, ncross=3, varmod.method="earth", glm=list(family=binomial))
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
> plotmo1(mod.earth.tita)
plotmo1(object=mod.earth.tita)
 plotmo grid:    pclass survived age sibsp  logage
             classthird        0  30     0 3.06991
> plotmo1(mod.earth.tita, type="link")
plotmo1(object=mod.earth.tita,type="link")
 plotmo grid:    pclass survived age sibsp  logage
             classthird        0  30     0 3.06991
> plotmo1(mod.earth.tita, type="class")
plotmo1(object=mod.earth.tita,type="class")
 plotmo grid:    pclass survived age sibsp  logage
             classthird        0  30     0 3.06991
> plotmo1(mod.earth.tita, level=.9, type="earth")
plotmo1(object=mod.earth.tita,level=0.9,type="earth")
 plotmo grid:    pclass survived age sibsp  logage
             classthird        0  30     0 3.06991
> 
> # check factor handling when factors are not ordered alphabetically
> tita.orgpclass <- etitanic[seq(1, nrow(etitanic), by=12), ]
> tita  <- get.tita()
> tita$logage <- NULL
> tita.orgpclass$parch <- NULL
> stopifnot(names(tita.orgpclass) == names(tita))
> a.tita.orgpclass <- earth(pclass~., degree=2, data=tita.orgpclass)
> a.tita           <- earth(pclass~., degree=2, data=tita)
> options(warn=2) # treat warnings as errors
> expect.err(try(plotmo(a.tita)), "Defaulting to nresponse=1, see above messages")

predict.earth[88,3]:
       class2 classthird      first
1   0.3179514  0.3141272 0.36792134
2   0.3179514  0.3141272 0.36792134
3   0.2400614  0.6754849 0.08445368
... 0.2180022  0.5645160 0.21748179
88  0.2400614  0.6754849 0.08445368

predict.earth returned multiple columns (see above) but nresponse is not specified
    Use the nresponse argument to specify a column.
         Example: nresponse=2
         Example: nresponse="classthird"

Error : (converted from warning) Defaulting to nresponse=1, see above messages
Got expected error from try(plotmo(a.tita))
> options(warn=1)
> # following two graphs should be identical
> plotmo1(a.tita.orgpclass, nresponse="1st",   all1=T, col.resp=3, type2="im")
plotmo1(object=a.tita.orgpclass,nresponse="1st",all1=T,col.resp=3,type2="im")
 plotmo grid:    survived  sex age sibsp
                        0 male  30     0
> plotmo1(a.tita,           nresponse="first", all1=T, col.resp=3, type2="im")
plotmo1(object=a.tita,nresponse="first",all1=T,col.resp=3,type2="im")
 plotmo grid:    survived  sex age sibsp
                        0 male  30     0
> # following two graphs should be identical
> plotmo1(a.tita.orgpclass, nresponse="2nd",    all1=T)
plotmo1(object=a.tita.orgpclass,nresponse="2nd",all1=T)
 plotmo grid:    survived  sex age sibsp
                        0 male  30     0
> plotmo1(a.tita,           nresponse="class2", all1=T)
plotmo1(object=a.tita,nresponse="class2",all1=T)
 plotmo grid:    survived  sex age sibsp
                        0 male  30     0
> 
> tita  <- get.tita()
> mod.earth.pclass <- earth(pclass~., data=tita, degree=2)
> options(warn=2) # treat warnings as errors
> expect.err(try(plotmo1(mod.earth.pclass)), "Defaulting to nresponse=1, see above messages")
plotmo1(object=mod.earth.pclass)

predict.earth[88,3]:
       class2 classthird     first
1   0.3197580  0.2991394 0.3811026
2   0.3197580  0.2991394 0.3811026
3   0.2490258  0.6472095 0.1037648
... 0.1984114  0.5220475 0.2795411
88  0.2490258  0.6472095 0.1037648

predict.earth returned multiple columns (see above) but nresponse is not specified
    Use the nresponse argument to specify a column.
         Example: nresponse=2
         Example: nresponse="classthird"

Error : (converted from warning) Defaulting to nresponse=1, see above messages
Got expected error from try(plotmo1(mod.earth.pclass))
> options(warn=1)
> plotmo1(mod.earth.pclass, nresponse="fi")
plotmo1(object=mod.earth.pclass,nresponse="fi")
 plotmo grid:    survived  sex age sibsp  logage
                        0 male  30     0 3.06991
> plotmo1(mod.earth.pclass, nresponse="first")
plotmo1(object=mod.earth.pclass,nresponse="first")
 plotmo grid:    survived  sex age sibsp  logage
                        0 male  30     0 3.06991
> plotmo1(mod.earth.pclass, nresponse=3)
plotmo1(object=mod.earth.pclass,nresponse=3)
 plotmo grid:    survived  sex age sibsp  logage
                        0 male  30     0 3.06991
> plotmo1(mod.earth.pclass, type="class")
plotmo1(object=mod.earth.pclass,type="class")
 plotmo grid:    survived  sex age sibsp  logage
                        0 male  30     0 3.06991
> plotmo1(mod.earth.pclass, nresponse=1,
+        type="class", grid.levels=list(sex="fem"),
+        smooth.col="indianred", smooth.lwd=2,
+        pt.col=as.numeric(tita$pclass)+1,
+        pt.pch=1)
plotmo1(object=mod.earth.pclass,nresponse=1,type="class",grid.levels=list(sex="fem"),smooth.col="indianred",smooth.lwd=2,pt.col=as.numeric(tita$pclass)+1,pt.pch=1)
 plotmo grid:    survived    sex age sibsp  logage
                        0 female  30     0 3.06991
> 
> # tita[,1] is pclass
> mod.earth.tita <- earth(tita[,-1], tita[,1], degree=2)
> options(warn=2) # treat warnings as errors
> expect.err(try(plotmo1(mod.earth.tita)), "Defaulting to nresponse=1, see above messages")
plotmo1(object=mod.earth.tita)

predict.earth[88,3]:
       class2 classthird     first
1   0.3197580  0.2991394 0.3811026
2   0.3197580  0.2991394 0.3811026
3   0.2490258  0.6472095 0.1037648
... 0.1984114  0.5220475 0.2795411
88  0.2490258  0.6472095 0.1037648

predict.earth returned multiple columns (see above) but nresponse is not specified
    Use the nresponse argument to specify a column.
         Example: nresponse=2
         Example: nresponse="classthird"

Error : (converted from warning) Defaulting to nresponse=1, see above messages
Got expected error from try(plotmo1(mod.earth.tita))
> options(warn=1)
> plotmo1(mod.earth.tita, nresponse="first")
plotmo1(object=mod.earth.tita,nresponse="first")
 plotmo grid:    survived  sex age sibsp  logage
                        0 male  30     0 3.06991
> plotmo1(mod.earth.tita, type="class")
plotmo1(object=mod.earth.tita,type="class")
 plotmo grid:    survived  sex age sibsp  logage
                        0 male  30     0 3.06991
> 
> mod.earth.pclass2 <- earth(pclass~., data=tita, degree=2, glm=list(family=binomial))
> # expect.err(try(plotmo1(mod.earth.pclass2)), "nresponse is not specified")
> plotmo1(mod.earth.pclass2, nresponse=3)
plotmo1(object=mod.earth.pclass2,nresponse=3)
 plotmo grid:    survived  sex age sibsp  logage
                        0 male  30     0 3.06991
> plotmo1(mod.earth.pclass2, type="link", nresponse=3)
plotmo1(object=mod.earth.pclass2,type="link",nresponse=3)
 plotmo grid:    survived  sex age sibsp  logage
                        0 male  30     0 3.06991
> plotmo1(mod.earth.pclass2, type="class")
plotmo1(object=mod.earth.pclass2,type="class")
 plotmo grid:    survived  sex age sibsp  logage
                        0 male  30     0 3.06991
> 
> # tita[,1] is pclass
> mod.earth.tita <- earth(tita[,-1], tita[,1], degree=2, glm=list(family=binomial))
> plotmo1(mod.earth.tita, nresponse=3)
plotmo1(object=mod.earth.tita,nresponse=3)
 plotmo grid:    survived  sex age sibsp  logage
                        0 male  30     0 3.06991
> plotmo1(mod.earth.tita, type="link", nresponse=3)
plotmo1(object=mod.earth.tita,type="link",nresponse=3)
 plotmo grid:    survived  sex age sibsp  logage
                        0 male  30     0 3.06991
> plotmo1(mod.earth.tita, type="class")
plotmo1(object=mod.earth.tita,type="class")
 plotmo grid:    survived  sex age sibsp  logage
                        0 male  30     0 3.06991
> 
> # plotmo vignette examples
> 
> # use a small set of variables for illustration
> printf("library(earth)\n")
library(earth)
> library(earth) # for ozone1 data
> data(ozone1)
> oz <- ozone1[, c("O3", "humidity", "temp", "ibt")]
> 
> lm.model.vignette <- lm(O3 ~ humidity + temp*ibt, data=oz) # linear model
> plotmo1(lm.model.vignette, pt.col="gray", nrug=-1)
plotmo1(object=lm.model.vignette,pt.col="gray",nrug=-1)
 plotmo grid:    humidity temp   ibt
                       64   62 167.5
> plotmo1(lm.model.vignette, level=.9)
plotmo1(object=lm.model.vignette,level=0.9)
 plotmo grid:    humidity temp   ibt
                       64   62 167.5
> 
> printf("library(mda)\n")
library(mda)
> library(mda)
Loading required package: class
Loaded mda 0.5-5

> mars.model.vignette1 <- mars(oz[,-1], oz[,1], degree=2)
> plotmo1(mars.model.vignette1)
plotmo1(object=mars.model.vignette1)
 plotmo grid:    humidity temp   ibt
                       64   62 167.5
> plotres1(mars.model.vignette1)
plotres1(object=mars.model.vignette1)
> mars.model.vignette2 <- mars(oz[,-1,drop=FALSE], oz[,1,drop=FALSE], degree=2)
> plotmo1(mars.model.vignette2)
plotmo1(object=mars.model.vignette2)
 plotmo grid:    humidity temp   ibt
                       64   62 167.5
> # TODO causes Error in lm.fit(object$x, y, singular.ok = FALSE) : (list) object cannot be coerced to type 'double'
> #      although still works
> #      the error is mars.to.earth try(hatvalues.lm.fit(lm.fit(object$x, y, singular.ok=FALSE)))
> plotres1(mars.model.vignette2, trace=1)
plotres1(object=mars.model.vignette2,trace=1)
stats::residuals(object=mars.object, type="response")
stats::fitted(object=mars.object)
got model response from getCall(object)$y
calling mars.to.earth (needed for the model selection plot)

training rsq 0.76
> 
> printf("library(rpart)\n")
library(rpart)
> library(rpart)                                          # rpart
> rpart.model.vignette <- rpart(O3 ~ ., data=oz)
> plotmo1(rpart.model.vignette, all2=TRUE)
plotmo1(object=rpart.model.vignette,all2=TRUE)
 plotmo grid:    humidity temp   ibt
                       64   62 167.5
> expect.err(try(plotmo1(rpart.model.vignette, level=.9)), "the level argument is not supported for \"rpart\" objects")
plotmo1(object=rpart.model.vignette,level=0.9)
Error : the level argument is not supported for "rpart" objects
Got expected error from try(plotmo1(rpart.model.vignette, level = 0.9))
> 
> # commented out because is slow and already tested in test.non.earth.R
> # printf("library(randomForest)\n")
> # library(randomForest)                                   # randomForest
> # rf.model.vignette <- randomForest(O3~., data=oz)
> # plotmo1(rf.model.vignette)
> # partialPlot(rf.model.vignette, oz, temp) # compare to partial-dependence plot
> 
> printf("library(gbm)\n")
library(gbm)
> library(gbm)                                            # gbm
Loaded gbm 2.2.2
This version of gbm is no longer under development. Consider transitioning to gbm3, https://github.com/gbm-developers/gbm3
> set.seed(2016)
> gbm.model.vignette <- gbm(O3~., data=oz, dist="gaussian", inter=2, n.trees=100)
> # commented out following because they always take the whole page
> # plot(gbm.model.vignette, i.var=2) # compare to partial-dependence plots
> # plot(gbm.model.vignette, i.var=c(2,3))
> set.seed(2016)
> plotmo1(gbm.model.vignette, caption="gbm.model.vignette")
plotmo1(object=gbm.model.vignette,caption="gbm.model.vignette")
 plotmo grid:    humidity temp   ibt
                       64   62 167.5
> 
> # commented out because is slow and already tested elsewhere
> # printf("library(mgcv)\n")
> # library(mgcv)                                           # gam
> # gam.model.vignette <- gam(O3 ~ s(humidity)+s(temp)+s(ibt)+s(temp,ibt), data=oz)
> # plotmo1(gam.model.vignette, level=.95, all2=TRUE)
> 
> printf("library(nnet)\n")
library(nnet)
> library(nnet)                                           # nnet
> set.seed(4)
> nnet.model.vignette <- nnet(O3~., data=scale(oz), size=2, decay=0.01, trace=FALSE)
> plotmo1(nnet.model.vignette, type="raw", all2=T)
plotmo1(object=nnet.model.vignette,type="raw",all2=T)
 plotmo grid:    humidity       temp        ibt
                0.2954793 0.01697621 0.08267399
> 
> printf("library(MASS)\n")
library(MASS)
> library(MASS)                                           # qda
> lcush <- data.frame(Type=as.numeric(Cushings$Type),log(Cushings[,1:2]))
> lcush <- lcush[1:21,]
> qda.model.vignette <- qda(Type~., data=lcush)
> plotmo1(qda.model.vignette, type="class", all2=TRUE,
+        type2="contour", ngrid2=100, contour.nlevels=2, contour.drawlabels=FALSE,
+        pt.col=as.numeric(lcush$Type)+1,
+        pt.pch=as.character(lcush$Type))
plotmo1(object=qda.model.vignette,type="class",all2=TRUE,type2="contour",ngrid2=100,contour.nlevels=2,contour.drawlabels=FALSE,pt.col=as.numeric(lcush$Type)+1,pt.pch=as.character(lcush$Type))
 plotmo grid:    Tetrahydrocortisone Pregnanetriol
                             2.04122     0.1823216
> 
> # miscellaneous other examples
> 
> tita <- get.tita()
> 
> mod.glm.sex <- glm(sex~., data=tita, family=binomial)
> plotmo1(mod.glm.sex, pt.col=as.numeric(tita$pclass)+1, do.par=2)
plotmo1(object=mod.glm.sex,pt.col=as.numeric(tita$pclass)+1,do.par=2)
 plotmo grid:    pclass survived age sibsp  logage
             classthird        0  30     0 3.06991
> empty.plot()
> 
> # glm binomial model: compare error band on plotmo versus manual plot
> 
> wasp <- read.table("darlingtonia.tab", header=T)
> iorder = order(wasp$leafHeight, wasp$visited)
> wasp = wasp[iorder,]
> iseq = seq(1, nrow(wasp), by = 2)
> dat = wasp[iseq, ]
> newdat = wasp[iseq+1, ]
> 
> print(summary(dat))
   leafHeight      visited      
 Min.   :14.0   Min.   :0.0000  
 1st Qu.:32.0   1st Qu.:0.0000  
 Median :43.0   Median :0.0000  
 Mean   :44.9   Mean   :0.1905  
 3rd Qu.:58.0   3rd Qu.:0.0000  
 Max.   :83.0   Max.   :1.0000  
> set.seed(2026)
> mod <- glm(visited ~ leafHeight, data=dat, family=binomial())
> print(summary(mod))

Call:
glm(formula = visited ~ leafHeight, family = binomial(), data = dat)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -4.71115    2.14347  -2.198   0.0280 *
leafHeight   0.06385    0.03652   1.748   0.0805 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 20.450  on 20  degrees of freedom
Residual deviance: 16.565  on 19  degrees of freedom
AIC: 20.565

Number of Fisher Scoring iterations: 5

> set.seed(2026)
> level95 = .95
> plotmo(mod, do.par=FALSE, trace=2,
+     level=level95, level.shade="lightgray",
+     # nrug=TRUE,
+     grid.col="lightgray",
+     pt.col=ifelse(dat$visited, "pink", "lightgreen"), jitter=0.1,
+     main="plotmo: visited vs leafHeight")
plotmo trace 2: plotmo(object=mod,
                       pt.col=ifelse(dat$visited,"pink","lightgreen"),
                       jitter=0.1, level=level95, grid.col="lightgray",
                       do.par=FALSE, trace=2, level.shade="lightgray",
                       main="plotmo: visited vs leafHeight")
--get.model.env for object with class glm
object call is glm(formula=visited~leafHeight, family=binomial(), data=dat)
using the environment saved in $terms of the glm model: R_GlobalEnv
--plotmo_prolog for glm object 'mod'
--plotmo_x for glm object

get.object.x:
object$x is NULL (and it has no colnames)

object call is glm(formula=visited~leafHeight, family=binomial(), data=dat)

get.x.from.model.frame:
formula(object) is visited ~ leafHeight
naked formula is the same
formula is valid, now looking for data for the model.frame
object$model is usable and has column names visited leafHeight
x=model.frame[,-1] is usable and has column name leafHeight
plotmo_x returned[21,1]:
    leafHeight
42          14
11          18
41          20
...         29
26          83

----Metadata: plotmo_predict with nresponse=NULL and newdata=NULL
plotmo_predict with NULL newdata (nrows=3), using plotmo_x to get the data
--plotmo_x for glm object

get.object.x:
object$x is NULL (and it has no colnames)

object call is glm(formula=visited~leafHeight, family=binomial(), data=dat)

get.x.from.model.frame:
formula(object) is visited ~ leafHeight
naked formula is the same
formula is valid, now looking for data for the model.frame
object$model is usable and has column names visited leafHeight
x=model.frame[,-1] is usable and has column name leafHeight
plotmo_x returned[21,1]:
    leafHeight
42          14
11          18
41          20
...         29
26          83
will use the above data instead of newdata=NULL for predict.glm
stats::predict(glm.object, data.frame[3,1], type="response")
predict returned[3,1] with no column names:
             
42 0.02151400
11 0.02760091
41 0.03124279
predict after processing with nresponse=NULL is [3,1] with no column names:
             
42 0.02151400
11 0.02760091
41 0.03124279

----Metadata: plotmo_fitted with nresponse=NULL
stats::fitted(object=glm.object)
fitted(object) returned[21,1] with no column names:
              
42  0.02151400
11  0.02760091
41  0.03124279
... 0.05418670
26  0.64292551
fitted(object) after processing with nresponse=NULL is [21,1] with no column names:
              
42  0.02151400
11  0.02760091
41  0.03124279
... 0.05418670
26  0.64292551

----Metadata: plotmo_y with nresponse=NULL
--plotmo_y with nresponse=NULL for glm object

get.object.y:
object$y is usable but without colnames so we will keep on searching

object call is glm(formula=visited~leafHeight, family=binomial(), data=dat)

get.y.from.model.frame:
formula(object) is visited ~ leafHeight
formula is valid, now looking for data for the model.frame
object$model is usable and has column names visited leafHeight
y=model.frame[,1] is usable and has column name visited
plotmo_y returned[21,1]:
    visited
42        0
11        0
41        0
...       0
26        0
plotmo_y after processing with nresponse=NULL is [21,1]:
    visited
42        0
11        0
41        0
...       0
26        0
converted nresponse=NA to nresponse=1
nresponse=1 (was NA) ncol(fitted) 1 ncol(predict) 1 ncol(y) 1

----Metadata: plotmo_y with nresponse=1
--plotmo_y with nresponse=1 for glm object

get.object.y:
object$y is usable but without colnames so we will keep on searching

object call is glm(formula=visited~leafHeight, family=binomial(), data=dat)

get.y.from.model.frame:
formula(object) is visited ~ leafHeight
formula is valid, now looking for data for the model.frame
object$model is usable and has column names visited leafHeight
y=model.frame[,1] is usable and has column name visited
got model response from object$model
plotmo_y returned[21,1]:
    visited
42        0
11        0
41        0
...       0
26        0
plotmo_y after processing with nresponse=1 is [21,1]:
    visited
1         0
2         0
3         0
...       0
21        0
got response name "visited" from yfull
resp.levs is NULL

----Metadata: done

number of x values: leafHeight 21

----plotmo_singles for glm object
singles: 1 leafHeight 

----plotmo_pairs for glm object
formula(object) returned visited ~ leafHeight
formula.vars "leafHeight"
term.labels "leafHeight"
plotmo_pairs_from_term_labels
term.labels: "leafHeight" "leafHeight"
pred.names:  "leafHeight"
considering leafHeight
considering leafHeight
no pairs

----Figuring out ylim
ylim c(-0.1, 1.1)    clip TRUE

--plot.degree1(draw.plot=TRUE)
degree1 plot1 (pmethod "plotmo") variable leafHeight
newdata[50,1]:
    leafHeight
1     14.00000
2     15.40816
3     16.81633
...   18.22449
50    83.00000
stats::predict(glm.object, data.frame[50,1], type="response")
predict returned[50,1] with no column names:
              
1   0.02151400
2   0.02349031
3   0.02564341
... 0.02798820
50  0.64292551
predict after processing with nresponse=1 is [50,1]:
       predict
1   0.02151400
2   0.02349031
3   0.02564341
... 0.02798820
50  0.64292551
plotmo_pint for glm object
plotmo: plot_degree1: get.degree1.data: plotmo_pint: plotmo.pint: plotmo.pint.glm:
get.quant_plotmo2(level 0.95 df 19)   level2 0.975   quant 2.09302
Warning: the level argument may not work correctly on glm objects built with weights
prediction intervals[50,2]:
        cint.lwr  cint.upr
1   0.0006797761 0.4154360
2   0.0008215046 0.4130813
3   0.0009922489 0.4108536
... 0.0011977732 0.4087640
50  0.1370955421 0.9532826
graphics::plot.default(x=c(14,15.41,16.8...), y=c(0.0215,0.0235...), type="n",
                       main="plotmo: visited vs leafHeight", xlab="", ylab="",
                       xaxt="s", yaxt="s", xlim=c(13.99,83.01),
                       ylim=c(-0.1,1.1))
Will shift and scale displayed points specified by pt.col: yshift -0 yscale 1
> 
> # wrong error bars (used prior to plotmo 3.7.0)
> pred = predict(mod, newdata=newdat, type='response', se.fit=TRUE)
> yhat = pred$fit
> wrong_se = pred$se.fit
> wrong_lwr = yhat - (2 * wrong_se)
> wrong_upr = yhat + (2 * wrong_se)
> points(newdat$leafHeight, wrong_lwr, pch=".", col="red")
> points(newdat$leafHeight, wrong_upr, pch=".",, col="red")
> 
> # correct error bars (fixed in plotmo 3.7.0)
> ilink = family(mod)$linkinv
> pred_link = predict(mod, newdata=newdat, type="link", se.fit=TRUE)
> q = qt((1 - level95) / 2, df = df.residual(mod), lower.tail = FALSE)
> printf("darlingtonia: level95 %.4f q %.4f\n", level95, q)
darlingtonia: level95 0.9500 q 2.0930
> lwr = ilink(pred_link$fit - (q * pred_link$se.fit))
> upr = ilink(pred_link$fit + (q * pred_link$se.fit))
> 
> points(newdat$leafHeight, lwr, pch="o")
> points(newdat$leafHeight, upr, pch="o")
> 
> legend(x="topleft",
+        legend=c("correct (plotmo 3.7.0)", "incorrect (prior to plotmo 3.7.0)"),
+        col=c(1,2),
+        pch=sprintf("%s%s", "o", "."))
> 
> plot(newdat$leafHeight, yhat, type="l", main="manual: visited vs leafHeight",
+      xlim=c(15,80), ylim=c(-.05,1.1))
> points(newdat$leafHeight, lwr, pch="o")
> points(newdat$leafHeight, upr, pch="o")
> 
> # tita[,4] is age, tita[,1] is pclass
> printf("library(lars)\n")
library(lars)
> library(lars)
Loaded lars 1.3

> set.seed(2015)
> xmat <- as.matrix(tita[,c(2,5,6)])
> mod.lars.xmat <- lars(xmat, tita[,4])
> par(mfrow=c(2,2))
> plot(mod.lars.xmat)
> plotmo1(mod.lars.xmat, nresponse=4, do.par=F)
plotmo1(object=mod.lars.xmat,nresponse=4,do.par=F)
 plotmo grid:    survived sibsp  logage
                        0     0 3.06991
> plotres(mod.lars.xmat, trace=0, nresponse=4)
> 
> if(0) { # TODO fails with R-3.4.2: object '.QP_qpgen2' not found
+     printf("library(cosso)\n")
+     library(cosso)
+     set.seed(2016)
+     cosso <- cosso(xmat,tita[,4],family="Gaussian")
+     # TODO tell maintainer of cosso that you have to do this
+     class(cosso) <- "cosso"
+     set.seed(2016)
+     plotmo1(cosso)
+     set.seed(2016)
+     plotres(cosso)
+ }
> # examples from James, Witten, et al. ISLR book
> # I tested all models in their scripts manually.
> # All worked except for exceptions below.
> 
> printf("library(pls)\n")
library(pls)
> library(pls)

Attaching package: 'pls'

The following object is masked from 'package:stats':

    loadings

> printf("library(ISLR)\n")
library(ISLR)
> library(ISLR)
> Hitters=na.omit(Hitters)
> 
> set.seed(1)
> x <- model.matrix(Salary~.,Hitters)[,-1]
> y <- Hitters$Salary
> train=sample(1:nrow(x), nrow(x)/2)
> pcr.fit1=pcr(Salary~., data=Hitters,subset=train,scale=TRUE, validation="CV")
> plotmo1(pcr.fit1, nresponse=10)
plotmo1(object=pcr.fit1,nresponse=10)
 plotmo grid:    AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun
                   394  102     8   50  44    36     6   1931   510     36
 CRuns CRBI CWalks League Division PutOuts Assists Errors NewLeague
   246  219    172      A        W     211      56      7         A
> 
> # set.seed(1)
> # x <- model.matrix(Salary~.,Hitters)[,-1]
> # y <- Hitters$Salary
> # train=sample(1:nrow(x), nrow(x)/2)
> # pcr.fit2=pcr(y~x,scale=TRUE,ncomp=7)
> # # TODO following gives Error: predictions returned the wrong length (got 263 but expected 50)
> # plotmo1(pcr.fit2, nresponse=5)
> 
> library(splines)
> fit.lm2=lm(wage~bs(age,knots=c(25,40,60)),data=Wage)
> par(mfrow=c(1,2),mar=c(4.5,4.5,1,1),oma=c(0,0,4,0))
> agelims=range(Wage$age)
> age.grid=seq(from=agelims[1],to=agelims[2])
> pred=predict(fit.lm2,newdata=list(age=age.grid),se=T)
> plot(Wage$age,Wage$wage,col="gray", ylim=c(0,320))
> lines(age.grid,pred$fit,lwd=2)
> lines(age.grid,pred$fit+2*pred$se,lty="dashed")
> lines(age.grid,pred$fit-2*pred$se,lty="dashed")
> fit.lm2=lm(wage~bs(age,knots=c(25,40,60)),data=Wage,model=F) # TODO delete
> plotmo1(fit.lm2, col.resp=2, do.par=F, level=.95, ylim=c(0,320),
+         nrug=TRUE, caption="fit.lm2", ylab="wage")
plotmo1(object=fit.lm2,col.resp=2,do.par=F,level=0.95,ylim=c(0,320),nrug=TRUE,ylab="wage",caption="fit.lm2")
> 
> fit.glm2 <- glm(I(wage>250)~poly(age,4),data=Wage,family=binomial)
> par(mfrow=c(1,2),mar=c(4.5,4.5,1,1),oma=c(0,0,4,0))
> agelims=range(Wage$age)
> age.grid=seq(from=agelims[1],to=agelims[2])
> # their plot
> preds=predict(fit.glm2,newdata=list(age=age.grid),se=T)
> pfit=exp(preds$fit)/(1+exp(preds$fit))
> se.bands.logit = cbind(preds$fit+2*preds$se.fit, preds$fit-2*preds$se.fit)
> se.bands = exp(se.bands.logit)/(1+exp(se.bands.logit))
> preds=predict(fit.glm2,newdata=list(age=age.grid),type="response",se=T)
> plot(Wage$age,I(Wage$wage>250),xlim=agelims,type="n",ylim=c(0,.2))
> points(jitter(Wage$age), I((Wage$wage>250)/5),cex=.5,pch="|",col="darkgrey")
> lines(age.grid,pfit,lwd=2, col="blue")
> matlines(age.grid,se.bands,lwd=1,col="blue",lty=3)
> # plotmo plot, side by side
> # TODO Warning: the level argument may not be properly supported on glm objects built with weights
> plotmo1(fit.glm2, level=.95, degree1.col="blue", ylim=c(0,.2), do.par=FALSE, nrug=-1, caption="fit.glm2", ylab="I(wage > 250)")
plotmo1(object=fit.glm2,level=0.95,degree1.col="blue",ylim=c(0,0.2),do.par=FALSE,nrug=-1,ylab="I(wage>250)",caption="fit.glm2")
Warning: the level argument may not work correctly on glm objects built with weights
> 
> # Test deparsing of the formula in plotmo.pairs.default
> # TODO Height is included in the plots even though formula says -Height
> Height2 <- trees$Height^2
> a <- lm(Volume~(Girth*Height2)-Height, data=trees, x=TRUE, model=FALSE)
> plotmo(a)
 plotmo grid:    Girth Height2 Height
                  12.9    5776     76
> 
> # test "the variable on the right side of the formula is a matrix or data.frame"
> # TODO would like to solve this problem
> 
> options(warn=2)
> data(gasoline, package="pls")
> earth.octane <- earth(octane ~ NIR, data=gasoline)
> print(summary(earth.octane)) # ok
Call: earth(formula=octane~NIR, data=gasoline)

                         coefficients
(Intercept)                 87.818970
h(NIR1016 nm- -0.050322)  -307.631441
h(NIR1036 nm- -0.060936)    83.025904
h(NIR1054 nm- -0.059068)   254.542458
h(NIR1134 nm-0.028475)      34.069219
h(0.484052-NIR1194 nm)     -45.522897
h(NIR1194 nm-0.484052)      50.623858
h(0.25499-NIR1208 nm)       81.506833
h(NIR1208 nm-0.25499)      -92.719551
h(NIR1686 nm-1.25012)       -7.936903
h(1.27324-NIR1690 nm)        3.531658

Selected 11 of 12 terms, and 8 of 401 predictors
Termination condition: RSq changed by less than 0.001 at 12 terms
Importance: NIR1208 nm, NIR1194 nm, NIR1134 nm, NIR1690 nm, NIR1016 nm, ...
Number of terms at each degree of interaction: 1 10 (additive model)
GCV 0.05120795    RSS 1.298122    GRSq 0.9784914    RSq 0.990602
> plotres(earth.octane) # ok
> expect.err(try(plotmo(earth.octane)), "the variable on the right side of the formula is a matrix or data.frame")
Error : (converted from warning) the variable on the right side of the formula is a matrix or data.frame
         plotmo often cannot process such variables
Got expected error from try(plotmo(earth.octane))
> options(warn=1)
> 
> # TODO May 2020 'ElemStatLearn' is not available (for R version 4.0.0)
> # library(ElemStatLearn)
> # x <- mixture.example$x
> # g <- mixture.example$y
> # lm.mixture.example <- lm(g ~ x)
> # options(warn=2)
> # expect.err(try(plotmo(lm.mixture.example)), "the variable on the right side of the formula is a matrix or data.frame")
> # options(warn=1)
> 
> # test variable names with $ are not supported
> 
> a <- earth(O3~ozone1$doy, data=ozone1)
> expect.err(try(plotmo(a)), "cannot get the original model predictors")
Warning: "$" in the formula is not supported by plotmo, will try to get the data elsewhere
         formula: ozone1$doy

Looked unsuccessfully for the original predictors in the following places:

(1) object$x: NULL

(2) model.frame: formula(object): "$" in formula is not allowed

(3) getCall(object)$x: NULL

Error : cannot get the original model predictors
Got expected error from try(plotmo(a))
> 
> a <- earth(O3~ozone1$doy + temp, data=ozone1)
> expect.err(try(plotmo(a)), "cannot get the original model predictors")
Warning: "$" in the formula is not supported by plotmo, will try to get the data elsewhere
         formula: ozone1$doy + temp

Looked unsuccessfully for the original predictors in the following places:

(1) object$x: NULL

(2) model.frame: formula(object): "$" in formula is not allowed

(3) getCall(object)$x: NULL

Error : cannot get the original model predictors
Got expected error from try(plotmo(a))
> 
> a <- lm(O3~ozone1$doy, data=ozone1)
> expect.err(try(plotmo(a)), "cannot get the original model predictors")
Warning: "$" in the formula is not supported by plotmo, will try to get the data elsewhere
         formula: ozone1$doy

Looked unsuccessfully for the original predictors in the following places:

(1) object$x: NULL

(2) model.frame: formula(object): "$" in formula is not allowed

(3) getCall(object)$x: NULL

Error : cannot get the original model predictors
Got expected error from try(plotmo(a))
> 
> a <- lm(O3~ozone1$doy + temp, data=ozone1)
> expect.err(try(plotmo(a)), "cannot get the original model predictors")
Warning: "$" in the formula is not supported by plotmo, will try to get the data elsewhere
         formula: ozone1$doy + temp

Looked unsuccessfully for the original predictors in the following places:

(1) object$x: NULL

(2) model.frame: formula(object): "$" in formula is not allowed

(3) getCall(object)$x: NULL

Error : cannot get the original model predictors
Got expected error from try(plotmo(a))
> 
> #--- test interaction of w1. and non w1 args -------------------------------------
> 
> par(mfrow=c(4,3), mar=c(3, 3, 4, 1), mgp=c(2, 0.6, 0))
> 
> mod78 <- earth(Volume ~ ., data = trees)
> par(mfrow=c(3,4), mar=c(3, 3, 3, 1), mgp=c(2, 0.6, 0))
> 
> # multiple which, earth model
> plotres(mod78, cex.main=1,
+         ylim=c(-.5, .8),    xlim=c(-2, 7),    col=2:3, do.par=FALSE,
+         w1.main=c("ylim=c(-.5, .8)\nxlim=c(-2, 7) col=2:3"))
> 
> # multiple which, earth model
> plotres(mod78, cex.main=.7,
+         w1.ylim=c(-.5, .8), w1.xlim=c(-2, 7), col=2:3, do.par=FALSE,
+         ylim=c(-10,10), xlim=c(-30, 100),
+         w1.main=c("w1.ylim=c(-.5, .8) w1.xlim=c(-2, 7)\nylim=c(-10,10), xlim=c(-30, 100)"))
> par(org.par)
> 
> par(mfrow=c(3,4), mar=c(3, 3, 3, 1), mgp=c(2, 0.6, 0))
> 
> # which=1, earth model
> 
> plotres(mod78, which=1, cex.main=.8,
+         col=2:3,
+         main="which=1, no other ylim args",
+         w1.main="which=1, no other ylim args")
> 
> plotres(mod78, which=1, cex.main=.8,
+         col=2:3, w1.ylim=c(.3,.98), w1.xlim=c(-2, 7),
+         main="w1.ylim=c(.3,.98)\nw1.xlim=c(-2, 7)")
> 
> plotres(mod78, which=1, cex.main=.8,
+         col=2:3, ylim=c(.3,.98),    xlim=c(-2, 7),
+         main="ylim=c(.3,.98)\nxlim=c(-2, 7)")  # ylim gets passed to modsel
> 
> plotres(mod78, which=1, cex.main=.75,
+         col=2:3, w1.ylim=c(.3,.98), ylim=c(-.5,.5),
+         w1.xlim=c(-2, 7), xlim=c(-90, 90),
+         main="w1.ylim=c(.3,.98), ylim=c(-.5,.5)\nw1.xlim=c(-2, 7), xlim=c(-90, 90)") # ignore ylim
> 
> # which=3, earth model
> plotres(mod78, which=3, cex.main=1,
+         col=2:3,
+         main="which=3, no other ylim args")
> 
> plotres(mod78, which=3, cex.main=1,
+         col=2:3, w1.ylim=c(.3,.98), w1.xlim=c(-2, 7),
+         main="w1.ylim=c(.3,.98)\nw1.xlim=c(-2, 7)") # not usual, ignore w1.ylim
> 
> plotres(mod78, which=3, cex.main=1,
+         col=2:3, ylim=c(-10,10), xlim=c(-90,90),
+          main="which=3, ylim=c(-10,10)\nxlim=c(-90,90)")
> 
> plotres(mod78, which=3, cex.main=1,
+         col=2:3, w1.ylim=c(.3,.98), ylim=c(-10,10), w1.xlim=c(-2, 7), xlim=c(-90,90),
+         main="w1.ylim=c(.3,.98) ylim=c(-10,10)\nw1.xlim=c(-2, 7), xlim=c(-90,90)")
> 
> par(org.par)
> 
> nullarg <- NULL
> expect.err(try(plotmo(nullarg)),   "argument 'nullarg' is NULL")
Error : argument 'nullarg' is NULL
Got expected error from try(plotmo(nullarg))
> expect.err(try(plotmo(NULL)),      "argument 'NULL' is NULL")
Error : argument 'NULL' is NULL
Got expected error from try(plotmo(NULL))
> expect.err(try(plotmo(0)),         "'0' is not an S3 model")
Error : '0' is not an S3 model
Got expected error from try(plotmo(0))
> expect.err(try(plotmo(list(1,2))), "'list(1, 2)' is a plain list, not an S3 model")
Error : 'list(1, 2)' is a plain list, not an S3 model
Got expected error from try(plotmo(list(1, 2)))
> expect.err(try(plotmo(list(1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0,1,2,3,4,5,6,7,8,0))),
+                "object is a plain list, not an S3 model")
Error : object is a plain list, not an S3 model
Got expected error from try(plotmo(list(1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 6,     7, 8, 0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 6, 7,     8, 0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 6, 7, 8,     0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 6, 7, 8, 0,     1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1,     2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2,     3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3,     4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4,     5, 6, 7, 8, 0)))
> 
> source("test.epilog.R")