Nothing
R/regplots.R
In procs: Recreates Some 'SAS®' Procedures in 'R'
Defines functions
remove_na_rows
get_reg_xlim
get_zero_scale
get_scale
draw_stats_box
collect_stats
render_dfbetas
render_dffits
render_residualboxplot
render_observedbypredicted
render_cooksd
render_rfplot
render_qqplot
render_residualhistogram
render_fitplot
render_rstudentbyleverage
render_rstudentbypredicted
render_residualbypredicted
render_residuals
render_diagnostics
regplot
Documented in
regplot
# Regplots Main Function --------------------------------------------------
# RESIDUALBYPREDICTED: Residuals vs. Predicted values.
# RSTUDENTBYPREDICTED: Externally Studentized Residuals (RStudent) vs. Predicted values.
# RSTUDENTBYLEVERAGE: Externally Studentized Residuals vs. Leverage.
# QQ (or QQPLOT): Normal Quantile-Quantile (Q-Q) plot of residuals.
# OBSERVEDBYPREDICTED: Dependent variable (Observed) vs. Predicted values.
# COOKSD: Cook’s D statistic vs. Observation number.
# RESIDUALHISTOGRAM: Histogram of residuals.
# RFPLOT: Residual-Fit (RF) spread plot.
# ALL: Requests all plots relevant to the current analysis and options.
# FITPLOT: Produces a scatter plot of the dependent variable against the regressor, including the fitted line and confidence/prediction bands. This is only available for models with a single regressor.
# RESIDUALS or RESIDUALPLOT: Produces a panel of residual plots against each independent variable in the model.
# PARTIAL (or PARTIALPLOT): Requests partial regression residual plots for each regressor. This requires the PARTIAL option to be specified in the MODEL statement.
# CRITERIONPANEL: When using model selection (like SELECTION=STEPWISE), this displays how criteria like AIC, SBC, or R-square change across different model sizes.
# panel = TRUE/FALSE
# type = diagnostics, residualbypredicted, RStudentbypredicted, rstudentbyleverage,
# qq/qqplot, observedbypredicted, cooksd, residualhistogram, rfplot, residuals/residualplot, fitplot
# stats = default aic sbc
# all = obs, edf, rsquare, sse, coefvar, bic, gmsep, pc, sp, parms, mse, adjrsquare, depmean, aic, cp, jp, sbc
#' @title Request Regression Plots
#' @description A function to request regression plots on a call to
#' \code{\link{proc_reg}}. The function allows you to specify the type of
#' regression plots to produce. It produces a combined diagnostics panel,
#' residuals dot plot, and a fit plot by default. You may also specify
#' individual plots by passing a vector of plot names
#' on the "type" parameter.
#' @details There are many types of regression plots. The plots have different
#' uses. Some of the plots help you assess normality of the data. Other plots
#' help you assess the quality of the model. You can select the type of plots
#' you want by passing a vector of plot names on the \code{type} parameter.
#'
#' Here is a list of possible plot types and a short description of each:
#' \itemize{
#' \item{\strong{diagnostics}: A fit diagnostics panel that contains 8 different
#' types of plots and a table of statistics.
#' }
#' \item{\strong{cooksd}: Cook’s D statistic vs. Observation number.
#' }
#' \item{\strong{dfbetas}: Displays the influence of each observation
#' for each coefficient in the model.
#' }
#' \item{\strong{dffits}: Displays the influence of each observation on fitted values.
#' }
#' \item{\strong{fitplot}: Produces a scatter plot of the dependent variable
#' against the regressor, including the fitted line and confidence/prediction bands.
#' The fitplot is only available for models with a single regressor.
#' }
#' \item{\strong{observedbypredicted}: Dependent variable (Observed) vs. Predicted values.
#' }
#' \item{\strong{qqplot}: Normal Quantile-Quantile (Q-Q) plot of residuals.
#' }
#' \item{\strong{residuals}: Produces a panel of residual plots against
#' each independent variable in the model.
#' }
#' \item{\strong{residualboxplot}: A box plot of residual values.
#' }
#' \item{\strong{residualbypredicted}: Residuals vs. Predicted values.
#' }
#' \item{\strong{residualhistogram}: Histogram of residuals,
#' with a normal and kernel curve overlay.
#' }
#' \item{\strong{rfplot}: Residual-Fit (RF) spread plot.
#' }
#' \item{\strong{rstudentbyleverage}: Externally Studentized Residuals vs. Leverage.
#' }
#' \item{\strong{rstudentbypredicted}: Externally Studentized Residuals (RStudent) vs. Predicted values.
#' }
#' }
#'
#' If possible, the statistics from the tabular report are used for the plots.
#' Otherwise, additional statistics functions are called to produce the
#' needed values.
#'
#' @section Plot Statistics:
#' The diagnostics panel and the fit plot each contain a small table of
#' statistics. This table is customizable. The following keywords may be
#' used to to customize the statistics table:
#' \itemize{
#' \item{\strong{adjrsq}: Adjusted R-square.
#' }
#' \item{\strong{aic}: Akaike's information criterion.
#' }
#' \item{\strong{coeffvar}: Coefficient of variation.
#' }
#' \item{\strong{depmean}: Mean of dependent.
#' }
#' \item{\strong{default}: A set of default statistics.
#' }
#' \item{\strong{edf}: Error degrees of freedom.
#' }
#' \item{\strong{mse}: Mean squared error.
#' }
#' \item{\strong{nobs}: Number of observations used.
#' }
#' \item{\strong{nparm}: Number of parameters in the model (including the intercept).
#' }
#' \item{\strong{rsquare}: The R-square statistic.
#' }
#' \item{\strong{sse}: Error sum of squares.
#' }
#' }
#' To use these keywords, pass them as a vector to the \code{stats}
#' parameter on the \code{regplot} function.
#' @section Labeling Outliers:
#' Some types of plots can be used to identify outliers in the data. For instance,
#' the Cook's D chart is excellent for such a task.
#' When identifying outliers, it is helpful to have them labelled on the chart.
#' The labels make it possible to trace the outliers back to the source
#' data.
#'
#' To get a basic row/observation label, set the \code{labels} parameter to TRUE.
#' If there is a column in the data that can be used to identify an individual
#' record, pass that column name on the \code{id} parameter. The values
#' from that column will then be used as outlier labels on those charts that
#' support labels.
#'
#' @section Additional Information:
#' Regression plots will be displayed on interactive reports only.
#' Plots are created as jpeg files, and stored in a temp directory. Those
#' temporary files are then referenced by the interactive report to display
#' the graphic.
#'
#' If desired, you may
#' output the report objects and pass to \code{\link{proc_print}}. To do this,
#' set \code{output = report} on the call to \code{\link{proc_reg}}, and pass
#' the entire list to \code{\link{proc_print}}.
#'
#' @section Discrepancies with SAS:
#' The histogram binning algorithm in R is different from the binning algorithm
#' in SAS®. R uses a "Sturges" algorithm, which more accurately represents
#' the distribution of the data. This algorithm may produce a different number
#' of bins and different height of bars than the corresponding SAS chart.
#'
#' Plots generated by \code{proc_reg} may also have some spacing, color,
#' and shape discrepancies with the corresponding SAS charts. The information
#' conveyed is expected to be similar.
#'
#' @param type The type(s) of plot to create. Multiple types should be passed
#' as a vector of strings. Valid values are "diagnostics", "cooksd", "dfbetas",
#' "dffits", "fitplot", "observedbypredicted", "qqplot", "residuals",
#' "residualboxplot", "residualbypredicted", "residualhistogram", "rfplot",
#' "rstudentbyleverage", and "rstudentbypredicted". The default value is
#' a vector with "diagnostics", "residuals", and "fitplot". The "diagnostics"
#' keyword produces a single combined chart with 8 different plots and a
#' selection of statistics in a small table. The statistics can be controlled
#' by the \code{stats} parameter. You may also pass the "all" keyword to produce
#' all charts available for the analysis.
#' @param panel Whether or not to display the diagnostics plots combined into
#' in a single panel. Default is TRUE. A value of FALSE will create
#' individual plots instead. This parameter is equivalent to the
#' "unpack" keyword in SAS.
#' @param stats The statistics to display on the diagnostics panel or fit plot.
#' Valid values
#' are: "adjrsq", "aic", "coeffvar", "depmean", "default", "edf", "mse", "nobs",
#' "nparm", "rsquare", and "sse". The
#' default value is "default", which produces the following statistics:
#' "nobs", "nparm", "edf", "mse", "rsquare", and "adjrsq". You may also pass
#' the value "none" if you do not want any statistics shown on the chart. In the
#' case of the diagnostics panel, if the statistics table is removed, it
#' will be replaced with a residual box plot.
#' @param label Whether or not to label outlier values. Valid values
#' are TRUE or FALSE. Default is FALSE. If TRUE, this option will assign
#' labels to outlier values on some charts. Charts that support labels are
#' as follows: "diagnostics", "cooksd", "dffits", "dfbetas", "rstudentbypredicted",
#' and "rstudentbyleverage".
#' @param id If the \code{label} parameter is TRUE, this parameter determines
#' which value is assigned to the label. By default, the row number will
#' be assigned. You may also assign a column name from the input dataset
#' to this parameter to use as the label value.
#' @seealso [proc_reg()]
#' @examples
#' library(procs)
#'
#' # Turn off printing for CRAN checks
#' # Set to TRUE to run in local environment
#' options("procs.print" = FALSE)
#'
#'
#' # Example 1: Regression statistics with default plots
#' res <- proc_reg(iris, model = "Sepal.Length = Petal.Length",
#' output = report,
#' plots = TRUE,
#' titles = "Iris Regression Statistics")
#'
#' # View results
#' res
#'
#' # Example 2: Regression statistics with custom plot request and by variable
#' res <- proc_reg(iris, model = "Sepal.Length = Petal.Length",
#' output = report,
#' by = Species,
#' plots = regplot(type = v(residualhistogram,
#' rstudentbypredicted, rstudentbyleverage),
#' label = TRUE),
#' titles = "Iris Regression Statistics")
#'
#' # View results
#' res
#'
#' # Example 3: Regression statistics with multiple models, same plot string
#' res <- proc_reg(iris, model = c("Sepal.Length = Petal.Length",
#' "Sepal.Length = Sepal.Width",
#' "Sepal.Length = Petal.Width"),
#' output = report,
#' plots = "diagnostics",
#' titles = "Iris Regression Statistics")
#'
#' # View results
#' res
#'
#' # Example 4: Regression statistics with multiple models, different plot strings
#' res <- proc_reg(iris, model = c("Sepal.Length = Petal.Length",
#' "Sepal.Length = Sepal.Width",
#' "Sepal.Length = Petal.Width"),
#' output = report,
#' plots = list("diagnostics",
#' "residualhistogram",
#' "fitplot"),
#' titles = "Iris Regression Statistics")
#'
#' # View results
#' res
#'
#' # Example 5: Regression statistics with multiple models, different plot functions
#' res <- proc_reg(iris, model = c("Sepal.Length = Petal.Length",
#' "Sepal.Length = Sepal.Width",
#' "Sepal.Length = Petal.Width"),
#' output = report,
#' plots = list(regplot(type = "diagnostics"),
#' regplot(type = "cooksd",
#' label = TRUE),
#' regplot(type = "fitplot",
#' stats = c("nobs", "mse", "rsquare"))),
#' titles = "Iris Regression Statistics")
#'
#' # View results
#' res
#'
#' # Example 6: Regression statistics with multiple models, influence charts
#' res <- proc_reg(iris, model = c("Sepal.Length = Petal.Length",
#' "Sepal.Length = Petal.Length Petal.Width",
#' "Sepal.Length = Petal.Length Petal.Width Sepal.Width"),
#' output = report,
#' plots = regplot(v(cooksd, dffits, dfbetas), label = TRUE),
#' titles = "Iris Regression Statistics")
#'
#' # View results
#' res
#' @export
regplot <- function(type = c("diagnostics", "residuals", "fitplot"), panel = TRUE,
stats = "default", label = FALSE, id = NULL) {
# Non-standard evaluation
type <- resolve_arg(type)
stats <- resolve_arg(stats)
id <- resolve_arg(id)
# Parameter Checks
vldvals <- c("diagnostics", "residuals", "fitplot", 'qqplot', 'rfplot',
"residualbypredicted", "rstudentbypredicted", "rstudentbyleverage",
"cooksd", "residualhistogram", "observedbypredicted", "residualboxplot",
"all", "dffits", "dfbetas")
if (any(!type %in% vldvals)) {
ivd <- type[!type %in% vldvals]
stop(paste0("Parameter value for 'type' invalid: ", paste0("'", ivd, "'", collapse = ", "),
"\nValid values are: ",
"'all', 'diagnostics', 'residuals', 'fitplot', 'qqplot', 'rfplot', ",
"'dffits', 'dfbetas', ",
"'residualbypredicted', 'rstudentbypredicted', 'rstudentbyleverage', ",
"'cooksd', 'residualhistogram', 'observedbypredicted', 'residualboxplot'."
))
}
if (!panel %in% c(TRUE, FALSE)) {
stop("Parameter value for 'panel' invalid. Valid values are TRUE or FALSE.")
}
# sse, bic, gmsep, pc, sp, aic, cp, jp, sbc
statlst <- c("adjrsq", "aic", "coeffvar", "depmean", "default", "edf", "mse",
"nobs", "none", "nparm", "rsquare", "sse")
if (!all(stats %in% statlst)) {
stop(paste0("Parameter value for 'stats' invalid. Valid values are ",
paste0("'", statlst, "'", collapse = ", "), "."))
}
# Create object
ret <- structure(list(), class = c("regplot", "list"))
ret$stats <- tolower(stats)
ret$panel <- panel
ret$type <- tolower(type)
ret$label <- label
ret$id <- id
return(ret)
}
# Plot Rendering ----------------------------------------------------------
#' @noRd
render_regplot <- function (dat, res, mdl, plt, alph) {
ret <- NULL
if (is.character(plt)) {
if (all(plt == "regplot")) {
plt <- regplot()
} else {
plt <- regplot(type = plt)
}
}
if ("regplot" %in% class(plt)) {
typs <- plt$type
# All logic
if (all(typs == "all")) {
typs <- c("residualhistogram", "residualbypredicted", "rstudentbypredicted",
"observedbypredicted", "cooksd", "rstudentbyleverage", 'qqplot', 'rfplot',
"residualboxplot", "diagnostics", "dffits", "dfbetas", "residuals", "fitplot")
}
# "Unpack" logic
if ("diagnostics" %in% typs & plt$panel == FALSE) {
pos <- match("diagnostics", typs)
vins <- c("residualbypredicted", "rstudentbypredicted", "rstudentbyleverage",
"qqplot", "observedbypredicted", "cooksd", "residualhistogram", "rfplot")
nv <- append(typs, vins, after = pos)
typs <- nv[-pos]
}
ret <- list()
# "diagnostics", "residuals", "fitplot", 'qqplot', 'rfplot',
# "residualbypredicted", "rstudentbypredicted", "rstudentbyleverage",
# "cooksd", "residualhistogram", "observedbypredicted"
for (tp in typs) {
if (tp == "diagnostics") {
ret[["diagnostics"]] <- render_diagnostics(dat, res, mdl, plt, alph)
} else if (tp == "residuals") { # Can be more than 1
resplts <- render_residuals(dat, res, mdl)
if ("plot_spec" %in% class(resplts)) {
ret[["residuals"]] <- resplts
} else {
for (idx in seq_len(length(resplts))) {
ret[[paste0("residuals", idx)]] <- resplts[[idx]]
}
}
} else if (tp == "fitplot") {
ret[["fitplot"]] <- render_fitplot(dat, res, mdl, plt, alph)
} else if (tp == "qqplot") {
ret[["qqplot"]] <- render_qqplot(dat, res, mdl)
} else if (tp == "rfplot") {
ret[["rfplot"]] <- render_rfplot(dat, res, mdl)
} else if (tp == "cooksd") {
ret[["cooksd"]] <- render_cooksd(dat, res, mdl, plt)
} else if (tp == "residualbypredicted") {
ret[["residualbypredicted"]] <- render_residualbypredicted(dat, res, mdl)
} else if (tp == "rstudentbypredicted") {
ret[["rstudentbypredicted"]] <- render_rstudentbypredicted(dat, res, mdl, plt)
} else if (tp == "rstudentbyleverage") {
ret[["rstudentbyleverage"]] <- render_rstudentbyleverage(dat, res, mdl, plt)
} else if (tp == "observedbypredicted") {
ret[["observedbypredicted"]] <- render_observedbypredicted(dat, res, mdl, plt)
} else if (tp == "residualhistogram") {
ret[["residualhistogram"]] <- render_residualhistogram(dat, res, mdl, plt)
} else if (tp == "residualboxplot") {
ret[["residualboxplot"]] <- render_residualboxplot(dat, res, mdl, plt)
} else if (tp == "dffits") {
ret[["dffits"]] <- render_dffits(dat, res, mdl, plt)
} else if (tp == "dfbetas") {
resplts <- render_dfbetas(dat, res, mdl, plt)
for (idx in seq_len(length(resplts))) {
ret[[paste0("dfbetas", idx)]] <- resplts[[idx]]
}
}
}
if (length(ret) == 0) {
ret <- NULL
}
}
return(ret)
}
# Diagnostics Panel -------------------------------------------------------
# Deal with "none" stats parameter. Adding extra box plot.
#' @import grDevices
#' @noRd
render_diagnostics <- function(dat, res, mdl, plt, alph) {
op <- par("mar")
om <- par("oma")
ofig <- par("fig")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 6 # Height in inches
wdi <- 6 # Width in inches
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Remove rows with NAs
datclean <- remove_na_rows(dat, dvr, ivr)
# Set margins
par(mfrow = c(3, 3),
mar = c(2, 2, 1, 0) + 0.1,
oma = c(0, 0, 3, .75))
#******************************
# 1) Residual By Predicted
#******************************
# Set margins
par(mar = c(4, 4, 0, 0) + 0.1)
# Prepare data
rdt <- res$OutputStatistics$RESID
pdt <- res$OutputStatistics$PREVAL
# Get scales
xscl <- get_scale(pdt, .05)
yscl <- get_zero_scale(rdt, .05)
# Generate plot
plot(pdt, rdt,
main = "",
xlab = "",
ylab = "",
pch = 1, # open circles
cex = 1.3,
col = "#05379B",
xlim = xscl,
ylim = yscl,
axes = FALSE
)
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# Generate zero line
abline(h = 0, col = "grey60", lwd = 1)
# Add labels
mtext("Predicted Value", side = 1, line = par("mar")[1] - 2, cex = .9)
mtext("Residual", side = 2, line = par("mar")[2] - 2, cex = .9)
# Frame
box(col = "grey70", lwd = 1)
#******************************
# 2) RStudent By Predicted
#******************************
# Prepare data
fit <- lm(mdl, data = dat)
rdt <- rstudent(fit)
pdt <- fitted(fit)
# Set margins
par(mar = c(4, 4, 0, 0) + 0.1)
# Get scales
xscl <- get_scale(pdt, .05)
yscl <- get_zero_scale(rdt, .1)
# Generate plot
plot(pdt, rdt,
main = "",
xlab = "",
ylab = "",
pch = 1, # open circles
cex = 1.3,
col = "#05379B",
xlim = xscl,
ylim = yscl,
axes = FALSE
)
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# Generate reference lines
abline(h = 2, col = "grey60", lwd = 1)
abline(h = -2, col = "grey60", lwd = 1)
# Add bottom label
mtext("Predicted Value", side = 1, line = par("mar")[1] - 2, cex = .9)
mtext("RStudent", side = 2, line = par("mar")[2] - 2, cex = .9)
# Frame
box(col = "grey70", lwd = 1)
#******************************
# 3) RStudent By Leverage
#******************************
# Prepare data
fit <- lm(mdl, data = dat)
rdt <- rstudent(fit)
ldt <- hatvalues(fit)
n <- length(ldt)
p <- length(coef(fit)) # includes intercept
# PROC REG-style cutoff
cut <- 2 * p / n
# Set margins
par(mar = c(4, 4, 0, 0) + 0.1)
# Get scales
xscl <- c(min(ldt) * .9, max(ldt, cut) * 1.10)
yscl <- get_zero_scale(rdt, .1)
# Define colors
nc <- "#05379B" # Normal color
# Generate plot
plot(ldt, rdt,
main = "",
xlab = "",
ylab = "",
pch = 1, # open circles
cex = 1.3,
col = nc,
xlim = xscl,
ylim = yscl,
axes = FALSE
)
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# Generate reference lines
abline(h = 2, col = "grey60", lwd = 1)
abline(h = -2, col = "grey60", lwd = 1)
abline(v = cut, col = "grey60", lwd = 1) # Cutoff
# Add labels
mtext("Leverage", side = 1, line = par("mar")[1] - 2, cex = .9)
mtext("RStudent", side = 2, line = par("mar")[2] - 2, cex = .9)
# Frame
box(col = "grey70", lwd = 1)
# Add labels
if (plt$label) {
# Identify points that need labels
lblids <- ifelse((rdt > 2 | rdt < -2) | ldt > cut, TRUE, FALSE)
if (!is.null(plt$id)) {
if (!plt$id %in% names(dat)) {
stop("Label ID variable '", plt$id, "' not found in input dataset.")
}
lbls <- dat[[plt$id]][lblids]
} else {
lbls <- seq(1, length(rdt))[lblids]
}
# Filter data
lrdt <- rdt[lblids]
lldt <- ldt[lblids]
# Assign labels
if (length(lbls) > 0) {
text(lldt, lrdt, labels = lbls, pos = 3)
}
}
#******************************
# 4) QQPlot
#******************************
# Prepare data
rdt <- res$OutputStatistics$RESID
# Set margins
par(mar = c(4, 4, 0, 0) + 0.1)
# Prepare Y data
ydat <- sort(rdt)
n <- length(ydat)
# SAS (Blom) plotting position formula: (i - 0.375) / (n + 0.25)
i <- 1:n
p_sas <- (i - 0.375) / (n + 0.25)
# Convert to theoretical normal quantiles
xdat <- qnorm(p_sas)
# Plot manually
plot(xdat, ydat,
main = "",
xlab = "",
ylab = "",
pch = 1,
col = "#05379B",
cex = 1.3,
axes = FALSE)
# Slope = Standard Deviation, Intercept = Mean
abline(a = mean(ydat), b = sd(ydat), col = "grey60")
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# X Axis label
mtext("Quantile", side = 1, line = par("mar")[1] - 2, cex = .9)
mtext("Residual", side = 2, line = par("mar")[2] - 2, cex = .9)
# Frame
box(col = "grey70", lwd = 1)
#******************************
# 5) Observed by Predicted
#******************************
# Prepare data
odt <- datclean[[dvr]]
pdt <- res$OutputStatistics$PREVAL
# Set margins
par(mar = c(4, 4, 0, 0) + 0.1)
# Get scales
xscl <- get_scale(pdt, 0)
yscl <- get_scale(odt, 0)
scl <- c(min(xscl[1], yscl[1]), max(xscl[2], yscl[2]))
# Generate plot
plot(pdt, odt,
main = "",
xlab = "",
ylab = "",
pch = 1, # open circles
cex = 1.3,
col = "#05379B",
xlim = scl,
ylim = scl,
axes = FALSE
)
# Add custom axes
axis(side = 1, col.ticks = "grey55", tck = -0.015, mgp = c(3, .5, 0))
axis(side = 2, las = 1, col.ticks = "grey55", tck = -0.015, mgp = c(3, .5, 0))
# Generate diagonal line
usr <- par("usr")
segments(usr[1], usr[3], usr[2], usr[4], col = "grey60")
# Add labels
mtext("Predicted Value", side = 1, line = par("mar")[1] - 2, cex = .9)
mtext(dvr, side = 2, line = par("mar")[1] - 2, cex = .9)
# Frame
box(col = "grey70", lwd = 1)
#******************************
# 6) Cook's D
#******************************
# Set margins
par(mar = c(4, 5, 0, 0) + 0.1)
# Calculate statistics
fit <- lm(mdl, data = dat)
cd <- cooks.distance(fit)
xcd <- seq_along(cd)
n <- length(cd)
p <- length(coef(fit)) # includes intercept
## Calculate cutoff
cutoff <- 4 / n
plot(xcd, cd,
type = "h", # vertical spikes
lwd = 1,
col = "#05379B",
xlab = "",
ylab = "",
main = "",
ylim = c(0, max(cd, cutoff) * 1.10),
axes = FALSE)
## Open circles at the top of each spike
points(seq_along(cd), cd,
pch = 1,
col = "#05379B", cex = 1.3)
## Reference line (PROC REG default)
abline(h = cutoff, col = "grey60", lwd = 1)
# Add custom axes
axis(side = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Generate zero line
abline(h = 0, col = "#05379B", lwd = 1)
# Add labels
mtext("Observation", side = 1, line = par("mar")[1] - 2, cex = .9)
mtext("Cook's D", side = 2, line = par("mar")[1] - 1, cex = .9)
# Frame
box(col = "grey70", lwd = 1)
# Add labels
if (plt$label) {
# Identify points that need labels
lblids <- ifelse(cd > cutoff, TRUE, FALSE)
if (!is.null(plt$id)) {
if (!plt$id %in% names(dat)) {
stop("Label ID variable '", plt$id, "' not found in input dataset.")
}
lbls <- dat[[plt$id]][lblids]
} else {
lbls <- seq(1, length(cd))[lblids]
}
# Filter data
lcd <- cd[lblids]
lxcd <- xcd[lblids]
# Filter colors
# lcol <- cols[cols != nc]
# Assign labels
if (length(lbls) > 0) {
text(lxcd, lcd, labels = lbls, pos = 3)
}
}
#******************************
# 7) Residual Histogram
#******************************
# Prepare data
rdt <- res$OutputStatistics$RESID
# Set margins
par(mar = c(4, 4, 0, 0) + 0.1)
# Calculate stats
n <- length(rdt)
mu <- mean(rdt)
sdx <- sd(rdt)
# Calculate breaks
brks <- pretty(range(rdt), n = nclass.Sturges(rdt),
min.n = 1, high.u.bias = 3)
# Histogram (Percent scale)
h <- hist(rdt,
breaks = brks,
plot = FALSE)
# Convert counts to percent
h$counts <- h$counts / sum(h$counts) * 100
# Use calculated breaks to get scale
scl <- get_reg_xlim(h$breaks, mu, sdx, pad_frac = 0, max_sigma = 4)
# Use calculated breaks to get scale
w <- diff(h$breaks)[1] # bin width
xmin <- min(h$breaks) - w
xmax <- max(h$breaks) + w
scl1 <- c(xmin, xmax)
# Use calculated break to get normal curve
grid <- seq(scl[1] , scl[2], length.out = 300)
y_norm_percent <- 100 * w * dnorm(grid, mean = mu, sd = sdx)
y_mx <- max(max(y_norm_percent), max(h$counts))
# Create plot using histogram values
plot(h,
col = "#CAD5E5", #"grey85",
border = "grey20",
main = "",
xlab = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.05),
axes = FALSE)
# Normal curve overlay (scaled to percent)
lines(grid, y_norm_percent, col = "steelblue4", lwd = 2)
# Add x axis label
mtext("Residual", side = 1, line = par("mar")[1] - 2, cex = .9)
mtext("Percent", side = 2, line = par("mar")[1] - 2, cex = .9)
# Calculate x axis ticks
df <- abs(h$mids[1] - h$mids[2])
te <- (df * seq(1, 10)) + h$mids[length(h$mids)]
be <- (-(df * seq(10, 1))) + h$mids[1]
df <- c(be, h$mids, te)
xtks <- df[df > scl[1] & df < scl[2]]
# Add custom axes
axis(side = 1, col.ticks = "grey55", at = xtks, labels = TRUE,
mgp = c(3, .5, 0), tck = -0.015, cex.axis = .75)
axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Frame
box(col = "grey70", lwd = 1)
#******************************
# 8) RF Plot
#******************************
# Prepare data
rdt <- res$OutputStatistics$RESID
fdt <- res$OutputStatistics$PREVAL
# dt <- dat[[ivr]]
# Get length of predicted values
n <- length(fdt)
## Proportion less (SAS-style)
p <- (seq_len(n) - 0.5) / n
## Left panel: Fit Mean
fit_centered <- sort(fdt - mean(fdt))
## Right panel: Residuals
res_sorted <- sort(rdt)
# Get y scale
smin <- min(fit_centered)
if (min(res_sorted) < smin) {
smin <- min(res_sorted)
}
smax <- max(fit_centered)
if (min(res_sorted) > smax) {
smax <- max(res_sorted)
}
# Finalize y scale
yscl <- c(smin, smax) * 1.05
# Set margins
par(mar = c(4, 4, 1.5, 0) + 0.1,
fig = c(.334, .53, 0, .334), new = TRUE) # Can place in specific spot
## Left: Fit Mean
plot(p, fit_centered,
pch = 1,
col = "#05379B",
xlab = "",
ylab = "",
ylim = yscl,
main = "Fit-Mean",
font.main = 1,
cex.main = 1.4,
axes = FALSE)
box(col = "grey70", lwd = 1)
# Add custom axes
axis(side = 1, col.ticks = "grey55", col = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", las = 1, col = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
## Plot layout
par(mar = c(4, .2, 1.5, 0) + 0.1,
fig = c(.53, .667, 0, .334), new = TRUE)
## Right: Residual
plot(p, res_sorted,
pch = 1,
col = "#05379B",
xlab = "",
ylab = "",
ylim = yscl,
main = "Residual",
font.main = 1,
cex.main = 1.4,
axes = FALSE)
box(col = "grey70", lwd = 1)
# Add custom axes
axis(side = 1, col.ticks = "grey55", col = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
## Add Labels
mtext(" Proportion Less",
side = 1, outer = TRUE,
line = par("oma")[1] - 2, cex = .9)
#******************************
# 9) Statistics or Boxplot
#******************************
if (!"none" %in% plt$stats) {
# Set margins and figure placement
par(mar = c(4, 5, 0, 2) + 0.1,
fig = c(.667, 1, 0, .334), new = TRUE)
# Initialize empty plot
# to create drawing area.
plot.new()
# Get needed statistics as named vector
stats <- collect_stats(res, plt)
if ("aic" %in% plt$stats) {
aicv <- round(extractAIC(fit)[2], 4)
stats <- append(stats, aicv)
names(stats)[length(stats)] <- "AIC"
}
# Set parameters
box_col = "white"
border_col = "grey70"
pad_y = 0.03 # vertical padding as fraction of y-range
pad_x = 0.03
# Get environment settings
usr <- par("usr") # c(x1, x2, y1, y2)
x_rng <- diff(usr[1:2])
y_rng <- diff(usr[3:4])
# x positions (user sets; relative to right edge)
x_left <- 0
x_right <- 1
# Determine line spacing
line_h <- strheight("Mg") # approx line height in user units
gap_h <- 0.8 * line_h # gap between lines
inner_h <- length(stats) * line_h + (length(stats) - 1) * gap_h
# box padding in user units
py <- pad_y * y_rng
px <- pad_x * x_rng
box_h <- inner_h + 2 * py
# center vertically in plot area
y_mid <- mean(usr[3:4])
y_bot <- y_mid - box_h / 2
y_top <- y_mid + box_h / 2
# draw box
rect(x_left, y_bot - (py/2), x_right, y_top, col = box_col, border = border_col)
# y positions for each row (top to bottom)
y_start <- y_top - py - line_h / 2
y_pos <- y_start - (seq_along(stats) - 1) * (line_h + gap_h)
# left column: names
text(x_left + px, y_pos,
labels = names(stats),
adj = c(0, 0.5), cex= 1.25)
# right column: values
text(x_right - px, y_pos,
labels = unname(stats),
adj = c(1, 0.5), cex= 1.25)
} else { # Boxplot
# Prepare data
rdt <- res$OutputStatistics$RESID
# Set margins
par(mar = c(4, 5, 0, 0) + 0.1,
fig = c(.667, 1, 0, .334), new = TRUE)
# Get scales
yscl <- get_zero_scale(rdt, .05)
# Calculate SAS compatible quantile
q <- quantile(rdt, probs = c(0, .25, .5, .75, 1), type = 2)
# Prepare to pass data
bp <- list(
stats = matrix(q, ncol = 1),
n = length(rdt),
conf = NULL,
out = numeric(0)
)
# Box plot
bxp(bp,
main = "",
ylab = "",
boxfill = "#CAD5E5",
border = "grey40",
outline = FALSE,
lty = 1,
ylim = yscl,
whiskcol = "#05379B",
staplecol = "#05379B",
medcol = "#05379B",
medlwd = 1,
boxwex = .65,
axes = FALSE)
# Mean diamond
points(1, mean(rdt),
pch = 5, # diamond
col = "#05379B",
cex = 1.4)
# Add custom axis
aval <- axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Label
mtext("Residual", side = 2, line = par("mar")[1] - 1, cex = .9)
# Frame
box(col = "grey70", lwd = 1)
}
#******************************
# Finish up
#******************************
# Set title
mtext(paste0("Fit Diagnostics for ", dvr), side = 3, cex = 1.2,
outer = TRUE, font = 2, line = par("oma")[3] - 2.25)
# Frame
box("outer", col = "grey70", lwd = 1)
# Restore margins
par(mar = op,
oma = om,
mfrow = c(1, 1),
fig = ofig,
new = FALSE)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
# Individual Plots --------------------------------------------------------
# What if there are multiple independent variables?
# Need to show separate plot for each independent variable, all in
# same panel.
#' @noRd
render_residuals <- function(dat, res, mdl) {
op <- par("mar")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Max per panel
pnlmx <- 6
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Set titles
ttl <- "Residuals for "
if (length(ivr) > 1) {
ttl <- "Residual by Regressors for "
}
# Determine number of panels and panel layout
if (length(ivr) == 1) {
pnls <- 1
rws <- c(1, 1)
} else if (length(ivr) == 2) {
pnls <- 1
rws <- c(1, 2)
hti <- 3.5
} else if (length(ivr) %in% c(3, 4)) {
pnls <- 1
rws <- c(2, 2)
} else if (length(ivr) %in% c(5, 6)) {
pnls <- 1
rws <- c(2, 3)
} else if (length(ivr) > pnlmx) {
pnls <- ceiling(length(ivr) / pnlmx)
rws <- c(2, 3)
}
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Set up object to return
ret <- list()
# Get residuals
rdt <- res$OutputStatistics$RESID
# Remove rows with NAs
dat <- remove_na_rows(dat, dvr, ivr)
# Split independent variables by panel
if (length(ivr) > pnlmx) {
fct <- ceiling(seq_along(ivr) / pnlmx)
ivrs <- split(ivr, fct)
} else {
ivrs <- list(ivr)
}
for (pnl in seq(1, pnls)) {
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Set up margins and grid
par(oma = c(0, 4, 2.5, .75) + 0.1,
mar = c(4, 1, 0, 0),
mfrow = rws)
# Flag for first chart
firstChart <- TRUE
# Chart counter
chtcnt <- 0
for (iv in ivrs[[pnl]]) {
# Increment counter
chtcnt <- chtcnt + 1
# Get data for independent variable
dt <- dat[[iv]]
# Get y scale
mx <- max(abs(rdt)) * 1.05
scl <- c(-mx, mx)
# Generate plot
plot(dt, rdt,
main = "",
xlab = "",
ylab = "",
pch = 1, # open circles
cex = 1.3,
col = "#05379B",
ylim = scl,
axes = FALSE
)
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
if (firstChart) {
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
}
# Put x axis label
mtext(iv, side = 1, line = par("mar")[1] - 1.5)
# Generate zero line
abline(h = 0, col = "grey60", lwd = 1)
# Frame
box(col = "grey70", lwd = 1)
# Set flag
if (chtcnt %% par("mfrow")[2] == 0) {
firstChart <- TRUE
} else {
firstChart <- FALSE
}
}
# Title
mtext(paste0(ttl, dvr), side = 3,
font = 2, cex = 1.2, outer = TRUE, line = par("oma")[3] - 2)
# Y Axis lable
mtext("Residual", side = 2, outer = TRUE, line = par("oma")[2] - 1.5)
# Draw line around entire chart
box("outer", col = "grey70", lwd = 1)
# Close device context
dev.off()
# Put plot in reporter plot object
ret[[length(ret) + 1]] <- create_plot(pth, height = hti, width = wdi)
}
# Restore margins
par(mar = op, mfrow = c(1, 1))
if (length(ret) == 1) {
ret <- ret[[1]]
}
return(ret)
}
render_residualbypredicted <- function(dat, res, mdl) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Prepare data
rdt <- res$OutputStatistics$RESID
pdt <- res$OutputStatistics$PREVAL
# Set margins
par(mar = c(4, 5, 2, .75) + 0.1)
# Get scales
xscl <- get_scale(pdt, .05)
yscl <- get_zero_scale(rdt, .05)
# Generate plot
plot(pdt, rdt,
main = paste0("Residual by Predicted for ", dvr),
xlab = "",
ylab = "Residual",
pch = 1, # open circles
cex = 1.3,
col = "#05379B",
xlim = xscl,
ylim = yscl,
axes = FALSE
)
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# Generate zero line
abline(h = 0, col = "grey60", lwd = 1)
# Add bottom label
mtext("Predicted Value", side = 1, line = par("mar")[1] - 2)
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
# Restore margins
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
render_rstudentbypredicted <- function(dat, res, mdl, plt) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Prepare data
fit <- lm(mdl, data = dat)
rdt <- rstudent(fit)
pdt <- fitted(fit)
# Set margins
par(mar = c(4, 5, 2, .75) + 0.1)
# Get scales
xscl <- get_scale(pdt, .05)
yscl <- get_zero_scale(rdt, .1)
# Generate plot
plot(pdt, rdt,
main = paste0("RStudent by Predicted for ", dvr),
xlab = "",
ylab = "RStudent",
pch = 1, # open circles
cex = 1.3,
col = "#05379B",
xlim = xscl,
ylim = yscl,
axes = FALSE
)
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# Generate reference lines
abline(h = 2, col = "grey60", lwd = 1)
abline(h = -2, col = "grey60", lwd = 1)
# Add bottom label
mtext("Predicted Value", side = 1, line = par("mar")[1] - 2)
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
# Add labels
if (plt$label) {
# Identify points that need labels
lblids <- ifelse(rdt > 2 | rdt < -2, TRUE, FALSE)
if (!is.null(plt$id)) {
if (!plt$id %in% names(dat)) {
stop("Label ID variable '", plt$id, "' not found in input dataset.")
}
lbls <- dat[[plt$id]][lblids]
} else {
lbls <- seq(1, length(rdt))[lblids]
}
# Filter data
lrdt <- rdt[lblids]
lpdt <- pdt[lblids]
# Filter colors
# lcol <- cols[cols != nc]
# Assign labels
if (length(lrdt) > 0) {
text(lpdt, lrdt, labels = lbls, pos = 3)
}
}
# Restore margins
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
render_rstudentbyleverage <- function(dat, res, mdl, plt) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Prepare data
fit <- lm(mdl, data = dat)
rdt <- rstudent(fit)
ldt <- hatvalues(fit)
n <- length(ldt)
p <- length(coef(fit)) # includes intercept
# PROC REG-style cutoff
cut <- 2 * p / n
# Set margins
par(mar = c(6, 5, 2, .75) + 0.1)
# Get scales
xscl <- c(min(ldt) * .9, max(ldt, cut) * 1.10)
yscl <- get_zero_scale(rdt, .1)
# Define colors
nc <- "#05379B" # Normal color
lvc <- "#01665E" # Leverage color
otc <- "#A23A2E" # Outlier color
olc <- "#543005" # Outlier and leverage
# Assign colors
cols <- rep(nc, length(ldt))
cols <- ifelse(ldt > cut, lvc, cols)
cols <- ifelse(rdt > 2 | rdt < -2, otc, cols)
cols <- ifelse((rdt > 2 | rdt < -2) & ldt > cut, olc, cols)
# Generate plot
plot(ldt, rdt,
main = paste0("Outlier and Leverage Diagnostics for ", dvr),
xlab = "",
ylab = "RStudent",
pch = 1, # open circles
cex = 1.3,
col = cols,
xlim = xscl,
ylim = yscl,
axes = FALSE
)
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# Generate reference lines
abline(h = 2, col = "grey60", lwd = 1)
abline(h = -2, col = "grey60", lwd = 1)
abline(v = cut, col = "grey60", lwd = 1) # Cutoff
# Add bottom label
mtext("Leverage", side = 1, line = par("mar")[1] - 4)
# Add labels
if (plt$label) {
# Identify points that need labels
lblids <- ifelse((rdt > 2 | rdt < -2) | ldt > cut, TRUE, FALSE)
if (!is.null(plt$id)) {
if (!plt$id %in% names(dat)) {
stop("Label ID variable '", plt$id, "' not found in input dataset.")
}
lbls <- dat[[plt$id]][lblids]
} else {
lbls <- seq(1, length(rdt))[lblids]
}
# Filter data
lrdt <- rdt[lblids]
lldt <- ldt[lblids]
# Filter colors
lcol <- cols[cols != nc]
# Assign labels
if (length(lbls) > 0) {
text(lldt, lrdt, labels = lbls, pos = 3, col = lcol)
}
}
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
# Legend (bottom, SAS-style)
legend("bottom",
legend = c("Outlier ",
"Leverage",
"Outlier and Leverage"),
pch = c(1, 1, 1),
pt.cex = 1.3,
col = c(otc,
lvc,
olc),
horiz = TRUE,
bty = "o",
box.col = "grey", # Grey border to match SAS
x.intersp = 1, # Spacing between group label and box
y.intersp = 0, # Spacing between content and borders
text.width = NA, # Compute label widths dynamically
inset = c(0, -0.23),
xpd = TRUE)
# Restore margins
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
# Kill if more than 1 independent variable.
# Kill stats if "none" parameter is passed.
#' @noRd
render_fitplot <- function(dat, res, mdl, plt, alph) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
if (length(ivr) == 1) {
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Set Confidence limits
alph1 <- 1 - alph
alph2 <- alph1 * 100
# Fit model
fit <- lm(mdl, data = dat)
# Prediction grid
xnew <- seq(min(dat[[ivr]], na.rm = TRUE), max(dat[[ivr]], na.rm = TRUE),
length.out = 200)
# Create data frame
df1 <- data.frame(xnew)
names(df1) <- ivr
conf <- predict(fit,
newdata = df1,
interval = "confidence", level = alph1)
pred <- predict(fit,
newdata = df1,
interval = "prediction", level = alph1)
# Set margins
if ("none" %in% plt$stats) {
par(mar = c(5, 4.5, 2, .75) + 0.1)
} else {
par(mar = c(5, 4.5, 2, 8.5) + 0.1)
}
# Plot
plot(dat[[ivr]], dat[[dvr]],
main = paste0("Fit Plot for ", dvr),
xlab = "",
ylab = dvr,
pch = 1,
cex = 1.23,
col = "#05379B",
ylim = range(pred),
axes = FALSE)
# X Axis label
mtext(ivr, side = 1, line = par("mar")[1] - 3.25)
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# 95% confidence band (shaded)
polygon(c(xnew, rev(xnew)),
c(conf[, "lwr"], rev(conf[, "upr"])),
col = adjustcolor("steelblue", alpha.f = 0.30),
border = NA)
# Fitted line
lines(xnew, conf[, "fit"],
col = "steelblue4",
lwd = 2)
# 95% prediction limits (dashed)
lines(xnew, pred[, "lwr"],
col = "steelblue3",
lty = 2)
lines(xnew, pred[, "upr"],
col = "steelblue3",
lty = 2)
# Legend (bottom, SAS-style)
legend("bottom",
legend = c("Fit ",
paste0(alph2, "% Confidence Limits "),
paste0(alph2, "% Prediction Limits")),
lwd = c(2, NA, 1),
lty = c(1, NA, 2),
pch = c(NA, 15, NA),
pt.cex = 2,
col = c("steelblue4",
adjustcolor("steelblue", 0.30),
"steelblue3"),
horiz = TRUE,
bty = "o",
box.col = "grey", # Grey border to match SAS
x.intersp = .7, # Spacing between group label and box
y.intersp = 0, # Spacing between content and borders
text.width = NA, # Compute label widths dynamically
inset = c(0, -0.19),
xpd = TRUE)
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
if (!"none" %in% plt$stats) {
# Gather statistics for box
sts <- collect_stats(res, plt)
if ("aic" %in% plt$stats) {
aicv <- round(extractAIC(fit)[2], 4)
sts <- append(sts, aicv)
names(sts)[length(sts)] <- "AIC"
}
# Create box
draw_stats_box(sts, x_frac = c(0.015, 0.3))
}
# Restore margins
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
} else {
ret <- NULL
}
return(ret)
}
#' @noRd
render_residualhistogram <- function(dat, res, mdl, plt) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Prepare data
rdt <- res$OutputStatistics$RESID
# dt <- dat[[ivr]]
# Set margins
par(mar = c(4, 5, 2, .75) + 0.1)
# Calculate stats
n <- length(rdt)
mu <- mean(rdt)
sdx <- sd(rdt)
# Calculate breaks - Closer to SAS algorithm
# brks <- pretty(range(rdt), n = nclass.Sturges(rdt),
# min.n = 1, high.u.bias = 3)
# Calculate breaks - Closer to SAS algorithm
brks <- pretty(range(rdt), n = nclass.Sturges(rdt),
min.n = 1, high.u.bias = 3)
# Histogram (Percent scale)
h <- hist(rdt,
breaks = brks,
plot = FALSE)
# Convert counts to percent
h$counts <- h$counts / sum(h$counts) * 100
# Use calculated breaks to get normal scale
w <- diff(h$breaks)[1] # bin width
xmin <- min(h$breaks) - w
xmax <- max(h$breaks) + w
scl1 <- c(xmin, xmax)
# Use calculated breaks to get scale
scl <- get_reg_xlim(h$breaks, mu, sdx, pad_frac = 0, max_sigma = 4)
# Use calculated break to get normal curve
grid <- seq(scl1[1] , scl1[2], length.out = 300)
y_norm_percent <- 100 * w * dnorm(grid, mean = mu, sd = sdx)
y_mx <- max(max(y_norm_percent), max(h$counts))
# Create plot using histogram values
plot(h,
col = "#CAD5E5", #"grey85",
border = "grey20",
main = paste0("Distribution of Residuals for ", dvr),
xlab = "",
ylab = "Percent",
xlim = scl,
ylim = c(0, y_mx * 1.05), # max(h$counts) * 1.05),
axes = FALSE)
# Normal curve overlay (scaled to percent)
# x <- seq(min(rdt) * 2, max(rdt) * 2, length.out = 300)
#
# y_norm <- dnorm(x, mean = mu, sd = sdx)
# y_norm <- y_norm / max(y_norm) * max(h$counts)
lines(grid, y_norm_percent, col = "steelblue4", lwd = 2)
# Kernel density overlay (scaled to percent)
# dens <- density(rdt)
#
# y_kern <- dens$y / max(dens$y) * max(h$counts)
# lines(dens$x, y_kern, col = "orangered2", lwd = 2)
# Bandwidth (Silverman / SAS)
bw <- 1.06 * sdx * n^(-1/5)
# Kernel density estimate
dens <- density(rdt, kernel = "gaussian", bw = bw, n = 512)
# Scale KDE to Percent axis
y_kernel_percent <- 100 * w * dens$y
yscl <- seq(50, length(y_kernel_percent) - 25)
# Overlay on histogram
lines(dens$x[yscl], y_kernel_percent[yscl], col = "orangered2", lwd = 2)
# Add x axis label
mtext("Residual", side = 1, line = par("mar")[1] - 2)
# Calculate x axis ticks
df <- abs(h$mids[1] - h$mids[2])
te <- (df * seq(1, 10)) + h$mids[length(h$mids)]
be <- (-(df * seq(10, 1))) + h$mids[1]
df <- c(be, h$mids, te)
xtks <- df[df > scl[1] & df < scl[2]]
# Add custom axes
axis(side = 1, col.ticks = "grey55", at = xtks, labels = TRUE,
mgp = c(3, .5, 0), tck = -0.015, cex.axis = .75)
axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
# Legend
mpos <- legend("topright",
legend = c("Normal", "Kernel"),
col = c("steelblue4", "orangered2"),
lwd = 2,
box.col = "grey80",
bty = "n",
cex = .9,
inset = c(.01, .01),
x.intersp = .5,
y.intersp = c(0.5, 1.8)
)
rect(mpos$rect$left, mpos$rect$top - (mpos$rect$h - .5),
mpos$rect$left + mpos$rect$w, mpos$rect$top,
border = "grey80"
)
# Restore margins
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
#' @noRd
render_qqplot <- function(dat, res, mdl) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Prepare data
rdt <- res$OutputStatistics$RESID
# Set margins
par(mar = c(4, 5, 2, .75) + 0.1)
# Prepare Y data
ydat <- sort(rdt)
n <- length(ydat)
# SAS (Blom) plotting position formula: (i - 0.375) / (n + 0.25)
i <- 1:n
p_sas <- (i - 0.375) / (n + 0.25)
# Convert to theoretical normal quantiles
xdat <- qnorm(p_sas)
# Plot manually
plot(xdat, ydat,
main = paste0("Q-Q Plot of Residuals for ", dvr),
xlab = "",
ylab = "Residual",
pch = 1,
col = "#05379B",
cex = 1.3,
axes = FALSE)
# Slope = Standard Deviation, Intercept = Mean
abline(a = mean(ydat), b = sd(ydat), col = "grey60")
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# X Axis label
mtext("Quantile", side = 1, line = par("mar")[1] - 2)
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
render_rfplot <- function(dat, res, mdl) {
# res <- proc_reg(cls,
# model = "Weight = Height",
# output = report,
# plots = regplot(type = "qqplot"),
# stats = p)
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Prepare data
rdt <- res$OutputStatistics$RESID
fdt <- res$OutputStatistics$PREVAL
# dt <- dat[[ivr]]
# Set margins
par(mar = c(5, 5, 2, .75) + 0.1)
## Proportion less (SAS-style)
n <- length(fdt)
p <- (seq_len(n) - 0.5) / n
## Left panel: Fit Mean
fit_centered <- sort(fdt - mean(fdt, na.rm = TRUE))
## Right panel: Residuals
res_sorted <- sort(rdt)
# Get y scale
smin <- min(fit_centered)
if (min(res_sorted, na.rm = TRUE) < smin) {
smin <- min(res_sorted, na.rm = TRUE)
}
smax <- max(fit_centered, na.rm = TRUE)
if (min(res_sorted, na.rm = TRUE) > smax) {
smax <- max(res_sorted, na.rm = TRUE)
}
# Finalize y scale
yscl <- c(smin, smax) * 1.05
## Plot layout
par(mfrow = c(1, 2),
mar = c(2, 0, 2, .75) + 0.1,
oma = c(2, 4, 2, .25))
## Left: Fit Mean
plot(p, fit_centered,
pch = 1,
col = "#05379B",
xlab = "",
ylab = "",
ylim = yscl,
main = "Fit-Mean",
font.main = 1,
cex.main = 1,
axes = FALSE)
box(col = "grey70", lwd = 1)
# Add custom axes
axis(side = 1, col.ticks = "grey55", col = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", las = 1, col = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
## Right: Residual
plot(p, res_sorted,
pch = 1,
col = "#05379B",
xlab = "",
ylab = "",
ylim = yscl,
main = "Residual",
font.main = 1,
cex.main = 1,
axes = FALSE)
box(col = "grey70", lwd = 1)
# Add custom axes
axis(side = 1, col.ticks = "grey55", col = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
## Overall title
mtext(paste0("Residual-Fit Spread Plot for ", dvr),
outer = TRUE, cex = 1.25, font = 2)
# X Axis Label
mtext("Proportion Less", side = 1, outer = TRUE, cex = 1.1)
# Frame
box("outer", col = "grey70", lwd = 1)
# Restore margins
par(mar = op, mfrow = c(1, 1), oma = c(0, 0, 0, 0))
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
}
#' @noRd
render_cooksd <- function(dat, res, mdl, plt) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Set margins
par(mar = c(4, 5, 2, .75) + 0.1)
# Calculate statistics
fit <- lm(mdl, data = dat)
cd <- cooks.distance(fit)
xcd <- seq_along(cd)
n <- length(cd)
p <- length(coef(fit)) # includes intercept
## Calculate cutoff
cutoff <- 4 / n
plot(xcd, cd,
type = "h", # vertical spikes
lwd = 1,
col = "#05379B",
xlab = "",
ylab = "Cook's D",
main = paste0("Cook's D for ", dvr),
ylim = c(0, max(cd, cutoff) * 1.10),
axes = FALSE)
## Open circles at the top of each spike
points(seq_along(cd), cd,
pch = 1,
col = "#05379B", cex = 1.25)
## Reference line (PROC REG default)
abline(h = cutoff, col = "grey60", lwd = 1)
# Add custom axes
axis(side = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Generate zero line
abline(h = 0, col = "#05379B", lwd = 1)
# X axis label
mtext("Observation", side = 1, line = par("mar")[1] - 2)
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
# Add labels
if (plt$label) {
# Identify points that need labels
lblids <- ifelse(cd > cutoff, TRUE, FALSE)
if (!is.null(plt$id)) {
if (!plt$id %in% names(dat)) {
stop("Label ID variable '", plt$id, "' not found in input dataset.")
}
lbls <- dat[[plt$id]][lblids]
} else {
lbls <- seq(1, length(cd))[lblids]
}
# Filter data
lcd <- cd[lblids]
lxcd <- xcd[lblids]
# Assign labels
if (length(lbls) > 0) {
text(lxcd, lcd, labels = lbls, pos = 3)
}
}
# Restore margins
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
render_observedbypredicted <- function(dat, res, mdl, plt) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 4.5 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Remove rows with NAs
dat <- remove_na_rows(dat, dvr, ivr)
# Prepare data
odt <- dat[[dvr]]
pdt <- res$OutputStatistics$PREVAL
# Set margins
par(mar = c(4, 5, 2, 1.05) + 0.1)
# Get scales
# Appears to take min of predicted and observed
# and max of predicted and observed, and uses
# those for both x and y scales. Scales have to be the same
# or the dots won't line up properly. Much trial and error
# to figure this out.
xscl <- get_scale(pdt, 0)
yscl <- get_scale(odt, 0)
scl <- c(min(xscl[1], yscl[1]), max(xscl[2], yscl[2]))
# Generate plot
plot(pdt, odt,
main = paste0("Observed by Predicted for ", dvr),
xlab = "",
ylab = dvr,
pch = 1, # open circles
cex = 1.3,
col = "#05379B",
xlim = scl,
ylim = scl,
axes = FALSE
)
# Add custom axes
axis(side = 1, col.ticks = "grey55", tck = -0.015, mgp = c(3, .5, 0))
axis(side = 2, las = 1, col.ticks = "grey55", tck = -0.015, mgp = c(3, .5, 0))
# Generate diagonal line
# abline() doesn't work at all. Not sure why. Have to use segments().
usr <- par("usr")
segments(usr[1], usr[3], usr[2], usr[4], col = "grey60")
# Add bottom label
mtext("Predicted Value", side = 1, line = par("mar")[1] - 2)
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
# Would be nice to get outlier labels in here,
# but I don't know how they are identified.
# There are no limits on this chart.
# Restore margins
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
render_residualboxplot <- function(dat, res, mdl, plt) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Prepare data
rdt <- res$OutputStatistics$RESID
pdt <- res$OutputStatistics$PREVAL
# Set margins
par(mar = c(.75, 5, 2, .75) + 0.1)
# Get scales
yscl <- get_zero_scale(rdt, .05)
# Calculate SAS compatible quantile
q <- quantile(rdt, probs = c(0, .25, .5, .75, 1), type = 2)
# Prepare to pass data
bp <- list(
stats = matrix(q, ncol = 1),
n = length(rdt),
conf = NULL,
out = numeric(0)
)
# Box plot
bxp(bp,
main = paste0("Residuals for ", dvr),
ylab = "Residual",
boxfill = "#CAD5E5",
border = "grey40",
outline = FALSE,
lty = 1,
ylim = yscl,
whiskcol = "#05379B",
staplecol = "#05379B",
medcol = "#05379B",
medlwd = 1,
boxwex = .65,
axes = FALSE)
# Mean diamond
points(1, mean(rdt, na.rm = TRUE),
pch = 5, # diamond
col = "#05379B",
cex = 1.3)
# Add custom axis
aval <- axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
## Optional: horizontal gridlines
abline(h = aval, col = "grey90", lwd = 1)
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
# Restore margins
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
render_dffits <- function(dat, res, mdl, plt) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Prepare data
fit <- lm(mdl, data = dat)
dff <- dffits(fit)
n <- length(dff)
p <- length(coef(fit)) # includes intercept
# Set margins
par(mar = c(4, 5, 2, .75) + 0.1)
## PROC REG cutoff for DFFITS
cutoff <- 2 * sqrt(p / n)
plot(seq_along(dff), dff,
type = "h", # vertical spikes
lwd = 1,
col = "#05379B",
xlab = "",
ylab = "DFFITS",
main = paste0("Influence Diagnostics for ", dvr),
ylim = range(c(dff, -cutoff, cutoff)) * 1.1,
axes = FALSE)
## Open circles at spike tips
points(seq_along(dff), dff,
pch = 1,
cex = 1.25,
col = "#05379B")
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
## Reference lines (PROC REG)
abline(h = 0, col = "#05379B", lwd = 1)
abline(h = cutoff, col = "grey60", lwd = 1)
abline(h = -cutoff, col = "grey60", lwd = 1)
# Add bottom label
mtext("Observation", side = 1, line = par("mar")[1] - 2)
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
# Add labels
if (plt$label) {
# Identify points that need labels
lblids1 <- ifelse(dff > cutoff, TRUE, FALSE)
lblids2 <- ifelse(dff < -cutoff, TRUE, FALSE)
if (!is.null(plt$id)) {
if (!plt$id %in% names(dat)) {
stop("Label ID variable '", plt$id, "' not found in input dataset.")
}
lbls1 <- dat[[plt$id]][lblids1]
lbls2 <- dat[[plt$id]][lblids2]
} else {
lbls1 <- seq(1, length(dff))[lblids1]
lbls2 <- seq(1, length(dff))[lblids2]
}
# Filter data
ly1 <- dff[lblids1]
lx1 <- seq_along(dff)[lblids1]
ly2 <- dff[lblids2]
lx2 <- seq_along(dff)[lblids2]
# Filter colors
# lcol <- cols[cols != nc]
# Assign labels
if (length(ly1) > 0) {
text(lx1, ly1, labels = lbls1, pos = 3)
}
if (length(ly2) > 0) {
text(lx2, ly2, labels = lbls2, pos = 1)
}
}
# Restore margins
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
render_dfbetas <- function(dat, res, mdl, plt) {
# Get current values
op <- par("mar")
om <- par("oma")
# Function can return multiple plots
ret <- list()
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Prepare data
fit <- lm(mdl, data = dat)
dfb <- dfbetas(fit) # matrix: n x p
n <- nrow(dfb)
p <- ncol(dfb)
# Get names for labels
param_names <- colnames(dfb)
# Baseline
chrts <- 1
rngs <- seq(1, p)
if (plt$panel) {
if (p == 1) {
mfrw <- c(1, 1)
} else if (p == 2) {
mfrw <- c(1, 2)
} else if (p > 2 & p < 5) {
mfrw <- c(2, 2)
} else if (p > 4 & p < 7) {
mfrw <- c(2, 3)
} else if (p > 6) {
chrts <- ceiling(p / 6)
mfrw <- c(2, 3)
}
# Split into a list of ranges per chart
if (chrts > 1) {
rngs <- split(rngs, ceiling(rngs/6))
} else {
rngs <- list(rngs)
}
} else {
# One chart per variable
chrts <- p
rngs <- split(rngs, rngs)
mfrw <- c(1, 1)
}
## PROC REG cutoff for DFBETAS
cutoff <- 2 / sqrt(n)
# Get yscale - Consistent Y scale on panel plot
if (plt$label) {
yscl <- c(min(dfb, -cutoff, na.rm = TRUE), max(dfb, cutoff, na.rm = TRUE)) * 1.2 # Give more room for labels
} else {
yscl <- c(min(dfb, -cutoff, na.rm = TRUE), max(dfb, cutoff, na.rm = TRUE)) * 1.1
}
# Get xscale
xscl <- c(0 - .01, nrow(dat) * 1.05)
for (i in seq_len(chrts)) {
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Output to image file
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Set margins
if (plt$panel) {
par(mfrow = mfrw,
mar = c(2, 0, 2, .75) + 0.1,
oma = c(2, 5, 2, 0) + .1)
} else {
par(mfrow = mfrw,
mar = c(2, 0, 0, .75) + 0.1,
oma = c(2, 5, 2, 0) + .1)
}
for (j in rngs[[i]]) {
# Get data for this chart
vals <- dfb[, j]
# Non-panel plot computes custom y scale for each variable
if (plt$panel == FALSE) {
yscl <- c(min(vals, -cutoff), max(vals, cutoff)) * 1.1
}
# Create empty plot
plot(seq_len(n), vals,
type = "n",
xlim = xscl,
ylim = yscl,
axes = FALSE)
# Determine grid position
grdPos <- par("mfg")
# X axis
xtcks <- axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# Y axis
if (grdPos[2] == 1) {
ytcks <- axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
}
# Grid lines
if (plt$panel) {
abline(h = ytcks, col = "grey90", lwd = 1)
abline(v = xtcks, col = "grey90", lwd = 1)
}
par(new = TRUE)
if (plt$panel) {
ttl <- param_names[j]
} else {
ttl <- ""
}
# Create actual data plot
plot(seq_len(n), vals,
type = "h",
font.main = 1,
lwd = 1,
col = "#05379B",
xlab = "",
ylab = "",
main = ttl,
xlim = xscl,
ylim = yscl, # range(c(vals, -cutoff, cutoff)) * 1.1,
axes = FALSE)
par(new = FALSE)
## Open circles at tips
if (plt$panel == FALSE) {
points(seq_len(n), vals,
pch = 1,
cex = 1.25,
col = "#05379B")
}
## Reference lines (PROC REG)
abline(h = 0, col = "#05379B", lwd = 1)
abline(h = cutoff, col = "grey60", lwd = 1)
abline(h = -cutoff, col = "grey60", lwd = 1)
# Add labels
if (plt$label) {
## Positive values
# Identify points that need labels
lblids <- ifelse(vals > cutoff, TRUE, FALSE)
if (!is.null(plt$id)) {
if (!plt$id %in% names(dat)) {
stop("Label ID variable '", plt$id, "' not found in input dataset.")
}
lbls <- dat[[plt$id]][lblids]
} else {
lbls <- seq(1, length(vals))[lblids]
}
# Filter data
xpos <- seq_len(n)[lblids]
ypos <- vals[lblids]
# Assign labels
if (length(xpos) > 0) {
text(xpos, ypos, labels = lbls, pos = 3)
}
## Negative values
# Identify points that need labels
lblids <- ifelse(vals < -cutoff, TRUE, FALSE)
if (!is.null(plt$id)) {
if (!plt$id %in% names(dat)) {
stop("Label ID variable '", plt$id, "' not found in input dataset.")
}
lbls <- dat[[plt$id]][lblids]
} else {
lbls <- seq(1, length(vals))[lblids]
}
# Filter data
xpos <- seq_len(n)[lblids]
ypos <- vals[lblids]
# Assign labels
if (length(xpos) > 0) {
text(xpos, ypos, labels = lbls, pos = 1)
}
}
box(col = "grey70", lwd = 1)
}
## Overall title
# Set title
if (plt$panel) {
mtext(paste0("Influence Diagnostics for ", dvr), side = 3, cex = 1.2,
outer = TRUE, font = 2, line = par("oma")[3] - 2.25)
} else {
vr <- param_names[j]
mtext(paste0("Influence on ", vr, " for ", dvr), side = 3, cex = 1.2,
outer = TRUE, font = 2, line = par("oma")[3] - 1.75)
}
# Labels
mtext("Observation", side = 1, outer = TRUE, line = par("oma")[1] - 2)
mtext("DFBETAS", side = 2, outer = TRUE, line = par("oma")[2] - 2)
# Frame
box("outer", col = "grey70", lwd = 1)
# Restore margins
par(mar = op, oma = om, mfrow = c(1, 1))
# Close device context
dev.off()
# Put plot in reporter plot object
ret[[length(ret) + 1]] <- create_plot(pth, height = hti, width = wdi)
}
return(ret)
}
# Utilities ---------------------------------------------------------------
#' @noRd
collect_stats <- function(res, plt) {
# mres <- proc_reg(cls, Weight ~ Height, stats = "clb", output = report)
ret <- c()
if ("default" %in% plt$stats) {
ret <- c(
"Observations" = res$NObs$NOBS[1],
"Parameters" = nrow(res$ParameterEstimates),
"Error DF" = res$ANOVA$DF[2],
"MSE" = roundup(res$ANOVA$MEANSQ[2], 2),
"R-Square" = roundup(res$FitStatistics$RSQ[1], 4),
"Adj R-Square" = roundup(res$FitStatistics$ADJRSQ[1], 4)
)
}
# Apparently got from AI and these are wrong
# obs, edf, rsquare, parms, mse, adjrsquare,
# coefvar, depmean
# sse, bic, gmsep, pc, sp, aic, cp, jp, sbc
# From SAS error message
# Documentation is not actually that good on how to do this.
# ADJRSQ, AIC, BIC, COEFFVAR, CP, DEFAULT, DEPMEAN, EDF, GMSEP, JP,
# MSE, NOBS, NPARM, PC, RSQUARE, SBC, SP, SSE
if ("nobs" %in% plt$stats) {
ret[length(ret) + 1] <- res$NObs$NOBS[1]
names(ret)[length(ret)] <- "Observations"
}
if ("nparm" %in% plt$stats) {
ret[length(ret) + 1] <- nrow(res$ParameterEstimates)
names(ret)[length(ret)] <- "Parameters"
}
if ("edf" %in% plt$stats) {
ret[length(ret) + 1] <- res$ANOVA$DF[2]
names(ret)[length(ret)] <- "Error DF"
}
if ("mse" %in% plt$stats) {
ret[length(ret) + 1] <- roundup(res$ANOVA$MEANSQ[2], 2)
names(ret)[length(ret)] <- "MSE"
}
if ("rsquare" %in% plt$stats) {
ret[length(ret) + 1] <- roundup(res$FitStatistics$RSQ[1], 4)
names(ret)[length(ret)] <- "R-Square"
}
if ("adjrsq" %in% plt$stats) {
ret[length(ret) + 1] <- roundup(res$FitStatistics$ADJRSQ[1], 4)
names(ret)[length(ret)] <- "Adj R-Square"
}
if ("depmean" %in% plt$stats) {
ret[length(ret) + 1] <- roundup(res$FitStatistics$DEPMEAN, 2)
names(ret)[length(ret)] <- "Dep Mean"
}
if ("coeffvar" %in% plt$stats) {
ret[length(ret) + 1] <- roundup(res$FitStatistics$COEFVAR, 4)
names(ret)[length(ret)] <- "Coeff Var"
}
if ("sse" %in% plt$stats) {
ret[length(ret) + 1] <- roundup(res$ANOVA$SUMSQ[2], 1)
names(ret)[length(ret)] <- "SSE"
}
return(ret)
}
#' @noRd
draw_stats_box <- function(stats,
x_frac = c(0.03, 0.30), # from right edge of plot (in x-range fractions)
cex = 0.85,
family = "sans",
box_col = "white",
border_col = "grey70",
pad_y = 0.01, # vertical padding as fraction of y-range
pad_x = 0.01) { # horizontal padding as fraction of x-range
stopifnot(is.numeric(stats), !is.null(names(stats)))
# allow drawing in the margin
op <- par(xpd = NA)
on.exit(par(op), add = TRUE)
usr <- par("usr") # c(x1, x2, y1, y2)
x_rng <- diff(usr[1:2])
y_rng <- diff(usr[3:4])
# x positions (user sets; relative to right edge)
x_left <- usr[2] + x_frac[1] * x_rng
x_right <- usr[2] + x_frac[2] * x_rng
# compute line height in user units (based on current cex)
old_cex <- par("cex")
old_family <- par("family")
par(cex = cex, family = family)
on.exit(par(cex = old_cex, family = old_family), add = TRUE)
line_h <- strheight("Mg") # approx line height in user units
gap_h <- 0.35 * line_h # gap between lines
inner_h <- length(stats) * line_h + (length(stats) - 1) * gap_h
# box padding in user units
py <- pad_y * y_rng
px <- pad_x * x_rng
box_h <- inner_h + 2 * py
# center vertically in plot area
y_mid <- mean(usr[3:4], na.rm = TRUE)
y_bot <- y_mid - box_h / 2
y_top <- y_mid + box_h / 2
# draw box
rect(x_left, y_bot, x_right, y_top, col = box_col, border = border_col)
# y positions for each row (top to bottom)
y_start <- y_top - py - line_h / 2
y_pos <- y_start - (seq_along(stats) - 1) * (line_h + gap_h)
# left column: names
text(x_left + px, y_pos,
labels = names(stats),
adj = c(0, 0.5))
# right column: values
text(x_right - px, y_pos,
labels = unname(stats),
adj = c(1, 0.5))
}
get_scale <- function(vct, pct = .1, h0 = NULL) {
rng <- range(vct, na.rm = TRUE)
if (!is.null(h0)) {
if (h0 < rng[1]) {
rng[1] <- h0
}
if (h0 > rng[2]) {
rng[2] <- h0
}
}
if (rng[1]>0){
mn <- rng[1] * (1 - pct)
} else{
mn <- rng[1] * (1 + pct)
}
if (rng[2]>0){
mx <- rng[2] * (1 + pct)
} else{
mx <- rng[2] * (1 - pct)
}
ret <- c(mn, mx)
return(ret)
}
get_zero_scale <- function(vct, pct = .1) {
rng <- range(vct, na.rm = TRUE)
mn <- abs(rng[1] * (1 + pct))
mx <- abs(rng[2] * (1 + pct))
tmx <- max(mn, mx, na.rm = TRUE)
ret <- c(-tmx, tmx)
return(ret)
}
# This function was mostly generated by AI and I don't fully understand it.
# But it does seem to work.
get_reg_xlim <- function(brks, mu, sdx,
pad_frac = 0.05,
normal_tail_frac = 0.001, # normal at edges <= this * peak
max_sigma = 6) {
# Get range of calculated breaks
rng <- range(brks, na.rm = TRUE)
# Add adding on bin range
pad0 <- pad_frac * diff(rng)
xlim <- c(rng[1] - pad0, rng[2] + pad0)
# Force symmetry around zero (expand, never shrink)
m <- max(abs(xlim), na.rm = TRUE)
xlim <- c(-m, m)
if (!is.finite(sdx) || sdx <= 0) return(xlim)
# Calculate peak of normal curve
peak <- stats::dnorm(mu, mean = mu, sd = sdx)
target <- normal_tail_frac * peak
# For symmetric limits, the "worst" edge is the one closer to mu.
# Ensure BOTH edges are far enough into the tails.
k <- 0
while ((stats::dnorm(xlim[1], mu, sdx) > target ||
stats::dnorm(xlim[2], mu, sdx) > target) &&
k < max_sigma) {
xlim <- xlim + c(-0.5 * sdx, 0.5 * sdx) # symmetric expansion
k <- k + 0.5
}
ret <- xlim - c(-0.5 * sdx, 0.5 * sdx) # xlim[seq(2, length(xlim) - 1)]
return(ret)
}
remove_na_rows <- function(dat, dvr, ivrs) {
# Determine rows with NAs
dvna <- is.na(dat[[dvr]])
ivna <- rep(FALSE, nrow(dat))
for (var in ivrs) {
ivna <- ivna | is.na(dat[[var]])
}
cvna <- dvna | ivna
# Strip any such rows
ret <- dat[!cvna, ]
return(ret)
}
Try the procs package in your browser
Any scripts or data that you put into this service are public.
procs documentation built on Nov. 2, 2026, 5:06 p.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 remove_na_rows get_reg_xlim get_zero_scale get_scale draw_stats_box collect_stats render_dfbetas render_dffits render_residualboxplot render_observedbypredicted render_cooksd render_rfplot render_qqplot render_residualhistogram render_fitplot render_rstudentbyleverage render_rstudentbypredicted render_residualbypredicted render_residuals render_diagnostics regplot
Documented in regplot
# Regplots Main Function --------------------------------------------------
# RESIDUALBYPREDICTED: Residuals vs. Predicted values.
# RSTUDENTBYPREDICTED: Externally Studentized Residuals (RStudent) vs. Predicted values.
# RSTUDENTBYLEVERAGE: Externally Studentized Residuals vs. Leverage.
# QQ (or QQPLOT): Normal Quantile-Quantile (Q-Q) plot of residuals.
# OBSERVEDBYPREDICTED: Dependent variable (Observed) vs. Predicted values.
# COOKSD: Cook’s D statistic vs. Observation number.
# RESIDUALHISTOGRAM: Histogram of residuals.
# RFPLOT: Residual-Fit (RF) spread plot.
# ALL: Requests all plots relevant to the current analysis and options.
# FITPLOT: Produces a scatter plot of the dependent variable against the regressor, including the fitted line and confidence/prediction bands. This is only available for models with a single regressor.
# RESIDUALS or RESIDUALPLOT: Produces a panel of residual plots against each independent variable in the model.
# PARTIAL (or PARTIALPLOT): Requests partial regression residual plots for each regressor. This requires the PARTIAL option to be specified in the MODEL statement.
# CRITERIONPANEL: When using model selection (like SELECTION=STEPWISE), this displays how criteria like AIC, SBC, or R-square change across different model sizes.
# panel = TRUE/FALSE
# type = diagnostics, residualbypredicted, RStudentbypredicted, rstudentbyleverage,
# qq/qqplot, observedbypredicted, cooksd, residualhistogram, rfplot, residuals/residualplot, fitplot
# stats = default aic sbc
# all = obs, edf, rsquare, sse, coefvar, bic, gmsep, pc, sp, parms, mse, adjrsquare, depmean, aic, cp, jp, sbc
#' @title Request Regression Plots
#' @description A function to request regression plots on a call to
#' \code{\link{proc_reg}}. The function allows you to specify the type of
#' regression plots to produce. It produces a combined diagnostics panel,
#' residuals dot plot, and a fit plot by default. You may also specify
#' individual plots by passing a vector of plot names
#' on the "type" parameter.
#' @details There are many types of regression plots. The plots have different
#' uses. Some of the plots help you assess normality of the data. Other plots
#' help you assess the quality of the model. You can select the type of plots
#' you want by passing a vector of plot names on the \code{type} parameter.
#'
#' Here is a list of possible plot types and a short description of each:
#' \itemize{
#' \item{\strong{diagnostics}: A fit diagnostics panel that contains 8 different
#' types of plots and a table of statistics.
#' }
#' \item{\strong{cooksd}: Cook’s D statistic vs. Observation number.
#' }
#' \item{\strong{dfbetas}: Displays the influence of each observation
#' for each coefficient in the model.
#' }
#' \item{\strong{dffits}: Displays the influence of each observation on fitted values.
#' }
#' \item{\strong{fitplot}: Produces a scatter plot of the dependent variable
#' against the regressor, including the fitted line and confidence/prediction bands.
#' The fitplot is only available for models with a single regressor.
#' }
#' \item{\strong{observedbypredicted}: Dependent variable (Observed) vs. Predicted values.
#' }
#' \item{\strong{qqplot}: Normal Quantile-Quantile (Q-Q) plot of residuals.
#' }
#' \item{\strong{residuals}: Produces a panel of residual plots against
#' each independent variable in the model.
#' }
#' \item{\strong{residualboxplot}: A box plot of residual values.
#' }
#' \item{\strong{residualbypredicted}: Residuals vs. Predicted values.
#' }
#' \item{\strong{residualhistogram}: Histogram of residuals,
#' with a normal and kernel curve overlay.
#' }
#' \item{\strong{rfplot}: Residual-Fit (RF) spread plot.
#' }
#' \item{\strong{rstudentbyleverage}: Externally Studentized Residuals vs. Leverage.
#' }
#' \item{\strong{rstudentbypredicted}: Externally Studentized Residuals (RStudent) vs. Predicted values.
#' }
#' }
#'
#' If possible, the statistics from the tabular report are used for the plots.
#' Otherwise, additional statistics functions are called to produce the
#' needed values.
#'
#' @section Plot Statistics:
#' The diagnostics panel and the fit plot each contain a small table of
#' statistics. This table is customizable. The following keywords may be
#' used to to customize the statistics table:
#' \itemize{
#' \item{\strong{adjrsq}: Adjusted R-square.
#' }
#' \item{\strong{aic}: Akaike's information criterion.
#' }
#' \item{\strong{coeffvar}: Coefficient of variation.
#' }
#' \item{\strong{depmean}: Mean of dependent.
#' }
#' \item{\strong{default}: A set of default statistics.
#' }
#' \item{\strong{edf}: Error degrees of freedom.
#' }
#' \item{\strong{mse}: Mean squared error.
#' }
#' \item{\strong{nobs}: Number of observations used.
#' }
#' \item{\strong{nparm}: Number of parameters in the model (including the intercept).
#' }
#' \item{\strong{rsquare}: The R-square statistic.
#' }
#' \item{\strong{sse}: Error sum of squares.
#' }
#' }
#' To use these keywords, pass them as a vector to the \code{stats}
#' parameter on the \code{regplot} function.
#' @section Labeling Outliers:
#' Some types of plots can be used to identify outliers in the data. For instance,
#' the Cook's D chart is excellent for such a task.
#' When identifying outliers, it is helpful to have them labelled on the chart.
#' The labels make it possible to trace the outliers back to the source
#' data.
#'
#' To get a basic row/observation label, set the \code{labels} parameter to TRUE.
#' If there is a column in the data that can be used to identify an individual
#' record, pass that column name on the \code{id} parameter. The values
#' from that column will then be used as outlier labels on those charts that
#' support labels.
#'
#' @section Additional Information:
#' Regression plots will be displayed on interactive reports only.
#' Plots are created as jpeg files, and stored in a temp directory. Those
#' temporary files are then referenced by the interactive report to display
#' the graphic.
#'
#' If desired, you may
#' output the report objects and pass to \code{\link{proc_print}}. To do this,
#' set \code{output = report} on the call to \code{\link{proc_reg}}, and pass
#' the entire list to \code{\link{proc_print}}.
#'
#' @section Discrepancies with SAS:
#' The histogram binning algorithm in R is different from the binning algorithm
#' in SAS®. R uses a "Sturges" algorithm, which more accurately represents
#' the distribution of the data. This algorithm may produce a different number
#' of bins and different height of bars than the corresponding SAS chart.
#'
#' Plots generated by \code{proc_reg} may also have some spacing, color,
#' and shape discrepancies with the corresponding SAS charts. The information
#' conveyed is expected to be similar.
#'
#' @param type The type(s) of plot to create. Multiple types should be passed
#' as a vector of strings. Valid values are "diagnostics", "cooksd", "dfbetas",
#' "dffits", "fitplot", "observedbypredicted", "qqplot", "residuals",
#' "residualboxplot", "residualbypredicted", "residualhistogram", "rfplot",
#' "rstudentbyleverage", and "rstudentbypredicted". The default value is
#' a vector with "diagnostics", "residuals", and "fitplot". The "diagnostics"
#' keyword produces a single combined chart with 8 different plots and a
#' selection of statistics in a small table. The statistics can be controlled
#' by the \code{stats} parameter. You may also pass the "all" keyword to produce
#' all charts available for the analysis.
#' @param panel Whether or not to display the diagnostics plots combined into
#' in a single panel. Default is TRUE. A value of FALSE will create
#' individual plots instead. This parameter is equivalent to the
#' "unpack" keyword in SAS.
#' @param stats The statistics to display on the diagnostics panel or fit plot.
#' Valid values
#' are: "adjrsq", "aic", "coeffvar", "depmean", "default", "edf", "mse", "nobs",
#' "nparm", "rsquare", and "sse". The
#' default value is "default", which produces the following statistics:
#' "nobs", "nparm", "edf", "mse", "rsquare", and "adjrsq". You may also pass
#' the value "none" if you do not want any statistics shown on the chart. In the
#' case of the diagnostics panel, if the statistics table is removed, it
#' will be replaced with a residual box plot.
#' @param label Whether or not to label outlier values. Valid values
#' are TRUE or FALSE. Default is FALSE. If TRUE, this option will assign
#' labels to outlier values on some charts. Charts that support labels are
#' as follows: "diagnostics", "cooksd", "dffits", "dfbetas", "rstudentbypredicted",
#' and "rstudentbyleverage".
#' @param id If the \code{label} parameter is TRUE, this parameter determines
#' which value is assigned to the label. By default, the row number will
#' be assigned. You may also assign a column name from the input dataset
#' to this parameter to use as the label value.
#' @seealso [proc_reg()]
#' @examples
#' library(procs)
#'
#' # Turn off printing for CRAN checks
#' # Set to TRUE to run in local environment
#' options("procs.print" = FALSE)
#'
#'
#' # Example 1: Regression statistics with default plots
#' res <- proc_reg(iris, model = "Sepal.Length = Petal.Length",
#' output = report,
#' plots = TRUE,
#' titles = "Iris Regression Statistics")
#'
#' # View results
#' res
#'
#' # Example 2: Regression statistics with custom plot request and by variable
#' res <- proc_reg(iris, model = "Sepal.Length = Petal.Length",
#' output = report,
#' by = Species,
#' plots = regplot(type = v(residualhistogram,
#' rstudentbypredicted, rstudentbyleverage),
#' label = TRUE),
#' titles = "Iris Regression Statistics")
#'
#' # View results
#' res
#'
#' # Example 3: Regression statistics with multiple models, same plot string
#' res <- proc_reg(iris, model = c("Sepal.Length = Petal.Length",
#' "Sepal.Length = Sepal.Width",
#' "Sepal.Length = Petal.Width"),
#' output = report,
#' plots = "diagnostics",
#' titles = "Iris Regression Statistics")
#'
#' # View results
#' res
#'
#' # Example 4: Regression statistics with multiple models, different plot strings
#' res <- proc_reg(iris, model = c("Sepal.Length = Petal.Length",
#' "Sepal.Length = Sepal.Width",
#' "Sepal.Length = Petal.Width"),
#' output = report,
#' plots = list("diagnostics",
#' "residualhistogram",
#' "fitplot"),
#' titles = "Iris Regression Statistics")
#'
#' # View results
#' res
#'
#' # Example 5: Regression statistics with multiple models, different plot functions
#' res <- proc_reg(iris, model = c("Sepal.Length = Petal.Length",
#' "Sepal.Length = Sepal.Width",
#' "Sepal.Length = Petal.Width"),
#' output = report,
#' plots = list(regplot(type = "diagnostics"),
#' regplot(type = "cooksd",
#' label = TRUE),
#' regplot(type = "fitplot",
#' stats = c("nobs", "mse", "rsquare"))),
#' titles = "Iris Regression Statistics")
#'
#' # View results
#' res
#'
#' # Example 6: Regression statistics with multiple models, influence charts
#' res <- proc_reg(iris, model = c("Sepal.Length = Petal.Length",
#' "Sepal.Length = Petal.Length Petal.Width",
#' "Sepal.Length = Petal.Length Petal.Width Sepal.Width"),
#' output = report,
#' plots = regplot(v(cooksd, dffits, dfbetas), label = TRUE),
#' titles = "Iris Regression Statistics")
#'
#' # View results
#' res
#' @export
regplot <- function(type = c("diagnostics", "residuals", "fitplot"), panel = TRUE,
stats = "default", label = FALSE, id = NULL) {
# Non-standard evaluation
type <- resolve_arg(type)
stats <- resolve_arg(stats)
id <- resolve_arg(id)
# Parameter Checks
vldvals <- c("diagnostics", "residuals", "fitplot", 'qqplot', 'rfplot',
"residualbypredicted", "rstudentbypredicted", "rstudentbyleverage",
"cooksd", "residualhistogram", "observedbypredicted", "residualboxplot",
"all", "dffits", "dfbetas")
if (any(!type %in% vldvals)) {
ivd <- type[!type %in% vldvals]
stop(paste0("Parameter value for 'type' invalid: ", paste0("'", ivd, "'", collapse = ", "),
"\nValid values are: ",
"'all', 'diagnostics', 'residuals', 'fitplot', 'qqplot', 'rfplot', ",
"'dffits', 'dfbetas', ",
"'residualbypredicted', 'rstudentbypredicted', 'rstudentbyleverage', ",
"'cooksd', 'residualhistogram', 'observedbypredicted', 'residualboxplot'."
))
}
if (!panel %in% c(TRUE, FALSE)) {
stop("Parameter value for 'panel' invalid. Valid values are TRUE or FALSE.")
}
# sse, bic, gmsep, pc, sp, aic, cp, jp, sbc
statlst <- c("adjrsq", "aic", "coeffvar", "depmean", "default", "edf", "mse",
"nobs", "none", "nparm", "rsquare", "sse")
if (!all(stats %in% statlst)) {
stop(paste0("Parameter value for 'stats' invalid. Valid values are ",
paste0("'", statlst, "'", collapse = ", "), "."))
}
# Create object
ret <- structure(list(), class = c("regplot", "list"))
ret$stats <- tolower(stats)
ret$panel <- panel
ret$type <- tolower(type)
ret$label <- label
ret$id <- id
return(ret)
}
# Plot Rendering ----------------------------------------------------------
#' @noRd
render_regplot <- function (dat, res, mdl, plt, alph) {
ret <- NULL
if (is.character(plt)) {
if (all(plt == "regplot")) {
plt <- regplot()
} else {
plt <- regplot(type = plt)
}
}
if ("regplot" %in% class(plt)) {
typs <- plt$type
# All logic
if (all(typs == "all")) {
typs <- c("residualhistogram", "residualbypredicted", "rstudentbypredicted",
"observedbypredicted", "cooksd", "rstudentbyleverage", 'qqplot', 'rfplot',
"residualboxplot", "diagnostics", "dffits", "dfbetas", "residuals", "fitplot")
}
# "Unpack" logic
if ("diagnostics" %in% typs & plt$panel == FALSE) {
pos <- match("diagnostics", typs)
vins <- c("residualbypredicted", "rstudentbypredicted", "rstudentbyleverage",
"qqplot", "observedbypredicted", "cooksd", "residualhistogram", "rfplot")
nv <- append(typs, vins, after = pos)
typs <- nv[-pos]
}
ret <- list()
# "diagnostics", "residuals", "fitplot", 'qqplot', 'rfplot',
# "residualbypredicted", "rstudentbypredicted", "rstudentbyleverage",
# "cooksd", "residualhistogram", "observedbypredicted"
for (tp in typs) {
if (tp == "diagnostics") {
ret[["diagnostics"]] <- render_diagnostics(dat, res, mdl, plt, alph)
} else if (tp == "residuals") { # Can be more than 1
resplts <- render_residuals(dat, res, mdl)
if ("plot_spec" %in% class(resplts)) {
ret[["residuals"]] <- resplts
} else {
for (idx in seq_len(length(resplts))) {
ret[[paste0("residuals", idx)]] <- resplts[[idx]]
}
}
} else if (tp == "fitplot") {
ret[["fitplot"]] <- render_fitplot(dat, res, mdl, plt, alph)
} else if (tp == "qqplot") {
ret[["qqplot"]] <- render_qqplot(dat, res, mdl)
} else if (tp == "rfplot") {
ret[["rfplot"]] <- render_rfplot(dat, res, mdl)
} else if (tp == "cooksd") {
ret[["cooksd"]] <- render_cooksd(dat, res, mdl, plt)
} else if (tp == "residualbypredicted") {
ret[["residualbypredicted"]] <- render_residualbypredicted(dat, res, mdl)
} else if (tp == "rstudentbypredicted") {
ret[["rstudentbypredicted"]] <- render_rstudentbypredicted(dat, res, mdl, plt)
} else if (tp == "rstudentbyleverage") {
ret[["rstudentbyleverage"]] <- render_rstudentbyleverage(dat, res, mdl, plt)
} else if (tp == "observedbypredicted") {
ret[["observedbypredicted"]] <- render_observedbypredicted(dat, res, mdl, plt)
} else if (tp == "residualhistogram") {
ret[["residualhistogram"]] <- render_residualhistogram(dat, res, mdl, plt)
} else if (tp == "residualboxplot") {
ret[["residualboxplot"]] <- render_residualboxplot(dat, res, mdl, plt)
} else if (tp == "dffits") {
ret[["dffits"]] <- render_dffits(dat, res, mdl, plt)
} else if (tp == "dfbetas") {
resplts <- render_dfbetas(dat, res, mdl, plt)
for (idx in seq_len(length(resplts))) {
ret[[paste0("dfbetas", idx)]] <- resplts[[idx]]
}
}
}
if (length(ret) == 0) {
ret <- NULL
}
}
return(ret)
}
# Diagnostics Panel -------------------------------------------------------
# Deal with "none" stats parameter. Adding extra box plot.
#' @import grDevices
#' @noRd
render_diagnostics <- function(dat, res, mdl, plt, alph) {
op <- par("mar")
om <- par("oma")
ofig <- par("fig")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 6 # Height in inches
wdi <- 6 # Width in inches
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Remove rows with NAs
datclean <- remove_na_rows(dat, dvr, ivr)
# Set margins
par(mfrow = c(3, 3),
mar = c(2, 2, 1, 0) + 0.1,
oma = c(0, 0, 3, .75))
#******************************
# 1) Residual By Predicted
#******************************
# Set margins
par(mar = c(4, 4, 0, 0) + 0.1)
# Prepare data
rdt <- res$OutputStatistics$RESID
pdt <- res$OutputStatistics$PREVAL
# Get scales
xscl <- get_scale(pdt, .05)
yscl <- get_zero_scale(rdt, .05)
# Generate plot
plot(pdt, rdt,
main = "",
xlab = "",
ylab = "",
pch = 1, # open circles
cex = 1.3,
col = "#05379B",
xlim = xscl,
ylim = yscl,
axes = FALSE
)
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# Generate zero line
abline(h = 0, col = "grey60", lwd = 1)
# Add labels
mtext("Predicted Value", side = 1, line = par("mar")[1] - 2, cex = .9)
mtext("Residual", side = 2, line = par("mar")[2] - 2, cex = .9)
# Frame
box(col = "grey70", lwd = 1)
#******************************
# 2) RStudent By Predicted
#******************************
# Prepare data
fit <- lm(mdl, data = dat)
rdt <- rstudent(fit)
pdt <- fitted(fit)
# Set margins
par(mar = c(4, 4, 0, 0) + 0.1)
# Get scales
xscl <- get_scale(pdt, .05)
yscl <- get_zero_scale(rdt, .1)
# Generate plot
plot(pdt, rdt,
main = "",
xlab = "",
ylab = "",
pch = 1, # open circles
cex = 1.3,
col = "#05379B",
xlim = xscl,
ylim = yscl,
axes = FALSE
)
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# Generate reference lines
abline(h = 2, col = "grey60", lwd = 1)
abline(h = -2, col = "grey60", lwd = 1)
# Add bottom label
mtext("Predicted Value", side = 1, line = par("mar")[1] - 2, cex = .9)
mtext("RStudent", side = 2, line = par("mar")[2] - 2, cex = .9)
# Frame
box(col = "grey70", lwd = 1)
#******************************
# 3) RStudent By Leverage
#******************************
# Prepare data
fit <- lm(mdl, data = dat)
rdt <- rstudent(fit)
ldt <- hatvalues(fit)
n <- length(ldt)
p <- length(coef(fit)) # includes intercept
# PROC REG-style cutoff
cut <- 2 * p / n
# Set margins
par(mar = c(4, 4, 0, 0) + 0.1)
# Get scales
xscl <- c(min(ldt) * .9, max(ldt, cut) * 1.10)
yscl <- get_zero_scale(rdt, .1)
# Define colors
nc <- "#05379B" # Normal color
# Generate plot
plot(ldt, rdt,
main = "",
xlab = "",
ylab = "",
pch = 1, # open circles
cex = 1.3,
col = nc,
xlim = xscl,
ylim = yscl,
axes = FALSE
)
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# Generate reference lines
abline(h = 2, col = "grey60", lwd = 1)
abline(h = -2, col = "grey60", lwd = 1)
abline(v = cut, col = "grey60", lwd = 1) # Cutoff
# Add labels
mtext("Leverage", side = 1, line = par("mar")[1] - 2, cex = .9)
mtext("RStudent", side = 2, line = par("mar")[2] - 2, cex = .9)
# Frame
box(col = "grey70", lwd = 1)
# Add labels
if (plt$label) {
# Identify points that need labels
lblids <- ifelse((rdt > 2 | rdt < -2) | ldt > cut, TRUE, FALSE)
if (!is.null(plt$id)) {
if (!plt$id %in% names(dat)) {
stop("Label ID variable '", plt$id, "' not found in input dataset.")
}
lbls <- dat[[plt$id]][lblids]
} else {
lbls <- seq(1, length(rdt))[lblids]
}
# Filter data
lrdt <- rdt[lblids]
lldt <- ldt[lblids]
# Assign labels
if (length(lbls) > 0) {
text(lldt, lrdt, labels = lbls, pos = 3)
}
}
#******************************
# 4) QQPlot
#******************************
# Prepare data
rdt <- res$OutputStatistics$RESID
# Set margins
par(mar = c(4, 4, 0, 0) + 0.1)
# Prepare Y data
ydat <- sort(rdt)
n <- length(ydat)
# SAS (Blom) plotting position formula: (i - 0.375) / (n + 0.25)
i <- 1:n
p_sas <- (i - 0.375) / (n + 0.25)
# Convert to theoretical normal quantiles
xdat <- qnorm(p_sas)
# Plot manually
plot(xdat, ydat,
main = "",
xlab = "",
ylab = "",
pch = 1,
col = "#05379B",
cex = 1.3,
axes = FALSE)
# Slope = Standard Deviation, Intercept = Mean
abline(a = mean(ydat), b = sd(ydat), col = "grey60")
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# X Axis label
mtext("Quantile", side = 1, line = par("mar")[1] - 2, cex = .9)
mtext("Residual", side = 2, line = par("mar")[2] - 2, cex = .9)
# Frame
box(col = "grey70", lwd = 1)
#******************************
# 5) Observed by Predicted
#******************************
# Prepare data
odt <- datclean[[dvr]]
pdt <- res$OutputStatistics$PREVAL
# Set margins
par(mar = c(4, 4, 0, 0) + 0.1)
# Get scales
xscl <- get_scale(pdt, 0)
yscl <- get_scale(odt, 0)
scl <- c(min(xscl[1], yscl[1]), max(xscl[2], yscl[2]))
# Generate plot
plot(pdt, odt,
main = "",
xlab = "",
ylab = "",
pch = 1, # open circles
cex = 1.3,
col = "#05379B",
xlim = scl,
ylim = scl,
axes = FALSE
)
# Add custom axes
axis(side = 1, col.ticks = "grey55", tck = -0.015, mgp = c(3, .5, 0))
axis(side = 2, las = 1, col.ticks = "grey55", tck = -0.015, mgp = c(3, .5, 0))
# Generate diagonal line
usr <- par("usr")
segments(usr[1], usr[3], usr[2], usr[4], col = "grey60")
# Add labels
mtext("Predicted Value", side = 1, line = par("mar")[1] - 2, cex = .9)
mtext(dvr, side = 2, line = par("mar")[1] - 2, cex = .9)
# Frame
box(col = "grey70", lwd = 1)
#******************************
# 6) Cook's D
#******************************
# Set margins
par(mar = c(4, 5, 0, 0) + 0.1)
# Calculate statistics
fit <- lm(mdl, data = dat)
cd <- cooks.distance(fit)
xcd <- seq_along(cd)
n <- length(cd)
p <- length(coef(fit)) # includes intercept
## Calculate cutoff
cutoff <- 4 / n
plot(xcd, cd,
type = "h", # vertical spikes
lwd = 1,
col = "#05379B",
xlab = "",
ylab = "",
main = "",
ylim = c(0, max(cd, cutoff) * 1.10),
axes = FALSE)
## Open circles at the top of each spike
points(seq_along(cd), cd,
pch = 1,
col = "#05379B", cex = 1.3)
## Reference line (PROC REG default)
abline(h = cutoff, col = "grey60", lwd = 1)
# Add custom axes
axis(side = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Generate zero line
abline(h = 0, col = "#05379B", lwd = 1)
# Add labels
mtext("Observation", side = 1, line = par("mar")[1] - 2, cex = .9)
mtext("Cook's D", side = 2, line = par("mar")[1] - 1, cex = .9)
# Frame
box(col = "grey70", lwd = 1)
# Add labels
if (plt$label) {
# Identify points that need labels
lblids <- ifelse(cd > cutoff, TRUE, FALSE)
if (!is.null(plt$id)) {
if (!plt$id %in% names(dat)) {
stop("Label ID variable '", plt$id, "' not found in input dataset.")
}
lbls <- dat[[plt$id]][lblids]
} else {
lbls <- seq(1, length(cd))[lblids]
}
# Filter data
lcd <- cd[lblids]
lxcd <- xcd[lblids]
# Filter colors
# lcol <- cols[cols != nc]
# Assign labels
if (length(lbls) > 0) {
text(lxcd, lcd, labels = lbls, pos = 3)
}
}
#******************************
# 7) Residual Histogram
#******************************
# Prepare data
rdt <- res$OutputStatistics$RESID
# Set margins
par(mar = c(4, 4, 0, 0) + 0.1)
# Calculate stats
n <- length(rdt)
mu <- mean(rdt)
sdx <- sd(rdt)
# Calculate breaks
brks <- pretty(range(rdt), n = nclass.Sturges(rdt),
min.n = 1, high.u.bias = 3)
# Histogram (Percent scale)
h <- hist(rdt,
breaks = brks,
plot = FALSE)
# Convert counts to percent
h$counts <- h$counts / sum(h$counts) * 100
# Use calculated breaks to get scale
scl <- get_reg_xlim(h$breaks, mu, sdx, pad_frac = 0, max_sigma = 4)
# Use calculated breaks to get scale
w <- diff(h$breaks)[1] # bin width
xmin <- min(h$breaks) - w
xmax <- max(h$breaks) + w
scl1 <- c(xmin, xmax)
# Use calculated break to get normal curve
grid <- seq(scl[1] , scl[2], length.out = 300)
y_norm_percent <- 100 * w * dnorm(grid, mean = mu, sd = sdx)
y_mx <- max(max(y_norm_percent), max(h$counts))
# Create plot using histogram values
plot(h,
col = "#CAD5E5", #"grey85",
border = "grey20",
main = "",
xlab = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.05),
axes = FALSE)
# Normal curve overlay (scaled to percent)
lines(grid, y_norm_percent, col = "steelblue4", lwd = 2)
# Add x axis label
mtext("Residual", side = 1, line = par("mar")[1] - 2, cex = .9)
mtext("Percent", side = 2, line = par("mar")[1] - 2, cex = .9)
# Calculate x axis ticks
df <- abs(h$mids[1] - h$mids[2])
te <- (df * seq(1, 10)) + h$mids[length(h$mids)]
be <- (-(df * seq(10, 1))) + h$mids[1]
df <- c(be, h$mids, te)
xtks <- df[df > scl[1] & df < scl[2]]
# Add custom axes
axis(side = 1, col.ticks = "grey55", at = xtks, labels = TRUE,
mgp = c(3, .5, 0), tck = -0.015, cex.axis = .75)
axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Frame
box(col = "grey70", lwd = 1)
#******************************
# 8) RF Plot
#******************************
# Prepare data
rdt <- res$OutputStatistics$RESID
fdt <- res$OutputStatistics$PREVAL
# dt <- dat[[ivr]]
# Get length of predicted values
n <- length(fdt)
## Proportion less (SAS-style)
p <- (seq_len(n) - 0.5) / n
## Left panel: Fit Mean
fit_centered <- sort(fdt - mean(fdt))
## Right panel: Residuals
res_sorted <- sort(rdt)
# Get y scale
smin <- min(fit_centered)
if (min(res_sorted) < smin) {
smin <- min(res_sorted)
}
smax <- max(fit_centered)
if (min(res_sorted) > smax) {
smax <- max(res_sorted)
}
# Finalize y scale
yscl <- c(smin, smax) * 1.05
# Set margins
par(mar = c(4, 4, 1.5, 0) + 0.1,
fig = c(.334, .53, 0, .334), new = TRUE) # Can place in specific spot
## Left: Fit Mean
plot(p, fit_centered,
pch = 1,
col = "#05379B",
xlab = "",
ylab = "",
ylim = yscl,
main = "Fit-Mean",
font.main = 1,
cex.main = 1.4,
axes = FALSE)
box(col = "grey70", lwd = 1)
# Add custom axes
axis(side = 1, col.ticks = "grey55", col = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", las = 1, col = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
## Plot layout
par(mar = c(4, .2, 1.5, 0) + 0.1,
fig = c(.53, .667, 0, .334), new = TRUE)
## Right: Residual
plot(p, res_sorted,
pch = 1,
col = "#05379B",
xlab = "",
ylab = "",
ylim = yscl,
main = "Residual",
font.main = 1,
cex.main = 1.4,
axes = FALSE)
box(col = "grey70", lwd = 1)
# Add custom axes
axis(side = 1, col.ticks = "grey55", col = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
## Add Labels
mtext(" Proportion Less",
side = 1, outer = TRUE,
line = par("oma")[1] - 2, cex = .9)
#******************************
# 9) Statistics or Boxplot
#******************************
if (!"none" %in% plt$stats) {
# Set margins and figure placement
par(mar = c(4, 5, 0, 2) + 0.1,
fig = c(.667, 1, 0, .334), new = TRUE)
# Initialize empty plot
# to create drawing area.
plot.new()
# Get needed statistics as named vector
stats <- collect_stats(res, plt)
if ("aic" %in% plt$stats) {
aicv <- round(extractAIC(fit)[2], 4)
stats <- append(stats, aicv)
names(stats)[length(stats)] <- "AIC"
}
# Set parameters
box_col = "white"
border_col = "grey70"
pad_y = 0.03 # vertical padding as fraction of y-range
pad_x = 0.03
# Get environment settings
usr <- par("usr") # c(x1, x2, y1, y2)
x_rng <- diff(usr[1:2])
y_rng <- diff(usr[3:4])
# x positions (user sets; relative to right edge)
x_left <- 0
x_right <- 1
# Determine line spacing
line_h <- strheight("Mg") # approx line height in user units
gap_h <- 0.8 * line_h # gap between lines
inner_h <- length(stats) * line_h + (length(stats) - 1) * gap_h
# box padding in user units
py <- pad_y * y_rng
px <- pad_x * x_rng
box_h <- inner_h + 2 * py
# center vertically in plot area
y_mid <- mean(usr[3:4])
y_bot <- y_mid - box_h / 2
y_top <- y_mid + box_h / 2
# draw box
rect(x_left, y_bot - (py/2), x_right, y_top, col = box_col, border = border_col)
# y positions for each row (top to bottom)
y_start <- y_top - py - line_h / 2
y_pos <- y_start - (seq_along(stats) - 1) * (line_h + gap_h)
# left column: names
text(x_left + px, y_pos,
labels = names(stats),
adj = c(0, 0.5), cex= 1.25)
# right column: values
text(x_right - px, y_pos,
labels = unname(stats),
adj = c(1, 0.5), cex= 1.25)
} else { # Boxplot
# Prepare data
rdt <- res$OutputStatistics$RESID
# Set margins
par(mar = c(4, 5, 0, 0) + 0.1,
fig = c(.667, 1, 0, .334), new = TRUE)
# Get scales
yscl <- get_zero_scale(rdt, .05)
# Calculate SAS compatible quantile
q <- quantile(rdt, probs = c(0, .25, .5, .75, 1), type = 2)
# Prepare to pass data
bp <- list(
stats = matrix(q, ncol = 1),
n = length(rdt),
conf = NULL,
out = numeric(0)
)
# Box plot
bxp(bp,
main = "",
ylab = "",
boxfill = "#CAD5E5",
border = "grey40",
outline = FALSE,
lty = 1,
ylim = yscl,
whiskcol = "#05379B",
staplecol = "#05379B",
medcol = "#05379B",
medlwd = 1,
boxwex = .65,
axes = FALSE)
# Mean diamond
points(1, mean(rdt),
pch = 5, # diamond
col = "#05379B",
cex = 1.4)
# Add custom axis
aval <- axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Label
mtext("Residual", side = 2, line = par("mar")[1] - 1, cex = .9)
# Frame
box(col = "grey70", lwd = 1)
}
#******************************
# Finish up
#******************************
# Set title
mtext(paste0("Fit Diagnostics for ", dvr), side = 3, cex = 1.2,
outer = TRUE, font = 2, line = par("oma")[3] - 2.25)
# Frame
box("outer", col = "grey70", lwd = 1)
# Restore margins
par(mar = op,
oma = om,
mfrow = c(1, 1),
fig = ofig,
new = FALSE)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
# Individual Plots --------------------------------------------------------
# What if there are multiple independent variables?
# Need to show separate plot for each independent variable, all in
# same panel.
#' @noRd
render_residuals <- function(dat, res, mdl) {
op <- par("mar")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Max per panel
pnlmx <- 6
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Set titles
ttl <- "Residuals for "
if (length(ivr) > 1) {
ttl <- "Residual by Regressors for "
}
# Determine number of panels and panel layout
if (length(ivr) == 1) {
pnls <- 1
rws <- c(1, 1)
} else if (length(ivr) == 2) {
pnls <- 1
rws <- c(1, 2)
hti <- 3.5
} else if (length(ivr) %in% c(3, 4)) {
pnls <- 1
rws <- c(2, 2)
} else if (length(ivr) %in% c(5, 6)) {
pnls <- 1
rws <- c(2, 3)
} else if (length(ivr) > pnlmx) {
pnls <- ceiling(length(ivr) / pnlmx)
rws <- c(2, 3)
}
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Set up object to return
ret <- list()
# Get residuals
rdt <- res$OutputStatistics$RESID
# Remove rows with NAs
dat <- remove_na_rows(dat, dvr, ivr)
# Split independent variables by panel
if (length(ivr) > pnlmx) {
fct <- ceiling(seq_along(ivr) / pnlmx)
ivrs <- split(ivr, fct)
} else {
ivrs <- list(ivr)
}
for (pnl in seq(1, pnls)) {
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Set up margins and grid
par(oma = c(0, 4, 2.5, .75) + 0.1,
mar = c(4, 1, 0, 0),
mfrow = rws)
# Flag for first chart
firstChart <- TRUE
# Chart counter
chtcnt <- 0
for (iv in ivrs[[pnl]]) {
# Increment counter
chtcnt <- chtcnt + 1
# Get data for independent variable
dt <- dat[[iv]]
# Get y scale
mx <- max(abs(rdt)) * 1.05
scl <- c(-mx, mx)
# Generate plot
plot(dt, rdt,
main = "",
xlab = "",
ylab = "",
pch = 1, # open circles
cex = 1.3,
col = "#05379B",
ylim = scl,
axes = FALSE
)
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
if (firstChart) {
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
}
# Put x axis label
mtext(iv, side = 1, line = par("mar")[1] - 1.5)
# Generate zero line
abline(h = 0, col = "grey60", lwd = 1)
# Frame
box(col = "grey70", lwd = 1)
# Set flag
if (chtcnt %% par("mfrow")[2] == 0) {
firstChart <- TRUE
} else {
firstChart <- FALSE
}
}
# Title
mtext(paste0(ttl, dvr), side = 3,
font = 2, cex = 1.2, outer = TRUE, line = par("oma")[3] - 2)
# Y Axis lable
mtext("Residual", side = 2, outer = TRUE, line = par("oma")[2] - 1.5)
# Draw line around entire chart
box("outer", col = "grey70", lwd = 1)
# Close device context
dev.off()
# Put plot in reporter plot object
ret[[length(ret) + 1]] <- create_plot(pth, height = hti, width = wdi)
}
# Restore margins
par(mar = op, mfrow = c(1, 1))
if (length(ret) == 1) {
ret <- ret[[1]]
}
return(ret)
}
render_residualbypredicted <- function(dat, res, mdl) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Prepare data
rdt <- res$OutputStatistics$RESID
pdt <- res$OutputStatistics$PREVAL
# Set margins
par(mar = c(4, 5, 2, .75) + 0.1)
# Get scales
xscl <- get_scale(pdt, .05)
yscl <- get_zero_scale(rdt, .05)
# Generate plot
plot(pdt, rdt,
main = paste0("Residual by Predicted for ", dvr),
xlab = "",
ylab = "Residual",
pch = 1, # open circles
cex = 1.3,
col = "#05379B",
xlim = xscl,
ylim = yscl,
axes = FALSE
)
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# Generate zero line
abline(h = 0, col = "grey60", lwd = 1)
# Add bottom label
mtext("Predicted Value", side = 1, line = par("mar")[1] - 2)
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
# Restore margins
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
render_rstudentbypredicted <- function(dat, res, mdl, plt) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Prepare data
fit <- lm(mdl, data = dat)
rdt <- rstudent(fit)
pdt <- fitted(fit)
# Set margins
par(mar = c(4, 5, 2, .75) + 0.1)
# Get scales
xscl <- get_scale(pdt, .05)
yscl <- get_zero_scale(rdt, .1)
# Generate plot
plot(pdt, rdt,
main = paste0("RStudent by Predicted for ", dvr),
xlab = "",
ylab = "RStudent",
pch = 1, # open circles
cex = 1.3,
col = "#05379B",
xlim = xscl,
ylim = yscl,
axes = FALSE
)
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# Generate reference lines
abline(h = 2, col = "grey60", lwd = 1)
abline(h = -2, col = "grey60", lwd = 1)
# Add bottom label
mtext("Predicted Value", side = 1, line = par("mar")[1] - 2)
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
# Add labels
if (plt$label) {
# Identify points that need labels
lblids <- ifelse(rdt > 2 | rdt < -2, TRUE, FALSE)
if (!is.null(plt$id)) {
if (!plt$id %in% names(dat)) {
stop("Label ID variable '", plt$id, "' not found in input dataset.")
}
lbls <- dat[[plt$id]][lblids]
} else {
lbls <- seq(1, length(rdt))[lblids]
}
# Filter data
lrdt <- rdt[lblids]
lpdt <- pdt[lblids]
# Filter colors
# lcol <- cols[cols != nc]
# Assign labels
if (length(lrdt) > 0) {
text(lpdt, lrdt, labels = lbls, pos = 3)
}
}
# Restore margins
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
render_rstudentbyleverage <- function(dat, res, mdl, plt) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Prepare data
fit <- lm(mdl, data = dat)
rdt <- rstudent(fit)
ldt <- hatvalues(fit)
n <- length(ldt)
p <- length(coef(fit)) # includes intercept
# PROC REG-style cutoff
cut <- 2 * p / n
# Set margins
par(mar = c(6, 5, 2, .75) + 0.1)
# Get scales
xscl <- c(min(ldt) * .9, max(ldt, cut) * 1.10)
yscl <- get_zero_scale(rdt, .1)
# Define colors
nc <- "#05379B" # Normal color
lvc <- "#01665E" # Leverage color
otc <- "#A23A2E" # Outlier color
olc <- "#543005" # Outlier and leverage
# Assign colors
cols <- rep(nc, length(ldt))
cols <- ifelse(ldt > cut, lvc, cols)
cols <- ifelse(rdt > 2 | rdt < -2, otc, cols)
cols <- ifelse((rdt > 2 | rdt < -2) & ldt > cut, olc, cols)
# Generate plot
plot(ldt, rdt,
main = paste0("Outlier and Leverage Diagnostics for ", dvr),
xlab = "",
ylab = "RStudent",
pch = 1, # open circles
cex = 1.3,
col = cols,
xlim = xscl,
ylim = yscl,
axes = FALSE
)
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# Generate reference lines
abline(h = 2, col = "grey60", lwd = 1)
abline(h = -2, col = "grey60", lwd = 1)
abline(v = cut, col = "grey60", lwd = 1) # Cutoff
# Add bottom label
mtext("Leverage", side = 1, line = par("mar")[1] - 4)
# Add labels
if (plt$label) {
# Identify points that need labels
lblids <- ifelse((rdt > 2 | rdt < -2) | ldt > cut, TRUE, FALSE)
if (!is.null(plt$id)) {
if (!plt$id %in% names(dat)) {
stop("Label ID variable '", plt$id, "' not found in input dataset.")
}
lbls <- dat[[plt$id]][lblids]
} else {
lbls <- seq(1, length(rdt))[lblids]
}
# Filter data
lrdt <- rdt[lblids]
lldt <- ldt[lblids]
# Filter colors
lcol <- cols[cols != nc]
# Assign labels
if (length(lbls) > 0) {
text(lldt, lrdt, labels = lbls, pos = 3, col = lcol)
}
}
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
# Legend (bottom, SAS-style)
legend("bottom",
legend = c("Outlier ",
"Leverage",
"Outlier and Leverage"),
pch = c(1, 1, 1),
pt.cex = 1.3,
col = c(otc,
lvc,
olc),
horiz = TRUE,
bty = "o",
box.col = "grey", # Grey border to match SAS
x.intersp = 1, # Spacing between group label and box
y.intersp = 0, # Spacing between content and borders
text.width = NA, # Compute label widths dynamically
inset = c(0, -0.23),
xpd = TRUE)
# Restore margins
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
# Kill if more than 1 independent variable.
# Kill stats if "none" parameter is passed.
#' @noRd
render_fitplot <- function(dat, res, mdl, plt, alph) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
if (length(ivr) == 1) {
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Set Confidence limits
alph1 <- 1 - alph
alph2 <- alph1 * 100
# Fit model
fit <- lm(mdl, data = dat)
# Prediction grid
xnew <- seq(min(dat[[ivr]], na.rm = TRUE), max(dat[[ivr]], na.rm = TRUE),
length.out = 200)
# Create data frame
df1 <- data.frame(xnew)
names(df1) <- ivr
conf <- predict(fit,
newdata = df1,
interval = "confidence", level = alph1)
pred <- predict(fit,
newdata = df1,
interval = "prediction", level = alph1)
# Set margins
if ("none" %in% plt$stats) {
par(mar = c(5, 4.5, 2, .75) + 0.1)
} else {
par(mar = c(5, 4.5, 2, 8.5) + 0.1)
}
# Plot
plot(dat[[ivr]], dat[[dvr]],
main = paste0("Fit Plot for ", dvr),
xlab = "",
ylab = dvr,
pch = 1,
cex = 1.23,
col = "#05379B",
ylim = range(pred),
axes = FALSE)
# X Axis label
mtext(ivr, side = 1, line = par("mar")[1] - 3.25)
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# 95% confidence band (shaded)
polygon(c(xnew, rev(xnew)),
c(conf[, "lwr"], rev(conf[, "upr"])),
col = adjustcolor("steelblue", alpha.f = 0.30),
border = NA)
# Fitted line
lines(xnew, conf[, "fit"],
col = "steelblue4",
lwd = 2)
# 95% prediction limits (dashed)
lines(xnew, pred[, "lwr"],
col = "steelblue3",
lty = 2)
lines(xnew, pred[, "upr"],
col = "steelblue3",
lty = 2)
# Legend (bottom, SAS-style)
legend("bottom",
legend = c("Fit ",
paste0(alph2, "% Confidence Limits "),
paste0(alph2, "% Prediction Limits")),
lwd = c(2, NA, 1),
lty = c(1, NA, 2),
pch = c(NA, 15, NA),
pt.cex = 2,
col = c("steelblue4",
adjustcolor("steelblue", 0.30),
"steelblue3"),
horiz = TRUE,
bty = "o",
box.col = "grey", # Grey border to match SAS
x.intersp = .7, # Spacing between group label and box
y.intersp = 0, # Spacing between content and borders
text.width = NA, # Compute label widths dynamically
inset = c(0, -0.19),
xpd = TRUE)
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
if (!"none" %in% plt$stats) {
# Gather statistics for box
sts <- collect_stats(res, plt)
if ("aic" %in% plt$stats) {
aicv <- round(extractAIC(fit)[2], 4)
sts <- append(sts, aicv)
names(sts)[length(sts)] <- "AIC"
}
# Create box
draw_stats_box(sts, x_frac = c(0.015, 0.3))
}
# Restore margins
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
} else {
ret <- NULL
}
return(ret)
}
#' @noRd
render_residualhistogram <- function(dat, res, mdl, plt) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Prepare data
rdt <- res$OutputStatistics$RESID
# dt <- dat[[ivr]]
# Set margins
par(mar = c(4, 5, 2, .75) + 0.1)
# Calculate stats
n <- length(rdt)
mu <- mean(rdt)
sdx <- sd(rdt)
# Calculate breaks - Closer to SAS algorithm
# brks <- pretty(range(rdt), n = nclass.Sturges(rdt),
# min.n = 1, high.u.bias = 3)
# Calculate breaks - Closer to SAS algorithm
brks <- pretty(range(rdt), n = nclass.Sturges(rdt),
min.n = 1, high.u.bias = 3)
# Histogram (Percent scale)
h <- hist(rdt,
breaks = brks,
plot = FALSE)
# Convert counts to percent
h$counts <- h$counts / sum(h$counts) * 100
# Use calculated breaks to get normal scale
w <- diff(h$breaks)[1] # bin width
xmin <- min(h$breaks) - w
xmax <- max(h$breaks) + w
scl1 <- c(xmin, xmax)
# Use calculated breaks to get scale
scl <- get_reg_xlim(h$breaks, mu, sdx, pad_frac = 0, max_sigma = 4)
# Use calculated break to get normal curve
grid <- seq(scl1[1] , scl1[2], length.out = 300)
y_norm_percent <- 100 * w * dnorm(grid, mean = mu, sd = sdx)
y_mx <- max(max(y_norm_percent), max(h$counts))
# Create plot using histogram values
plot(h,
col = "#CAD5E5", #"grey85",
border = "grey20",
main = paste0("Distribution of Residuals for ", dvr),
xlab = "",
ylab = "Percent",
xlim = scl,
ylim = c(0, y_mx * 1.05), # max(h$counts) * 1.05),
axes = FALSE)
# Normal curve overlay (scaled to percent)
# x <- seq(min(rdt) * 2, max(rdt) * 2, length.out = 300)
#
# y_norm <- dnorm(x, mean = mu, sd = sdx)
# y_norm <- y_norm / max(y_norm) * max(h$counts)
lines(grid, y_norm_percent, col = "steelblue4", lwd = 2)
# Kernel density overlay (scaled to percent)
# dens <- density(rdt)
#
# y_kern <- dens$y / max(dens$y) * max(h$counts)
# lines(dens$x, y_kern, col = "orangered2", lwd = 2)
# Bandwidth (Silverman / SAS)
bw <- 1.06 * sdx * n^(-1/5)
# Kernel density estimate
dens <- density(rdt, kernel = "gaussian", bw = bw, n = 512)
# Scale KDE to Percent axis
y_kernel_percent <- 100 * w * dens$y
yscl <- seq(50, length(y_kernel_percent) - 25)
# Overlay on histogram
lines(dens$x[yscl], y_kernel_percent[yscl], col = "orangered2", lwd = 2)
# Add x axis label
mtext("Residual", side = 1, line = par("mar")[1] - 2)
# Calculate x axis ticks
df <- abs(h$mids[1] - h$mids[2])
te <- (df * seq(1, 10)) + h$mids[length(h$mids)]
be <- (-(df * seq(10, 1))) + h$mids[1]
df <- c(be, h$mids, te)
xtks <- df[df > scl[1] & df < scl[2]]
# Add custom axes
axis(side = 1, col.ticks = "grey55", at = xtks, labels = TRUE,
mgp = c(3, .5, 0), tck = -0.015, cex.axis = .75)
axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
# Legend
mpos <- legend("topright",
legend = c("Normal", "Kernel"),
col = c("steelblue4", "orangered2"),
lwd = 2,
box.col = "grey80",
bty = "n",
cex = .9,
inset = c(.01, .01),
x.intersp = .5,
y.intersp = c(0.5, 1.8)
)
rect(mpos$rect$left, mpos$rect$top - (mpos$rect$h - .5),
mpos$rect$left + mpos$rect$w, mpos$rect$top,
border = "grey80"
)
# Restore margins
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
#' @noRd
render_qqplot <- function(dat, res, mdl) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Prepare data
rdt <- res$OutputStatistics$RESID
# Set margins
par(mar = c(4, 5, 2, .75) + 0.1)
# Prepare Y data
ydat <- sort(rdt)
n <- length(ydat)
# SAS (Blom) plotting position formula: (i - 0.375) / (n + 0.25)
i <- 1:n
p_sas <- (i - 0.375) / (n + 0.25)
# Convert to theoretical normal quantiles
xdat <- qnorm(p_sas)
# Plot manually
plot(xdat, ydat,
main = paste0("Q-Q Plot of Residuals for ", dvr),
xlab = "",
ylab = "Residual",
pch = 1,
col = "#05379B",
cex = 1.3,
axes = FALSE)
# Slope = Standard Deviation, Intercept = Mean
abline(a = mean(ydat), b = sd(ydat), col = "grey60")
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# X Axis label
mtext("Quantile", side = 1, line = par("mar")[1] - 2)
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
render_rfplot <- function(dat, res, mdl) {
# res <- proc_reg(cls,
# model = "Weight = Height",
# output = report,
# plots = regplot(type = "qqplot"),
# stats = p)
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Prepare data
rdt <- res$OutputStatistics$RESID
fdt <- res$OutputStatistics$PREVAL
# dt <- dat[[ivr]]
# Set margins
par(mar = c(5, 5, 2, .75) + 0.1)
## Proportion less (SAS-style)
n <- length(fdt)
p <- (seq_len(n) - 0.5) / n
## Left panel: Fit Mean
fit_centered <- sort(fdt - mean(fdt, na.rm = TRUE))
## Right panel: Residuals
res_sorted <- sort(rdt)
# Get y scale
smin <- min(fit_centered)
if (min(res_sorted, na.rm = TRUE) < smin) {
smin <- min(res_sorted, na.rm = TRUE)
}
smax <- max(fit_centered, na.rm = TRUE)
if (min(res_sorted, na.rm = TRUE) > smax) {
smax <- max(res_sorted, na.rm = TRUE)
}
# Finalize y scale
yscl <- c(smin, smax) * 1.05
## Plot layout
par(mfrow = c(1, 2),
mar = c(2, 0, 2, .75) + 0.1,
oma = c(2, 4, 2, .25))
## Left: Fit Mean
plot(p, fit_centered,
pch = 1,
col = "#05379B",
xlab = "",
ylab = "",
ylim = yscl,
main = "Fit-Mean",
font.main = 1,
cex.main = 1,
axes = FALSE)
box(col = "grey70", lwd = 1)
# Add custom axes
axis(side = 1, col.ticks = "grey55", col = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", las = 1, col = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
## Right: Residual
plot(p, res_sorted,
pch = 1,
col = "#05379B",
xlab = "",
ylab = "",
ylim = yscl,
main = "Residual",
font.main = 1,
cex.main = 1,
axes = FALSE)
box(col = "grey70", lwd = 1)
# Add custom axes
axis(side = 1, col.ticks = "grey55", col = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
## Overall title
mtext(paste0("Residual-Fit Spread Plot for ", dvr),
outer = TRUE, cex = 1.25, font = 2)
# X Axis Label
mtext("Proportion Less", side = 1, outer = TRUE, cex = 1.1)
# Frame
box("outer", col = "grey70", lwd = 1)
# Restore margins
par(mar = op, mfrow = c(1, 1), oma = c(0, 0, 0, 0))
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
}
#' @noRd
render_cooksd <- function(dat, res, mdl, plt) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Set margins
par(mar = c(4, 5, 2, .75) + 0.1)
# Calculate statistics
fit <- lm(mdl, data = dat)
cd <- cooks.distance(fit)
xcd <- seq_along(cd)
n <- length(cd)
p <- length(coef(fit)) # includes intercept
## Calculate cutoff
cutoff <- 4 / n
plot(xcd, cd,
type = "h", # vertical spikes
lwd = 1,
col = "#05379B",
xlab = "",
ylab = "Cook's D",
main = paste0("Cook's D for ", dvr),
ylim = c(0, max(cd, cutoff) * 1.10),
axes = FALSE)
## Open circles at the top of each spike
points(seq_along(cd), cd,
pch = 1,
col = "#05379B", cex = 1.25)
## Reference line (PROC REG default)
abline(h = cutoff, col = "grey60", lwd = 1)
# Add custom axes
axis(side = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Generate zero line
abline(h = 0, col = "#05379B", lwd = 1)
# X axis label
mtext("Observation", side = 1, line = par("mar")[1] - 2)
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
# Add labels
if (plt$label) {
# Identify points that need labels
lblids <- ifelse(cd > cutoff, TRUE, FALSE)
if (!is.null(plt$id)) {
if (!plt$id %in% names(dat)) {
stop("Label ID variable '", plt$id, "' not found in input dataset.")
}
lbls <- dat[[plt$id]][lblids]
} else {
lbls <- seq(1, length(cd))[lblids]
}
# Filter data
lcd <- cd[lblids]
lxcd <- xcd[lblids]
# Assign labels
if (length(lbls) > 0) {
text(lxcd, lcd, labels = lbls, pos = 3)
}
}
# Restore margins
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
render_observedbypredicted <- function(dat, res, mdl, plt) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 4.5 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Remove rows with NAs
dat <- remove_na_rows(dat, dvr, ivr)
# Prepare data
odt <- dat[[dvr]]
pdt <- res$OutputStatistics$PREVAL
# Set margins
par(mar = c(4, 5, 2, 1.05) + 0.1)
# Get scales
# Appears to take min of predicted and observed
# and max of predicted and observed, and uses
# those for both x and y scales. Scales have to be the same
# or the dots won't line up properly. Much trial and error
# to figure this out.
xscl <- get_scale(pdt, 0)
yscl <- get_scale(odt, 0)
scl <- c(min(xscl[1], yscl[1]), max(xscl[2], yscl[2]))
# Generate plot
plot(pdt, odt,
main = paste0("Observed by Predicted for ", dvr),
xlab = "",
ylab = dvr,
pch = 1, # open circles
cex = 1.3,
col = "#05379B",
xlim = scl,
ylim = scl,
axes = FALSE
)
# Add custom axes
axis(side = 1, col.ticks = "grey55", tck = -0.015, mgp = c(3, .5, 0))
axis(side = 2, las = 1, col.ticks = "grey55", tck = -0.015, mgp = c(3, .5, 0))
# Generate diagonal line
# abline() doesn't work at all. Not sure why. Have to use segments().
usr <- par("usr")
segments(usr[1], usr[3], usr[2], usr[4], col = "grey60")
# Add bottom label
mtext("Predicted Value", side = 1, line = par("mar")[1] - 2)
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
# Would be nice to get outlier labels in here,
# but I don't know how they are identified.
# There are no limits on this chart.
# Restore margins
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
render_residualboxplot <- function(dat, res, mdl, plt) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
# All output types accept jpeg
# So start with that
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Prepare data
rdt <- res$OutputStatistics$RESID
pdt <- res$OutputStatistics$PREVAL
# Set margins
par(mar = c(.75, 5, 2, .75) + 0.1)
# Get scales
yscl <- get_zero_scale(rdt, .05)
# Calculate SAS compatible quantile
q <- quantile(rdt, probs = c(0, .25, .5, .75, 1), type = 2)
# Prepare to pass data
bp <- list(
stats = matrix(q, ncol = 1),
n = length(rdt),
conf = NULL,
out = numeric(0)
)
# Box plot
bxp(bp,
main = paste0("Residuals for ", dvr),
ylab = "Residual",
boxfill = "#CAD5E5",
border = "grey40",
outline = FALSE,
lty = 1,
ylim = yscl,
whiskcol = "#05379B",
staplecol = "#05379B",
medcol = "#05379B",
medlwd = 1,
boxwex = .65,
axes = FALSE)
# Mean diamond
points(1, mean(rdt, na.rm = TRUE),
pch = 5, # diamond
col = "#05379B",
cex = 1.3)
# Add custom axis
aval <- axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
## Optional: horizontal gridlines
abline(h = aval, col = "grey90", lwd = 1)
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
# Restore margins
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
render_dffits <- function(dat, res, mdl, plt) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Prepare data
fit <- lm(mdl, data = dat)
dff <- dffits(fit)
n <- length(dff)
p <- length(coef(fit)) # includes intercept
# Set margins
par(mar = c(4, 5, 2, .75) + 0.1)
## PROC REG cutoff for DFFITS
cutoff <- 2 * sqrt(p / n)
plot(seq_along(dff), dff,
type = "h", # vertical spikes
lwd = 1,
col = "#05379B",
xlab = "",
ylab = "DFFITS",
main = paste0("Influence Diagnostics for ", dvr),
ylim = range(c(dff, -cutoff, cutoff)) * 1.1,
axes = FALSE)
## Open circles at spike tips
points(seq_along(dff), dff,
pch = 1,
cex = 1.25,
col = "#05379B")
# Add custom axes
axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
## Reference lines (PROC REG)
abline(h = 0, col = "#05379B", lwd = 1)
abline(h = cutoff, col = "grey60", lwd = 1)
abline(h = -cutoff, col = "grey60", lwd = 1)
# Add bottom label
mtext("Observation", side = 1, line = par("mar")[1] - 2)
# Frame
box(col = "grey70", lwd = 1)
box("figure", col = "grey70", lwd = 1)
# Add labels
if (plt$label) {
# Identify points that need labels
lblids1 <- ifelse(dff > cutoff, TRUE, FALSE)
lblids2 <- ifelse(dff < -cutoff, TRUE, FALSE)
if (!is.null(plt$id)) {
if (!plt$id %in% names(dat)) {
stop("Label ID variable '", plt$id, "' not found in input dataset.")
}
lbls1 <- dat[[plt$id]][lblids1]
lbls2 <- dat[[plt$id]][lblids2]
} else {
lbls1 <- seq(1, length(dff))[lblids1]
lbls2 <- seq(1, length(dff))[lblids2]
}
# Filter data
ly1 <- dff[lblids1]
lx1 <- seq_along(dff)[lblids1]
ly2 <- dff[lblids2]
lx2 <- seq_along(dff)[lblids2]
# Filter colors
# lcol <- cols[cols != nc]
# Assign labels
if (length(ly1) > 0) {
text(lx1, ly1, labels = lbls1, pos = 3)
}
if (length(ly2) > 0) {
text(lx2, ly2, labels = lbls2, pos = 1)
}
}
# Restore margins
par(mar = op)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
render_dfbetas <- function(dat, res, mdl, plt) {
# Get current values
op <- par("mar")
om <- par("oma")
# Function can return multiple plots
ret <- list()
# Standard height and width
hti <- 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Get analysis variables
vrs <- get_vars(mdl)
dvr <- vrs$dvar
ivr <- vrs$ivar
# Prepare data
fit <- lm(mdl, data = dat)
dfb <- dfbetas(fit) # matrix: n x p
n <- nrow(dfb)
p <- ncol(dfb)
# Get names for labels
param_names <- colnames(dfb)
# Baseline
chrts <- 1
rngs <- seq(1, p)
if (plt$panel) {
if (p == 1) {
mfrw <- c(1, 1)
} else if (p == 2) {
mfrw <- c(1, 2)
} else if (p > 2 & p < 5) {
mfrw <- c(2, 2)
} else if (p > 4 & p < 7) {
mfrw <- c(2, 3)
} else if (p > 6) {
chrts <- ceiling(p / 6)
mfrw <- c(2, 3)
}
# Split into a list of ranges per chart
if (chrts > 1) {
rngs <- split(rngs, ceiling(rngs/6))
} else {
rngs <- list(rngs)
}
} else {
# One chart per variable
chrts <- p
rngs <- split(rngs, rngs)
mfrw <- c(1, 1)
}
## PROC REG cutoff for DFBETAS
cutoff <- 2 / sqrt(n)
# Get yscale - Consistent Y scale on panel plot
if (plt$label) {
yscl <- c(min(dfb, -cutoff, na.rm = TRUE), max(dfb, cutoff, na.rm = TRUE)) * 1.2 # Give more room for labels
} else {
yscl <- c(min(dfb, -cutoff, na.rm = TRUE), max(dfb, cutoff, na.rm = TRUE)) * 1.1
}
# Get xscale
xscl <- c(0 - .01, nrow(dat) * 1.05)
for (i in seq_len(chrts)) {
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Output to image file
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Set margins
if (plt$panel) {
par(mfrow = mfrw,
mar = c(2, 0, 2, .75) + 0.1,
oma = c(2, 5, 2, 0) + .1)
} else {
par(mfrow = mfrw,
mar = c(2, 0, 0, .75) + 0.1,
oma = c(2, 5, 2, 0) + .1)
}
for (j in rngs[[i]]) {
# Get data for this chart
vals <- dfb[, j]
# Non-panel plot computes custom y scale for each variable
if (plt$panel == FALSE) {
yscl <- c(min(vals, -cutoff), max(vals, cutoff)) * 1.1
}
# Create empty plot
plot(seq_len(n), vals,
type = "n",
xlim = xscl,
ylim = yscl,
axes = FALSE)
# Determine grid position
grdPos <- par("mfg")
# X axis
xtcks <- axis(side = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
# Y axis
if (grdPos[2] == 1) {
ytcks <- axis(side = 2, las = 1, col.ticks = "grey55", mgp = c(3, .5, 0), tck = -0.015)
}
# Grid lines
if (plt$panel) {
abline(h = ytcks, col = "grey90", lwd = 1)
abline(v = xtcks, col = "grey90", lwd = 1)
}
par(new = TRUE)
if (plt$panel) {
ttl <- param_names[j]
} else {
ttl <- ""
}
# Create actual data plot
plot(seq_len(n), vals,
type = "h",
font.main = 1,
lwd = 1,
col = "#05379B",
xlab = "",
ylab = "",
main = ttl,
xlim = xscl,
ylim = yscl, # range(c(vals, -cutoff, cutoff)) * 1.1,
axes = FALSE)
par(new = FALSE)
## Open circles at tips
if (plt$panel == FALSE) {
points(seq_len(n), vals,
pch = 1,
cex = 1.25,
col = "#05379B")
}
## Reference lines (PROC REG)
abline(h = 0, col = "#05379B", lwd = 1)
abline(h = cutoff, col = "grey60", lwd = 1)
abline(h = -cutoff, col = "grey60", lwd = 1)
# Add labels
if (plt$label) {
## Positive values
# Identify points that need labels
lblids <- ifelse(vals > cutoff, TRUE, FALSE)
if (!is.null(plt$id)) {
if (!plt$id %in% names(dat)) {
stop("Label ID variable '", plt$id, "' not found in input dataset.")
}
lbls <- dat[[plt$id]][lblids]
} else {
lbls <- seq(1, length(vals))[lblids]
}
# Filter data
xpos <- seq_len(n)[lblids]
ypos <- vals[lblids]
# Assign labels
if (length(xpos) > 0) {
text(xpos, ypos, labels = lbls, pos = 3)
}
## Negative values
# Identify points that need labels
lblids <- ifelse(vals < -cutoff, TRUE, FALSE)
if (!is.null(plt$id)) {
if (!plt$id %in% names(dat)) {
stop("Label ID variable '", plt$id, "' not found in input dataset.")
}
lbls <- dat[[plt$id]][lblids]
} else {
lbls <- seq(1, length(vals))[lblids]
}
# Filter data
xpos <- seq_len(n)[lblids]
ypos <- vals[lblids]
# Assign labels
if (length(xpos) > 0) {
text(xpos, ypos, labels = lbls, pos = 1)
}
}
box(col = "grey70", lwd = 1)
}
## Overall title
# Set title
if (plt$panel) {
mtext(paste0("Influence Diagnostics for ", dvr), side = 3, cex = 1.2,
outer = TRUE, font = 2, line = par("oma")[3] - 2.25)
} else {
vr <- param_names[j]
mtext(paste0("Influence on ", vr, " for ", dvr), side = 3, cex = 1.2,
outer = TRUE, font = 2, line = par("oma")[3] - 1.75)
}
# Labels
mtext("Observation", side = 1, outer = TRUE, line = par("oma")[1] - 2)
mtext("DFBETAS", side = 2, outer = TRUE, line = par("oma")[2] - 2)
# Frame
box("outer", col = "grey70", lwd = 1)
# Restore margins
par(mar = op, oma = om, mfrow = c(1, 1))
# Close device context
dev.off()
# Put plot in reporter plot object
ret[[length(ret) + 1]] <- create_plot(pth, height = hti, width = wdi)
}
return(ret)
}
# Utilities ---------------------------------------------------------------
#' @noRd
collect_stats <- function(res, plt) {
# mres <- proc_reg(cls, Weight ~ Height, stats = "clb", output = report)
ret <- c()
if ("default" %in% plt$stats) {
ret <- c(
"Observations" = res$NObs$NOBS[1],
"Parameters" = nrow(res$ParameterEstimates),
"Error DF" = res$ANOVA$DF[2],
"MSE" = roundup(res$ANOVA$MEANSQ[2], 2),
"R-Square" = roundup(res$FitStatistics$RSQ[1], 4),
"Adj R-Square" = roundup(res$FitStatistics$ADJRSQ[1], 4)
)
}
# Apparently got from AI and these are wrong
# obs, edf, rsquare, parms, mse, adjrsquare,
# coefvar, depmean
# sse, bic, gmsep, pc, sp, aic, cp, jp, sbc
# From SAS error message
# Documentation is not actually that good on how to do this.
# ADJRSQ, AIC, BIC, COEFFVAR, CP, DEFAULT, DEPMEAN, EDF, GMSEP, JP,
# MSE, NOBS, NPARM, PC, RSQUARE, SBC, SP, SSE
if ("nobs" %in% plt$stats) {
ret[length(ret) + 1] <- res$NObs$NOBS[1]
names(ret)[length(ret)] <- "Observations"
}
if ("nparm" %in% plt$stats) {
ret[length(ret) + 1] <- nrow(res$ParameterEstimates)
names(ret)[length(ret)] <- "Parameters"
}
if ("edf" %in% plt$stats) {
ret[length(ret) + 1] <- res$ANOVA$DF[2]
names(ret)[length(ret)] <- "Error DF"
}
if ("mse" %in% plt$stats) {
ret[length(ret) + 1] <- roundup(res$ANOVA$MEANSQ[2], 2)
names(ret)[length(ret)] <- "MSE"
}
if ("rsquare" %in% plt$stats) {
ret[length(ret) + 1] <- roundup(res$FitStatistics$RSQ[1], 4)
names(ret)[length(ret)] <- "R-Square"
}
if ("adjrsq" %in% plt$stats) {
ret[length(ret) + 1] <- roundup(res$FitStatistics$ADJRSQ[1], 4)
names(ret)[length(ret)] <- "Adj R-Square"
}
if ("depmean" %in% plt$stats) {
ret[length(ret) + 1] <- roundup(res$FitStatistics$DEPMEAN, 2)
names(ret)[length(ret)] <- "Dep Mean"
}
if ("coeffvar" %in% plt$stats) {
ret[length(ret) + 1] <- roundup(res$FitStatistics$COEFVAR, 4)
names(ret)[length(ret)] <- "Coeff Var"
}
if ("sse" %in% plt$stats) {
ret[length(ret) + 1] <- roundup(res$ANOVA$SUMSQ[2], 1)
names(ret)[length(ret)] <- "SSE"
}
return(ret)
}
#' @noRd
draw_stats_box <- function(stats,
x_frac = c(0.03, 0.30), # from right edge of plot (in x-range fractions)
cex = 0.85,
family = "sans",
box_col = "white",
border_col = "grey70",
pad_y = 0.01, # vertical padding as fraction of y-range
pad_x = 0.01) { # horizontal padding as fraction of x-range
stopifnot(is.numeric(stats), !is.null(names(stats)))
# allow drawing in the margin
op <- par(xpd = NA)
on.exit(par(op), add = TRUE)
usr <- par("usr") # c(x1, x2, y1, y2)
x_rng <- diff(usr[1:2])
y_rng <- diff(usr[3:4])
# x positions (user sets; relative to right edge)
x_left <- usr[2] + x_frac[1] * x_rng
x_right <- usr[2] + x_frac[2] * x_rng
# compute line height in user units (based on current cex)
old_cex <- par("cex")
old_family <- par("family")
par(cex = cex, family = family)
on.exit(par(cex = old_cex, family = old_family), add = TRUE)
line_h <- strheight("Mg") # approx line height in user units
gap_h <- 0.35 * line_h # gap between lines
inner_h <- length(stats) * line_h + (length(stats) - 1) * gap_h
# box padding in user units
py <- pad_y * y_rng
px <- pad_x * x_rng
box_h <- inner_h + 2 * py
# center vertically in plot area
y_mid <- mean(usr[3:4], na.rm = TRUE)
y_bot <- y_mid - box_h / 2
y_top <- y_mid + box_h / 2
# draw box
rect(x_left, y_bot, x_right, y_top, col = box_col, border = border_col)
# y positions for each row (top to bottom)
y_start <- y_top - py - line_h / 2
y_pos <- y_start - (seq_along(stats) - 1) * (line_h + gap_h)
# left column: names
text(x_left + px, y_pos,
labels = names(stats),
adj = c(0, 0.5))
# right column: values
text(x_right - px, y_pos,
labels = unname(stats),
adj = c(1, 0.5))
}
get_scale <- function(vct, pct = .1, h0 = NULL) {
rng <- range(vct, na.rm = TRUE)
if (!is.null(h0)) {
if (h0 < rng[1]) {
rng[1] <- h0
}
if (h0 > rng[2]) {
rng[2] <- h0
}
}
if (rng[1]>0){
mn <- rng[1] * (1 - pct)
} else{
mn <- rng[1] * (1 + pct)
}
if (rng[2]>0){
mx <- rng[2] * (1 + pct)
} else{
mx <- rng[2] * (1 - pct)
}
ret <- c(mn, mx)
return(ret)
}
get_zero_scale <- function(vct, pct = .1) {
rng <- range(vct, na.rm = TRUE)
mn <- abs(rng[1] * (1 + pct))
mx <- abs(rng[2] * (1 + pct))
tmx <- max(mn, mx, na.rm = TRUE)
ret <- c(-tmx, tmx)
return(ret)
}
# This function was mostly generated by AI and I don't fully understand it.
# But it does seem to work.
get_reg_xlim <- function(brks, mu, sdx,
pad_frac = 0.05,
normal_tail_frac = 0.001, # normal at edges <= this * peak
max_sigma = 6) {
# Get range of calculated breaks
rng <- range(brks, na.rm = TRUE)
# Add adding on bin range
pad0 <- pad_frac * diff(rng)
xlim <- c(rng[1] - pad0, rng[2] + pad0)
# Force symmetry around zero (expand, never shrink)
m <- max(abs(xlim), na.rm = TRUE)
xlim <- c(-m, m)
if (!is.finite(sdx) || sdx <= 0) return(xlim)
# Calculate peak of normal curve
peak <- stats::dnorm(mu, mean = mu, sd = sdx)
target <- normal_tail_frac * peak
# For symmetric limits, the "worst" edge is the one closer to mu.
# Ensure BOTH edges are far enough into the tails.
k <- 0
while ((stats::dnorm(xlim[1], mu, sdx) > target ||
stats::dnorm(xlim[2], mu, sdx) > target) &&
k < max_sigma) {
xlim <- xlim + c(-0.5 * sdx, 0.5 * sdx) # symmetric expansion
k <- k + 0.5
}
ret <- xlim - c(-0.5 * sdx, 0.5 * sdx) # xlim[seq(2, length(xlim) - 1)]
return(ret)
}
remove_na_rows <- function(dat, dvr, ivrs) {
# Determine rows with NAs
dvna <- is.na(dat[[dvr]])
ivna <- rep(FALSE, nrow(dat))
for (var in ivrs) {
ivna <- ivna | is.na(dat[[var]])
}
cvna <- dvna | ivna
# Strip any such rows
ret <- dat[!cvna, ]
return(ret)
}
Try the procs package in your browser
Any scripts or data that you put into this service are public.
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.