R/plot.R
In sahirbhatnagar/sail: Sparse Additive Interaction Learning
Defines functions
plotInter
plotMain
plotSailCoef
Documented in
plotInter
plotMain
######################################
# R Source code file for plotting functions
# plotSailCoef is called by plot.sail and is not exported
# plotMain and plotInter are exported
# Author: Sahir Bhatnagar
# Created: 2016
# Updated: April 9, 2018
#####################################
#' @importFrom grDevices hcl
plotSailCoef <- function(coefs, lambda, group, df, dev, vnames, environ,
alpha = 1, legend.loc, label = FALSE, log.l = TRUE,
norm = FALSE, ...) {
# browser()
if (alpha < 0 | alpha > 1) {
warning("alpha must be in the range [0,1]. Setting alpha = 1")
alpha <- 1
}
if (norm) { # not implemented for now
# Y <- predict(x, type = "norm")
# index <- Y
# approx.f = 1
# if (any(x$group == 0))
# Y <- Y[-1, ]
# nonzero <- which(apply(abs(Y), 1, sum) != 0)
# Y <- Y[nonzero, ]
# g <- 1:nrow(Y)
} else {
if (length(dim(coefs)) == 3) {
beta <- matrix(coefs[, -1, , drop = FALSE], ncol = dim(coefs)[3])
} else {
beta <- coefs
}
penalized <- which(group != 0)
nonzero <- which(apply(abs(beta), 1, sum) != 0)
ind <- intersect(penalized, nonzero)
Y <- as.matrix(beta[ind, , drop = FALSE])
g <- as.numeric(as.factor(group[ind]))
}
p <- nrow(Y)
l <- lambda
n.g <- max(g)
if (log.l) {
l <- log(l)
index <- l
approx.f <- 0
xlab <- expression(log(lambda))
} else {
xlab <- expression(lambda)
index <- lambda
approx.f <- 0
}
ylims <- if (!missing(environ)) range(Y, environ) else range(Y)
plot.args <- list(
x = l, y = 1:length(l), ylim = ylims,
xlab = xlab, ylab = "", type = "n",
xlim = rev(range(l)),
# las = 1,
cex.lab = 1.5,
cex.axis = 1.5,
cex = 1.5,
# bty = "n",
# mai=c(1,1,0.1,0.2),
# tcl = -0.5,
# omi = c(0.2,1,0.2,0.2),
family = "serif"
)
new.args <- list(...)
if (length(new.args)) {
new.plot.args <- new.args[names(new.args) %in% c(
names(par()),
names(formals(plot.default))
)]
plot.args[names(new.plot.args)] <- new.plot.args
}
do.call("plot", plot.args)
if (plot.args$ylab == "") {
ylab <- if (norm) {
expression("||" * hat(theta) * "||")
} else {
expression(hat(theta))
}
mtext(ylab, 2, 3.5, las = 1, adj = 0, cex = 2)
}
abline(h = 0, lwd = 0.8, col = "gray")
cols <- hcl(
h = seq(15, 375, len = max(4, n.g + 1)), l = 60,
c = 150, alpha = alpha
)
cols <- if (n.g == 2) cols[c(1, 3)] else cols[1:n.g]
line.args <- list(
col = cols, lwd = 1 + 2 * exp(-p / 20),
lty = 1, pch = ""
)
if (length(new.args)) {
line.args[names(new.args)] <- new.args
}
line.args$x <- l
line.args$y <- t(Y)
line.args$col <- line.args$col[g]
line.args$lty <- rep(line.args$lty, length.out = max(g))
line.args$lty <- line.args$lty[g]
do.call("matlines", line.args)
if (!missing(environ)) lines(l, environ, lwd = line.args$lwd)
if (!missing(legend.loc)) {
legend.args <- list(
col = cols, lwd = line.args$lwd,
lty = line.args$lty, legend = vnames
)
if (length(new.args)) {
new.legend.args <- new.args[names(new.args) %in%
names(formals(legend))]
legend.args[names(new.legend.args)] <- new.legend.args
}
legend.args$x <- legend.loc
do.call("legend", legend.args)
}
if (label) {
ypos <- Y[, ncol(Y)]
text(min(l), ypos, names(ypos), xpd = NA, adj = c(
0,
NA
))
}
atdf <- pretty(index)
prettydf <- stats::approx(
x = index, y = df, xout = atdf, rule = 2,
method = "constant", f = approx.f
)$y
axis(3,
at = atdf, labels = prettydf, mgp = c(3, .3, 0),
tcl = NA,
cex.axis = 1.2
)
}
#' @title Plot Estimated Component Smooth Functions for Main Effects
#' @description Takes a fitted sail object produced by \code{sail()} or
#' \code{cv.sail()$sail.fit} and plots the component smooth function for a
#' pre-specified variable at a given value of lambda and on the scale of the
#' linear predictor. Currently only implemented for \code{type="gaussian"}
#' @param object a fitted \code{sail} object as produced by \code{sail()} or
#' \code{cv.sail()$sail.fit}
#' @param x original data supplied to the original call to \code{\link{sail}}
#' @param xvar a character corresponding to the predictor to be plotted. Only
#' one variable name should be supplied, if more than one is supplied, only
#' the first element will be plotted. This variable name must be in
#' \code{colnames(x)}.
#' @param s a single value of the penalty parameter \code{lambda} at which
#' coefficients will be extracted via the \code{coef} method for objects of
#' class \code{"sail"}. If more than one is supplied, only the first one will
#' be used.
#' @param f.truth true function. Only used for simulation purposes when the
#' truth is known. The function takes as a input a numeric vector
#' corresponding the \code{xvar} column in \code{x} of length \code{nrow(x)}.
#' A second line will be plotted for the truth and a legend is added to the
#' plot.
#' @param col color of the line. The first element corresponds to the color used
#' for the estimated function and the second element is for the true function
#' (if \code{f.truth} is specified). Default: c("#D55E00", "#009E73")
#' @param legend.position position of the legend. Only used when \code{f.truth}
#' is specified. Default: 'bottomleft'. Can be a single keyword from the list
#' "bottomright", "bottom", "bottomleft", "left", "topleft", "top",
#' "topright", "right" and "center". This places the legend on the inside of
#' the plot frame at the given location. Partial argument matching is used.
#' @param rug adds a rug representation (1-d plot) of the data to the plot, logical. Default: TRUE.
#' @param ... other graphical paramters passed to \code{plot}.
#' @return A plot is produced and nothing is returned
#' @details The linear predictor \eqn{basis(xvar) * \beta_xvar} is
#' plotted against \code{xvar}, where \code{basis} is the expansion provided
#' in the original call to \code{sail}.
#' @examples
#' f.basis <- function(i) splines::bs(i, degree = 3)
#' # Parallel
#' library(doParallel)
#' cl <- makeCluster(2)
#' registerDoParallel(cl)
#' cvfit <- cv.sail(x = sailsim$x, y = sailsim$y, e = sailsim$e,
#' parallel = TRUE, nlambda = 10,
#' maxit = 100, basis = f.basis,
#' nfolds = 3, dfmax = 10)
#' stopCluster(cl)
#' # plot cv-error curve
#' plot(cvfit)
#' # non-zero estimated coefficients at lambda.min
#' predict(cvfit, type = "nonzero", s="lambda.min")
#' # plot main effect for X4 with a line for the truth also
#' plotMain(cvfit$sail.fit, x = sailsim$x, xvar = "X4",
#' s = cvfit$lambda.min, f.truth = sailsim$f4)
#' @seealso \code{\link{coef.sail}} \code{\link{predict.sail}}, \code{\link[graphics]{rug}}
#' @rdname plotMain
#' @importFrom graphics abline axis legend lines mtext par plot.default segments text
#' @importFrom stats coef
#' @export
plotMain <- function(object, x, xvar, s, f.truth, col = c("#D55E00", "#009E73"),
legend.position = "bottomleft", rug = TRUE, ...) {
# browser()
if (length(xvar) > 1) {
xvar <- xvar[[1]]
warning("More than 1 xvar provided. Only first element will be plotted.")
}
if (length(s) > 1) {
s <- s[[1]]
warning("More than 1 s value provided. Only first element will be used for the estimated coefficients.")
}
ind <- object$group == which(object$vnames == xvar)
allCoefs <- coef(object, s = s)
a0 <- allCoefs[1, ]
betas <- as.matrix(allCoefs[object$main.effect.names[ind], , drop = FALSE])
design.mat <- object$design[, object$main.effect.names[ind], drop = FALSE]
originalX <- x[, unique(object$group[ind])]
# f.hat <- drop(a0 + design.mat %*% betas)
f.hat <- drop(design.mat %*% betas)
if (!missing(f.truth)) {
seqs <- seq(range(originalX)[1], range(originalX)[2], length.out = 100)
f.truth.eval <- f.truth(seqs)
ylims <- range(f.truth.eval, f.hat)
} else {
ylims <- range(f.hat)
}
plot.args <- list(
x = originalX[order(originalX)],
y = f.hat[order(originalX)],
ylim = ylims,
xlab = xvar,
ylab = sprintf("f(%s)", xvar),
type = "n",
# xlim = rev(range(l)),
# las = 1,
cex.lab = 1.5,
cex.axis = 1.5,
cex = 1.5,
# bty = "n",
# mai=c(1,1,0.1,0.2),
# tcl = -0.5,
# omi = c(0.2,1,0.2,0.2),
family = "serif"
)
new.args <- list(...)
if (length(new.args)) {
new.plot.args <- new.args[names(new.args) %in% c(
names(par()),
names(formals(plot.default))
)]
plot.args[names(new.plot.args)] <- new.plot.args
}
do.call("plot", plot.args)
abline(h = 0, lwd = 1, col = "gray")
lines(originalX[order(originalX)], f.hat[order(originalX)], col = col[1], lwd = 3)
if (rug) graphics::rug(originalX, side = 1)
if (!missing(f.truth)) {
lines(seqs[order(seqs)], f.truth.eval[order(seqs)], col = col[2], lwd = 3)
}
if (!missing(f.truth)) {
legend(legend.position,
c("Estimated", "Truth"),
col = col, cex = 1, bty = "n", lwd = 3
)
}
}
#' @title Plot Interaction Effects from sail object
#' @description Takes a fitted sail object produced by \code{sail()} or
#' \code{cv.sail()$sail.fit} and plots a \code{\link[graphics]{persp}} for a
#' pre-specified variable at a given value of lambda and on the scale of the
#' linear predictor. Currently only implemented for \code{type="gaussian"}
#' @inheritParams plotMain
#' @param f.truth true function. Only used for simulation purposes when the
#' truth is known. The function takes as a input two numeric vectors e.g.
#' \code{f(x,e)} corresponding the \code{xvar} column in \code{x} of length
#' \code{nrow(x)} and the exposure variable contained in the \code{sail
#' object}. A second \code{persp} will be plotted for the truth
#' @param interation.only if \code{TRUE} only the interaction part is used to
#' calculate the linear predictor, i.e., \eqn{linear predictor = E * f(X) *
#' interaction_effects}. If \code{FALSE}, then \eqn{linear predictor = E *
#' \beta_E + f(X) * interaction_effects + E * f(X) * interaction_effects}.
#' Default: TRUE
#' @param truthonly only plot the truth. \code{f.truth} must be specified if
#' this argument is set to \code{TRUE}. Default: FALSE
#' @param npoints number of points in the grid to calculate the perspective
#' plot. Default: 30
#' @param title_z title for the plot, Default: ''
#' @param xlab character for xlabel. if missing, variable name is used
#' @param ylab character for ylabel. if missing, variable name is used
#' @param zlab character for zlabel. if missing, variable name is used
#' @param ... currently ignored
#' @return A plot is produced and nothing is returned
#' @examples
#' f.basis <- function(i) splines::bs(i, degree = 3)
#' # Parallel
#' library(doParallel)
#' cl <- makeCluster(2)
#' registerDoParallel(cl)
#' cvfit <- cv.sail(x = sailsim$x, y = sailsim$y, e = sailsim$e,
#' parallel = TRUE, nlambda = 10,
#' maxit = 100, basis = f.basis,
#' nfolds = 3, dfmax = 10)
#' stopCluster(cl)
#' # plot cv-error curve
#' plot(cvfit)
#' # non-zero estimated coefficients at lambda.min
#' predict(cvfit, type = "nonzero", s="lambda.min")
#' # plot interaction effect for X4 and the true interaction effect also
#' plotInter(cvfit$sail.fit, x = sailsim$x, xvar = "X3",
#' f.truth = sailsim$f4.inter,
#' s = cvfit$lambda.min,
#' title_z = "Estimated")
#' @seealso \code{\link[graphics]{persp}} \code{\link{coef.sail}}
#' \code{\link{predict.sail}}, \code{\link[graphics]{rug}}
#' @rdname plotInter
#' @importFrom graphics abline axis legend lines mtext par plot.default segments text
#' @importFrom stats coef
#' @export
plotInter <- function(object, x, xvar, s, f.truth, interation.only = TRUE, truthonly = FALSE,
npoints = 30, col = c("#56B4E9", "#D55E00"), title_z = "", xlab, ylab, zlab,
...) {
# browser()
if (length(xvar) > 1) {
xvar <- xvar[[1]]
warning("More than 1 xvar provided. Only first element will be plotted.")
}
if (length(s) > 1) {
s <- s[[1]]
warning("More than 1 s value provided. Only first element will be used for the estimated coefficients.")
}
ind <- object$group == which(object$vnames == xvar)
allCoefs <- coef(object, s = s)
a0 <- allCoefs[1, ]
betas <- as.matrix(allCoefs[object$main.effect.names[ind], , drop = FALSE])
alphas <- as.matrix(allCoefs[object$interaction.names[ind], , drop = FALSE])
betaE <- as.matrix(allCoefs["E", , drop = FALSE])
design.mat.main <- object$design[, object$main.effect.names[ind], drop = FALSE]
design.mat.int <- object$design[, object$interaction.names[ind], drop = FALSE]
originalE <- object$design[, "E", drop = FALSE] # this is the centered E
originalX <- x[, unique(object$group[ind])]
# f.hat <- drop(a0 + design.mat %*% betas)
# f.hat <- drop(originalE %*% betaE + design.mat.main %*% betas + design.mat.int %*% alphas)
# all.equal((e * standardize(splines::bs(x, df = object$df, degree = object$degree))$x),
# sweep(standardize(splines::bs(x, df = object$df, degree = object$degree))$x, 1, e, FUN = "*"))
x <- seq(range(originalX)[1], range(originalX)[2], length.out = npoints)
e <- seq(range(originalE)[1], range(originalE)[2], length.out = npoints)
# show interaction effect only for simulation
if (interation.only) {
f.est <- function(X, E) {
# E * as.vector(betaE) + standardize(splines::bs(X, df = object$df, degree = object$degree))$x %*% betas +
(drop(standardize(E, center = object$center.e)$x) * standardize(object$basis(X), center = object$center.e)$x) %*% alphas
# (E * standardize(object$basis(X), center = FALSE)$x) %*% alphas
}
} else {
f.est <- function(X, E) {
standardize(E, center = object$center.e)$x * as.vector(betaE) +
standardize(object$basis(X), center = object$center.e)$x %*% betas +
(drop(standardize(E, center = object$center.e)$x) * standardize(object$basis(X), center = object$center.e)$x) %*% alphas
}
}
# f.truth <- function(x, e) { e * DT$f4(x) }
z.est <- outer(x, e, f.est)
if (!missing(f.truth)) {
z.truth <- outer(x, e, f.truth)
z_range <- c(min(z.est, z.truth), max(z.est, z.truth))
} else {
z_range <- c(min(z.est), max(z.est))
}
# c(bottom, left, top, right)
if (truthonly) {
oldpar <- par(mfrow = c(1, 1), tcl = -0.5, family = "serif", omi = c(0.2, 0.2, 0, 0),
mai = c(0., 0.2, 0.4, 0.), bg = "white")
on.exit(par(oldpar))
graphics::persp(x, e, z.truth,
zlim = z_range,
theta = 30, phi = 30,
ltheta = 120, expand = 0.5,
r = 2, shade = 0.3, axes = TRUE, scale = TRUE, box = T,
nticks = 5,
# ticktype="detailed",
col = col[2],
cex.lab = 3,
cex.main = 3,
xlab = if(missing(xlab)) sprintf("f(%s)", xvar) else xlab,
ylab = if(missing(ylab)) "X_E" else ylab,
zlab = if(missing(zlab)) "Y" else zlab,
main = "Truth"
)
} else if (!missing(f.truth)) {
oldpar <- par(mfrow = c(1, 2), tcl = -0.5, family = "serif", omi = c(0.2, 0.2, 0, 0),
mai = c(0., 0.8, 0.6, 0.), bg = "white")
on.exit(par(oldpar))
graphics::persp(x, e, z.truth,
zlim = z_range,
theta = 30, phi = 30,
ltheta = 120, expand = 0.5,
r = 2, shade = 0.3, axes = TRUE, scale = TRUE, box = T,
nticks = 5,
# ticktype="detailed",
col = col[2],
cex.lab = 3,
cex.main = 3,
xlab = if(missing(xlab)) sprintf("f(%s)", xvar) else xlab,
ylab = if(missing(ylab)) "X_E" else ylab,
zlab = if(missing(zlab)) "Y" else zlab, main = "Truth"
)
graphics::persp(x, e, z.est,
theta = 30, phi = 30,
ltheta = 120, expand = 0.5,
r = 2, shade = 0.3, axes = TRUE, scale = TRUE, box = T,
nticks = 5,
# zlim = z_range,
cex.lab = 3,
cex.main = 3,
# ticktype="detailed",
col = col[1],
xlab = if(missing(xlab)) sprintf("f(%s)", xvar) else xlab,
ylab = if(missing(ylab)) "X_E" else ylab,
zlab = if(missing(zlab)) "Y" else zlab,
main = title_z
)
} else {
oldpar <- par(mfrow = c(1, 1), tcl = -0.5, family = "serif", omi = c(0.2, 0.2, 0, 0),
mai = c(0., 0.2, 0.4, 0.), bg = "white")
on.exit(par(oldpar))
graphics::persp(x, e, z.est,
cex.lab = 3,
cex.main = 3,
theta = 30, phi = 30,
ltheta = 120, expand = 0.5,
r = 2, shade = 0.3, axes = TRUE, scale = TRUE, box = T,
nticks = 5,
zlim = z_range,
col = col[1],
# xlab=sprintf("f(x_%s)",as.numeric(gsub("X","",4))),
xlab = if(missing(xlab)) sprintf("f(%s)", xvar) else xlab,
ylab = if(missing(ylab)) "X_E" else ylab,
zlab = if(missing(zlab)) "Y" else zlab,
# main=sprintf("Estimated Interaction Effect for %s",xvar)
main = title_z
)
}
}
sahirbhatnagar/sail documentation built on July 17, 2021, 5:10 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions plotInter plotMain plotSailCoef
Documented in plotInter plotMain
######################################
# R Source code file for plotting functions
# plotSailCoef is called by plot.sail and is not exported
# plotMain and plotInter are exported
# Author: Sahir Bhatnagar
# Created: 2016
# Updated: April 9, 2018
#####################################
#' @importFrom grDevices hcl
plotSailCoef <- function(coefs, lambda, group, df, dev, vnames, environ,
alpha = 1, legend.loc, label = FALSE, log.l = TRUE,
norm = FALSE, ...) {
# browser()
if (alpha < 0 | alpha > 1) {
warning("alpha must be in the range [0,1]. Setting alpha = 1")
alpha <- 1
}
if (norm) { # not implemented for now
# Y <- predict(x, type = "norm")
# index <- Y
# approx.f = 1
# if (any(x$group == 0))
# Y <- Y[-1, ]
# nonzero <- which(apply(abs(Y), 1, sum) != 0)
# Y <- Y[nonzero, ]
# g <- 1:nrow(Y)
} else {
if (length(dim(coefs)) == 3) {
beta <- matrix(coefs[, -1, , drop = FALSE], ncol = dim(coefs)[3])
} else {
beta <- coefs
}
penalized <- which(group != 0)
nonzero <- which(apply(abs(beta), 1, sum) != 0)
ind <- intersect(penalized, nonzero)
Y <- as.matrix(beta[ind, , drop = FALSE])
g <- as.numeric(as.factor(group[ind]))
}
p <- nrow(Y)
l <- lambda
n.g <- max(g)
if (log.l) {
l <- log(l)
index <- l
approx.f <- 0
xlab <- expression(log(lambda))
} else {
xlab <- expression(lambda)
index <- lambda
approx.f <- 0
}
ylims <- if (!missing(environ)) range(Y, environ) else range(Y)
plot.args <- list(
x = l, y = 1:length(l), ylim = ylims,
xlab = xlab, ylab = "", type = "n",
xlim = rev(range(l)),
# las = 1,
cex.lab = 1.5,
cex.axis = 1.5,
cex = 1.5,
# bty = "n",
# mai=c(1,1,0.1,0.2),
# tcl = -0.5,
# omi = c(0.2,1,0.2,0.2),
family = "serif"
)
new.args <- list(...)
if (length(new.args)) {
new.plot.args <- new.args[names(new.args) %in% c(
names(par()),
names(formals(plot.default))
)]
plot.args[names(new.plot.args)] <- new.plot.args
}
do.call("plot", plot.args)
if (plot.args$ylab == "") {
ylab <- if (norm) {
expression("||" * hat(theta) * "||")
} else {
expression(hat(theta))
}
mtext(ylab, 2, 3.5, las = 1, adj = 0, cex = 2)
}
abline(h = 0, lwd = 0.8, col = "gray")
cols <- hcl(
h = seq(15, 375, len = max(4, n.g + 1)), l = 60,
c = 150, alpha = alpha
)
cols <- if (n.g == 2) cols[c(1, 3)] else cols[1:n.g]
line.args <- list(
col = cols, lwd = 1 + 2 * exp(-p / 20),
lty = 1, pch = ""
)
if (length(new.args)) {
line.args[names(new.args)] <- new.args
}
line.args$x <- l
line.args$y <- t(Y)
line.args$col <- line.args$col[g]
line.args$lty <- rep(line.args$lty, length.out = max(g))
line.args$lty <- line.args$lty[g]
do.call("matlines", line.args)
if (!missing(environ)) lines(l, environ, lwd = line.args$lwd)
if (!missing(legend.loc)) {
legend.args <- list(
col = cols, lwd = line.args$lwd,
lty = line.args$lty, legend = vnames
)
if (length(new.args)) {
new.legend.args <- new.args[names(new.args) %in%
names(formals(legend))]
legend.args[names(new.legend.args)] <- new.legend.args
}
legend.args$x <- legend.loc
do.call("legend", legend.args)
}
if (label) {
ypos <- Y[, ncol(Y)]
text(min(l), ypos, names(ypos), xpd = NA, adj = c(
0,
NA
))
}
atdf <- pretty(index)
prettydf <- stats::approx(
x = index, y = df, xout = atdf, rule = 2,
method = "constant", f = approx.f
)$y
axis(3,
at = atdf, labels = prettydf, mgp = c(3, .3, 0),
tcl = NA,
cex.axis = 1.2
)
}
#' @title Plot Estimated Component Smooth Functions for Main Effects
#' @description Takes a fitted sail object produced by \code{sail()} or
#' \code{cv.sail()$sail.fit} and plots the component smooth function for a
#' pre-specified variable at a given value of lambda and on the scale of the
#' linear predictor. Currently only implemented for \code{type="gaussian"}
#' @param object a fitted \code{sail} object as produced by \code{sail()} or
#' \code{cv.sail()$sail.fit}
#' @param x original data supplied to the original call to \code{\link{sail}}
#' @param xvar a character corresponding to the predictor to be plotted. Only
#' one variable name should be supplied, if more than one is supplied, only
#' the first element will be plotted. This variable name must be in
#' \code{colnames(x)}.
#' @param s a single value of the penalty parameter \code{lambda} at which
#' coefficients will be extracted via the \code{coef} method for objects of
#' class \code{"sail"}. If more than one is supplied, only the first one will
#' be used.
#' @param f.truth true function. Only used for simulation purposes when the
#' truth is known. The function takes as a input a numeric vector
#' corresponding the \code{xvar} column in \code{x} of length \code{nrow(x)}.
#' A second line will be plotted for the truth and a legend is added to the
#' plot.
#' @param col color of the line. The first element corresponds to the color used
#' for the estimated function and the second element is for the true function
#' (if \code{f.truth} is specified). Default: c("#D55E00", "#009E73")
#' @param legend.position position of the legend. Only used when \code{f.truth}
#' is specified. Default: 'bottomleft'. Can be a single keyword from the list
#' "bottomright", "bottom", "bottomleft", "left", "topleft", "top",
#' "topright", "right" and "center". This places the legend on the inside of
#' the plot frame at the given location. Partial argument matching is used.
#' @param rug adds a rug representation (1-d plot) of the data to the plot, logical. Default: TRUE.
#' @param ... other graphical paramters passed to \code{plot}.
#' @return A plot is produced and nothing is returned
#' @details The linear predictor \eqn{basis(xvar) * \beta_xvar} is
#' plotted against \code{xvar}, where \code{basis} is the expansion provided
#' in the original call to \code{sail}.
#' @examples
#' f.basis <- function(i) splines::bs(i, degree = 3)
#' # Parallel
#' library(doParallel)
#' cl <- makeCluster(2)
#' registerDoParallel(cl)
#' cvfit <- cv.sail(x = sailsim$x, y = sailsim$y, e = sailsim$e,
#' parallel = TRUE, nlambda = 10,
#' maxit = 100, basis = f.basis,
#' nfolds = 3, dfmax = 10)
#' stopCluster(cl)
#' # plot cv-error curve
#' plot(cvfit)
#' # non-zero estimated coefficients at lambda.min
#' predict(cvfit, type = "nonzero", s="lambda.min")
#' # plot main effect for X4 with a line for the truth also
#' plotMain(cvfit$sail.fit, x = sailsim$x, xvar = "X4",
#' s = cvfit$lambda.min, f.truth = sailsim$f4)
#' @seealso \code{\link{coef.sail}} \code{\link{predict.sail}}, \code{\link[graphics]{rug}}
#' @rdname plotMain
#' @importFrom graphics abline axis legend lines mtext par plot.default segments text
#' @importFrom stats coef
#' @export
plotMain <- function(object, x, xvar, s, f.truth, col = c("#D55E00", "#009E73"),
legend.position = "bottomleft", rug = TRUE, ...) {
# browser()
if (length(xvar) > 1) {
xvar <- xvar[[1]]
warning("More than 1 xvar provided. Only first element will be plotted.")
}
if (length(s) > 1) {
s <- s[[1]]
warning("More than 1 s value provided. Only first element will be used for the estimated coefficients.")
}
ind <- object$group == which(object$vnames == xvar)
allCoefs <- coef(object, s = s)
a0 <- allCoefs[1, ]
betas <- as.matrix(allCoefs[object$main.effect.names[ind], , drop = FALSE])
design.mat <- object$design[, object$main.effect.names[ind], drop = FALSE]
originalX <- x[, unique(object$group[ind])]
# f.hat <- drop(a0 + design.mat %*% betas)
f.hat <- drop(design.mat %*% betas)
if (!missing(f.truth)) {
seqs <- seq(range(originalX)[1], range(originalX)[2], length.out = 100)
f.truth.eval <- f.truth(seqs)
ylims <- range(f.truth.eval, f.hat)
} else {
ylims <- range(f.hat)
}
plot.args <- list(
x = originalX[order(originalX)],
y = f.hat[order(originalX)],
ylim = ylims,
xlab = xvar,
ylab = sprintf("f(%s)", xvar),
type = "n",
# xlim = rev(range(l)),
# las = 1,
cex.lab = 1.5,
cex.axis = 1.5,
cex = 1.5,
# bty = "n",
# mai=c(1,1,0.1,0.2),
# tcl = -0.5,
# omi = c(0.2,1,0.2,0.2),
family = "serif"
)
new.args <- list(...)
if (length(new.args)) {
new.plot.args <- new.args[names(new.args) %in% c(
names(par()),
names(formals(plot.default))
)]
plot.args[names(new.plot.args)] <- new.plot.args
}
do.call("plot", plot.args)
abline(h = 0, lwd = 1, col = "gray")
lines(originalX[order(originalX)], f.hat[order(originalX)], col = col[1], lwd = 3)
if (rug) graphics::rug(originalX, side = 1)
if (!missing(f.truth)) {
lines(seqs[order(seqs)], f.truth.eval[order(seqs)], col = col[2], lwd = 3)
}
if (!missing(f.truth)) {
legend(legend.position,
c("Estimated", "Truth"),
col = col, cex = 1, bty = "n", lwd = 3
)
}
}
#' @title Plot Interaction Effects from sail object
#' @description Takes a fitted sail object produced by \code{sail()} or
#' \code{cv.sail()$sail.fit} and plots a \code{\link[graphics]{persp}} for a
#' pre-specified variable at a given value of lambda and on the scale of the
#' linear predictor. Currently only implemented for \code{type="gaussian"}
#' @inheritParams plotMain
#' @param f.truth true function. Only used for simulation purposes when the
#' truth is known. The function takes as a input two numeric vectors e.g.
#' \code{f(x,e)} corresponding the \code{xvar} column in \code{x} of length
#' \code{nrow(x)} and the exposure variable contained in the \code{sail
#' object}. A second \code{persp} will be plotted for the truth
#' @param interation.only if \code{TRUE} only the interaction part is used to
#' calculate the linear predictor, i.e., \eqn{linear predictor = E * f(X) *
#' interaction_effects}. If \code{FALSE}, then \eqn{linear predictor = E *
#' \beta_E + f(X) * interaction_effects + E * f(X) * interaction_effects}.
#' Default: TRUE
#' @param truthonly only plot the truth. \code{f.truth} must be specified if
#' this argument is set to \code{TRUE}. Default: FALSE
#' @param npoints number of points in the grid to calculate the perspective
#' plot. Default: 30
#' @param title_z title for the plot, Default: ''
#' @param xlab character for xlabel. if missing, variable name is used
#' @param ylab character for ylabel. if missing, variable name is used
#' @param zlab character for zlabel. if missing, variable name is used
#' @param ... currently ignored
#' @return A plot is produced and nothing is returned
#' @examples
#' f.basis <- function(i) splines::bs(i, degree = 3)
#' # Parallel
#' library(doParallel)
#' cl <- makeCluster(2)
#' registerDoParallel(cl)
#' cvfit <- cv.sail(x = sailsim$x, y = sailsim$y, e = sailsim$e,
#' parallel = TRUE, nlambda = 10,
#' maxit = 100, basis = f.basis,
#' nfolds = 3, dfmax = 10)
#' stopCluster(cl)
#' # plot cv-error curve
#' plot(cvfit)
#' # non-zero estimated coefficients at lambda.min
#' predict(cvfit, type = "nonzero", s="lambda.min")
#' # plot interaction effect for X4 and the true interaction effect also
#' plotInter(cvfit$sail.fit, x = sailsim$x, xvar = "X3",
#' f.truth = sailsim$f4.inter,
#' s = cvfit$lambda.min,
#' title_z = "Estimated")
#' @seealso \code{\link[graphics]{persp}} \code{\link{coef.sail}}
#' \code{\link{predict.sail}}, \code{\link[graphics]{rug}}
#' @rdname plotInter
#' @importFrom graphics abline axis legend lines mtext par plot.default segments text
#' @importFrom stats coef
#' @export
plotInter <- function(object, x, xvar, s, f.truth, interation.only = TRUE, truthonly = FALSE,
npoints = 30, col = c("#56B4E9", "#D55E00"), title_z = "", xlab, ylab, zlab,
...) {
# browser()
if (length(xvar) > 1) {
xvar <- xvar[[1]]
warning("More than 1 xvar provided. Only first element will be plotted.")
}
if (length(s) > 1) {
s <- s[[1]]
warning("More than 1 s value provided. Only first element will be used for the estimated coefficients.")
}
ind <- object$group == which(object$vnames == xvar)
allCoefs <- coef(object, s = s)
a0 <- allCoefs[1, ]
betas <- as.matrix(allCoefs[object$main.effect.names[ind], , drop = FALSE])
alphas <- as.matrix(allCoefs[object$interaction.names[ind], , drop = FALSE])
betaE <- as.matrix(allCoefs["E", , drop = FALSE])
design.mat.main <- object$design[, object$main.effect.names[ind], drop = FALSE]
design.mat.int <- object$design[, object$interaction.names[ind], drop = FALSE]
originalE <- object$design[, "E", drop = FALSE] # this is the centered E
originalX <- x[, unique(object$group[ind])]
# f.hat <- drop(a0 + design.mat %*% betas)
# f.hat <- drop(originalE %*% betaE + design.mat.main %*% betas + design.mat.int %*% alphas)
# all.equal((e * standardize(splines::bs(x, df = object$df, degree = object$degree))$x),
# sweep(standardize(splines::bs(x, df = object$df, degree = object$degree))$x, 1, e, FUN = "*"))
x <- seq(range(originalX)[1], range(originalX)[2], length.out = npoints)
e <- seq(range(originalE)[1], range(originalE)[2], length.out = npoints)
# show interaction effect only for simulation
if (interation.only) {
f.est <- function(X, E) {
# E * as.vector(betaE) + standardize(splines::bs(X, df = object$df, degree = object$degree))$x %*% betas +
(drop(standardize(E, center = object$center.e)$x) * standardize(object$basis(X), center = object$center.e)$x) %*% alphas
# (E * standardize(object$basis(X), center = FALSE)$x) %*% alphas
}
} else {
f.est <- function(X, E) {
standardize(E, center = object$center.e)$x * as.vector(betaE) +
standardize(object$basis(X), center = object$center.e)$x %*% betas +
(drop(standardize(E, center = object$center.e)$x) * standardize(object$basis(X), center = object$center.e)$x) %*% alphas
}
}
# f.truth <- function(x, e) { e * DT$f4(x) }
z.est <- outer(x, e, f.est)
if (!missing(f.truth)) {
z.truth <- outer(x, e, f.truth)
z_range <- c(min(z.est, z.truth), max(z.est, z.truth))
} else {
z_range <- c(min(z.est), max(z.est))
}
# c(bottom, left, top, right)
if (truthonly) {
oldpar <- par(mfrow = c(1, 1), tcl = -0.5, family = "serif", omi = c(0.2, 0.2, 0, 0),
mai = c(0., 0.2, 0.4, 0.), bg = "white")
on.exit(par(oldpar))
graphics::persp(x, e, z.truth,
zlim = z_range,
theta = 30, phi = 30,
ltheta = 120, expand = 0.5,
r = 2, shade = 0.3, axes = TRUE, scale = TRUE, box = T,
nticks = 5,
# ticktype="detailed",
col = col[2],
cex.lab = 3,
cex.main = 3,
xlab = if(missing(xlab)) sprintf("f(%s)", xvar) else xlab,
ylab = if(missing(ylab)) "X_E" else ylab,
zlab = if(missing(zlab)) "Y" else zlab,
main = "Truth"
)
} else if (!missing(f.truth)) {
oldpar <- par(mfrow = c(1, 2), tcl = -0.5, family = "serif", omi = c(0.2, 0.2, 0, 0),
mai = c(0., 0.8, 0.6, 0.), bg = "white")
on.exit(par(oldpar))
graphics::persp(x, e, z.truth,
zlim = z_range,
theta = 30, phi = 30,
ltheta = 120, expand = 0.5,
r = 2, shade = 0.3, axes = TRUE, scale = TRUE, box = T,
nticks = 5,
# ticktype="detailed",
col = col[2],
cex.lab = 3,
cex.main = 3,
xlab = if(missing(xlab)) sprintf("f(%s)", xvar) else xlab,
ylab = if(missing(ylab)) "X_E" else ylab,
zlab = if(missing(zlab)) "Y" else zlab, main = "Truth"
)
graphics::persp(x, e, z.est,
theta = 30, phi = 30,
ltheta = 120, expand = 0.5,
r = 2, shade = 0.3, axes = TRUE, scale = TRUE, box = T,
nticks = 5,
# zlim = z_range,
cex.lab = 3,
cex.main = 3,
# ticktype="detailed",
col = col[1],
xlab = if(missing(xlab)) sprintf("f(%s)", xvar) else xlab,
ylab = if(missing(ylab)) "X_E" else ylab,
zlab = if(missing(zlab)) "Y" else zlab,
main = title_z
)
} else {
oldpar <- par(mfrow = c(1, 1), tcl = -0.5, family = "serif", omi = c(0.2, 0.2, 0, 0),
mai = c(0., 0.2, 0.4, 0.), bg = "white")
on.exit(par(oldpar))
graphics::persp(x, e, z.est,
cex.lab = 3,
cex.main = 3,
theta = 30, phi = 30,
ltheta = 120, expand = 0.5,
r = 2, shade = 0.3, axes = TRUE, scale = TRUE, box = T,
nticks = 5,
zlim = z_range,
col = col[1],
# xlab=sprintf("f(x_%s)",as.numeric(gsub("X","",4))),
xlab = if(missing(xlab)) sprintf("f(%s)", xvar) else xlab,
ylab = if(missing(ylab)) "X_E" else ylab,
zlab = if(missing(zlab)) "Y" else zlab,
# main=sprintf("Estimated Interaction Effect for %s",xvar)
main = title_z
)
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.