Nothing
# earth.varmod.R: build variance models for estimating prediction intervals
#
# TODO Extend the coverage table (print_inconf_tab) to show percentages in
# lower and upper intervals, so the user can check for asymmetry of the
# residuals.
#
# TODO Add QQ plot for prediction intervals, a "PIQ" plot
#
# TODO Consider making the code automatically detect non-monotonicity and
# issuing a warning. Probably only possible for univariate models.
#
# TODO Could maybe prevent "Error in numericDeriv" by internally passing
# an explicit derivative in the call to nls.
VARMOD.METHODS <- c("const", "power", "power0",
"lm", "rlm", "earth", "gam",
"x.lm", "x.rlm", "x.earth", "x.gam")
TRACE.VARMOD <- .3
TRACE.VARMOD.DETAILS <- .31 # will also cause plotting
# varmod returns a "varmod" object.
# y is the observed response (it is a n x 1 matrix)
varmod <- function(parent,
method, exponent, conv, clamp, minspan,
trace, x, y, model.var, ...)
{
UseMethod("varmod")
}
varmod.earth <- function(parent,
method, exponent, conv, clamp, minspan,
trace, x, y, model.var, ...)
{
check.classname(parent, substitute(parent), "earth")
varmod_internal(parent,
method, exponent, conv, clamp, minspan,
trace, x, y, model.var, ...)
}
varmod.default <- function(parent,
method, exponent, conv, clamp, minspan,
trace, x, y, model.var, ...)
{
warning0("varmod.default: varmods are not supported for \"",
class(parent)[1], "\" objects\nContinuing anyway")
# TODO this won't work: varmod_internal assumes an earth parent model in some places
varmod_internal(parent,
method, exponent, conv, clamp, minspan,
trace, x, y, model.var, ...)
}
varmod_internal <- function(parent,
method, exponent=1, conv=1, clamp=.1,
minspan=-5,
trace=0, parent.x=NULL, parent.y=NULL,
model.var, ...)
{
# The following constant "lambda" was an argument to earth but I removed
# it and hardcoded it here for simplicity in the earth interface.
# We use lambda to transform the squared residuals as follows:
# transformed.resids = squared.resids ^ (lambda / 2)
# So with lambda=1, we transform to absolute residuals, and if
# lambda=2, then there is no transform. We call the transformed
# residuals the abs.resids in the code (which is actually the correct
# nomenclature only when lambda is 1). See also get.resids.name.
lambda <- 1
# likewise, rmethod is hardcoded here instead of being an arg to earth
rmethod <- "hc12" # TODO doesn't match documentation "Variance models in earth"
trace <- as.numeric(check.numeric.scalar(trace, logical.ok=TRUE))
check.lambda.arg(lambda)
check.exponent.arg(exponent, method)
check.conv.arg(conv)
check.clamp.arg(clamp)
stopifnot(!is.null(parent.x))
stopifnot(is.matrix(parent.x))
stopifnot(!is.null(parent.y))
stopifnot(is.matrix(parent.y))
if(NCOL(parent.y) != 1)
stop0("variance models are not supported for multiple response models")
if(trace >= TRACE.VARMOD) {
printf(
"\nvarmod method=\"%s\" rmethod=\"%s\" lambda=%g exponent=%g conv=%g clamp=%g minspan=%g:\n",
method, rmethod, lambda, exponent, conv, clamp, minspan)
if(trace == TRACE.VARMOD.DETAILS) {
oldpar <- par(no.readonly=TRUE)
on.exit(par(oldpar))
par(mfrow=c(2, 3), mar=c(3, 3, 3, 1), mgp=c(1.5, .5, 0))
}
}
n <- nrow(parent.x)
df <- length(parent$selected.terms)
leverages <- parent$leverages
stopifnot(!is.null(leverages))
leverages[leverages > .9] <- .9 # prevent any residual from being too influential
correction <- switch(match.choices(rmethod,
c("hc0", "hc1", "hc2", "hc3", "hc12"), "varmod.rmethod"),
hc0 = 1,
hc1 = n / (n - df),
hc2 = 1 / (1 - leverages),
hc3 = 1 / (1 - leverages)^2,
hc12 = n / ((n - df) * (1 - leverages)))
squared.resids <- (parent.y - parent.predict(parent))^2
squared.resids <- correction * squared.resids + model.var
abs.resids <- squared.resids ^ (lambda / 2) # by default lambda=1, so this takes sqrt
temp <- iterate.residmod(parent, abs.resids,
method, exponent, lambda, conv, clamp, minspan,
trace, parent.x, parent.y, ...)
residmod <- temp$residmod
converged <- temp$converged
iters <- temp$iters
if(trace >= TRACE.VARMOD)
printf("\n")
varmod <- NULL
varmod$call <- make.call.generic(match.call(), "varmod")
varmod$parent <- parent
varmod$method <- method
varmod$package <- which.package(method)
varmod$exponent <- exponent
varmod$lambda <- lambda
varmod$rmethod <- rmethod
varmod$converged <- temp$converged
varmod$iters <- temp$iters
varmod$residmod <- residmod
varmod$min.sd <- get.min.sd(residmod, lambda, clamp)
varmod$model.var <- model.var
varmod$abs.resids <- abs.resids # transformed residuals (actually only abs when lambda is 1)
varmod$parent.x <- parent.x
varmod$parent.y <- parent.y
class(varmod) <- "varmod"
varmod$iter.rsq <- get.iter.rsq(varmod, abs.resids)
varmod$iter.stderr <- get.iter.stderr(varmod, trace)
attr(varmod, ".Environment") <- get.model.env(residmod, "varmod", trace)
varmod
}
get.iter.rsq <- function(object, abs.resids) # return NULL if can't get rsq
{
check.classname(object, substitute(object), "varmod")
if(object$method == "const")
return(NULL)
fitted <- fitted(object$residmod)
weights <- weights(object$residmod)
if(is.null(fitted) || is.null(weights))
return(NULL)
get.weighted.rsq(abs.resids, fitted, weights)
}
get.iter.stderr <- function(object, trace) # return if NULL if can't get stderr
{
check.classname(object, substitute(object), "varmod")
residmod <- object$residmod
if(class(residmod)[1] %in% c("lm", "rlm", "nls")) {
coef <- summary(residmod)$coefficients
stopifnot(!is.null(coef[,"Std. Error"]))
coef[,"Std. Error"]
} else if(class(residmod)[1] == "Gam") { # package gam version 1.15 or higher
coef <- coefficients(summary.glm(residmod))
stopifnot(!is.null(coef[,"Std. Error"]))
coef[,"Std. Error"]
} else if(class(residmod)[1] == "gam" && object$package == "gam") { # package gam version less than 1.15
coef <- coefficients(summary.glm(residmod))
stopifnot(!is.null(coef[,"Std. Error"]))
coef[,"Std. Error"]
} else if(class(residmod)[1] == "gam" && object$package == "mgcv") {
# only the stderr for the intercept is available
len.coef <- length(coefficients(residmod))
std.err <- repl(NA, len.coef)
std.err[1] <- summary(residmod)$p.table[1, "Std. Error"] # se of intercept
std.err
} else if(class(residmod)[1] == "earth") {
trace2(trace, "--get.iter.stderr\n")
y <- plotmo::plotmo_y(residmod, nresponse=1, trace)$y
bx <- model.matrix(residmod)
coef <- summary(lm(y~bx))$coefficients
coef[,"Std. Error"]
} else
NULL
}
# iteratively reweighted least squares
issued.singularities.warning.global <- FALSE
iterate.residmod <- function(parent, abs.resids,
method, exponent, lambda, conv, clamp, minspan,
trace, parent.x, parent.y, ...)
{
varmod <- NULL
max.iter <- 50
weights <- rep(1, nrow(parent.y))
residmod <- NULL
# following is needed because we want to issue singular warning at most once
assignInMyNamespace("issued.singularities.warning.global", FALSE)
# we always build the trace tab but only print it
# if tracing is enabled or convergence failed
trace.tab <- NULL
for(iter in 1:max.iter) {
residmod <- get.residmod(method, exponent, minspan, parent.x, parent.y,
abs.resids, weights, trace, iter, parent, residmod, ...)
coef.change <- get.coef.change(method, iter, residmod)
# fill in enough of varmod for predict.varmod in get.residmod.weights
varmod$parent <- parent
varmod$method <- method
varmod$exponent <- exponent
varmod$lambda <- lambda
varmod$residmod <- residmod
varmod$min.sd <- get.min.sd(residmod, lambda, clamp)
class(varmod) <- "varmod"
weights <- get.residmod.weights(varmod, iter, trace)
trace.tab <- update_trace_tab(trace.tab, trace, iter, max.iter,
residmod, varmod, coef, coef.change,
parent.y, weights, exponent)
if(residmod.converged(coef.change, conv, iter, max.iter, method))
break
}
converged <- residmod.converged(coef.change, conv, iter, max.iter, method)
if(trace >= TRACE.VARMOD) {
# || anyNA(coef(residmod)) || (conv >= 0 && !converged))
if(trace == TRACE.VARMOD.DETAILS && inherits(residmod, "nls"))
printf("\n")
print(trace.tab[1:iter,], row.names=FALSE, digits=2)
}
if(!converged)
warnf("varmod did not converge after %d iters (%s), final coefchange %.1f%%",
iter, get.non.convergence.reason(trace.tab), mean(abs(coef.change)))
list(residmod=residmod, converged=converged, iters=iter)
}
blank.plot <- function(main=NULL)
{
plot(0, 0, col=0, bty="n", xlab="", ylab="",
xaxt="n", yaxt="n", main=main)
}
trace.ncoef.global <- 0
update_trace_tab <- function(trace.tab, trace, iter, max.iter, residmod, varmod,
coef, coef.change, parent.y, weights, exponent)
{
if(trace == TRACE.VARMOD.DETAILS) {
if(inherits(residmod, "nls"))
printf("\n")
plotmo::plotmo(varmod, type="abs.residual",
do.par=FALSE, degree1=1, degree2=0, trace=-1,
pt.col=1, degree1.col=2, degree1.lwd=3,
main="residmod first predictor")
}
coef <- coef(residmod)
if(iter == 1) {
# following needed because nbr of earth coefs can change
assignInMyNamespace("trace.ncoef.global", length(coef))
trace.tab <- as.data.frame(matrix(NA, nrow=max.iter, ncol=3+trace.ncoef.global))
colnames(trace.tab) <- c(" iter", "weight.ratio", "coefchange%",
fix.coef.names(names(coef), colnames(parent.y), exponent))
}
trace.tab[iter,] <-
c(iter, max(weights) / min(weights),
mean(abs(coef.change)),
c(coef, repl(NA, trace.ncoef.global))[1:trace.ncoef.global])
trace.tab
}
# For debugging non-convergence. In simulation about .4% of runs
# did not converge, mostly "oscillating", nearly all with sample
# sizes of less than 100.
get.non.convergence.reason <- function(trace.tab)
{
nrow <- nrow(trace.tab)
if(nrow < 7)
return <- "short tab" # should never return this
coefchange <- trace.tab[,"coefchange%"]
c6 <- coefchange[nrow-6]
c5 <- coefchange[nrow-5]
c4 <- coefchange[nrow-4]
c3 <- coefchange[nrow-3]
c2 <- coefchange[nrow-2]
c1 <- coefchange[nrow-1]
c0 <- coefchange[nrow-0]
# Example for oscillating:
# iter weight.ratio coefchange% (Intercept) y
# 27 57 big c3 3.7 0.00972 0.089
# 28 63 small c2 3.4 0.00921 0.091
# 29 57 big c1 3.5 0.00971 0.089
# 30 62 small c0 3.3 0.00922 0.091
if(c0 < c1 && c1 > c2 && c2 < c3 && c3 > c0)
return("oscillating-lo") # oscillating, last iter lower than prev
else if(c0 > c1 && c1 < c2 && c2 > c3 && c3 < c0)
return("oscillating-hi") # oscillating, last iter higher than prev
if(c2 > c1 && c1 > c0) { # only last two converged
reason <- "converging2"
if(c3 > c2) { # only last three converged
reason <- "converging3"
if(c4 > c3) { # only last four converged
reason <- "converging4"
if(c5 > c4) { # only last five converged
reason <- "converging5"
if(c6 > c5) # last six converged
reason <- "converging6"
}
}
}
return(reason)
}
if(c2 < c1 && c1 < c0) { # last two diverged
reason <- "diverging2"
if(c3 < c2) { # last three diverged
reason <- "diverging3"
if(c4 < c3) { # last four diverged
reason <- "diverging4"
if(c5 < c4) { # last five diverged
reason <- "diverging5"
if(c6 < c5) # last six diverged
reason <- "diverging6"
}
}
}
return(reason)
}
"non-monotonic" # can find no other pattern
}
get.residmod <- function(method, exponent, minspan, parent.x, parent.y,
abs.resids, weights, trace, iter, parent, prev.residmod, ...)
{
switch(method,
const = residmod.const(parent.x, parent.y, abs.resids,
weights, trace),
power = residmod.power(exponent, parent.x, parent.y, abs.resids,
weights, trace, parent, prev.residmod, iter, ...),
power0 = residmod.power0(exponent, parent.x, parent.y, abs.resids,
weights, trace, parent, prev.residmod, iter, ...),
lm = residmod.lm(exponent, parent.x, parent.y, abs.resids,
weights, trace, parent, ...),
rlm = residmod.rlm(exponent, parent.x, parent.y, abs.resids,
weights, trace, parent, ...),
earth = residmod.earth(exponent, minspan, parent.x, parent.y, abs.resids,
weights, trace, parent, ...),
gam = residmod.gam(exponent, parent.x, parent.y, abs.resids,
weights, trace, parent, iter, ...),
x.lm = residmod.x.lm(exponent, parent.x, parent.y, abs.resids,
weights, trace, ...),
x.rlm = residmod.x.rlm(exponent, parent.x, parent.y, abs.resids,
weights, trace, ...),
x.earth = residmod.x.earth(exponent, minspan, parent.x, parent.y, abs.resids,
weights, trace, ...),
x.gam = residmod.x.gam(exponent, parent.x, parent.y, abs.resids,
weights, trace, iter, ...),
stop0("illegal varmod.method \"", method, "\""))
}
prev.coef.global <- NULL
get.coef.change <- function(method, iter, residmod) # returns percents for each coef
{
coef <- coef(residmod)
# lm sometimes returns 2nd coef as NA if resids are all the same
# TODO this should be an error? --- but that halts simulation tests
if(anyNA(coef)) {
if(!issued.singularities.warning.global) {
warning0("singularities in residual model (coefs ",
paste.collapse(coef), ")")
assignInMyNamespace("issued.singularities.warning.global", TRUE)
}
coef[is.na(coef)] <- 0
}
if(iter == 1)
assignInMyNamespace("prev.coef.global", coef)
if(method %in% c("earth", "x.earth")) # see comments in residmod.converged
if(length(prev.coef.global) != length(coef))
return(9999)
coef.change <- abs(coef - prev.coef.global)
prev <- abs(prev.coef.global)
assignInMyNamespace("prev.coef.global", coef)
# Divide by absolute value of previous coefficients.
# But ensure no divide by near zero, by downweighting extremely small
# coefs, thus preventing them from completely dominating the mean
# change if the rest are large.
# The 1e-8 prevents noise floor coefficients from preventing convergence.
min <- max(.01 * max(prev), 1e-8)
prev[prev < min] <- min
100 * coef.change / prev # a percentage for each coef
}
residmod.converged <- function(coef.change, conv, iter, max.iter, method)
{
if(conv < 0)
iter >= abs(conv)
else {
method == "const" ||
# TODO Since the earth basis funcs can change, looking at changes
# in the coefs can't be used to determine convergence.
# So for now, always do 1 iter for earth residual models.
(method %in% c("earth", "x.earth")) ||
# TODO following will sometimes create an intercept only model so unused
# (method %in% c("earth", "x.earth") && iter >= 2) ||
iter > 1 && mean(abs(coef.change)) < conv
}
}
draw.residmod.weights <- function(w, main, min=NA, max=NA, median=NA) # for debugging
{
plot(w, type="l", main=main, ylim=c(0, max(w, if(is.na(max)) 0 else max)))
if(!is.na(min)) {
abline(h=min, col=2)
abline(h=max, col=2)
abline(h=median, col=2)
}
else
legend("topright", sprint("max/min %.0f", max(w) / min(w)))
lines(w) # replot over other annotations
}
# clamp to prevent extreme weights after squaring and inverse in get.residmod.weights
clamp.se <- function(se, iter, trace)
{
median <- median(se)
max.ratio <- 5 # 5 seems ok with (limited) simulation studies
min <- median / max.ratio
max <- max.ratio * median
if(trace == TRACE.VARMOD.DETAILS)
draw.residmod.weights(se, main=sprint("iter %d: se", iter), min, max, median)
se[se < min] <- min
se[se > max] <- max
se
}
# The variance for a regression on absolute residuals is proportional to
# the square of the regression model predicted value (Carrol and Ruppert
# book Section 3.3.3 and Table 3.3).
get.residmod.weights <- function(object, iter=0, trace=0)
{
check.classname(object, substitute(object), "varmod")
# square to convert se to variance, inverse to convert variance to weight
weights <- 1 / clamp.se(predict.varmod(object, type="se"), iter, trace)^2
# normalization is not strictly necessary, may help numerics
weights <- weights / mean(weights)
if(trace == TRACE.VARMOD.DETAILS)
draw.residmod.weights(weights, "weights")
weights
}
# We calculate lamba.factor.global only when necessary because the
# calculation can be slow. Hence we need the following global variables.
lamba.global <- lamba.factor.global <- -999
update_lambda_factor <- function(lambda, trace)
{
approx.equal <- function(x, y)
{
# allow for limited precision in doubles, also allows .33 for 1/3
abs(x - y) < 1e-2
}
#--- update_lambda_factor starts here ---
if(lambda != lamba.global) {
assignInMyNamespace("lamba.global", lambda)
# some values have been precalculated
if(approx.equal(lambda, 2))
assignInMyNamespace("lamba.factor.global", 1)
# sqrt(pi / 2) = 1.2533, ratio mean dev to stddev, Geary 1935
else if(approx.equal(lambda, 1))
assignInMyNamespace("lamba.factor.global", sqrt(pi / 2))
# (residuals^2)^(1/3) is approx normal by the Wilson-Hilferty
# transform, although the left tail will still be short
else if(approx.equal(lambda, 2/3))
assignInMyNamespace("lamba.factor.global", 1.2464)
else
{
rnorm(1) # seems to be necessary to make .Random.seed available
old.seed <- .Random.seed
set.seed(1) # for reproducibility
# 1e6 below could be bigger but then slow
assignInMyNamespace("lamba.factor.global",
1 / mean(rnorm(1e6)^2 ^ (lambda/2)))
set.seed(old.seed)
}
if(trace >= TRACE.VARMOD)
printf("lambda %g lamba.factor %g\n", lambda, lamba.factor.global)
}
}
# scale a prediction by the residmod back to a standard deviation
to.sd <- function(abs.resids, lambda, trace=0)
{
update_lambda_factor(lambda, trace)
# pmax is necessary to prevent e.g. sqrt of neg prediction from residmod
(lamba.factor.global * pmax(abs.resids, 0)) ^ (1 / lambda)
}
get.min.sd <- function(residmod, lambda, clamp=.1)
{
predict <- predict(residmod)
predict <- predict[predict > 0]
stopifnot(length(predict) > 0)
stopifnot(clamp >= 0, clamp <= 1)
clamp * mean(to.sd(predict, lambda, 0))
}
check.lambda.arg <- function(lambda)
{
check.numeric.scalar(lambda)
if(lambda < .25 || lambda > 2)
stop0("lambda=", lambda, " but it should be between 0.25 and 2")
}
# TRUE if estimation of variance depends only on the fitted response (not on x)
method.uses.fitted.response <- function(method)
{
method %in% c("power", "power0", "lm", "rlm", "earth", "gam")
}
check.exponent.arg <- function(exponent, method)
{
check.numeric.scalar(exponent)
# TODO following restriction could be lifted but currently only partially implemented
if(exponent != 1 && !method.uses.fitted.response(method))
stop0("varmod.exponent argument is not allowed with method=\"", method, "\"\n",
"(varmod.exponent is only allowed for varmod.methods that depend only ",
"on the fitted response)")
if(exponent < .1 || exponent > 5)
stop0("varmod.exponent=", exponent, " but it should be between .1 and 5")
}
check.conv.arg <- function(conv)
{
err <- function(conv)
stop0("varmod.conv=", conv,
" but it should be a negative integer ",
"or a percent between 0 and 100")
check.numeric.scalar(conv)
if(conv < 0) {
if(floor(conv) != conv) # conv is negative, check that it is an integer
err(conv)
} else if(conv == 0 || conv > 100)
err(conv)
}
check.clamp.arg <- function(clamp)
{
check.numeric.scalar(clamp)
if(clamp < 0 || clamp > 1)
stop0("varmod.clamp=", clamp, " but it should be between 0 and 1")
}
residmod.const <- function(parent.x, parent.y, abs.resids, weights, trace)
{
# Predictions can be handled in a simple consistent way in
# residmod.predict if instead of calculating the variance directly
# here, we achieve the same result by building an intercept-only model
# which always predicts mean(abs.resids).
#
# The conversion to a dataframe is necessary if the user later calls
# plot(parent$varmod$residmod) or plotmo(parent$varmod$residmod).
# Note that plotmo will call predict.varmod via predict.earth.
data <- data.frame(abs.resids, parent.x)
colnames(data) <- c("abs.resids", colnames(parent.x))
lm(abs.resids~1, data=data, weights=weights, y=TRUE)
}
apply.exponent <- function(yhat, exponent)
{
check.vec(yhat, "yhat")
# exponents of neg numbers are allowed only for integer exponents
if(floor(exponent) != exponent) {
check.that.most.are.positive(
yhat, "parent.fit", sprint("exponent=%g", exponent), "nonpositive")
yhat[yhat < 0] <- 0 # don't want to take say sqrt of a neg number
}
yhat ^ exponent
}
nls.wrapper <- function(form, data, start, weights, abs.resids, trace)
{
# We use algorithm="port" below because the default algorithm more often causes
# "Error in numericDeriv: Missing value or an infinity produced"
# Also, on test data we sometimes need more iterations than the default 50
mod <- nls(formula=form,
data=data, start=start, weights=weights,
trace=(trace == TRACE.VARMOD.DETAILS),
algorithm="port", control=list(maxiter=100))
# make model data available for plotmo and plotres
mod$x <- data[,-1,drop=FALSE]
mod$y <- abs.resids
# nls doesn't save the terms, so call$formula can confuse plotmo and plotres
mod$call <- NULL
mod
}
estimate.power.start.values <- function(prev.residmod, abs.resids, data, weights, trace, iter)
{
if(is.null(prev.residmod)) { # first iteration in iterate.residmod?
# use a linear model to estimate the start values
lm <- lm(abs.resids~., data=data, weights=weights)
coefs <- coef(lm)
if(trace == TRACE.VARMOD.DETAILS) {
plotmo::plotmo(lm, pt.col=2, do.par=FALSE, trace=-1,
main=sprint("iter 1: lm for start vals\nvarmod.method=\"power\""))
plot(lm, which=1)
blank.plot()
}
start <- list(coefs[1], coefs[2], exponent=1)
if(trace >= TRACE.VARMOD)
printf(
"\n start: (Intercept)=%.2g coef=%.2g exponent=%.2g\n\n",
start[[1]], start[[2]], start[[3]])
} else { # not first iteration
# use previous model values as starting values
coefs <- coef(prev.residmod)
stopifnot(length(coefs) == 3)
start <- list(coefs[1], coefs[2], coefs[3])
}
names(start) <- c("(Intercept)", "coef", "exponent")
if(trace == TRACE.VARMOD.DETAILS)
cat(sprint("iter %d RSS: ", iter), names(start), "\n")
start
}
residmod.power <- function(exponent, parent.x, parent.y, abs.resids,
weights, trace, parent, prev.residmod, iter, ...)
{
if(exponent != 1) # TODO allow this?
stop0("the exponent argument is not allowed with varmod.method=\"power\"")
parent.fit <- parent.predict(parent)
check.that.most.are.positive(
parent.fit, "parent.predict(parent)", "varmod.method=\"power\"", "nonpositive")
parent.fit[parent.fit < 0] <- 0 # force negative values to zero
form <- abs.resids~`(Intercept)` + coef * RHS^exponent
data <- data.frame(abs.resids, apply.exponent(parent.fit, exponent))
colnames(data) <- c("abs.resids", "RHS")
start <- estimate.power.start.values(prev.residmod, abs.resids,
data, weights, trace, iter)
nls.wrapper(form, data, start, weights, abs.resids, trace)
}
estimate.power0.start.values <- function(prev.residmod, abs.resids,
data, weights, trace, iter)
{
if(is.null(prev.residmod)) { # first iteration in iterate.residmod?
# use a linear model to estimate the start values
lm <- lm(abs.resids~.-1, data=data, weights=weights)
coefs <- coef(lm)
if(trace == TRACE.VARMOD.DETAILS) {
plotmo::plotmo(lm, pt.col=2, do.par=FALSE, trace=-1,
main=sprint("iter 1: lm for start vals\nvarmod.method=\"power0\""))
plot(lm, which=1)
blank.plot()
}
start <- list(coefs[1], exponent=1)
if(trace >= TRACE.VARMOD)
printf(
"\n start: coef=%.2g exponent=%.2g\n\n",
start[[1]], start[[2]])
} else { # not first iteration
# use previous model values as starting values
coefs <- coef(prev.residmod)
stopifnot(length(coefs) == 2)
start <- list(coefs[1], coefs[2])
}
names(start) <- c("coef", "exponent")
if(trace == TRACE.VARMOD.DETAILS)
cat(sprint("iter %d RSS: ", iter), names(start), "\n")
start
}
residmod.power0 <- function(exponent, parent.x, parent.y, abs.resids,
weights, trace, parent, prev.residmod, iter, ...)
{
if(exponent != 1) # TODO allow this?
stop0("the exponent argument is not allowed with varmod.method=\"power0\"")
parent.fit <- parent.predict(parent)
check.that.most.are.positive(
parent.fit, "parent.predict(parent)", "varmod.method=\"power0\"", "nonpositive")
parent.fit[parent.fit < 0] <- 0 # force negative values to zero
data <- data.frame(abs.resids, apply.exponent(parent.fit, exponent))
colnames(data) <- c("abs.resids", "RHS")
start <- estimate.power0.start.values(prev.residmod, abs.resids,
data, weights, trace, iter)
nls.wrapper(abs.resids~coef * RHS^exponent,
data, start, weights, abs.resids, trace)
}
residmod.lm <- function(exponent, parent.x, parent.y, abs.resids,
weights, trace, parent, ...)
{
parent.fit <- parent.predict(parent)
data <- data.frame(abs.resids, apply.exponent(parent.fit, exponent))
# we use RHS instead of colnames(parent.y) because we have applied exponent
colnames(data) <- c("abs.resids", "RHS")
lm(abs.resids~., data=data, weights=weights, y=TRUE)
}
residmod.rlm <- function(exponent, parent.x, parent.y, abs.resids,
weights, trace, parent, ...)
{
parent.fit <- parent.predict(parent)
data <- data.frame(abs.resids, apply.exponent(parent.fit, exponent))
# we use RHS instead of colnames(parent.y) because we have applied exponent
colnames(data) <- c("abs.resids", "RHS")
mod <- MASS::rlm(abs.resids~., data=data, weights=weights, method="MM")
# make model data available for plotmo and plotres
mod$data <- data
mod
}
residmod.earth <- function(exponent, minspan, parent.x, parent.y, abs.resids,
weights, trace, parent, ...)
{
parent.fit <- parent.predict(parent)
data <- data.frame(abs.resids, apply.exponent(parent.fit, exponent))
colnames(data) <- c("abs.resids", "RHS")
earth(abs.resids~., data=data, weights=weights,
keepxy=TRUE, trace=trace, minspan=minspan, ...)
}
please.load.gam.package <- function()
{
stop0("please load either the gam or mgcv package before using varmod.method=\"gam\"")
}
# Do we use the gam function in the gam or the mgcv package?
# Note that library(gam) has to be used before calling gam::gam, else the
# wrong "s" function is invoked. This is because requireNamespace(gam)
# doesn't work there, even if we use gam::s when invoking gam.
# But CRAN check disallows library(gam) in the code (as from Jan 2015).
# So we have to ask the user to manually load the package if it is not loaded.
which.gam.package.is.loaded <- function()
{
gam.package.loaded <- "package:gam" %in% search()
mgcv.package.loaded <- "package:mgcv" %in% search()
if(mgcv.package.loaded && gam.package.loaded) {
# prevent downstream confusing error messages
stop0("varmod.method=\"gam\" is not allowed when both the ",
"'gam' and 'mgcv' packages are loaded")
}
if(gam.package.loaded)
return("gam")
if(mgcv.package.loaded)
return("mgcv")
please.load.gam.package()
}
which.package <- function(method)
{
if(method %in% c("gam", "x.gam")) {
if("package:gam" %in% search())
return("gam")
if("package:mgcv" %in% search())
return("mgcv")
}
NULL
}
residmod.gam.aux <- function(form, data, weights, trace, iter)
{
package.name <- which.gam.package.is.loaded()
if(package.name == "gam") {
if(trace >= TRACE.VARMOD && iter==1)
printf("using the gam function from the 'gam' package\n")
residmod <- gam::gam(formula=form, data=data, weights=weights)
# We don't use x=TRUE else the x has colnames like s(x) which
# confuses things later. But we do save the data for plotmo.
residmod$data <- data
} else if(package.name == "mgcv") {
if(trace >= TRACE.VARMOD && iter==1)
printf("using the gam function from the 'mgcv' package\n")
residmod <- mgcv::gam(formula=form, data=data, weights=weights)
residmod$data <- data # for later access by plotmo etc.
} else
please.load.gam.package()
residmod
}
residmod.gam <- function(exponent, parent.x, parent.y, abs.resids,
weights, trace, parent, iter, ...)
{
form <- abs.resids ~ s(RHS)
parent.fit <- parent.predict(parent)
RHS <- apply.exponent(parent.fit, exponent)
data <- data.frame(abs.resids, RHS)
colnames(data) <- c("abs.resids", "RHS")
residmod.gam.aux(form, data, weights, trace, iter)
}
residmod.x.lm <- function(exponent, parent.x, parent.y, abs.resids,
weights, trace, ...)
{
if(exponent != 1)
stop0("the exponent argument is not allowed with varmod.method=\"x.lm\"")
data <- data.frame(abs.resids, parent.x)
colnames(data) <- c("abs.resids", colnames(parent.x))
lm(abs.resids~., data=data, weights=weights, y=TRUE)
}
residmod.x.rlm <- function(exponent, parent.x, parent.y, abs.resids,
weights, trace, ...)
{
if(exponent != 1)
stop0("the exponent argument is not allowed with varmod.method=\"x.rlm\"")
data <- data.frame(abs.resids, parent.x)
colnames(data) <- c("abs.resids", colnames(parent.x))
mod <- MASS::rlm(abs.resids~., data=data, weights=weights,
method="MM", y.ret=TRUE)
# make model data available for plotmo and plotres
mod$y <- abs.resids
mod
}
residmod.x.earth <- function(exponent, minspan, parent.x, parent.y, abs.resids,
weights, trace, ...)
{
if(exponent != 1)
stop0("the exponent argument is not allowed with varmod.method=\"x.earth\"")
data <- data.frame(abs.resids, parent.x)
colnames(data) <- c("abs.resids", colnames(parent.x))
earth(abs.resids~., data=data, weights=weights,
keepxy=TRUE, trace=trace, minspan=minspan, ...)
}
residmod.x.gam <- function(exponent, parent.x, parent.y, abs.resids,
weights, trace, iter, ...)
{
if(exponent != 1)
stop0("the exponent argument is not allowed with varmod.method=\"x.gam\"")
if(ncol(parent.x) != 1)
stop0("varmod.method=\"x.gam\" is not allowed when x has more than one column")
form <- abs.resids ~ s(RHS)
RHS <- parent.x[,1]
data <- data.frame(abs.resids=abs.resids, RHS=RHS)
colnames(data) <- c("abs.resids", "RHS")
residmod.gam.aux(form, data, weights, trace, iter)
}
get.quant <- function(level) # e.g for level=.95 return 1.96
{
check.level.arg(level, zero.ok=FALSE)
stopifnot(level > 0, level < 1)
level <- 1 - (1 - level) / 2 # .95 becomes .975
qnorm(level) # .975 becomes 1.96
}
predict_se <- function(object, newdata)
{
to.sd(predict_abs_residual(object, newdata), object$lambda)
}
predict_abs_residual <- function(object, newdata)
{
# unfortunately needed to get model formulas to work for some varmod methods
hack.colnames <- function(newdata, method)
{
if(!is.null(newdata) &&
method %in% c("power", "power0", "lm", "rlm", "earth", "gam", "x.gam")) {
if(NCOL(newdata) != 1) {
stop0("predict.varmod: NCOL(newdata) must be 1 ",
"when method=\"", method, "\" (implementation restriction)")
}
newdata <- as.data.frame(newdata)
colnames(newdata) <- "RHS"
}
newdata
}
if(is.null(newdata))
abs.resid <- predict(object$residmod)
else if(method.uses.fitted.response(object$method)) {
parent.fit <- parent.predict(object$parent, newdata=newdata)
parent.fit <- apply.exponent(parent.fit, object$exponent)
parent.fit <- data.frame(parent.fit)
parent.fit <- hack.colnames(parent.fit, object$method)
if(object$method %in% c("power", "power0"))
parent.fit[parent.fit < 0] <- 0 # force negative values to zero
abs.resid <- predict(object$residmod, newdata=parent.fit)
stopifnot(length(abs.resid) == NROW(parent.fit))
} else {
newdata <- hack.colnames(newdata, object$method)
abs.resid <- predict(object$residmod, newdata=newdata)
stopifnot(length(abs.resid) == NROW(newdata))
}
abs.resid <- as.vector(abs.resid)
# clamp at object$min.sd
min.abs.resid <- (object$min.sd ^ object$lambda) / lamba.factor.global
pmax(abs.resid, min.abs.resid)
}
predict_pint <- function(object, newdata, level) # newdata allowed
{
se <- predict_se(object, newdata)
parent.fit <- parent.predict(object$parent, newdata)
stopifnot(length(parent.fit) == length(se))
quant <- get.quant(level)
data.frame(fit = parent.fit,
lwr = parent.fit - quant * se,
upr = parent.fit + quant * se)
}
predict_cint <- function(object, newdata, level)
{
if(!is.null(newdata))
stop0("predict.varmod: newdata is not allowed with interval=\"cint\"")
parent.fit <- parent.predict(object$parent, newdata)
se <- sqrt(object$model.var)
stopifnot(length(se) == length(parent.fit))
quant <- get.quant(level)
data.frame(fit = parent.fit,
lwr = parent.fit - quant * se,
upr = parent.fit + quant * se)
}
predict.varmod <- function(
object = stop("no 'object' argument"),
newdata = NULL,
type = c("pint", "cint", "se", "abs.residual"),
level = .95,
trace = FALSE, # unused but needed for plotmo
...)
{
check.level95 <- function(level)
{
check.level.arg(level, zero.ok=TRUE)
if(level != .95)
stop0("predict.varmod: the level argument is not allowed with type=\"",
type, "\"")
}
check.classname(object, substitute(object), "varmod")
warn.if.dots(...)
switch(match.arg1(type, "type"),
pint = predict_pint(object, newdata, level),
cint = predict_cint(object, newdata, level),
se = {
check.level95(level)
predict_se(object, newdata)
},
abs.residual = {
check.level95(level)
predict_abs_residual(object, newdata)
})
}
# Example: if digits=3, then "%.*f" becomes "%.3f"
# Needed because R printf doesn't support * in printf formats
# and we need it to make the digits arg work in printfs
dot.star.to.digits <- function(s, digits)
{
check.integer.scalar(digits, min=1)
stopifnot(floor(digits) == digits)
stopifnot(digits > 0, digits < 20)
gsub(".*", sprint(".%d", digits), s, fixed=TRUE)
}
# Example:
# fix.coef.names(coef.names=h(y-123), resp.name="y", exponent=.5)
# returns
# h(sqrt(y)-123) # the func knows that the special case of exponent=.5 is sqrt
fix.coef.names <- function(coef.names, resp.name, exponent)
{
if(length(coef.names) == 1)
return(coef.names) # do nothing if intercept only model
stopifnot(length(resp.name) == 1)
stopifnot(exponent > 0)
new.resp.name <-
if(exponent > .33 && exponent < .34)
sprint("cbrt(%s)", resp.name)
else if(exponent == .5)
sprint("sqrt(%s)", resp.name)
else if(exponent == 1)
resp.name
else if(exponent == 2)
sprint("sq(%s)", resp.name)
else
sprint("%s^%.3g", resp.name, exponent)
coef.names <- gsub("RHS", resp.name, coef.names, fixed=TRUE)
if(exponent == 1)
coef.names
else {
# TODO revisit, will fail if resp.name is substring of a token in
# coef.names or if resp.name="h" and coef.names="h(h-12)"
gsub(resp.name, new.resp.name, coef.names, fixed=TRUE)
}
}
# restore original exponent, it doesn't get scaled like the other coefficients
restore.exponent <- function(coef, org.coef, method)
{
if(method == "power")
coef[3] <- org.coef[3] # exponent is in coef[3]
else if(method == "power0")
coef[2] <- org.coef[2] # exponent is in coef[2]
coef
}
coef.varmod <- function(object, as.sd=TRUE, ...)
{
warn.if.dots(...)
coef <- coef(object$residmod)
if(is.null(coef))
stop0("coef.varmod: cannot get coefficients for \"",
class(object$residmod)[1], "\" residmod")
as.sd <- check.boolean(as.sd)
if(as.sd) {
org.coef <- coef
negs <- coef < 0
coef <- to.sd(abs(coef), object$lambda)
coef[negs] <- -coef[negs]
coef <- restore.exponent(coef, org.coef, object$method)
}
names(coef) <- fix.coef.names(names(coef),
colnames(object$parent.y), object$exponent)
coef
}
VARMOD.COEF.TAB.STYLES <- c("standard", "unit")
get.varmod.coef.tab <- function(
object,
style = VARMOD.COEF.TAB.STYLES)
{
style <- match.arg1(style, "style")
coef <- coef.varmod(object, as.sd=TRUE)
# if style="unit", normalize coef if possible
unit <- 1
if(style == "unit") {
# choose which coef will be the unit
if(length(coef) == 1) # method == "const"?
unit <- abs(coef[1])
else
unit <- abs(coef[2])
if(unit < 1e-3) {
warning0("coef=", unit, " is very small, forcing style=\"standard\"\n")
style <- "standard"
unit <- 1
}
}
org.coef <- coef
coef <- coef / unit
coef <- restore.exponent(coef, org.coef, object$method)
# get stderr
NAs <- repl(NA, length(coef))
stderr <- NAs
if(!is.null(object$iter.stderr)) {
org.stderr <- object$iter.stderr
stderr <- to.sd(org.stderr, object$lambda) / unit
stderr <- restore.exponent(stderr, org.stderr, object$method)
}
abs.coef <- abs(coef)
stderr.percent <- 100 * stderr / abs.coef
coef.names <- names(coef)
coef.names <- gsub("\`", "", coef.names, fixed=TRUE) # remove backquotes added by lm etc.
if(style == "unit") {
coef.tab <- data.frame(c(coef, unit), c(stderr, NA), c(stderr.percent, NA))
rownames(coef.tab) <- c(coef.names, "unit")
} else {
coef.tab <- data.frame(coef, stderr, stderr.percent)
rownames(coef.tab) <- coef.names
}
if(object$method %in% c("earth", "x.earth")) {
order <- reorder.earth(object$residmod, decomp="anova")
coef.tab <- coef.tab[order,]
}
colnames(coef.tab) <- c("coefficients", "iter.stderr", "iter.stderr.percent")
coef.tab
}
get.interval.tab <- function(object, level)
{
level <- check.level.arg(level, zero.ok=FALSE)
predict <- predict.varmod(object, type="pint", level=level)
interval <- predict$upr - predict$lwr
interval <- interval[order(interval)]
tab <- data.frame(
" ", mean(interval),
" ", interval[1],
" ", interval[length(interval)],
" ", interval[length(interval)] / interval[1])
colnames(tab) <- c(
" ", "mean",
" ", "smallest",
" ", "largest",
" ", "ratio")
rownames(tab) <- sprint("%g%% prediction interval", 100*level)
tab
}
percent.inconf <- function(object, level, parent.y, newdata)
{
predict <- predict.varmod(object, newdata, type="pint", level=level)
inconf <- parent.y >= predict$lwr & parent.y <= predict$upr
100 * sum(inconf) / length(inconf)
}
print_inconf_tab <- function(object, parent.y, newdata)
{
if(NCOL(parent.y) != 1) {
warning0("multiple response model: the table is for the first response")
parent.y <- parent.y[,1]
}
stopifnot(is.numeric(parent.y) || is.logical(parent.y))
inconf68 <- percent.inconf(object, .68, parent.y, newdata)
inconf80 <- percent.inconf(object, .80, parent.y, newdata)
inconf90 <- percent.inconf(object, .90, parent.y, newdata)
inconf95 <- percent.inconf(object, .95, parent.y, newdata)
# .5 below adjusts for rounding in printf %.0f
lt <- function(x, level) if(x < level-.5) "<" else " "
tab <- data.frame(
" ", sprint("%.0f%s ", inconf68, lt(inconf68, 68)),
" ", sprint("%.0f%s ", inconf80, lt(inconf80, 80)),
" ", sprint("%.0f%s ", inconf90, lt(inconf90, 90)),
" ", sprint("%.0f%s ", inconf95, lt(inconf95, 95)))
colnames(tab) <- c(
" ", "68% ",
" ", "80% ",
" ", "90% ",
" ", "95% ")
if(is.null(newdata))
rowname <- "response values in prediction interval"
else
rowname <- "newdata in prediction interval"
rownames(tab) <- rowname
print(tab)
# return value is the table but not in string form
tab <- data.frame(inconf68, inconf80, inconf90, inconf95)
colnames(tab) <- c("68%", "80%", "90%", "95%")
rownames(tab) <- rowname
tab
}
print.varmod <- function(
x = stop("no 'x' argument"), # x is a varmod object
level = .95, # use 0 to not print the interval tabs
style = "standard", # one of VARMOD.COEF.TAB.STYLES
digits = 2,
newdata = NULL,
...)
{
check.classname(x, substitute(x), "varmod")
object <- x # minimize confusion with x, the regression input matrix
remove(x) # not necessary but prevents mistakes later
warn.if.dots(...)
if(!is.null(newdata)) { # if newdata, print just the inconf table
object$inconf.tab <-
print_inconf_tab(object,
plotmo::plotmo_response(object$parent, newdata, trace=0, ...),
newdata)
return(invisible(object))
}
printf("method \"%s\"", object$method)
if(!is.null(object$package))
printf(" (%s package)", object$package)
space <- if(object$exponent != 1 || object$lambda != 1) "" else " "
if(object$exponent != 1)
printf("%s exponent %.3f", space, object$exponent)
if(object$lambda != 1)
printf("%s lambda %g", space, object$lambda)
printf("%s min.sd %.3g", space, object$min.sd)
if(!is.null(object$iter.rsq)) {
printf("%s iter.rsq %.3f", space, object$iter.rsq)
# TODO prints too many digits
# printf(dot.star.to.digits(", iter.rsq %.*f", digits+1), object$iter.rsq)
}
# coef tab
printf("\n\nstddev of predictions%s:\n",
if(style == "unit") " (scaled by unit)" else "")
tab <- object$coef.tab
if(is.null(tab)) { # needed if did not come here via summary.varmod
tab <- get.varmod.coef.tab(object, style)
object$coef.tab <- tab # for return value of this function
}
tab$coefficients <- zapsmall(tab$coefficients, digits+1)
# sprint below so print "NA" not "<NA>"
tab$iter.stderr <- sprint("%g", zapsmall(tab$iter.stderr, digits))
# convert iter.stderr.percent to character and print "big" if appropriate
tab$iter.stderr.percent.as.char <- sprint("%.0f", tab$iter.stderr.percent)
tab$iter.stderr.percent.as.char[tab$iter.stderr.percent >= 1e3] <- "big"
tab$iter.stderr.percent <- NULL
colnames(tab) <- c("coefficients", "iter.stderr", "iter.stderr%")
print(tab, digits=digits)
# interval and inconf tabs
level <- check.level.arg(level, zero.ok=TRUE)
if(is.specified(level)) {
stopifnot(level == object$level)
printf("\n")
tab <- object$interval.tab
if(is.null(tab)) {
tab <- get.interval.tab(object, level)
object$interval.tab <- tab # for return value of this function
}
print(tab, digits=digits)
printf("\n")
object$inconf.tab <- print_inconf_tab(object, object$parent.y, newdata=NULL)
}
invisible(object)
}
print.summary.varmod <- function(
x = stop("no 'x' argument"), # x is a summary.varmod object
level = x$level,
style = x$style,
digits = x$digits,
newdata = x$newdata,
...)
{
check.classname(x, substitute(x), "varmod")
warn.if.dots(...)
if(is.null(level))
level <- .95
if(is.null(style))
style <- "standard"
if(is.null(digits))
digits <- 2
if(!is.null(newdata)) # if newdata, print just the inconf table
print.varmod(x, level, style, digits, newdata)
else {
printcall("Parent model: ", x$parent$call)
printf("\n")
print.varmod(x, level, style, digits)
printf("\nRegression submodel (%s):\n", get.resids.name(x))
if(!(class(x$residmod)[1] %in% c("lm", "x.lm")))
printf("\n")
print(x$residmod, digits=digits)
}
invisible(x)
}
summary.varmod <- function(
object = stop("no 'object' argument"),
level = .95,
style = "standard", # one of VARMOD.COEF.TAB.STYLES
digits = 2,
newdata = NULL,
...)
{
check.classname(object, substitute(object), "varmod")
warn.if.dots(...)
object$level <- level # pass level on to print.summary.varmod
object$style <- style # ditto
object$digits <- digits # ditto
object$newdata <- newdata # ditto
object$coef.tab <- get.varmod.coef.tab(object, style)
object$interval.tab <- get.interval.tab(object, level)
# TODO add inconf table here too
class(object) <- c("summary.varmod", "varmod")
object
}
get.resids.name <- function(object)
{
if(object$lambda == 1)
sprint("Abs Residuals")
else if(object$lambda == 2)
sprint("Squared Residuals")
else
sprint("(Squared Residuals) ^ %.2g", object$lambda / 2)
}
get.varmod.ylab <- function(object, as.sd)
{
sprint("ParentMod %s", if(as.sd) "StdDev" else get.resids.name(object))
}
min_sd_line <- function(object, min.sd.col, lwd) # draw horizontal line at min.sd
{
if(is.specified(min.sd.col)) {
# TODO need to apply lambda exponent here?
abline(h=object$min.sd / lamba.factor.global, col=min.sd.col, lty=2, lwd=lwd)
}
}
sd.axis <- function(object) # draw righthand axis in standard deviation scale
{
# for righthand axis
sd <- lamba.factor.global * object$abs.resids ^ (1 / object$lambda)
pretty.sd <- pretty(range(sd))
axis(side=4, at=pretty.sd / lamba.factor.global, labels=pretty.sd, srt=90)
# the line setting depends on the axis margin lines
mtext(get.varmod.ylab(object, as.sd=TRUE), side=4, cex=par("cex"),
line=if(par("mgp")[1] < 1.8) 1.4 else 1.8)
}
plot.varmod <- function(
x = stop("no 'x' argument"),
which = 1:4,
do.par = NULL,
info = FALSE,
cex = NULL,
caption = NULL,
line.col = 2,
min.sd.col = line.col,
trace = 0,
...) # unused, for compat with the generic
{
check.classname(x, substitute(x), "varmod")
object <- x # minimize confusion with x, the regression input matrix
remove(x) # needed else plotmo.x gets this x instead of the x matrix
warn.if.dots(...)
trace <- as.numeric(check.integer.scalar(trace, logical.ok=TRUE))
info <- check.boolean(info)
check.index(which, "which", 1:4)
do.par <- check.do.par(do.par, length(which)) # do.par is 0, 1, or 2
# prepare caption --- we need it now for do.par() but
# can only display it later after at least one plot
stopifnot.string(caption, allow.empty=TRUE, null.ok=TRUE)
if(length(which) > 1 && do.par && is.null(caption)) # auto caption?
caption <- sprint("Variance Model method=\"%s\"\nParentMod: %s",
object$method, strip.deparse(object$parent$call))
main <- dota("main", ...)
if(do.par) {
oldpar <- par(no.readonly=TRUE)
do.par(nfigs=length(which), caption=caption, main1=main,
xlab1=NULL, ylab1=NULL, trace=trace,
nlines.in.main=2, def.right.mar=3,
def.font.main=1, # for compat with lm.plot
...)
if(do.par == 1)
on.exit(par(oldpar), add=TRUE)
} else { # do.par=FALSE
oldpar <- do.par.dots(..., trace=trace)
if(length(oldpar))
on.exit(do.call(par, oldpar), add=TRUE)
}
if(is.null(cex))
cex <- pt.cex(length(object$parent.y))
if(is.specified(main))
main <- repl(main, 4) # recycle for up to 4 plots
ylim <- fix.lim(c(min(object$abs.resids, 0), max(object$abs.resids)))
parent.fit <- parent.predict(object$parent, newdata=NULL)
order <- order(parent.fit)
smooth.col <- if(info) 2 else 0
lwd <- 1
for(iwhich in seq_along(which)) {
if(which[iwhich] == 1) { #--- fitted vs parent fitted ---
plot(parent.fit[order], object$abs.resids[order],
main=if(!is.specified(main))
sprint("%s vs Fitted", get.resids.name(object))
else
main[iwhich],
ylim=ylim, pch=20, cex=cex, xlab="Fitted",
ylab=get.varmod.ylab(object, as.sd=FALSE))
min_sd_line(object, min.sd.col, lwd) # horizontal line at min.sd
sd.axis(object) # right hand axis in stddev scale
# fitted values of residual model
fit <- predict.varmod(object, type="abs.residual")
lines(parent.fit[order], fit[order], col=line.col, lwd=lwd)
if(info) {
# lowess smooth
smooth <- lowess(parent.fit[order],
object$abs.resids[order], f=.5)
lines(smooth$x, smooth$y, col=smooth.col, lwd=1)
}
} else if(which[iwhich] == 2) { #--- fitted vs parent first pred ---
plotmo::plotmo(object, type="abs.residual",
ylim=ylim, degree1=1, degree2=0, do.par=FALSE,
trace=if(trace==0) -1 else trace,
pt.col=1, pt.cex=cex, degree1.col=line.col,
degree1.lwd=lwd, smooth.col=smooth.col,
ylab=get.varmod.ylab(object, as.sd=FALSE),
main=if(!is.specified(main))
sprint("%s vs First Predictor", get.resids.name(object))
else
main[iwhich])
min_sd_line(object, min.sd.col, lwd) # horizontal line at min.sd
sd.axis(object) # right hand axis in stddev scale
} else if(which[iwhich] == 3) { #--- residual plot ---
plotmo::plotres(object$residmod, which=3,
do.par=FALSE, center=FALSE,
xlab=get.varmod.ylab(object, as.sd=FALSE),
ylab="VarMod Residuals", info=info)
} else if(which[iwhich] == 4) { #--- model selection graph ---
if(class(object$residmod)[1] == "earth")
plot.earth(object$residmod, which=1, do.par=FALSE,
main=if(!is.specified(main)) "VarMod Model Selection"
else main[iwhich])
} else
stopf("plot.varmod: illegal value %g in 'which' argument",
which[iwhich])
}
draw.caption(caption, ...)
invisible()
}
# This func exists because when predicting for variance calculations with earth-glm
# models, we want to predict using the earth model itself (not the glm submodel).
# Therefore for earth models, we force type="earth".
parent.predict <- function(parent, newdata=NULL)
{
stopifnot(!is.null(parent))
type <- if(inherits(parent, "earth"))
"earth" # ignore glm submodel of earth model, if any
else
plotmo::plotmo_type(parent, trace, "varmod")
parent.fit <- predict(parent, newdata=newdata, type=type)
check.vec(parent.fit, "parent.fit")
stopifnot(!is.null(dim(parent.fit))) # check parent.fit is a matrix or dataframe
parent.fit[,1]
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.