Nothing
R/ttestplots.R
In procs: Recreates Some 'SAS®' Procedures in 'R'
Defines functions
get_sas_bins
boxplot_stats2
boxplot_stats1
render_agreement
render_profiles
render_tqqplot2
render_tqqplot1
render_interval2
render_interval1
render_boxplot2
render_boxplot1
render_histogram2
render_histogram1
render_summary2
render_summary1
ttestplot
Documented in
ttestplot
# TTestplots Main Function --------------------------------------------------
# AGREEMENT, ALL, BOX, HISTOGRAM, INTERVAL, NONE, PROFILES, QQ, SUMMARY
# agreement, boxplot, histogram, interval, profiles, qqplot, summary
#' @title Request T-Test Plots
#' @description A function to request T-test plots on a call to
#' \code{\link{proc_ttest}}. The function allows you to specify the type of
#' T-Test plots to produce. By default, it produces plots appropriate
#' to the analysis being performed.
#' You may also request specific plots
#' by passing a vector of plot names to the \code{type} parameter.
#' @details Plots produced by \code{\link{proc_ttest}} can change depending on
#' the type of analysis being performed.
#'
#' For one sample T-Tests, the function creates single variable
#' charts. The "summary" and
#' "qqplot" charts are generated by default. For this type of analysis,
#' it is helpful to set the "showh0" parameter
#' to TRUE so you can see the null hypothesis on the chart.
#'
#' For two sample T-Tests, the function produces two variable charts. For example,
#' on a two sample T-Test, the "histogram" chart will show histograms for
#' each class value on the same plot. Likewise, the "boxplot" and "qqplot" are
#' separated for each variable. The "interval" plot is split into "Satterthwaite" and
#' "Pooled" intervals. Note that the "h0" option also produces a null
#' hypothesis line for the "interval" plot, if requested.
#'
#' For a paired analysis, the "profiles" and "agreement" charts are available
#' in addition to the "histogram", "boxplot", "interval", and "qqplot" charts.
#' Paired analysis also produces single variable charts.
#'
#' To see all plots available for the requested analysis, pass the keyword
#' "all" to the \code{type} parameter.
#'
#' Any requested
#' 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_freq}}, and pass
#' the entire list to \code{\link{proc_print}}.
#' @section Plots:
#' The \code{type} parameter allows you to request several types of T-Test
#' plots. Below are the types of plots that are supported. The list shows
#' the plot type keyword needed to request the plot, and a brief description:
#' \itemize{
#' \item{\strong{agreement}: An agreement plot for the input data. Only
#' available for paired comparisons.
#' }
#' \item{\strong{boxplot/box}: Displays a box and whisker plot for group comparisons,
#' including confidence intervals (if appropriate).
#' }
#' \item{\strong{histogram/hist}: A histogram with normal curve and kernel density
#' overlays.
#' }
#' \item{\strong{interval}: Visualizes the confidence intervals for means.
#' }
#' \item{\strong{profiles}: A line plot mapping responses
#' between analysis variables. Available only for paired analysis.
#' }
#' \item{\strong{qqplot/qq}: A quantile-quantile plot used to assess the assumption
#' of normality.
#' }
#' \item{\strong{summary}: Combines the histogram and boxplot charts
#' onto the same panel. For two-sample T-Tests, both the histogram and boxplots
#' are separated for each class value.
#' }
#' }
#' The above plots may be requested as a vector of keywords to the \code{plots}
#' parameter on \code{\link{proc_ttest}}, or as vector of values
#' to the \code{type} parameter on \code{ttestplots}.
#'
#' Note that, when passed as a vector of keywords, only the requested plots
#' will be produced. That is to say, the "only" keyword in SAS is implied
#' at all times for \code{\link{proc_ttest}}.
#'
#' @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_ttest} 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 "agreement", "boxplot", "histogram",
#' "interval", "profiles", "qqplot", and "summary". The default value is "default", which
#' will produce default plots appropriate for the analysis. You may use the
#' \code{\link[common]{v}} function to pass the plot keywords unquoted.
#' @param panel Whether or not to display the summary plot combined into
#' in a single panel. Default is TRUE. A value of FALSE will create
#' individual "histogram" and "boxplot" charts instead. This parameter is
#' equivalent to the "unpack" keyword in SAS.
#' @param showh0 Whether or not to show the null hypothesis on the chart.
#' Default is FALSE. If this parameter is TRUE, a vertical line will be created
#' on some charts to show the value specified for the "h0=" option on
#' \code{\link{proc_ttest}}.
#' @param label Whether or not to label outlier values. Valid values
#' are TRUE or FALSE. Default is TRUE. If TRUE, this option will assign
#' labels to outlier values on the "summary" and "boxplot" charts.
#' @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_ttest()]
#' @examples
#' library(procs)
#'
#' # Turn off printing for CRAN checks
#' options("procs.print" = FALSE)
#'
#' # View sample data
#' sleep
#' # extra group ID
#' # 1 0.7 1 1
#' # 2 -1.6 1 2
#' # 3 -0.2 1 3
#' # 4 -1.2 1 4
#' # 5 -0.1 1 5
#' # 6 3.4 1 6
#' # 7 3.7 1 7
#' # 8 0.8 1 8
#' # 9 0.0 1 9
#' # 10 2.0 1 10
#' # 11 1.9 2 1
#' # 12 0.8 2 2
#' # 13 1.1 2 3
#' # 14 0.1 2 4
#' # 15 -0.1 2 5
#' # 16 4.4 2 6
#' # 17 5.5 2 7
#' # 18 1.6 2 8
#' # 19 4.6 2 9
#' # 20 3.4 2 10
#'
#'
#' # Example 1: plots = TRUE
#' proc_ttest(sleep,
#' var = "extra",
#' class = "group",
#' plots = TRUE)
#'
#' # Example 2: plots = "all"
#' proc_ttest(sleep,
#' var = "extra",
#' class = "group",
#' plots = "all")
#'
#' # Example 3: plots = {vector}
#' proc_ttest(sleep,
#' var = "extra",
#' class = "group",
#' plots = c("histogram", "boxplot", "interval"))
#'
#' # Example 4: plots = ttestplot()
#' proc_ttest(sleep,
#' var = "extra",
#' class = "group",
#' plots = ttestplot("summary", panel = FALSE))
#'
#' @export
ttestplot <- function(type = "default", panel = TRUE, showh0 = FALSE, label = TRUE, id = NULL) {
# Non-standard evaluation
type <- resolve_arg(type)
# agreement, boxplot, histogram, interval, profiles, qqplot, summary
# Parameter Checks
vldvals <- c("agreement", "boxplot", "box", "histogram", 'hist', 'interval', 'profiles',
"qqplot", "qq", "summary", "default", "all")
if (any(!type %in% vldvals)) {
ivd <- type[!type %in% vldvals]
stop(paste0("Parameter value for 'type' invalid: ", paste0("'", ivd, "'", collapse = ", "),
"\nValid values are: ",
"'agreement', 'boxplot', 'box', 'histogram', 'hist', 'interval', 'profiles', ",
"'qqplot', 'qq', 'summary', 'default', 'all'."
))
}
if (!panel %in% c(TRUE, FALSE)) {
stop("Parameter value for 'panel' invalid. Valid values are TRUE or FALSE.")
}
if (!showh0 %in% c(TRUE, FALSE)) {
stop("Parameter value for 'showh0' invalid. Valid values are TRUE or FALSE.")
}
# Create object
ret <- structure(list(), class = c("ttestplot", "list"))
ret$type <- tolower(type)
ret$panel <- panel
ret$showh0 <- showh0
ret$label <- label
ret$id <- id
ret$alph <- 0.95 # default for now
ret$h0 <- 0 # default for now
return(ret)
}
# Plot Rendering ----------------------------------------------------------
#' @noRd
render_ttestplot <- function (dat, var, plt, class, res) {
ret <- NULL
if ("ttestplot" %in% class(plt)) {
if (all(plt$type == "default")) {
if (!is.null(plt$varlbl)) {
plt$type <- c("summary", "profiles", "agreement", "qqplot")
} else {
plt$type <- c("summary", "qqplot")
}
}
if (all(plt$type == "all")) {
# agreement, boxplot, histogram, interval, profiles, qqplot, summary
if (!is.null(plt$varlbl)) {
plt$type <- c("summary", "histogram", "boxplot", "interval", "profiles", "agreement", "qqplot")
} else {
plt$type <- c("summary", "histogram", "boxplot", "interval", "qqplot")
}
}
typs <- plt$type
# "Unpack" logic
if ("summary" %in% typs & plt$panel == FALSE) {
# Find and replace
pos <- match("summary", typs)
vins <- c("histogram", "boxplot")
nv <- append(typs, vins, after = pos)
typs <- nv[-pos]
}
ret <- list()
# agreement, boxplot, histogram, interval, profiles, qqplot, summary
for (tp in typs) {
if (tp == "summary") {
if (is.null(class)) {
ret[["summary"]] <- render_summary1(dat, var, plt, res)
} else {
ret[["summary"]] <- render_summary2(dat, var, plt, class, res)
}
} else if (tp == "histogram" || tp == "hist") {
if (is.null(class)) {
ret[["histogram"]] <- render_histogram1(dat, var, plt, res)
} else {
ret[["histogram"]] <- render_histogram2(dat, var, plt, class, res)
}
} else if (tp == "boxplot" || tp == "box") {
if (is.null(class)) {
ret[["boxplot"]] <- render_boxplot1(dat, var, plt, res)
} else {
ret[["boxplot"]] <- render_boxplot2(dat, var, plt, class, res)
}
} else if (tp == "qqplot" || tp == "qq") {
if (is.null(class)) {
ret[["qqplot"]] <- render_tqqplot1(dat, var, plt, res)
} else {
ret[["qqplot"]] <- render_tqqplot2(dat, var, plt, class, res)
}
} else if (tp == "interval") {
if (is.null(class)) {
ret[["interval"]] <- render_interval1(dat, var, plt, res)
} else {
ret[["interval"]] <- render_interval2(dat, var, plt, res)
}
} else if (tp == "profiles") {
ret[["profiles"]] <- render_profiles(dat, var, plt, res)
} else if (tp == "agreement") {
ret[["agreement"]] <- render_agreement(dat, var, plt, res)
}
}
if (length(ret) == 0) {
ret <- NULL
}
}
return(ret)
}
# Summary Plot ------------------------------------------------------------
# unpack: TRUE or FALSE
#' @noRd
render_summary1 <- function(dat, var, plt, res) {
op <- par("mar")
of <- par("fig")
# 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
alph <- (1 - plt$alph) * 100 # Percentage
alpha <- plt$alph # Actual value
# 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 <- dat[[var]]
wgt <- if (!is.null(plt$weight)) dat[[plt$weight]] else NULL
freq <- if (!is.null(plt$freq)) dat[[plt$freq]] else NULL
if (!is.null(freq)) {
freq <- floor(freq)
idx <- !is.na(rdt) & !is.na(freq) & freq > 0
if (!is.null(wgt))
idx <- idx & !is.na(wgt) & wgt > 0
rdt_hist <- rep(rdt[idx], times = freq[idx])
wgt_hist <- if (!is.null(wgt)) rep(wgt[idx], times = freq[idx]) else NULL
} else {
rdt_hist <- rdt
wgt_hist <- wgt
}
# dt <- dat[[ivr]]
# Assign labels
xlbl <- var
tlbl <- paste(" ", var)
if (!is.null(plt$varlbl)) {
xlbl <- "Difference"
tlbl <- paste0(" Difference: ", plt$varlbl)
}
#**********************************
#* Histogram
#**********************************
# Set margins
par(mar = c(0, 5, 3, .75) + 0.1,
fig = c(0, 1, .3, 1))
# Calculate stats
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
clm_tbl <- res[[grep("ConfLimits$", names(res))[length(grep("ConfLimits$", names(res)))]]]
n <- stat_tbl$N[1]
# mu <- stat_tbl$MEAN[1]
# sdx <- stat_tbl$STD[1]
mu <- mean(rdt, na.rm = TRUE)
sdx <- sd(rdt, na.rm = TRUE)
ci <- c(clm_tbl$LCLM[1], clm_tbl$UCLM[1])
# Calculate breaks - Closer to SAS algorithm
brks <- pretty(range(rdt, na.rm = TRUE), n = nclass.Sturges(rdt),
min.n = 1, high.u.bias = 3)
# Calculate breaks
# brks <- get_sas_bins(rdt_hist)
# Calculate breaks and y scale.
h <- hist(rdt_hist,
breaks = brks,
right = FALSE,
include.lowest = TRUE,
plot = FALSE)
h$counts <- h$counts / sum(h$counts, na.rm = TRUE) * 100
# Use calculated breaks to get scale (1.5 bin-widths)
w <- diff(h$breaks)[1] # bin width
xmin <- min(h$breaks, na.rm = TRUE) - 1.5 * w
xmax <- max(h$breaks, na.rm = TRUE) + 1.5 * w
scl <- c(xmin, xmax)
# Consider h0 if needed
if (plt$showh0) {
if (plt$h0 < scl[1]) {
scl[1] <- plt$h0
}
if (plt$h0 > scl[2]) {
scl[2] <- plt$h0
}
}
# Use calculated break to get normal curve
curve_lo <- min(h$breaks, na.rm = TRUE) - w
curve_hi <- max(h$breaks, na.rm = TRUE) + w
grid <- seq(curve_lo, curve_hi, length.out = 300)
y_norm_percent <- 100 * w * dnorm(grid, mean = mu, sd = sdx)
y_mx <- max(max(y_norm_percent, na.rm = TRUE), max(h$counts, na.rm = TRUE), na.rm = TRUE)
# Initial chart to draw vertical lines
hist(Inf,
breaks = brks,
right = FALSE,
include.lowest = TRUE,
main = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.025),
plot = TRUE,
axes = FALSE)
# Get axis tick marks
aval <- pretty(scl, n = 3)
# Orientation lines
abline(v = aval, col = "grey90", lwd = 1)
# Create plot using histogram values
plot(h,
col = "#CAD5E5", #"grey85",
border = "grey20",
main = "",
xlab = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.025),
axes = FALSE,
add = TRUE)
# Title
mtext(paste0("Distribution of", tlbl), side = 3,
line = par("mar")[3] - 1.5,
font = 2, cex = 1.25)
# Subtitle
mtext(paste0("With ", alph, "% Confidence Interval for Mean"), side = 3,
line = par("mar")[3] - 2.5,
font = 1)
# Y Label
mtext("Percent", side = 2, line = par("mar")[2] - 2)
# Normal curve overlay (scaled to percent)
lines(grid, y_norm_percent, col = "steelblue4", lwd = 2)
# Kernel density overlay (scaled to percent)
# Bandwidth (Silverman / SAS) - SAS uses the distinct (non-freq-expanded)
# observation count for the bandwidth. Kish's effective n when weighted.
if (!is.null(wgt_hist)&is.null(freq)) {
n_eff <- (sum(wgt_hist, na.rm = TRUE))^2 / sum(wgt_hist^2, na.rm = TRUE)
} else {
n_eff <- length(rdt)
}
bw <- 1.06 * sdx * n_eff^(-1/5)
# Kernel density estimate
kernel_lo <- min(h$breaks, na.rm = TRUE) - 1.25 * w
kernel_hi <- max(h$breaks, na.rm = TRUE) + 1.25 * w
dens <- density(rdt_hist, kernel = "gaussian", bw = bw, n = 512,
from = kernel_lo, to = kernel_hi, na.rm = TRUE)
# Scale KDE to Percent axis
y_kernel_percent <- 100 * w * dens$y
# Overlay density on histogram
lines(dens$x, y_kernel_percent, col = "orangered2", lwd = 2)
# Add custom Y axis
axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Frame
box(col = "grey70", lwd = 1)
box("outer", 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)
)
# Custom Legend border
rect(mpos$rect$left, mpos$rect$top - (mpos$rect$h - .5),
mpos$rect$left + mpos$rect$w, mpos$rect$top,
border = "grey80"
)
#**********************************
#* Boxplot
#**********************************
par(mar = c(4, 5, 0, .75) + 0.1,
fig = c(0, 1, 0, .3), new = TRUE)
# Get data
dt <- dat[[var]]
# Get scales
xscl <- get_scale(dt, .05) # Not used?
# Calculations
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
clm_tbl <- res[[grep("ConfLimits$", names(res))[length(grep("ConfLimits$", names(res)))]]]
n <- stat_tbl$N[1]
mu <- mean(dt, na.rm = TRUE)
sdx <- stat_tbl$STD[1]
ci <- c(clm_tbl$LCLM[1], clm_tbl$UCLM[1])
## Draw empty plot first (for layering)
plot(dt, rep(1, length(dt)),
type = "n",
xlab = "",
ylab = "",
yaxt = "n",
xlim = scl,
main = "",
axes = FALSE) #With 95% Confidence Interval for Mean")
# X Label
mtext(xlbl, side = 1,
line = par("mar")[1] - 2,
font = 1)
## Shaded CI band and make sure infinite values are out of boundry
usr <- par("usr")
if (is.infinite(ci[1])) ci[1] <- usr[1] - abs(usr[1])*10
if (is.infinite(ci[2])) ci[2] <- usr[2] + abs(usr[2])*10
rect(ci[1], usr[3], ci[2], usr[4],
col = "#B3D2D0",
border = NA)
# Draw axis (SAS-style: fewer, rounder ticks)
aval <- axis(side = 1, at = pretty(scl, n = 3), las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Orientation lines
abline(v = aval, col = "grey90", lwd = 1)
# Show h0
if (plt$showh0) {
abline(v = plt$h0, col = "grey70", lwd = 2)
}
# Get boxplot stats
bp <- boxplot_stats1(dt, if (!is.null(plt$freq)) dat[[plt$freq]] else NULL)
# Boxplot (horizontal)
bxp(bp,
horizontal = TRUE,
add = TRUE,
boxfill = "#CAD5E5",
border = "grey40",
lty = 1,
whiskcol = "#05379B",
staplecol = "#05379B",
medcol = "#05379B",
medlwd = 1,
boxwex = .65,
outline = FALSE,
axes = FALSE) # Already drawn above
## Mean diamond
points(mu, 1,
pch = 5,
col = "#05379B",
lwd = 1,
cex = 1.25)
# Display outliers
if (length(bp$out) > 0) {
points(bp$out, rep(1, length(bp$out)),
pch = 1,
col = "grey20",
lwd = 1.6,
cex = 1.25)
# Add labels
if (plt$label) {
if (is.null(plt$id)) {
olbls <- bp$onm
} else {
olbls <- dat[[plt$id]][as.integer(bp$onm)]
}
if (length(olbls) > 0) {
text(bp$out, rep(1, length(bp$out)),
labels = olbls,
cex = .9,
pos = 1)
}
}
}
# Legend location - Move to left side if needed
lgnd <- "topright"
pt <- (bp$stats[4] - xscl[1])/diff(xscl)
if (pt > .75) {
lgnd <- "topleft"
}
# Create legend
if (plt$showh0) {
mpos <- legend(lgnd,
legend = c(paste0(alph, "% Confidence"), "Null Value"),
pch = c(15, NA),
fill = c("#B3D2D0", NA),
col = c(NA, "grey70"),
lwd = c(NA, 2),
border = c("grey60", NA),
inset = c(.01, .02),
x.intersp = .5,
y.intersp = c(0.5, 1.5, .5),
cex = .8,
bty = "n")
rect(mpos$rect$left, mpos$rect$top - (mpos$rect$h - .1),
mpos$rect$left + mpos$rect$w, mpos$rect$top,
border = "grey80"
)
} else {
legend(lgnd,
legend = paste0(alph, "% Confidence"),
fill = "#B3D2D0",
border = "grey60",
box.col = "grey80",
inset = c(.01, .02),
x.intersp = .5,
y.intersp = c(0.5, 1.8),
cex = .9,
bty = "o")
}
# Frame
box(col = "grey70", lwd = 1)
#**********************************
#* Clean up
#**********************************
# Restore margins
par(mar = op, fig = of)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
#' @noRd
render_summary2 <- function(dat, var, plt, class, res) {
op <- par("mar")
om <- par("oma")
of <- par("fig")
# 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
alph <- (1 - plt$alph) * 100 # Percentage
alpha <- plt$alph # Actual value
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Convert to data frame just in case it is a tibble
dat <- as.data.frame(dat)
if (!is.null(plt$id)) {
lbls <- dat[[plt$id]]
} else {
lbls <- NULL
}
# Add rownumbers
rownames(dat) <- seq(1, nrow(dat))
dat <- sort(dat, by = class)
# Get analysis variables
cvls <- unique(dat[[class]])
cvls <- cvls[!is.na(cvls)]
vl1 <- cvls[1]
vl2 <- cvls[2]
# Prepare data
dt1 <- dat[dat[[class]] == vl1, var]
dt2 <- dat[dat[[class]] == vl2, var]
wgt1 <- if (!is.null(plt$weight)) dat[dat[[class]] == vl1, plt$weight] else NULL
wgt2 <- if (!is.null(plt$weight)) dat[dat[[class]] == vl2, plt$weight] else NULL
freq1 <- if (!is.null(plt$freq)) dat[dat[[class]] == vl1, plt$freq] else NULL
freq2 <- if (!is.null(plt$freq)) dat[dat[[class]] == vl2, plt$freq] else NULL
# Apply freq-expansion to mirror SAS hist behavior
if (!is.null(freq1)) {
freq1 <- floor(freq1)
idx1 <- !is.na(dt1) & !is.na(freq1) & freq1 > 0
if (!is.null(wgt1))
idx1 <- idx1 & !is.na(wgt1) & wgt1 > 0
dt1_hist <- rep(dt1[idx1], times = freq1[idx1])
wgt1_hist <- if (!is.null(wgt1)) rep(wgt1[idx1], times = freq1[idx1]) else NULL
} else {
dt1_hist <- dt1
wgt1_hist <- wgt1
}
if (!is.null(freq2)) {
freq2 <- floor(freq2)
idx2 <- !is.na(dt2) & !is.na(freq2) & freq2 > 0
if (!is.null(wgt2))
idx2 <- idx2 & !is.na(wgt2) & wgt2 > 0
dt2_hist <- rep(dt2[idx2], times = freq2[idx2])
wgt2_hist <- if (!is.null(wgt2)) rep(wgt2[idx2], times = freq2[idx2]) else NULL
} else {
dt2_hist <- dt2
wgt2_hist <- wgt2
}
# Get number of lines needed to show labels
minlns <- get_line_count(as.character(cvls)) + 1
lmg <- 3
if (minlns > 2) {
lmg <- minlns + 1
}
# Set margins
par(oma = c(5, 1, 2, .75) + 0.1,
mar = c(.1, lmg, 0, 0) + 0.1)
#******************************
#* Calculate scales
#******************************
# Overall scale
scl <- range(dat[[var]], na.rm = TRUE)
scl <- c(scl[1] * .875, scl[2] * 1.125)
# Calculate breaks - Histogram 1
brks1 <- pretty(range(dt1, na.rm = TRUE), n = nclass.Sturges(dt1),
min.n = 1, high.u.bias = 3)
# Calculate breaks - Histogram 1
# brks1 <- get_sas_bins(dt1_hist)
# Calculate breaks - Histogram 2
brks2 <- pretty(range(dt2, na.rm = TRUE), n = nclass.Sturges(dt2),
min.n = 1, high.u.bias = 3)
# Calculate breaks - Histogram 2
# brks2 <- get_sas_bins(dt2_hist)
# Scale for Histogram 1 (SAS-style: 1.5 bin-widths of padding)
w1 <- diff(brks1)[1] # bin width
xmin1 <- min(brks1, na.rm = TRUE) - 1.5 * w1
xmax1 <- max(brks1, na.rm = TRUE) + 1.5 * w1
scl1 <- c(xmin1, xmax1)
# Scale for Histogram 2 (SAS-style: 1.5 bin-widths of padding)
w2 <- diff(brks2)[1] # bin width
xmin2 <- min(brks2, na.rm = TRUE) - 1.5 * w2
xmax2 <- max(brks2, na.rm = TRUE) + 1.5 * w2
scl2 <- c(xmin2, xmax2)
# Adjust overall scale if needed by histogram 1
if (scl1[1] < scl[1]) {
scl[1] <- scl1[1]
}
if (scl1[2] > scl[2]) {
scl[2] <- scl1[2]
}
# Adjust overall scale if needed by histogram 2
if (scl2[1] < scl[1]) {
scl[1] <- scl2[1]
}
if (scl2[2] > scl[2]) {
scl[2] <- scl2[2]
}
#******************************
#* Histogram 1
#******************************
par(fig = c(0, 1, .6, 1))
# Calculate stats
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
clm_tbl <- res[[grep("ConfLimits$", names(res))[length(grep("ConfLimits$", names(res)))]]]
n <- stat_tbl$N[as.character(stat_tbl$CLASS) == as.character(vl1)]
mu <- mean(dt1, na.rm = TRUE)
sdx <- sd(dt1, na.rm = TRUE)
# Histogram (Percent scale) - SAS uses left-closed bins [a, b)
h <- hist(dt1_hist,
breaks = brks1,
right = FALSE,
include.lowest = TRUE,
plot = FALSE)
h$counts <- h$counts / sum(h$counts, na.rm = TRUE) * 100
# Use calculated break to get normal curve
grid <- seq(scl1[1] , scl1[2], length.out = 300)
y_norm_percent <- 100 * w1 * dnorm(grid, mean = mu, sd = sdx)
y_mx <- max(max(y_norm_percent, na.rm = TRUE), max(h$counts, na.rm = TRUE), na.rm = TRUE)
# Initial chart to draw vertical orientation lines underneath bars
hist(Inf,
breaks = brks1,
right = FALSE,
include.lowest = TRUE,
main = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.025),
plot = TRUE,
axes = FALSE)
# Get axis tick marks and draw orientation lines
aval <- pretty(scl, n = 3)
abline(v = aval, col = "grey90", lwd = 1)
# Create plot using histogram values
plot(h,
col = "#CAD5E5", #"grey85",
border = "grey20",
main = "",
xlab = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.025), # max(h$counts) * 1.05),
axes = FALSE,
add = TRUE)
# Normal curve overlay (scaled to percent)
lines(grid, y_norm_percent, col = "steelblue4", lwd = 2)
# Kernel density overlay (scaled to percent)
# Bandwidth (Silverman / SAS) - SAS uses distinct (non-freq-expanded) n;
# Kish's effective n when only weight is supplied.
if (!is.null(wgt1_hist) & is.null(freq1)) {
n_eff <- (sum(wgt1_hist, na.rm = TRUE))^2 / sum(wgt1_hist^2, na.rm = TRUE)
} else {
n_eff <- length(dt1)
}
bw <- 1.06 * sdx * n_eff^(-1/5)
# Kernel density estimate
dens <- density(dt1_hist, kernel = "gaussian", bw = bw, n = 512, na.rm = TRUE)
# Scale KDE to Percent axis
y_kernel_percent <- 100 * w1 * dens$y
yscl <- seq(50, length(y_kernel_percent) - 25)
# Overlay on histogram
lines(dens$x, y_kernel_percent, col = "orangered2", lwd = 2)
# Add custom axes
axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Add y label
mtext("Percent", side = 2, line = par()$mar[2] - .75)
# Add class label
legend("topleft", legend = vl1, bty = "n",
x.intersp = 0, # Spacing between group label and box
y.intersp = 0, # Spacing between content and borders
cex = .9
)
# Frame
box(col = "grey70", lwd = 1)
#******************************
#* Histogram 2
#******************************
par(fig = c(0, 1, .2, .6), new = TRUE)
# Calculate stats
n <- stat_tbl$N[as.character(stat_tbl$CLASS) == as.character(vl2)]
mu <- mean(dt2, na.rm = TRUE)
sdx <- sd(dt2, na.rm = TRUE)
# Histogram (Percent scale) - SAS uses left-closed bins [a, b)
h <- hist(dt2_hist,
breaks = brks2,
right = FALSE,
include.lowest = TRUE,
plot = FALSE)
h$counts <- h$counts / sum(h$counts, na.rm = TRUE) * 100
# Use calculated break to get normal curve
grid <- seq(scl2[1] , scl2[2], length.out = 300)
y_norm_percent <- 100 * w2 * dnorm(grid, mean = mu, sd = sdx)
y_mx <- max(max(y_norm_percent, na.rm = TRUE),
max(h$counts, na.rm = TRUE), na.rm = TRUE)
# Initial chart to draw vertical orientation lines underneath bars
hist(Inf,
breaks = brks2,
right = FALSE,
include.lowest = TRUE,
main = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.025),
plot = TRUE,
axes = FALSE)
# Get axis tick marks and draw orientation lines
aval <- pretty(scl, n = 3)
abline(v = aval, col = "grey90", lwd = 1)
# Create plot using histogram values
plot(h,
col = "#CAD5E5", #"grey85",
border = "grey20",
main = "",
xlab = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.025), # max(h$counts) * 1.05),
axes = FALSE,
add = TRUE)
# Normal curve overlay (scaled to percent)
lines(grid, y_norm_percent, col = "steelblue4", lwd = 2)
# Kernel density overlay (scaled to percent)
# Bandwidth (Silverman / SAS) - distinct n; Kish's effective n when weight only
if (!is.null(wgt2_hist) & is.null(freq2)) {
n_eff <- (sum(wgt2_hist, na.rm = TRUE))^2 / sum(wgt2_hist^2, na.rm = TRUE)
} else {
n_eff <- length(dt2)
}
bw <- 1.06 * sdx * n_eff^(-1/5) # 1.06
# Kernel density estimate
dens <- density(dt2_hist, kernel = "gaussian", bw = bw, n = 512, na.rm = TRUE)
# Scale KDE to Percent axis
y_kernel_percent <- 100 * w2 * dens$y
# Overlay on histogram
lines(dens$x, y_kernel_percent, col = "orangered2", lwd = 2)
# Add custom axes
axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Add y label
mtext("Percent", side = 2, line = par()$mar[2] - .75)
# Add class label
legend("topleft", legend = vl2, bty = "n",
x.intersp = 0, # Spacing between group label and box
y.intersp = 0, # Spacing between content and borders
cex = .9
)
# Frame
box(col = "grey70", lwd = 1)
#**********************************
#* Boxplot
#**********************************
par(fig = c(0, 1, 0, .2), new = TRUE)
# Prepare data
dt <- dat[[var]]
rwnms <- rownames(dat)
nms1 <- rwnms[dat[[class]] == vl1]
nms2 <- rwnms[dat[[class]] == vl2]
# Get x scale
# xscl <- get_scale(dt, .05)
# Calculate
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
clm_tbl <- res[[grep("ConfLimits$", names(res))[length(grep("ConfLimits$", names(res)))]]]
n1 <- stat_tbl$N[as.character(stat_tbl$CLASS) == as.character(vl1)]
mu1 <- stat_tbl$MEAN[as.character(stat_tbl$CLASS) == as.character(vl1)]
sdx1 <- stat_tbl$STD[as.character(stat_tbl$CLASS) == as.character(vl1)]
n2 <- stat_tbl$N[as.character(stat_tbl$CLASS) == as.character(vl2)]
mu2 <- stat_tbl$MEAN[as.character(stat_tbl$CLASS) == as.character(vl2)]
sdx2 <- stat_tbl$STD[as.character(stat_tbl$CLASS) == as.character(vl2)]
# Draw empty plot first, so we can put ablines() behind the boxes
plot(dt, rep(1, length(dt)),
ylim = c(.5, 2.5),
xlim = scl,
type = "n",
xlab = "",
ylab = "",
yaxt = "n",
main = "",
axes = FALSE) #With 95% Confidence Interval for Mean")
# Y Label
mtext(class, side = 2,
line = par()$mar[2] - .75,
font = 1)
# Draw X axis
aval <- axis(side = 1, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# # Draw Y axis
axis(side = 2, las = 1, col.ticks = "grey55", labels = cvls,
at = c(2, 1), mgp = c(3, .5, 0), tck = -0.015)
# Orientation lines
abline(v = aval, col = "grey90", lwd = 1)
# Get boxplot stats for 2 plots
freq1 <- if (!is.null(plt$freq)) dat[dat[[class]] == vl1, plt$freq] else NULL
freq2 <- if (!is.null(plt$freq)) dat[dat[[class]] == vl2, plt$freq] else NULL
bp <- boxplot_stats2(dt1, dt2, nms1, nms2, cvls, freq1, freq2)
# Boxplot (horizontal)
bxp(bp,
horizontal = TRUE,
add = TRUE,
boxfill = "#CAD5E5",
border = "grey40",
lty = 1,
whiskcol = "#05379B",
staplecol = "#05379B",
medcol = "#05379B",
medlwd = 1,
show.names = TRUE,
boxwex = .425,
outline = FALSE,
axes = FALSE) # Already drawn above
## Mean diamond
points(c(mu1, mu2), c(2, 1),
pch = 5,
col = "#05379B",
lwd = 1,
cex = 1.25)
# Outliers for first class value
if (length(bp$out1) > 0) {
points(bp$out1, rep(2, length(bp$out1)),
pch = 1,
col = "grey20",
lwd = 1.6,
cex = 1.25)
# Add labels
if (plt$label) {
if (is.null(plt$id)) {
olbls1 <- bp$onm1
} else {
olbls1 <- lbls[as.integer(bp$onm1)]
}
if (length(olbls1) > 0) {
text(bp$out1, rep(2, length(bp$out1)),
labels = olbls1,
cex = .9,
pos = 1)
}
}
}
# Outliers for second class value
if (length(bp$out2) > 0) {
# Plot outlier points
points(bp$out2, rep(1, length(bp$out2)),
pch = 1,
col = "grey20",
lwd = 1.6,
cex = 1.25)
# Add labels
if (plt$label) {
if (is.null(plt$id)) {
olbls2 <- bp$onm2
} else {
olbls2 <- lbls[as.integer(bp$onm2)]
}
if (length(olbls2) > 0) {
text(bp$out2, rep(1, length(bp$out2)),
labels = olbls2,
cex = .9,
pos = 3)
}
}
}
# Frame
box(col = "grey70", lwd = 1)
#**********************************
#* Clean up
#**********************************
# Title
mtext(paste0("Distribution of ", var), side = 3,
line = par("oma")[3] - 1.75,
font = 2, cex = 1.25, outer = TRUE)
# X Label
mtext(var, side = 1,
line = par("oma")[1] - 3.5,
font = 1, outer = TRUE)
box("outer", col = "grey70", lwd = 1)
# Legend
mpos <- legend("bottom",
legend = c("Normal", "Kernel"),
col = c("steelblue4", "orangered2"),
horiz = TRUE,
lwd = 2,
box.col = "grey80",
bty = "o",
cex = .9,
xpd = NA,
inset = c(0, -1),
x.intersp = .5,
y.intersp = c(0.3, .3) # c(0.5, 1.8)
)
# Restore margins
par(mar = op, fig = of, oma = om)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
# Individual Plots --------------------------------------------------------
# agreement, boxplot, histogram, interval, profiles, qqplot, summary
# Curves for normal and density most correct on this one.
#' @noRd
render_histogram1 <- function(dat, var, plt, res) {
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 <- dat[[var]]
# rdt <- rdt[!is.na(rdt)]
wgt <- if (!is.null(plt$weight)) dat[[plt$weight]] else NULL
freq <- if (!is.null(plt$freq)) dat[[plt$freq]] else NULL
if (!is.null(freq)) {
freq <- floor(freq)
idx <- !is.na(rdt) & !is.na(freq) & freq > 0
if (!is.null(wgt))
idx <- idx & !is.na(wgt) & wgt > 0
rdt_hist <- rep(rdt[idx], times = freq[idx])
wgt_hist <- if (!is.null(wgt)) rep(wgt[idx], times = freq[idx]) else NULL
} else {
rdt_hist <- rdt
wgt_hist <- wgt
}
# Set margins
par(mar = c(4, 5, 2, .75) + 0.1)
# Assign labels
xlbl <- var
tlbl <- paste(" ", var)
if (!is.null(plt$varlbl)) {
xlbl <- "Difference"
tlbl <- paste0(" Difference: ", plt$varlbl)
}
# Calculate stats
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
clm_tbl <- res[[grep("ConfLimits$", names(res))[length(grep("ConfLimits$", names(res)))]]]
n <- stat_tbl$N[1]
mu <- mean(rdt, na.rm = TRUE)
sdx <- sd(rdt, na.rm = TRUE)
ci <- c(clm_tbl$LCLM[1], clm_tbl$UCLM[1])
# Calculate breaks - Closer to SAS algorithm
brks <- pretty(range(rdt, na.rm = TRUE), n = nclass.Sturges(rdt),
min.n = 1, high.u.bias = 3)
# Calculate breaks
# brks <- get_sas_bins(rdt_hist)
# Histogram (Percent scale) - SAS uses left-closed bins [a, b)
h <- hist(rdt_hist,
breaks = brks,
right = FALSE,
include.lowest = TRUE,
plot = FALSE)
h$counts <- h$counts / sum(h$counts, na.rm = TRUE) * 100
# Use calculated breaks to get scale (1.5 bin-widths)
w <- diff(h$breaks)[1] # bin width
xmin <- min(h$breaks, na.rm = TRUE) - 1.5 * w
xmax <- max(h$breaks, na.rm = TRUE) + 1.5 * w
scl <- c(xmin, xmax)
curve_lo <- min(h$breaks, na.rm = TRUE) - w
curve_hi <- max(h$breaks, na.rm = TRUE) + w
kernel_lo <- min(h$breaks, na.rm = TRUE) - 1.4 * w
kernel_hi <- max(h$breaks, na.rm = TRUE) + 1.4 * w
# Bandwidth (Silverman / SAS) - distinct n; Kish's effective n when weight only
if (!is.null(wgt_hist) & is.null(freq)) {
n_eff <- (sum(wgt_hist, na.rm = TRUE))^2 / sum(wgt_hist^2, na.rm = TRUE)
} else {
n_eff <- length(rdt)
}
bw <- 1.06 * sdx * n_eff^(-1/5)
# Kernel density estimate
dens <- density(rdt_hist, kernel = "gaussian", bw = bw, n = 512,
from = kernel_lo, to = kernel_hi, na.rm = TRUE)
# Use calculated break to get normal curve
grid <- seq(curve_lo, curve_hi, length.out = 300)
y_norm_percent <- 100 * w * dnorm(grid, mean = mu, sd = sdx)
y_kernel_percent <- 100 * w * dens$y
y_mx <- max(max(y_norm_percent, na.rm = TRUE), max(h$counts, na.rm = TRUE), na.rm = TRUE)
# Initial chart to draw vertical orientation lines underneath bars
hist(Inf, # Putting values here was drawing a single bar, and messing up the plot
breaks = brks,
right = FALSE,
include.lowest = TRUE,
main = paste0("Distribution of", tlbl),
xlab = "",
ylab = "Percent",
xlim = scl,
ylim = c(0, y_mx * 1.025),
plot = TRUE,
axes = FALSE)
# Get axis tick marks and draw orientation lines
aval <- pretty(scl, n = 3)
abline(v = aval, col = "grey90", lwd = 1)
# Create plot using histogram values
plot(h,
col = "#CAD5E5", #"grey85",
border = "grey20",
main = "",
xlab = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.025),
axes = FALSE,
add = TRUE)
# Draw normal curve line
lines(grid, y_norm_percent, col = "steelblue4", lwd = 2)
# Overlay on histogram
lines(dens$x, y_kernel_percent, col = "orangered2", lwd = 2)
# Add x axis label
mtext(xlbl, side = 1, line = par("mar")[1] - 2)
# Add custom axes
axis(side = 1, at = pretty(scl, n = 3), col.ticks = "grey55",
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_histogram2 <- function(dat, var, plt, class, res) {
op <- par("mar")
om <- par("oma")
# 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
# Convert to data frame
dat <- as.data.frame(dat)
dat <- sort(dat, by = class)
# Get analysis variables
cvls <- unique(dat[[class]])
vl1 <- cvls[1]
vl2 <- cvls[2]
# Prepare data
dt1 <- dat[dat[[class]] == vl1, var]
dt2 <- dat[dat[[class]] == vl2, var]
wgt1 <- if (!is.null(plt$weight)) dat[dat[[class]] == vl1, plt$weight] else NULL
wgt2 <- if (!is.null(plt$weight)) dat[dat[[class]] == vl2, plt$weight] else NULL
freq1 <- if (!is.null(plt$freq)) dat[dat[[class]] == vl1, plt$freq] else NULL
freq2 <- if (!is.null(plt$freq)) dat[dat[[class]] == vl2, plt$freq] else NULL
# Apply freq-expansion to mirror SAS hist behavior
if (!is.null(freq1)) {
freq1 <- floor(freq1)
idx1 <- !is.na(dt1) & !is.na(freq1) & freq1 > 0
if (!is.null(wgt1))
idx1 <- idx1 & !is.na(wgt1) & wgt1 > 0
dt1_hist <- rep(dt1[idx1], times = freq1[idx1])
wgt1_hist <- if (!is.null(wgt1)) rep(wgt1[idx1], times = freq1[idx1]) else NULL
} else {
dt1_hist <- dt1
wgt1_hist <- wgt1
}
if (!is.null(freq2)) {
freq2 <- floor(freq2)
idx2 <- !is.na(dt2) & !is.na(freq2) & freq2 > 0
if (!is.null(wgt2))
idx2 <- idx2 & !is.na(wgt2) & wgt2 > 0
dt2_hist <- rep(dt2[idx2], times = freq2[idx2])
wgt2_hist <- if (!is.null(wgt2)) rep(wgt2[idx2], times = freq2[idx2]) else NULL
} else {
dt2_hist <- dt2
wgt2_hist <- wgt2
}
# Output to image file
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Set margins
par(oma = c(5, 1, 2, .75) + 0.1,
mar = c(0, 3, 0, 0) + 0.1,
mfrow = c(2, 1))
#******************************
#* Calculate scales
#******************************
# Overall scale
scl <- range(dat[[var]], na.rm = TRUE)
scl <- c(scl[1] * .875, scl[2] * 1.125)
# Calculate breaks - Histogram 1
brks1 <- pretty(range(dt1, na.rm = TRUE), n = nclass.Sturges(dt1),
min.n = 1, high.u.bias = 3)
# Calculate breaks - Histogram 2
brks2 <- pretty(range(dt2, na.rm = TRUE), n = nclass.Sturges(dt2),
min.n = 1, high.u.bias = 3)
# Calculate breaks - Histogram 1
# brks1 <- get_sas_bins(dt1_hist)
# Calculate breaks - Histogram 2
# brks2 <- get_sas_bins(dt2_hist)
# Scale for Histogram 1 (SAS-style: 1.5 bin-widths of padding)
w1 <- diff(brks1)[1] # bin width
xmin1 <- min(brks1, na.rm = TRUE) - 1.5 * w1
xmax1 <- max(brks1, na.rm = TRUE) + 1.5 * w1
scl1 <- c(xmin1, xmax1)
# Scale for Histogram 2 (SAS-style: 1.5 bin-widths of padding)
w2 <- diff(brks2)[1] # bin width
xmin2 <- min(brks2, na.rm = TRUE) - 1.5 * w2
xmax2 <- max(brks2, na.rm = TRUE) + 1.5 * w2
scl2 <- c(xmin2, xmax2)
# Adjust overall scale if needed by histogram 1
if (scl1[1] < scl[1]) {
scl[1] <- scl1[1]
}
if (scl1[2] > scl[2]) {
scl[2] <- scl1[2]
}
# Adjust overall scale if needed by histogram 2
if (scl2[1] < scl[1]) {
scl[1] <- scl2[1]
}
if (scl2[2] > scl[2]) {
scl[2] <- scl2[2]
}
#******************************
#* Histogram 1
#******************************
# Calculate stats
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
n <- stat_tbl$N[as.character(stat_tbl$CLASS) == as.character(vl1)]
mu <- mean(dt1, na.rm = TRUE)
sdx <- sd(dt1, na.rm = TRUE)
# Histogram (Percent scale) - SAS uses left-closed bins [a, b)
h <- hist(dt1_hist,
breaks = brks1,
right = FALSE,
include.lowest = TRUE,
plot = FALSE)
h$counts <- h$counts / sum(h$counts, na.rm = TRUE) * 100
# Use calculated break to get normal curve
grid <- seq(scl1[1] , scl1[2], length.out = 300)
y_norm_percent <- 100 * w1 * dnorm(grid, mean = mu, sd = sdx)
y_mx <- max(max(y_norm_percent, na.rm = TRUE), max(h$counts, na.rm = TRUE), na.rm = TRUE)
# Initial chart to draw vertical orientation lines underneath bars
hist(Inf,
breaks = brks1,
right = FALSE,
include.lowest = TRUE,
main = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.025),
plot = TRUE,
axes = FALSE)
# Get axis tick marks and draw orientation lines
aval <- pretty(scl, n = 3)
abline(v = aval, col = "grey90", lwd = 1)
# Create plot using histogram values
plot(h,
col = "#CAD5E5", #"grey85",
border = "grey20",
main = "",
xlab = "",
ylab = "Percent",
xlim = scl,
ylim = c(0, y_mx * 1.025), # max(h$counts) * 1.05),
axes = FALSE,
add = TRUE)
# Normal curve overlay (scaled to percent)
lines(grid, y_norm_percent, col = "steelblue4", lwd = 2)
# Bandwidth (Silverman / SAS) - distinct n; Kish's effective n when weight only
if (!is.null(wgt1_hist) & is.null(freq1)) {
n_eff <- (sum(wgt1_hist, na.rm = TRUE))^2 / sum(wgt1_hist^2, na.rm = TRUE)
} else {
n_eff <- length(dt1)
}
bw <- 1.06 * sdx * n_eff^(-1/5)
# Kernel density estimate
dens <- density(dt1_hist, kernel = "gaussian", bw = bw, n = 512, na.rm = TRUE)
# Scale KDE to Percent axis
y_kernel_percent <- 100 * w1 * dens$y
yscl <- seq(50, length(y_kernel_percent) - 25)
# Overlay on histogram
lines(dens$x, y_kernel_percent, col = "orangered2", lwd = 2)
# Add custom axes
axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Add y label
mtext("Percent", side = 2, line = 2.5)
# Add class label
legend("topleft", legend = vl1, bty = "n",
x.intersp = 0, # Spacing between group label and box
y.intersp = 0, # Spacing between content and borders
cex = .9
)
# Frame
box(col = "grey70", lwd = 1)
#******************************
#* Histogram 2
#******************************
# Calculate stats
n <- stat_tbl$N[as.character(stat_tbl$CLASS) == as.character(vl2)]
mu <- mean(dt2, na.rm = TRUE)
sdx <- sd(dt2, na.rm = TRUE)
# Histogram (Percent scale) - SAS uses left-closed bins [a, b)
h <- hist(dt2_hist,
breaks = brks2,
right = FALSE,
include.lowest = TRUE,
plot = FALSE)
h$counts <- h$counts / sum(h$counts, na.rm = TRUE) * 100
# Use calculated break to get normal curve
grid <- seq(scl2[1] , scl2[2], length.out = 300)
y_norm_percent <- 100 * w2 * dnorm(grid, mean = mu, sd = sdx)
y_mx <- max(max(y_norm_percent, na.rm = TRUE),
max(h$counts, na.rm = TRUE), na.rm = TRUE)
# Initial chart to draw vertical orientation lines underneath bars
hist(Inf,
breaks = brks2,
right = FALSE,
include.lowest = TRUE,
main = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.025),
plot = TRUE,
axes = FALSE)
# Get axis tick marks and draw orientation lines
aval <- pretty(scl, n = 3)
abline(v = aval, col = "grey90", lwd = 1)
# Create plot using histogram values
plot(h,
col = "#CAD5E5", #"grey85",
border = "grey20",
main = "",
xlab = "",
ylab = "Percent",
xlim = scl,
ylim = c(0, y_mx * 1.025), # max(h$counts) * 1.05),
axes = FALSE,
add = TRUE)
# Normal curve overlay (scaled to percent)
lines(grid, y_norm_percent, col = "steelblue4", lwd = 2)
# Bandwidth (Silverman / SAS) - distinct n; Kish's effective n when weight only
if (!is.null(wgt2_hist) & is.null(freq2)) {
n_eff <- (sum(wgt2_hist, na.rm = TRUE))^2 / sum(wgt2_hist^2, na.rm = TRUE)
} else {
n_eff <- length(dt2)
}
bw <- 1.06 * sdx * n_eff^(-1/5) # 1.06
# Kernel density estimate
dens <- density(dt2_hist, kernel = "gaussian", bw = bw, n = 512, na.rm = TRUE)
# Scale KDE to Percent axis
y_kernel_percent <- 100 * w2 * dens$y
# Overlay on histogram
lines(dens$x, y_kernel_percent, col = "orangered2", lwd = 2)
# Add custom axes (SAS-style: fewer, rounder ticks on x)
axis(side = 1, at = pretty(scl, n = 3), col.ticks = "grey55",
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)
# Add y label
mtext("Percent", side = 2, line = 2.5)
# Add class label
legend("topleft", legend = vl2, bty = "n",
x.intersp = 0, # Spacing between group label and box
y.intersp = 0, # Spacing between content and borders
cex = .9
)
# Frame
box(col = "grey70", lwd = 1)
#******************************
#* Clean up
#******************************
box("outer", col = "grey70", lwd = 1)
# Add titles
mtext(paste0("Distribution of ", var), side = 3, line = par("oma")[3] - 1.75,
outer = TRUE, font = 2, cex = 1.25)
# Add x axis label
mtext(var, side = 1, line = par("oma")[1] - 3.5, outer = TRUE)
# Add legend
mpos <- legend("bottom",
legend = c("Normal ",
"Kernel"),
col = c("steelblue4", "orangered2"),
lwd = 2,
cex = .9,
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, -.38),
xpd = NA)
# Restore margins
par(mar = op, mfrow = c(1, 1), oma = om)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
# Legend moves to left if mean too far to the right
#' @noRd
render_boxplot1 <- function(dat, var, plt, res) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 2.25 # 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
alph <- (1 - plt$alph) * 100 # Percentage
alpha <- plt$alph # Actual value
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Prepare data
dt <- dat[[var]]
# Set margins
par(mar = c(3, 1, 3, .75) + 0.1)
# Assign labels
xlbl <- var
tlbl <- paste(" ", var)
if (!is.null(plt$varlbl)) {
xlbl <- "Difference"
tlbl <- paste0(" Difference: ", plt$varlbl)
}
# Get scales
if (plt$showh0) {
xscl <- get_scale(dt, .05, plt$h0)
} else {
xscl <- get_scale(dt, .05)
}
# Get boxplot stats
bp <- boxplot_stats1(dt, if (!is.null(plt$freq)) dat[[plt$freq]] else NULL)
if (length(bp$out) > 0) {
if (min(bp$out, na.rm = TRUE) < xscl[1]) {
xscl[1] <- min(bp$out, na.rm = TRUE) * 1.05
}
if (max(bp$out, na.rm = TRUE) > xscl[2]) {
xscl[2] <- max(bp$out, na.rm = TRUE) * 1.05
}
}
# Get stats
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
clm_tbl <- res[[grep("ConfLimits$", names(res))[length(grep("ConfLimits$", names(res)))]]]
n <- stat_tbl$N[1]
mu <- stat_tbl$MEAN[1]
sdx <- stat_tbl$STD[1]
ci <- c(clm_tbl$LCLM[1], clm_tbl$UCLM[1])
## Draw empty plot first (for layering)
plot(dt, rep(1, length(dt)),
type = "n",
xlab = "",
ylab = "",
yaxt = "n",
main = "",
xlim = xscl,
axes = FALSE) #With 95% Confidence Interval for Mean")
# Title
mtext(paste0("Distribution of", tlbl), side = 3,
line = par("mar")[3] - 1.5,
font = 2, cex = 1.25)
# Subtitle
mtext(paste0("With ", alph, "% Confidence Interval for Mean"), side = 3,
line = par("mar")[3] - 2.5,
font = 1)
# X Label
mtext(xlbl, side = 1,
line = par("mar")[1] - 1.5,
font = 1)
## Shaded CI band
usr <- par("usr")
if (is.infinite(ci[1])) ci[1] <- usr[1] - abs(usr[1])*10
if (is.infinite(ci[2])) ci[2] <- usr[2] + abs(usr[2])*10
rect(ci[1], usr[3], ci[2], usr[4],
col = "#B3D2D0",
border = NA)
# Draw axis
aval <- axis(side = 1, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Orientation lines
abline(v = aval, col = "grey90", lwd = 1)
# Show h0
if (plt$showh0) {
abline(v = plt$h0, col = "grey70", lwd = 2)
}
# Boxplot (horizontal)
bxp(bp,
horizontal = TRUE,
add = TRUE,
boxfill = "#CAD5E5",
border = "grey40",
lty = 1,
whiskcol = "#05379B",
staplecol = "#05379B",
medcol = "#05379B",
medlwd = 1,
boxwex = .65,
outline = FALSE,
axes = FALSE) # Already drawn above
## Mean diamond
points(mu, 1,
pch = 5,
col = "#05379B",
lwd = 1,
cex = 1.25)
# Display outliers
if (length(bp$out) > 0) {
points(bp$out, rep(1, length(bp$out)),
pch = 1,
col = "grey20",
lwd = 1.6,
cex = 1.25)
# Add labels
if (plt$label) {
if (is.null(plt$id)) {
olbls <- bp$onm
} else {
olbls <- dat[[plt$id]][as.integer(bp$onm)]
}
if (length(olbls) > 0) {
text(bp$out, rep(1, length(bp$out)),
labels = olbls,
cex = .9,
pos = 1)
}
}
}
# Legend location - Move to left side if needed
lgnd <- "topright"
pt <- (bp$stats[4] - xscl[1])/diff(xscl)
if (pt > .75) {
lgnd <- "topleft"
}
# Create legend
if (plt$showh0) {
if ( plt$h0 > mu) {
lgnd <- "topleft"
}
legend(lgnd,
legend = c(paste0(alph, "% Confidence"), "Null Value"),
pch = c(15, NA),
fill = c("#B3D2D0", NA),
col = c(NA, "grey70"),
lwd = c(NA, 2),
border = c("grey60", NA),
box.col = "grey80",
inset = c(.01, .02),
x.intersp = .5,
y.intersp = c(0.5, 1.8),
cex = .9,
bty = "o")
} else {
legend(lgnd,
legend = paste0(alph, "% Confidence"),
fill = "#B3D2D0",
border = "grey60",
box.col = "grey80",
inset = c(.01, .02),
x.intersp = .5,
y.intersp = c(0.5, 1.8),
cex = .9,
bty = "o")
}
# 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)
}
# Need to deal with outliers
#' @noRd
render_boxplot2 <- function(dat, var, plt, class, res) {
# browser()
op <- par("mar")
om <- par("oma")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 3.25 # 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
alph <- (1 - plt$alph) * 100 # Percentage
alpha <- plt$alph # Actual value
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Convert to data frame just in case it is a tibble
dat <- as.data.frame(dat)
if (!is.null(plt$id)) {
lbls <- dat[[plt$id]]
} else {
lbls <- NULL
}
# Add rownumbers
rownames(dat) <- seq(1, nrow(dat))
# Get filter variable and values
fvr <- "" # Filter variable
fvls <- c() # Filter values
if (!is.null(class)) {
fvr <- class
dat <- sort(dat, by = fvr)
fvls <- unique(dat[[class]])
fvls <- fvls[!is.na(fvls)]
} else {
stop("Something is wrong")
}
# Prepare data
dt <- dat[[var]]
rwnms <- rownames(dat)
dt1 <- dat[dat[[fvr]] == fvls[1], var]
nms1 <- rwnms[dat[[fvr]] == fvls[1]]
dt2 <- dat[dat[[fvr]] == fvls[2], var]
nms2 <- rwnms[dat[[fvr]] == fvls[2]]
dat1 <- dat[dat[[fvr]] == fvls[1], ]
dat2 <- dat[dat[[fvr]] == fvls[2], ]
# Get number of lines needed to show labels
minlns <- get_line_count(as.character(fvls)) + 1
lmg <- 4
if (minlns > 2) {
lmg <- minlns + 2
}
# Set margins
par(mar = c(4, lmg, 2, .75) + 0.1)
# Get x scale
# xscl <- get_scale(dt, .05)
if (plt$showh0) {
xscl1 <- get_scale(dt1, .05, plt$h0)
} else {
xscl1 <- get_scale(dt1, .05)
}
if (plt$showh0) {
xscl2 <- get_scale(dt2, .05, plt$h0)
} else {
xscl2 <- get_scale(dt2, .05)
}
# Get analysis variables
cvls <- unique(dat[[class]])
cvls <- cvls[!is.na(cvls)]
vl1 <- cvls[1]
vl2 <- cvls[2]
# Get combined scale
xscl <- c()
xscl[1] <- min(xscl1[1], xscl2[1], na.rm = TRUE)
xscl[2] <- max(xscl1[2], xscl2[2], na.rm = TRUE)
# Get boxplot stats for 2 plots
freq1 <- if (!is.null(plt$freq)) dat[dat[[fvr]] == fvls[1], plt$freq] else NULL
freq2 <- if (!is.null(plt$freq)) dat[dat[[fvr]] == fvls[2], plt$freq] else NULL
bp <- boxplot_stats2(dt1, dt2, nms1, nms2, fvls, freq1, freq2)
if (length(bp$out1) > 0) {
if (min(bp$out1, na.rm = TRUE) < xscl[1]) {
xscl[1] <- min(bp$out1, na.rm = TRUE) * 1.05
}
if (max(bp$out1, na.rm = TRUE) > xscl[2]) {
xscl[2] <- max(bp$out1, na.rm = TRUE) * 1.05
}
}
if (length(bp$out2) > 0) {
if (min(bp$out2, na.rm = TRUE) < xscl[1]) {
xscl[1] <- min(bp$out2, na.rm = TRUE) * 1.05
}
if (max(bp$out2, na.rm = TRUE) > xscl[2]) {
xscl[2] <- max(bp$out2, na.rm = TRUE) * 1.05
}
}
# Calculate
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
n1 <- stat_tbl$N[as.character(stat_tbl$CLASS) == as.character(vl1)]
mu1 <- stat_tbl$MEAN[as.character(stat_tbl$CLASS) == as.character(vl1)]
sdx1 <- stat_tbl$STD[as.character(stat_tbl$CLASS) == as.character(vl1)]
n2 <- stat_tbl$N[as.character(stat_tbl$CLASS) == as.character(vl2)]
mu2 <- stat_tbl$MEAN[as.character(stat_tbl$CLASS) == as.character(vl2)]
sdx2 <- stat_tbl$STD[as.character(stat_tbl$CLASS) == as.character(vl2)]
# Draw empty plot first, so we can put ablines() behind the boxes
plot(dt, rep(1, length(dt)),
xlim = xscl,
ylim = c(.5, 2.5),
type = "n",
xlab = "",
ylab = "",
yaxt = "n",
main = "",
axes = FALSE) #With 95% Confidence Interval for Mean")
# Title
mtext(paste0("Distribution of ", var), side = 3,
line = par("mar")[3] - 1.5,
font = 2, cex = 1.25)
# X Label
mtext(var, side = 1,
line = par("mar")[1] - 2,
font = 1)
# Y Label
mtext(fvr, side = 2,
line = par("mar")[2] - 2,
font = 1)
# Draw X axis
aval <- axis(side = 1, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# # Draw Y axis
axis(side = 2, las = 1, col.ticks = "grey55", labels = fvls,
at = c(2, 1), mgp = c(3, .5, 0), tck = -0.015)
# Orientation lines
abline(v = aval, col = "grey90", lwd = 1)
# Boxplot (horizontal)
bxp(bp,
horizontal = TRUE,
add = TRUE,
boxfill = "#CAD5E5",
border = "grey40",
lty = 1,
whiskcol = "#05379B",
staplecol = "#05379B",
medcol = "#05379B",
medlwd = 1,
show.names = TRUE,
boxwex = .425,
outline = FALSE,
axes = FALSE) # Already drawn above
## Mean diamond
points(c(mu1, mu2), c(2, 1),
pch = 5,
col = "#05379B",
lwd = 1,
cex = 1.25)
# Outliers for first class value
if (length(bp$out1) > 0) {
points(bp$out1, rep(2, length(bp$out1)),
pch = 1,
col = "grey20",
lwd = 1.6,
cex = 1.25)
# Add labels
if (plt$label) {
if (is.null(plt$id)) {
olbls1 <- bp$onm1
} else {
olbls1 <- lbls[as.integer(bp$onm1)]
}
if (length(olbls1) > 0) {
text(bp$out1, rep(2, length(bp$out1)),
labels = olbls1,
cex = .9,
pos = 1)
}
}
}
# Outliers for second class value
if (length(bp$out2) > 0) {
# Plot outlier points
points(bp$out2, rep(1, length(bp$out2)),
pch = 1,
col = "grey20",
lwd = 1.6,
cex = 1.25)
# browser()
# Add labels
if (plt$label) {
if (is.null(plt$id)) {
olbls2 <- bp$onm2
} else {
olbls2 <- lbls[as.integer(bp$onm2)]
}
if (length(olbls2) > 0) {
text(bp$out2, rep(1, length(bp$out2)),
labels = olbls2,
cex = .9,
pos = 1)
}
}
}
# Frame
box(col = "grey70", lwd = 1)
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 <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
# type: pergroup or period
#' @noRd
render_interval1 <- function(dat, var, plt, res) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 2.75 # 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
alph <- (1 - plt$alph) * 100 # Percentage
alpha <- plt$alph # Actual value
# 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")
# Prepare data
dt <- dat[[var]]
# Set margins
par(mar = c(3, 1, 3, .75) + 0.1)
# Calculate basic parameters
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
clm_tbl <- res[[grep("ConfLimits$", names(res))[length(grep("ConfLimits$", names(res)))]]]
n <- stat_tbl$N[1]
mu <- stat_tbl$MEAN[1]
sdx <- stat_tbl$STD[1]
ci <- c(clm_tbl$LCLM[1], clm_tbl$UCLM[1])
# Assign labels
xlbl <- var
tlbl <- paste0(" ", var)
if (!is.null(plt$varlbl)) {
xlbl <- "Difference"
tlbl <- paste0(" Difference: ", plt$varlbl)
}
# Get scales
xscl <- ci
if (is.infinite(xscl[1])) {xscl[1] <- mu-1.5*(xscl[2]-mu)}
if (is.infinite(xscl[2])) {xscl[2] <- mu+1.5*(mu-xscl[1])}
if (plt$showh0) {
xscl <- get_scale(xscl, .001, plt$h0)
} else {
xscl <- get_scale(xscl, .001)
}
# Create plot
plot(ci, rep(1, length(ci)),
type = "n",
yaxt = "n",
ylab = "",
xlab = "",
xlim = xscl,
main = "",
axes = FALSE)
# Title
mtext(paste0("Mean of ", tlbl), side = 3,
line = par("mar")[3] - 1.5,
font = 2, cex = 1.25)
# Subtitle
mtext(paste0("With ", alph, "% Confidence Interval"), side = 3,
line = par("mar")[3] - 2.5,
font = 1)
# X Label
mtext(xlbl, side = 1,
line = par("mar")[1] - 1.5,
font = 1)
# Draw axis
aval <- axis(side = 1, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
## CI line
if (is.infinite(ci[2])) {
# Upper is Inf: Draw arrow pointing right to the edge of the plot
arrows(mu, 1, xscl[2]/1.001, 1, col = "#05379B", lwd = 2, length = 0.1)
segments(ci[1], 1, mu, 1, col = "firebrick3", lwd = 2)
} else if (is.infinite(ci[1])) {
# Lower is -Inf: Draw arrow pointing left to the edge of the plot
arrows(ci[2], 1, xscl[1]/0.999, 1, col = "#05379B", lwd = 2, length = 0.1)
segments(ci[2], 1, mu, 1, col = "firebrick3", lwd = 2)
} else {
# Finite: Draw normal segment
segments(ci[1], 1, ci[2], 1, col = "firebrick3", lwd = 2)
}
## CI end caps
if (is.finite(ci[1])) segments(ci[1], 0.95, ci[1], 1.05, col = "firebrick3", lwd = 2)
if (is.finite(ci[2])) segments(ci[2], 0.95, ci[2], 1.05, col = "firebrick3", lwd = 2)
## Mean diamond
points(mu, 1,
pch = 5,
col = "#05379B",
cex = 1.25)
# Show h0 line
if (plt$showh0) {
abline(v = plt$h0, col = "grey70", lwd = 2)
if (plt$h0 > mu) {
lgnd <- "topleft"
} else {
lgnd <- "topright"
}
legend(lgnd,
horiz = TRUE,
legend = c("Mean", "Null Value"),
pch = c(5, NA),
col = c("#05379B", "grey70"),
lwd = c(NA, 2),
box.col = "grey80",
inset = c(.01, .02),
x.intersp = .5,
y.intersp = c(0.5, 1.8),
pt.cex = 1,
cex = .9,
bty = "o")
} else {
# Create legend
legend("topright",
legend = "Mean",
pch = 5,
col = "#05379B",
box.col = "grey80",
inset = c(.01, .02),
x.intersp = .5,
y.intersp = c(0.5, 1.8),
cex = .9,
pt.cex = 1,
bty = "o")
}
# 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_interval2 <- function(dat, var, plt, res) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 2.75 # 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
alph <- (1 - plt$alph) * 100 # Percentage
alpha <- plt$alph # Actual value
# 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")
# Prepare data
dt <- res[[grep("ConfLimits$", names(res))[length(grep("ConfLimits$", names(res)))]]]
# Set margins
par(mar = c(4, 1, 3, .75) + 0.1)
# Get class label
clslbl <- paste0("(", dt$CLASS[1], " - ", dt$CLASS[2], ")")
# Get labels for means
labels <- c("Satterthwaite", "Pooled")
# Pull off data values
mean_diff <- c(dt$MEAN[4], dt$MEAN[3])
lower_ci <- c(dt$LCLM[4], dt$LCLM[3])
upper_ci <- c(dt$UCLM[4], dt$UCLM[3])
# Create y and x values
y <- c(2, 1) # vertical positions
ci <- c(lower_ci, upper_ci)
# Get scales - seed from finite CI bounds plus the means so one-sided tests
# (where all upper_ci or all lower_ci are Inf) produce a readable axis range.
xscl <- c(min(lower_ci, na.rm = TRUE), max(upper_ci, na.rm = TRUE))
mean_mid <- mean(mean_diff, na.rm = TRUE)
if (is.infinite(xscl[1])) { xscl[1] <- mean_mid - 1.5 * (xscl[2] - mean_mid) }
if (is.infinite(xscl[2])) { xscl[2] <- mean_mid + 1.5 * (mean_mid - xscl[1]) }
xscl <- c(min(c(xscl, mean_diff), na.rm = TRUE), max(c(xscl, mean_diff, na.rm = TRUE), na.rm = TRUE))
if (plt$showh0) {
xscl <- get_scale(xscl, .001, plt$h0)
} else {
xscl <- get_scale(xscl, .001)
}
# Create base plot
plot(mean_diff, y,
xlim = xscl,
ylim = c(0.5, 2.5),
xlab = "",
ylab = "",
yaxt = "n",
pch = NA,
main = "",
axes = FALSE)
# Title
mtext(paste("Mean of ", var, " Difference ", clslbl), side = 3,
line = 1.5, cex = 1.25, font = 2)
## Subtitle
sides_lbl <- if (any(is.infinite(upper_ci))) "Upper " else if (any(is.infinite(lower_ci))) "Lower " else ""
mtext(paste0("With ", alph, "% ", sides_lbl, "Confidence Intervals"), side = 3,
line = .5, cex = 1)
# X axis label
mtext("Difference", side = 1, line = 2)
# Draw axis
aval <- axis(side = 1, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
for (k in seq_along(y)) {
if (is.infinite(upper_ci[k])) {
arrows(mean_diff[k], y[k], xscl[2] / 1.001, y[k], col = "#05379B", lwd = 3, length = 0.1)
segments(lower_ci[k], y[k], mean_diff[k], y[k], col = "brown3", lwd = 3)
} else if (is.infinite(lower_ci[k])) {
arrows(mean_diff[k], y[k], xscl[1] / 0.999, y[k], col = "#05379B", lwd = 3, length = 0.1)
segments(upper_ci[k], y[k], mean_diff[k], y[k], col = "brown3", lwd = 3)
} else {
segments(lower_ci[k], y[k], upper_ci[k], y[k], col = "brown3", lwd = 3)
}
if (is.finite(lower_ci[k])) segments(lower_ci[k], y[k] - 0.08, lower_ci[k], y[k] + 0.08, col = "brown3", lwd = 3)
if (is.finite(upper_ci[k])) segments(upper_ci[k], y[k] - 0.08, upper_ci[k], y[k] + 0.08, col = "brown3", lwd = 3)
}
## Mean diamonds
points(mean_diff, y,
pch = 5, # diamond
col = "#05379B",
cex = 1.25)
## Labels above lines
text(mean_diff, y + 0.15, labels,
cex = 0.9)
if (plt$showh0) {
abline(v = plt$h0, col = "grey70", lwd = 2)
if (plt$h0 > mean(ci, na.rm = TRUE)) {
lgnd <- "topleft"
} else {
lgnd <- "topright"
}
legend(lgnd,
horiz = TRUE,
legend = c("Mean", "Null Value"),
pch = c(5, NA),
col = c("#05379B", "grey70"),
lwd = c(NA, 2),
box.col = "grey80",
inset = c(.01, .02),
x.intersp = .5,
y.intersp = c(0.5, 1.8),
pt.cex = 1,
cex = .9,
bty = "o")
} else {
legend("topright",
legend = "Mean",
pch = 5,
col = "#05379B",
box.col = "grey80",
inset = c(.01, .02),
x.intersp = .5,
y.intersp = c(0.5, 1.8),
cex = .9,
bty = "o")
}
# 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)
}
#' @noRd
render_tqqplot1 <- function(dat, var, plt, res) {
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")
# Extract data
rdt <- dat[[var]]
freq <- if (!is.null(plt$freq)) dat[[plt$freq]] else NULL
if (!is.null(freq)) {
freq <- floor(freq)
idx <- !is.na(rdt) & !is.na(freq) & freq > 0
rdt <- rep(rdt[idx], times = freq[idx])
}
# Assign label
ylbl <- var
tlbl <- paste(" ", var)
if (!is.null(plt$varlbl)) {
ylbl <- "Difference"
tlbl <- paste0(": ", plt$varlbl)
}
# 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", tlbl),
xlab = "",
ylab = ylbl,
pch = 1,
col = "#05379B",
cex = 1.3,
axes = FALSE)
# Get stats table
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
# Slope = Standard Deviation, Intercept = Mean
if (!is.null(plt$weight)) {
# For some reason, the qqplot trend line uses the unweighted values.
# So have to calculate them here.
ovar <- dat[[var]]
owgt <- dat[[plt$weight]]
adjmean <- mean(ovar, na.rm = TRUE)
adjstd <- sd(ovar, na.rm = TRUE)
abline(a = adjmean, b = adjstd, col = "grey60")
} else {
abline(a = stat_tbl$MEAN[1], b = stat_tbl$STD[1], 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_tqqplot1_back <- function(dat, var, 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
# rdt <- dat[[var]]
#
# # Assign label
# ylbl <- var
# tlbl <- paste(" ", var)
# if (!is.null(plt$varlbl)) {
# ylbl <- "Difference"
# tlbl <- paste0(": ", plt$varlbl)
# }
#
# # Set margins
# par(mar = c(4, 5, 2, .75) + 0.1)
#
# # Get y scale
# # mx <- max(abs(rdt)) * 1.125 #1.0025
# # scl <- c(-mx, mx)
#
# # Draw plot
# qqnorm(rdt,
# main = paste0("Q-Q Plot of", tlbl),
# xlab = "",
# ylab = ylbl,
# pch = 1, # open circles
# col = "#05379B",
# cex = 1.3,
# # asp = 1,
# # ylim = scl,
# # asp = 1,
# axes = FALSE)
#
# # Add diagonal line
# # This line is a problem. Right now just putting it diagonal
# # between the corners. All qqline() functions are worse.
# usr <- par("usr")
# segments(usr[1], usr[3], usr[2], usr[4], col = "grey60")
# # qqline(rdt, col = "grey60", lwd = 1, qtype = 3)
#
# # 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)
#
#
# }
#' @noRd
render_tqqplot2 <- function(dat, var, plt, class, res) {
op <- par("mar")
om <- par("oma")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 3.25 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Prepare data
dat <- as.data.frame(dat)
dat <- sort(dat, by = class)
# Subset
cvls <- unique(dat[[class]])
dt1 <- dat[dat[[class]] == cvls[1], var]
dt2 <- dat[dat[[class]] == cvls[2], var]
freq1 <- if (!is.null(plt$freq)) dat[dat[[class]] == cvls[1], plt$freq] else NULL
freq2 <- if (!is.null(plt$freq)) dat[dat[[class]] == cvls[2], plt$freq] else NULL
if (!is.null(freq1)) {
freq1 <- floor(freq1)
idx1 <- !is.na(dt1) & !is.na(freq1) & freq1 > 0
dt1 <- rep(dt1[idx1], times = freq1[idx1])
}
if (!is.null(freq2)) {
freq2 <- floor(freq2)
idx2 <- !is.na(dt2) & !is.na(freq2) & freq2 > 0
dt2 <- rep(dt2[idx2], times = freq2[idx2])
}
# Output to image file
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Set margins
par(oma = c(0, 0, 2, .75) + 0.1,
mar = c(4, 4, 0, .0) + 0.1,
mfrow = c(1, 2))
#*******************************
#* QQPlot 1
#*******************************
# Prepare Y data
ydat <- sort(dt1)
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,
axes = FALSE)
# Get stats table
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
# Slope = Standard Deviation, Intercept = Mean
if (!is.null(plt$weight)) {
# For some reason, the qqplot trend line uses the unweighted values.
# So have to calculate them here.
ovar <- dt1
adjmean <- mean(ovar, na.rm = TRUE)
adjstd <- sd(ovar, na.rm = TRUE)
abline(a = adjmean, b = adjstd, col = "grey60")
} else {
abline(a = stat_tbl$MEAN[as.character(stat_tbl$CLASS) == as.character(cvls[1])],
b = stat_tbl$STD[as.character(stat_tbl$CLASS) == as.character(cvls[1])], 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)
# Y Axis label
mtext(var, side = 2, line = par("mar")[2] - 1.5)
# Add class label
legend("topleft", legend = cvls[1], bty = "n",
x.intersp = 0, # Spacing between group label and box
y.intersp = 0, # Spacing between content and borders
cex = .9
)
# Frame
box(col = "grey70", lwd = 1)
#*******************************
#* QQPlot 2
#*******************************
# Prepare Y data
ydat <- sort(dt2)
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,
axes = FALSE)
# Slope = Standard Deviation, Intercept = Mean
if (!is.null(plt$weight)) {
# For some reason, the qqplot trend line uses the unweighted values.
# So have to calculate them here.
ovar <- dt2
adjmean <- mean(ovar, na.rm = TRUE)
adjstd <- sd(ovar, na.rm = TRUE)
abline(a = adjmean, b = adjstd, col = "grey60")
} else {
abline(a = stat_tbl$MEAN[as.character(stat_tbl$CLASS) == as.character(cvls[2])],
b = stat_tbl$STD[as.character(stat_tbl$CLASS) == as.character(cvls[2])], 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)
# Y Axis label
mtext(var, side = 2, line = par("mar")[2] - 1.5)
# Add class label
legend("topleft", legend = cvls[2], bty = "n",
x.intersp = 0, # Spacing between group label and box
y.intersp = 0, # Spacing between content and borders
cex = .9
)
# Frame
box(col = "grey70", lwd = 1)
#*******************************
#* Clean up
#*******************************
# Title
mtext(paste0("Q-Q Plots of ", var), side = 3, outer = TRUE,
line = par("oma")[3] - 1.75, font = 2, cex = 1.25)
box("outer", col = "grey70", lwd = 1)
par(mar = op, oma = om, mfrow = c(1, 1))
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
#' @noRd
render_profiles <- function(dat, var, plt, res) {
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
alph <- (1 - plt$alph) * 100 # Percentage
alpha <- plt$alph # Actual value
# 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")
# Prepare data
splt <- trimws(strsplit(plt$varlbl, "-", fixed = TRUE)[[1]])
v1 <- splt[1]
v2 <- splt[2]
vrs <- c(v1, v2)
dt <- as.matrix(dat[ , vrs ])
## X positions (categorical)
x <- c(1, 2)
par(mar = c(3, 3, 3, 3) + 0.1)
## Empty plot frame
plot(x, range(dt, na.rm = TRUE),
type = "n",
xaxt = "n",
xlab = "",
ylab = "",
main = paste0("Paired Profiles for (", v1, ", ", v2, ")"),
axes = FALSE)
# X Axis
axis(1, at = x, labels = c(v1, v2))
# Y Axis
axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
axis(side = 4, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
## Individual subject profiles (thin lines)
for (i in seq_len(nrow(dt))) {
lines(x, dt[i, ],
col = rgb(0.2, 0.4, 0.8, 0.45), # light blue
lwd = 1)
}
## Mean profile (bold line)
mean_profile <- colMeans(dt)
lines(x, mean_profile,
col = "orangered2",
lwd = 3)
## Mean legend
legend("top",
legend = "Mean",
col = "orangered2",
lwd = 3,
box.col = "grey80",
inset = c(.01, .02),
x.intersp = .5,
y.intersp = c(0.5, 1.8),
cex = .9,
bty = "o")
# 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)
}
#' @noRd
render_agreement <- function(dat, var, plt, res) {
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
alph <- (1 - plt$alph) * 100 # Percentage
alpha <- plt$alph # Actual value
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Prepare data
splt <- trimws(strsplit(plt$varlbl, "-", fixed = TRUE)[[1]])
v1 <- splt[1]
v2 <- splt[2]
vrs <- c(v1, v2)
# Get vectors
dt1 <- dat[[v1]]
dt2 <- dat[[v2]]
## Means
m1 <- mean(dt1, na.rm = TRUE)
m2 <- mean(dt2, na.rm = TRUE)
## Axis limits
lims <- range(c(dt1, dt2), na.rm = TRUE)
# Set margins
par(mar = c(4, 10, 2, 6) + 0.1)
## Scatter plot
plot(dt1, dt2,
xlim = lims,
ylim = lims,
pch = 1,
cex = 1.25,
col = "#05379B",
xlab = "",
ylab = v2,
main = paste0("Agreement of ", v2, " and ", v1),
axes = FALSE,
asp = 1) # equal scaling (critical)
# X label
mtext(v1, side = 1, line = par("mar")[1] - 2)
# Axes
axis(side = 1, las = 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)
## Identity (45-degree) line
usr <- par("usr")
segments(usr[1], usr[3], usr[2], usr[4],
col = "grey60", lwd = 1)
## Mean point
points(m1, m2,
pch = 21,
bg = "white",
col = "#05379B",
lwd = 2,
cex = 2.5)
## Legend
legend("right",
legend = "Mean",
pch = 21,
pt.bg = "white",
col = "#05379B",
pt.lwd = 2,
# lwd = 3,
box.col = "grey80",
inset = c(-.2, 0),
x.intersp = 1,
y.intersp = c(0.5, 1.8),
pt.cex = 2.5,
xpd = NA,
bty = "o")
# 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)
}
# Utilities ---------------------------------------------------------------
# Boxplot.stats(data1) can prepare the boxplot data more easily,
# but has no parameter to use SAS-style quantiles. So need to calculate
# everything manually.
#' @noRd
boxplot_stats1 <- function(dt1, freq = NULL) {
# Apply freq expansion if specified, tracking original row indices
if (!is.null(freq)) {
f <- floor(freq)
valid <- !is.na(dt1) & !is.na(f) & f > 0
row_map <- rep(which(valid), times = f[valid])
dt1 <- rep(dt1[valid], times = f[valid])
} else {
row_map <- seq_along(dt1)
}
# Calculate quantile
q1 <- quantile(dt1, probs = c(0, .25, .5, .75, 1), type = 2, na.rm = TRUE)
# Calculate Inter-Quartile Range
iqr1 <- IQR(dt1, type = 2, na.rm = TRUE) * 1.5
# Calculate outlier limits
lm1 <- c(q1[2] - iqr1, q1[4] + iqr1)
# Determine if there are outliers
ol1 <- dt1 <= lm1[1] | dt1 >= lm1[2]
ol1 <- ifelse(is.na(ol1), FALSE, ol1)
out1 <- dt1[ol1]
# Change min and max to remove outliers
onm1 <- c()
if (length(out1) > 0) {
mn1 <- min(dt1[dt1 > lm1[1]], na.rm = TRUE)
mx1 <- max(dt1[dt1 < lm1[2]], na.rm = TRUE)
q1[1] <- mn1
q1[5] <- mx1
onm1 <- row_map[seq(1, length(dt1))[ol1]] # Obs within class: This matches SAS
# onm1 <- nms1[ol1] # Original obs #: This seems more correct
}
# Get length
lndt <- length(dt1)
# Prepare to pass data
bp <- list(
stats = matrix(q1, ncol = 1),
n = lndt,
conf = NULL,
out = out1,
onm = onm1
)
return(bp)
}
# Boxplot.stats(data1) can prepare the boxplot data more easily,
# but has no parameter to use SAS-style quantiles. So need to calculate
# everything manually.
#' @noRd
boxplot_stats2 <- function(dt1, dt2, nms1, nms2, fvls, freq1, freq2) {
# browser()
if (!is.null(freq1)) {
f1 <- floor(freq1)
valid1 <- !is.na(dt1) & !is.na(f1) & f1 > 0
nms1 <- rep(nms1[valid1], times = f1[valid1])
dt1 <- rep(dt1[valid1], times = f1[valid1])
}
if (!is.null(freq2)) {
f2 <- floor(freq2)
valid2 <- !is.na(dt2) & !is.na(f2) & f2 > 0
nms2 <- rep(nms2[valid2], times = f2[valid2])
dt2 <- rep(dt2[valid2], times = f2[valid2])
}
# Calculate quantile
q1 <- quantile(dt1, probs = c(0, .25, .5, .75, 1), type = 2, na.rm = TRUE)
q2 <- quantile(dt2, probs = c(0, .25, .5, .75, 1), type = 2, na.rm = TRUE)
# Calculate Inter-Quartile Range
iqr1 <- IQR(dt1, type = 2, na.rm = TRUE) * 1.5
iqr2 <- IQR(dt2, type = 2, na.rm = TRUE) * 1.5
# Calculate outlier limits
lm1 <- c(q1[2] - iqr1, q1[4] + iqr1)
lm2 <- c(q2[2] - iqr2, q2[4] + iqr2)
# Determine if there are outliers
ol1 <- dt1 <= lm1[1] | dt1 >= lm1[2]
ol1 <- ifelse(is.na(ol1), FALSE, ol1)
out1 <- dt1[ol1]
ol2 <- dt2 <= lm2[1] | dt2 >= lm2[2]
ol2 <- ifelse(is.na(ol2), FALSE, ol2)
out2 <- dt2[ol2]
# Change min and max to remove outliers
onm1 <- c()
if (length(out1) > 0) {
mn1 <- min(dt1[dt1 > lm1[1]], na.rm = TRUE)
mx1 <- max(dt1[dt1 < lm1[2]], na.rm = TRUE)
q1[1] <- mn1
q1[5] <- mx1
onm1 <- nms1[ol1] # Original obs #: This seems more correct
}
# Change min and max to remove outliers
onm2 <- c()
if (length(out2) > 0) {
mn2 <- min(dt1[dt2 > lm2[1]], na.rm = TRUE)
mx2 <- max(dt1[dt2 < lm2[2]], na.rm = TRUE)
q2[1] <- mn2
q2[5] <- mx2
onm2 <- nms2[ol2] # Original obs #: This seems more correct
}
# Get length
lndt <- length(dt1) + length(dt2)
# Prepare to pass data
bp <- list(
stats = cbind(q2, q1),
n = cbind(lndt, lndt),
conf = NULL,
out = numeric(0),
names = fvls,
"out1" = out1,
"out2" = out2,
"onm1" = onm1,
"onm2" = onm2
)
return(bp)
}
#' @noRd
get_sas_bins <- function(x) {
x <- x[!is.na(x)]
rng <- range(x, na.rm = TRUE)
n <- length(x)
n_bins <- ceiling((2 * n)^(1/3))
raw_w <- diff(rng) / n_bins
if (!is.finite(raw_w) || raw_w <= 0) {
return(pretty(rng, n = n_bins, min.n = 1))
}
mag <- 10^floor(log10(raw_w))
cand <- c(1, 2, 2.5, 3, 4, 5) * mag
# SAS rounds DOWN to the nearest nice candidate
le <- cand[cand <= raw_w]
nice_w <- if (length(le) > 0) max(le, na.rm = TRUE) else min(cand, na.rm = TRUE)
# SAS picks edge-aligned vs midpoint-centered bins based on whether
# X_min already lies on the nice_w grid.
m <- floor(rng[1] / nice_w) * nice_w
on_grid <- isTRUE(all.equal(rng[1] - m, 0))
if (on_grid) {
lo <- m - nice_w / 2
} else {
anchor <- 10 * mag
lo <- floor(rng[1] / anchor) * anchor
while (lo + nice_w <= rng[1]) lo <- lo + nice_w
while (lo > rng[1]) lo <- lo - nice_w
}
hi <- lo
while (hi < rng[2]) hi <- hi + nice_w
seq(lo, hi, by = nice_w)
}
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 get_sas_bins boxplot_stats2 boxplot_stats1 render_agreement render_profiles render_tqqplot2 render_tqqplot1 render_interval2 render_interval1 render_boxplot2 render_boxplot1 render_histogram2 render_histogram1 render_summary2 render_summary1 ttestplot
Documented in ttestplot
# TTestplots Main Function --------------------------------------------------
# AGREEMENT, ALL, BOX, HISTOGRAM, INTERVAL, NONE, PROFILES, QQ, SUMMARY
# agreement, boxplot, histogram, interval, profiles, qqplot, summary
#' @title Request T-Test Plots
#' @description A function to request T-test plots on a call to
#' \code{\link{proc_ttest}}. The function allows you to specify the type of
#' T-Test plots to produce. By default, it produces plots appropriate
#' to the analysis being performed.
#' You may also request specific plots
#' by passing a vector of plot names to the \code{type} parameter.
#' @details Plots produced by \code{\link{proc_ttest}} can change depending on
#' the type of analysis being performed.
#'
#' For one sample T-Tests, the function creates single variable
#' charts. The "summary" and
#' "qqplot" charts are generated by default. For this type of analysis,
#' it is helpful to set the "showh0" parameter
#' to TRUE so you can see the null hypothesis on the chart.
#'
#' For two sample T-Tests, the function produces two variable charts. For example,
#' on a two sample T-Test, the "histogram" chart will show histograms for
#' each class value on the same plot. Likewise, the "boxplot" and "qqplot" are
#' separated for each variable. The "interval" plot is split into "Satterthwaite" and
#' "Pooled" intervals. Note that the "h0" option also produces a null
#' hypothesis line for the "interval" plot, if requested.
#'
#' For a paired analysis, the "profiles" and "agreement" charts are available
#' in addition to the "histogram", "boxplot", "interval", and "qqplot" charts.
#' Paired analysis also produces single variable charts.
#'
#' To see all plots available for the requested analysis, pass the keyword
#' "all" to the \code{type} parameter.
#'
#' Any requested
#' 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_freq}}, and pass
#' the entire list to \code{\link{proc_print}}.
#' @section Plots:
#' The \code{type} parameter allows you to request several types of T-Test
#' plots. Below are the types of plots that are supported. The list shows
#' the plot type keyword needed to request the plot, and a brief description:
#' \itemize{
#' \item{\strong{agreement}: An agreement plot for the input data. Only
#' available for paired comparisons.
#' }
#' \item{\strong{boxplot/box}: Displays a box and whisker plot for group comparisons,
#' including confidence intervals (if appropriate).
#' }
#' \item{\strong{histogram/hist}: A histogram with normal curve and kernel density
#' overlays.
#' }
#' \item{\strong{interval}: Visualizes the confidence intervals for means.
#' }
#' \item{\strong{profiles}: A line plot mapping responses
#' between analysis variables. Available only for paired analysis.
#' }
#' \item{\strong{qqplot/qq}: A quantile-quantile plot used to assess the assumption
#' of normality.
#' }
#' \item{\strong{summary}: Combines the histogram and boxplot charts
#' onto the same panel. For two-sample T-Tests, both the histogram and boxplots
#' are separated for each class value.
#' }
#' }
#' The above plots may be requested as a vector of keywords to the \code{plots}
#' parameter on \code{\link{proc_ttest}}, or as vector of values
#' to the \code{type} parameter on \code{ttestplots}.
#'
#' Note that, when passed as a vector of keywords, only the requested plots
#' will be produced. That is to say, the "only" keyword in SAS is implied
#' at all times for \code{\link{proc_ttest}}.
#'
#' @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_ttest} 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 "agreement", "boxplot", "histogram",
#' "interval", "profiles", "qqplot", and "summary". The default value is "default", which
#' will produce default plots appropriate for the analysis. You may use the
#' \code{\link[common]{v}} function to pass the plot keywords unquoted.
#' @param panel Whether or not to display the summary plot combined into
#' in a single panel. Default is TRUE. A value of FALSE will create
#' individual "histogram" and "boxplot" charts instead. This parameter is
#' equivalent to the "unpack" keyword in SAS.
#' @param showh0 Whether or not to show the null hypothesis on the chart.
#' Default is FALSE. If this parameter is TRUE, a vertical line will be created
#' on some charts to show the value specified for the "h0=" option on
#' \code{\link{proc_ttest}}.
#' @param label Whether or not to label outlier values. Valid values
#' are TRUE or FALSE. Default is TRUE. If TRUE, this option will assign
#' labels to outlier values on the "summary" and "boxplot" charts.
#' @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_ttest()]
#' @examples
#' library(procs)
#'
#' # Turn off printing for CRAN checks
#' options("procs.print" = FALSE)
#'
#' # View sample data
#' sleep
#' # extra group ID
#' # 1 0.7 1 1
#' # 2 -1.6 1 2
#' # 3 -0.2 1 3
#' # 4 -1.2 1 4
#' # 5 -0.1 1 5
#' # 6 3.4 1 6
#' # 7 3.7 1 7
#' # 8 0.8 1 8
#' # 9 0.0 1 9
#' # 10 2.0 1 10
#' # 11 1.9 2 1
#' # 12 0.8 2 2
#' # 13 1.1 2 3
#' # 14 0.1 2 4
#' # 15 -0.1 2 5
#' # 16 4.4 2 6
#' # 17 5.5 2 7
#' # 18 1.6 2 8
#' # 19 4.6 2 9
#' # 20 3.4 2 10
#'
#'
#' # Example 1: plots = TRUE
#' proc_ttest(sleep,
#' var = "extra",
#' class = "group",
#' plots = TRUE)
#'
#' # Example 2: plots = "all"
#' proc_ttest(sleep,
#' var = "extra",
#' class = "group",
#' plots = "all")
#'
#' # Example 3: plots = {vector}
#' proc_ttest(sleep,
#' var = "extra",
#' class = "group",
#' plots = c("histogram", "boxplot", "interval"))
#'
#' # Example 4: plots = ttestplot()
#' proc_ttest(sleep,
#' var = "extra",
#' class = "group",
#' plots = ttestplot("summary", panel = FALSE))
#'
#' @export
ttestplot <- function(type = "default", panel = TRUE, showh0 = FALSE, label = TRUE, id = NULL) {
# Non-standard evaluation
type <- resolve_arg(type)
# agreement, boxplot, histogram, interval, profiles, qqplot, summary
# Parameter Checks
vldvals <- c("agreement", "boxplot", "box", "histogram", 'hist', 'interval', 'profiles',
"qqplot", "qq", "summary", "default", "all")
if (any(!type %in% vldvals)) {
ivd <- type[!type %in% vldvals]
stop(paste0("Parameter value for 'type' invalid: ", paste0("'", ivd, "'", collapse = ", "),
"\nValid values are: ",
"'agreement', 'boxplot', 'box', 'histogram', 'hist', 'interval', 'profiles', ",
"'qqplot', 'qq', 'summary', 'default', 'all'."
))
}
if (!panel %in% c(TRUE, FALSE)) {
stop("Parameter value for 'panel' invalid. Valid values are TRUE or FALSE.")
}
if (!showh0 %in% c(TRUE, FALSE)) {
stop("Parameter value for 'showh0' invalid. Valid values are TRUE or FALSE.")
}
# Create object
ret <- structure(list(), class = c("ttestplot", "list"))
ret$type <- tolower(type)
ret$panel <- panel
ret$showh0 <- showh0
ret$label <- label
ret$id <- id
ret$alph <- 0.95 # default for now
ret$h0 <- 0 # default for now
return(ret)
}
# Plot Rendering ----------------------------------------------------------
#' @noRd
render_ttestplot <- function (dat, var, plt, class, res) {
ret <- NULL
if ("ttestplot" %in% class(plt)) {
if (all(plt$type == "default")) {
if (!is.null(plt$varlbl)) {
plt$type <- c("summary", "profiles", "agreement", "qqplot")
} else {
plt$type <- c("summary", "qqplot")
}
}
if (all(plt$type == "all")) {
# agreement, boxplot, histogram, interval, profiles, qqplot, summary
if (!is.null(plt$varlbl)) {
plt$type <- c("summary", "histogram", "boxplot", "interval", "profiles", "agreement", "qqplot")
} else {
plt$type <- c("summary", "histogram", "boxplot", "interval", "qqplot")
}
}
typs <- plt$type
# "Unpack" logic
if ("summary" %in% typs & plt$panel == FALSE) {
# Find and replace
pos <- match("summary", typs)
vins <- c("histogram", "boxplot")
nv <- append(typs, vins, after = pos)
typs <- nv[-pos]
}
ret <- list()
# agreement, boxplot, histogram, interval, profiles, qqplot, summary
for (tp in typs) {
if (tp == "summary") {
if (is.null(class)) {
ret[["summary"]] <- render_summary1(dat, var, plt, res)
} else {
ret[["summary"]] <- render_summary2(dat, var, plt, class, res)
}
} else if (tp == "histogram" || tp == "hist") {
if (is.null(class)) {
ret[["histogram"]] <- render_histogram1(dat, var, plt, res)
} else {
ret[["histogram"]] <- render_histogram2(dat, var, plt, class, res)
}
} else if (tp == "boxplot" || tp == "box") {
if (is.null(class)) {
ret[["boxplot"]] <- render_boxplot1(dat, var, plt, res)
} else {
ret[["boxplot"]] <- render_boxplot2(dat, var, plt, class, res)
}
} else if (tp == "qqplot" || tp == "qq") {
if (is.null(class)) {
ret[["qqplot"]] <- render_tqqplot1(dat, var, plt, res)
} else {
ret[["qqplot"]] <- render_tqqplot2(dat, var, plt, class, res)
}
} else if (tp == "interval") {
if (is.null(class)) {
ret[["interval"]] <- render_interval1(dat, var, plt, res)
} else {
ret[["interval"]] <- render_interval2(dat, var, plt, res)
}
} else if (tp == "profiles") {
ret[["profiles"]] <- render_profiles(dat, var, plt, res)
} else if (tp == "agreement") {
ret[["agreement"]] <- render_agreement(dat, var, plt, res)
}
}
if (length(ret) == 0) {
ret <- NULL
}
}
return(ret)
}
# Summary Plot ------------------------------------------------------------
# unpack: TRUE or FALSE
#' @noRd
render_summary1 <- function(dat, var, plt, res) {
op <- par("mar")
of <- par("fig")
# 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
alph <- (1 - plt$alph) * 100 # Percentage
alpha <- plt$alph # Actual value
# 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 <- dat[[var]]
wgt <- if (!is.null(plt$weight)) dat[[plt$weight]] else NULL
freq <- if (!is.null(plt$freq)) dat[[plt$freq]] else NULL
if (!is.null(freq)) {
freq <- floor(freq)
idx <- !is.na(rdt) & !is.na(freq) & freq > 0
if (!is.null(wgt))
idx <- idx & !is.na(wgt) & wgt > 0
rdt_hist <- rep(rdt[idx], times = freq[idx])
wgt_hist <- if (!is.null(wgt)) rep(wgt[idx], times = freq[idx]) else NULL
} else {
rdt_hist <- rdt
wgt_hist <- wgt
}
# dt <- dat[[ivr]]
# Assign labels
xlbl <- var
tlbl <- paste(" ", var)
if (!is.null(plt$varlbl)) {
xlbl <- "Difference"
tlbl <- paste0(" Difference: ", plt$varlbl)
}
#**********************************
#* Histogram
#**********************************
# Set margins
par(mar = c(0, 5, 3, .75) + 0.1,
fig = c(0, 1, .3, 1))
# Calculate stats
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
clm_tbl <- res[[grep("ConfLimits$", names(res))[length(grep("ConfLimits$", names(res)))]]]
n <- stat_tbl$N[1]
# mu <- stat_tbl$MEAN[1]
# sdx <- stat_tbl$STD[1]
mu <- mean(rdt, na.rm = TRUE)
sdx <- sd(rdt, na.rm = TRUE)
ci <- c(clm_tbl$LCLM[1], clm_tbl$UCLM[1])
# Calculate breaks - Closer to SAS algorithm
brks <- pretty(range(rdt, na.rm = TRUE), n = nclass.Sturges(rdt),
min.n = 1, high.u.bias = 3)
# Calculate breaks
# brks <- get_sas_bins(rdt_hist)
# Calculate breaks and y scale.
h <- hist(rdt_hist,
breaks = brks,
right = FALSE,
include.lowest = TRUE,
plot = FALSE)
h$counts <- h$counts / sum(h$counts, na.rm = TRUE) * 100
# Use calculated breaks to get scale (1.5 bin-widths)
w <- diff(h$breaks)[1] # bin width
xmin <- min(h$breaks, na.rm = TRUE) - 1.5 * w
xmax <- max(h$breaks, na.rm = TRUE) + 1.5 * w
scl <- c(xmin, xmax)
# Consider h0 if needed
if (plt$showh0) {
if (plt$h0 < scl[1]) {
scl[1] <- plt$h0
}
if (plt$h0 > scl[2]) {
scl[2] <- plt$h0
}
}
# Use calculated break to get normal curve
curve_lo <- min(h$breaks, na.rm = TRUE) - w
curve_hi <- max(h$breaks, na.rm = TRUE) + w
grid <- seq(curve_lo, curve_hi, length.out = 300)
y_norm_percent <- 100 * w * dnorm(grid, mean = mu, sd = sdx)
y_mx <- max(max(y_norm_percent, na.rm = TRUE), max(h$counts, na.rm = TRUE), na.rm = TRUE)
# Initial chart to draw vertical lines
hist(Inf,
breaks = brks,
right = FALSE,
include.lowest = TRUE,
main = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.025),
plot = TRUE,
axes = FALSE)
# Get axis tick marks
aval <- pretty(scl, n = 3)
# Orientation lines
abline(v = aval, col = "grey90", lwd = 1)
# Create plot using histogram values
plot(h,
col = "#CAD5E5", #"grey85",
border = "grey20",
main = "",
xlab = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.025),
axes = FALSE,
add = TRUE)
# Title
mtext(paste0("Distribution of", tlbl), side = 3,
line = par("mar")[3] - 1.5,
font = 2, cex = 1.25)
# Subtitle
mtext(paste0("With ", alph, "% Confidence Interval for Mean"), side = 3,
line = par("mar")[3] - 2.5,
font = 1)
# Y Label
mtext("Percent", side = 2, line = par("mar")[2] - 2)
# Normal curve overlay (scaled to percent)
lines(grid, y_norm_percent, col = "steelblue4", lwd = 2)
# Kernel density overlay (scaled to percent)
# Bandwidth (Silverman / SAS) - SAS uses the distinct (non-freq-expanded)
# observation count for the bandwidth. Kish's effective n when weighted.
if (!is.null(wgt_hist)&is.null(freq)) {
n_eff <- (sum(wgt_hist, na.rm = TRUE))^2 / sum(wgt_hist^2, na.rm = TRUE)
} else {
n_eff <- length(rdt)
}
bw <- 1.06 * sdx * n_eff^(-1/5)
# Kernel density estimate
kernel_lo <- min(h$breaks, na.rm = TRUE) - 1.25 * w
kernel_hi <- max(h$breaks, na.rm = TRUE) + 1.25 * w
dens <- density(rdt_hist, kernel = "gaussian", bw = bw, n = 512,
from = kernel_lo, to = kernel_hi, na.rm = TRUE)
# Scale KDE to Percent axis
y_kernel_percent <- 100 * w * dens$y
# Overlay density on histogram
lines(dens$x, y_kernel_percent, col = "orangered2", lwd = 2)
# Add custom Y axis
axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Frame
box(col = "grey70", lwd = 1)
box("outer", 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)
)
# Custom Legend border
rect(mpos$rect$left, mpos$rect$top - (mpos$rect$h - .5),
mpos$rect$left + mpos$rect$w, mpos$rect$top,
border = "grey80"
)
#**********************************
#* Boxplot
#**********************************
par(mar = c(4, 5, 0, .75) + 0.1,
fig = c(0, 1, 0, .3), new = TRUE)
# Get data
dt <- dat[[var]]
# Get scales
xscl <- get_scale(dt, .05) # Not used?
# Calculations
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
clm_tbl <- res[[grep("ConfLimits$", names(res))[length(grep("ConfLimits$", names(res)))]]]
n <- stat_tbl$N[1]
mu <- mean(dt, na.rm = TRUE)
sdx <- stat_tbl$STD[1]
ci <- c(clm_tbl$LCLM[1], clm_tbl$UCLM[1])
## Draw empty plot first (for layering)
plot(dt, rep(1, length(dt)),
type = "n",
xlab = "",
ylab = "",
yaxt = "n",
xlim = scl,
main = "",
axes = FALSE) #With 95% Confidence Interval for Mean")
# X Label
mtext(xlbl, side = 1,
line = par("mar")[1] - 2,
font = 1)
## Shaded CI band and make sure infinite values are out of boundry
usr <- par("usr")
if (is.infinite(ci[1])) ci[1] <- usr[1] - abs(usr[1])*10
if (is.infinite(ci[2])) ci[2] <- usr[2] + abs(usr[2])*10
rect(ci[1], usr[3], ci[2], usr[4],
col = "#B3D2D0",
border = NA)
# Draw axis (SAS-style: fewer, rounder ticks)
aval <- axis(side = 1, at = pretty(scl, n = 3), las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Orientation lines
abline(v = aval, col = "grey90", lwd = 1)
# Show h0
if (plt$showh0) {
abline(v = plt$h0, col = "grey70", lwd = 2)
}
# Get boxplot stats
bp <- boxplot_stats1(dt, if (!is.null(plt$freq)) dat[[plt$freq]] else NULL)
# Boxplot (horizontal)
bxp(bp,
horizontal = TRUE,
add = TRUE,
boxfill = "#CAD5E5",
border = "grey40",
lty = 1,
whiskcol = "#05379B",
staplecol = "#05379B",
medcol = "#05379B",
medlwd = 1,
boxwex = .65,
outline = FALSE,
axes = FALSE) # Already drawn above
## Mean diamond
points(mu, 1,
pch = 5,
col = "#05379B",
lwd = 1,
cex = 1.25)
# Display outliers
if (length(bp$out) > 0) {
points(bp$out, rep(1, length(bp$out)),
pch = 1,
col = "grey20",
lwd = 1.6,
cex = 1.25)
# Add labels
if (plt$label) {
if (is.null(plt$id)) {
olbls <- bp$onm
} else {
olbls <- dat[[plt$id]][as.integer(bp$onm)]
}
if (length(olbls) > 0) {
text(bp$out, rep(1, length(bp$out)),
labels = olbls,
cex = .9,
pos = 1)
}
}
}
# Legend location - Move to left side if needed
lgnd <- "topright"
pt <- (bp$stats[4] - xscl[1])/diff(xscl)
if (pt > .75) {
lgnd <- "topleft"
}
# Create legend
if (plt$showh0) {
mpos <- legend(lgnd,
legend = c(paste0(alph, "% Confidence"), "Null Value"),
pch = c(15, NA),
fill = c("#B3D2D0", NA),
col = c(NA, "grey70"),
lwd = c(NA, 2),
border = c("grey60", NA),
inset = c(.01, .02),
x.intersp = .5,
y.intersp = c(0.5, 1.5, .5),
cex = .8,
bty = "n")
rect(mpos$rect$left, mpos$rect$top - (mpos$rect$h - .1),
mpos$rect$left + mpos$rect$w, mpos$rect$top,
border = "grey80"
)
} else {
legend(lgnd,
legend = paste0(alph, "% Confidence"),
fill = "#B3D2D0",
border = "grey60",
box.col = "grey80",
inset = c(.01, .02),
x.intersp = .5,
y.intersp = c(0.5, 1.8),
cex = .9,
bty = "o")
}
# Frame
box(col = "grey70", lwd = 1)
#**********************************
#* Clean up
#**********************************
# Restore margins
par(mar = op, fig = of)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
#' @noRd
render_summary2 <- function(dat, var, plt, class, res) {
op <- par("mar")
om <- par("oma")
of <- par("fig")
# 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
alph <- (1 - plt$alph) * 100 # Percentage
alpha <- plt$alph # Actual value
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Convert to data frame just in case it is a tibble
dat <- as.data.frame(dat)
if (!is.null(plt$id)) {
lbls <- dat[[plt$id]]
} else {
lbls <- NULL
}
# Add rownumbers
rownames(dat) <- seq(1, nrow(dat))
dat <- sort(dat, by = class)
# Get analysis variables
cvls <- unique(dat[[class]])
cvls <- cvls[!is.na(cvls)]
vl1 <- cvls[1]
vl2 <- cvls[2]
# Prepare data
dt1 <- dat[dat[[class]] == vl1, var]
dt2 <- dat[dat[[class]] == vl2, var]
wgt1 <- if (!is.null(plt$weight)) dat[dat[[class]] == vl1, plt$weight] else NULL
wgt2 <- if (!is.null(plt$weight)) dat[dat[[class]] == vl2, plt$weight] else NULL
freq1 <- if (!is.null(plt$freq)) dat[dat[[class]] == vl1, plt$freq] else NULL
freq2 <- if (!is.null(plt$freq)) dat[dat[[class]] == vl2, plt$freq] else NULL
# Apply freq-expansion to mirror SAS hist behavior
if (!is.null(freq1)) {
freq1 <- floor(freq1)
idx1 <- !is.na(dt1) & !is.na(freq1) & freq1 > 0
if (!is.null(wgt1))
idx1 <- idx1 & !is.na(wgt1) & wgt1 > 0
dt1_hist <- rep(dt1[idx1], times = freq1[idx1])
wgt1_hist <- if (!is.null(wgt1)) rep(wgt1[idx1], times = freq1[idx1]) else NULL
} else {
dt1_hist <- dt1
wgt1_hist <- wgt1
}
if (!is.null(freq2)) {
freq2 <- floor(freq2)
idx2 <- !is.na(dt2) & !is.na(freq2) & freq2 > 0
if (!is.null(wgt2))
idx2 <- idx2 & !is.na(wgt2) & wgt2 > 0
dt2_hist <- rep(dt2[idx2], times = freq2[idx2])
wgt2_hist <- if (!is.null(wgt2)) rep(wgt2[idx2], times = freq2[idx2]) else NULL
} else {
dt2_hist <- dt2
wgt2_hist <- wgt2
}
# Get number of lines needed to show labels
minlns <- get_line_count(as.character(cvls)) + 1
lmg <- 3
if (minlns > 2) {
lmg <- minlns + 1
}
# Set margins
par(oma = c(5, 1, 2, .75) + 0.1,
mar = c(.1, lmg, 0, 0) + 0.1)
#******************************
#* Calculate scales
#******************************
# Overall scale
scl <- range(dat[[var]], na.rm = TRUE)
scl <- c(scl[1] * .875, scl[2] * 1.125)
# Calculate breaks - Histogram 1
brks1 <- pretty(range(dt1, na.rm = TRUE), n = nclass.Sturges(dt1),
min.n = 1, high.u.bias = 3)
# Calculate breaks - Histogram 1
# brks1 <- get_sas_bins(dt1_hist)
# Calculate breaks - Histogram 2
brks2 <- pretty(range(dt2, na.rm = TRUE), n = nclass.Sturges(dt2),
min.n = 1, high.u.bias = 3)
# Calculate breaks - Histogram 2
# brks2 <- get_sas_bins(dt2_hist)
# Scale for Histogram 1 (SAS-style: 1.5 bin-widths of padding)
w1 <- diff(brks1)[1] # bin width
xmin1 <- min(brks1, na.rm = TRUE) - 1.5 * w1
xmax1 <- max(brks1, na.rm = TRUE) + 1.5 * w1
scl1 <- c(xmin1, xmax1)
# Scale for Histogram 2 (SAS-style: 1.5 bin-widths of padding)
w2 <- diff(brks2)[1] # bin width
xmin2 <- min(brks2, na.rm = TRUE) - 1.5 * w2
xmax2 <- max(brks2, na.rm = TRUE) + 1.5 * w2
scl2 <- c(xmin2, xmax2)
# Adjust overall scale if needed by histogram 1
if (scl1[1] < scl[1]) {
scl[1] <- scl1[1]
}
if (scl1[2] > scl[2]) {
scl[2] <- scl1[2]
}
# Adjust overall scale if needed by histogram 2
if (scl2[1] < scl[1]) {
scl[1] <- scl2[1]
}
if (scl2[2] > scl[2]) {
scl[2] <- scl2[2]
}
#******************************
#* Histogram 1
#******************************
par(fig = c(0, 1, .6, 1))
# Calculate stats
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
clm_tbl <- res[[grep("ConfLimits$", names(res))[length(grep("ConfLimits$", names(res)))]]]
n <- stat_tbl$N[as.character(stat_tbl$CLASS) == as.character(vl1)]
mu <- mean(dt1, na.rm = TRUE)
sdx <- sd(dt1, na.rm = TRUE)
# Histogram (Percent scale) - SAS uses left-closed bins [a, b)
h <- hist(dt1_hist,
breaks = brks1,
right = FALSE,
include.lowest = TRUE,
plot = FALSE)
h$counts <- h$counts / sum(h$counts, na.rm = TRUE) * 100
# Use calculated break to get normal curve
grid <- seq(scl1[1] , scl1[2], length.out = 300)
y_norm_percent <- 100 * w1 * dnorm(grid, mean = mu, sd = sdx)
y_mx <- max(max(y_norm_percent, na.rm = TRUE), max(h$counts, na.rm = TRUE), na.rm = TRUE)
# Initial chart to draw vertical orientation lines underneath bars
hist(Inf,
breaks = brks1,
right = FALSE,
include.lowest = TRUE,
main = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.025),
plot = TRUE,
axes = FALSE)
# Get axis tick marks and draw orientation lines
aval <- pretty(scl, n = 3)
abline(v = aval, col = "grey90", lwd = 1)
# Create plot using histogram values
plot(h,
col = "#CAD5E5", #"grey85",
border = "grey20",
main = "",
xlab = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.025), # max(h$counts) * 1.05),
axes = FALSE,
add = TRUE)
# Normal curve overlay (scaled to percent)
lines(grid, y_norm_percent, col = "steelblue4", lwd = 2)
# Kernel density overlay (scaled to percent)
# Bandwidth (Silverman / SAS) - SAS uses distinct (non-freq-expanded) n;
# Kish's effective n when only weight is supplied.
if (!is.null(wgt1_hist) & is.null(freq1)) {
n_eff <- (sum(wgt1_hist, na.rm = TRUE))^2 / sum(wgt1_hist^2, na.rm = TRUE)
} else {
n_eff <- length(dt1)
}
bw <- 1.06 * sdx * n_eff^(-1/5)
# Kernel density estimate
dens <- density(dt1_hist, kernel = "gaussian", bw = bw, n = 512, na.rm = TRUE)
# Scale KDE to Percent axis
y_kernel_percent <- 100 * w1 * dens$y
yscl <- seq(50, length(y_kernel_percent) - 25)
# Overlay on histogram
lines(dens$x, y_kernel_percent, col = "orangered2", lwd = 2)
# Add custom axes
axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Add y label
mtext("Percent", side = 2, line = par()$mar[2] - .75)
# Add class label
legend("topleft", legend = vl1, bty = "n",
x.intersp = 0, # Spacing between group label and box
y.intersp = 0, # Spacing between content and borders
cex = .9
)
# Frame
box(col = "grey70", lwd = 1)
#******************************
#* Histogram 2
#******************************
par(fig = c(0, 1, .2, .6), new = TRUE)
# Calculate stats
n <- stat_tbl$N[as.character(stat_tbl$CLASS) == as.character(vl2)]
mu <- mean(dt2, na.rm = TRUE)
sdx <- sd(dt2, na.rm = TRUE)
# Histogram (Percent scale) - SAS uses left-closed bins [a, b)
h <- hist(dt2_hist,
breaks = brks2,
right = FALSE,
include.lowest = TRUE,
plot = FALSE)
h$counts <- h$counts / sum(h$counts, na.rm = TRUE) * 100
# Use calculated break to get normal curve
grid <- seq(scl2[1] , scl2[2], length.out = 300)
y_norm_percent <- 100 * w2 * dnorm(grid, mean = mu, sd = sdx)
y_mx <- max(max(y_norm_percent, na.rm = TRUE),
max(h$counts, na.rm = TRUE), na.rm = TRUE)
# Initial chart to draw vertical orientation lines underneath bars
hist(Inf,
breaks = brks2,
right = FALSE,
include.lowest = TRUE,
main = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.025),
plot = TRUE,
axes = FALSE)
# Get axis tick marks and draw orientation lines
aval <- pretty(scl, n = 3)
abline(v = aval, col = "grey90", lwd = 1)
# Create plot using histogram values
plot(h,
col = "#CAD5E5", #"grey85",
border = "grey20",
main = "",
xlab = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.025), # max(h$counts) * 1.05),
axes = FALSE,
add = TRUE)
# Normal curve overlay (scaled to percent)
lines(grid, y_norm_percent, col = "steelblue4", lwd = 2)
# Kernel density overlay (scaled to percent)
# Bandwidth (Silverman / SAS) - distinct n; Kish's effective n when weight only
if (!is.null(wgt2_hist) & is.null(freq2)) {
n_eff <- (sum(wgt2_hist, na.rm = TRUE))^2 / sum(wgt2_hist^2, na.rm = TRUE)
} else {
n_eff <- length(dt2)
}
bw <- 1.06 * sdx * n_eff^(-1/5) # 1.06
# Kernel density estimate
dens <- density(dt2_hist, kernel = "gaussian", bw = bw, n = 512, na.rm = TRUE)
# Scale KDE to Percent axis
y_kernel_percent <- 100 * w2 * dens$y
# Overlay on histogram
lines(dens$x, y_kernel_percent, col = "orangered2", lwd = 2)
# Add custom axes
axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Add y label
mtext("Percent", side = 2, line = par()$mar[2] - .75)
# Add class label
legend("topleft", legend = vl2, bty = "n",
x.intersp = 0, # Spacing between group label and box
y.intersp = 0, # Spacing between content and borders
cex = .9
)
# Frame
box(col = "grey70", lwd = 1)
#**********************************
#* Boxplot
#**********************************
par(fig = c(0, 1, 0, .2), new = TRUE)
# Prepare data
dt <- dat[[var]]
rwnms <- rownames(dat)
nms1 <- rwnms[dat[[class]] == vl1]
nms2 <- rwnms[dat[[class]] == vl2]
# Get x scale
# xscl <- get_scale(dt, .05)
# Calculate
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
clm_tbl <- res[[grep("ConfLimits$", names(res))[length(grep("ConfLimits$", names(res)))]]]
n1 <- stat_tbl$N[as.character(stat_tbl$CLASS) == as.character(vl1)]
mu1 <- stat_tbl$MEAN[as.character(stat_tbl$CLASS) == as.character(vl1)]
sdx1 <- stat_tbl$STD[as.character(stat_tbl$CLASS) == as.character(vl1)]
n2 <- stat_tbl$N[as.character(stat_tbl$CLASS) == as.character(vl2)]
mu2 <- stat_tbl$MEAN[as.character(stat_tbl$CLASS) == as.character(vl2)]
sdx2 <- stat_tbl$STD[as.character(stat_tbl$CLASS) == as.character(vl2)]
# Draw empty plot first, so we can put ablines() behind the boxes
plot(dt, rep(1, length(dt)),
ylim = c(.5, 2.5),
xlim = scl,
type = "n",
xlab = "",
ylab = "",
yaxt = "n",
main = "",
axes = FALSE) #With 95% Confidence Interval for Mean")
# Y Label
mtext(class, side = 2,
line = par()$mar[2] - .75,
font = 1)
# Draw X axis
aval <- axis(side = 1, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# # Draw Y axis
axis(side = 2, las = 1, col.ticks = "grey55", labels = cvls,
at = c(2, 1), mgp = c(3, .5, 0), tck = -0.015)
# Orientation lines
abline(v = aval, col = "grey90", lwd = 1)
# Get boxplot stats for 2 plots
freq1 <- if (!is.null(plt$freq)) dat[dat[[class]] == vl1, plt$freq] else NULL
freq2 <- if (!is.null(plt$freq)) dat[dat[[class]] == vl2, plt$freq] else NULL
bp <- boxplot_stats2(dt1, dt2, nms1, nms2, cvls, freq1, freq2)
# Boxplot (horizontal)
bxp(bp,
horizontal = TRUE,
add = TRUE,
boxfill = "#CAD5E5",
border = "grey40",
lty = 1,
whiskcol = "#05379B",
staplecol = "#05379B",
medcol = "#05379B",
medlwd = 1,
show.names = TRUE,
boxwex = .425,
outline = FALSE,
axes = FALSE) # Already drawn above
## Mean diamond
points(c(mu1, mu2), c(2, 1),
pch = 5,
col = "#05379B",
lwd = 1,
cex = 1.25)
# Outliers for first class value
if (length(bp$out1) > 0) {
points(bp$out1, rep(2, length(bp$out1)),
pch = 1,
col = "grey20",
lwd = 1.6,
cex = 1.25)
# Add labels
if (plt$label) {
if (is.null(plt$id)) {
olbls1 <- bp$onm1
} else {
olbls1 <- lbls[as.integer(bp$onm1)]
}
if (length(olbls1) > 0) {
text(bp$out1, rep(2, length(bp$out1)),
labels = olbls1,
cex = .9,
pos = 1)
}
}
}
# Outliers for second class value
if (length(bp$out2) > 0) {
# Plot outlier points
points(bp$out2, rep(1, length(bp$out2)),
pch = 1,
col = "grey20",
lwd = 1.6,
cex = 1.25)
# Add labels
if (plt$label) {
if (is.null(plt$id)) {
olbls2 <- bp$onm2
} else {
olbls2 <- lbls[as.integer(bp$onm2)]
}
if (length(olbls2) > 0) {
text(bp$out2, rep(1, length(bp$out2)),
labels = olbls2,
cex = .9,
pos = 3)
}
}
}
# Frame
box(col = "grey70", lwd = 1)
#**********************************
#* Clean up
#**********************************
# Title
mtext(paste0("Distribution of ", var), side = 3,
line = par("oma")[3] - 1.75,
font = 2, cex = 1.25, outer = TRUE)
# X Label
mtext(var, side = 1,
line = par("oma")[1] - 3.5,
font = 1, outer = TRUE)
box("outer", col = "grey70", lwd = 1)
# Legend
mpos <- legend("bottom",
legend = c("Normal", "Kernel"),
col = c("steelblue4", "orangered2"),
horiz = TRUE,
lwd = 2,
box.col = "grey80",
bty = "o",
cex = .9,
xpd = NA,
inset = c(0, -1),
x.intersp = .5,
y.intersp = c(0.3, .3) # c(0.5, 1.8)
)
# Restore margins
par(mar = op, fig = of, oma = om)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
# Individual Plots --------------------------------------------------------
# agreement, boxplot, histogram, interval, profiles, qqplot, summary
# Curves for normal and density most correct on this one.
#' @noRd
render_histogram1 <- function(dat, var, plt, res) {
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 <- dat[[var]]
# rdt <- rdt[!is.na(rdt)]
wgt <- if (!is.null(plt$weight)) dat[[plt$weight]] else NULL
freq <- if (!is.null(plt$freq)) dat[[plt$freq]] else NULL
if (!is.null(freq)) {
freq <- floor(freq)
idx <- !is.na(rdt) & !is.na(freq) & freq > 0
if (!is.null(wgt))
idx <- idx & !is.na(wgt) & wgt > 0
rdt_hist <- rep(rdt[idx], times = freq[idx])
wgt_hist <- if (!is.null(wgt)) rep(wgt[idx], times = freq[idx]) else NULL
} else {
rdt_hist <- rdt
wgt_hist <- wgt
}
# Set margins
par(mar = c(4, 5, 2, .75) + 0.1)
# Assign labels
xlbl <- var
tlbl <- paste(" ", var)
if (!is.null(plt$varlbl)) {
xlbl <- "Difference"
tlbl <- paste0(" Difference: ", plt$varlbl)
}
# Calculate stats
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
clm_tbl <- res[[grep("ConfLimits$", names(res))[length(grep("ConfLimits$", names(res)))]]]
n <- stat_tbl$N[1]
mu <- mean(rdt, na.rm = TRUE)
sdx <- sd(rdt, na.rm = TRUE)
ci <- c(clm_tbl$LCLM[1], clm_tbl$UCLM[1])
# Calculate breaks - Closer to SAS algorithm
brks <- pretty(range(rdt, na.rm = TRUE), n = nclass.Sturges(rdt),
min.n = 1, high.u.bias = 3)
# Calculate breaks
# brks <- get_sas_bins(rdt_hist)
# Histogram (Percent scale) - SAS uses left-closed bins [a, b)
h <- hist(rdt_hist,
breaks = brks,
right = FALSE,
include.lowest = TRUE,
plot = FALSE)
h$counts <- h$counts / sum(h$counts, na.rm = TRUE) * 100
# Use calculated breaks to get scale (1.5 bin-widths)
w <- diff(h$breaks)[1] # bin width
xmin <- min(h$breaks, na.rm = TRUE) - 1.5 * w
xmax <- max(h$breaks, na.rm = TRUE) + 1.5 * w
scl <- c(xmin, xmax)
curve_lo <- min(h$breaks, na.rm = TRUE) - w
curve_hi <- max(h$breaks, na.rm = TRUE) + w
kernel_lo <- min(h$breaks, na.rm = TRUE) - 1.4 * w
kernel_hi <- max(h$breaks, na.rm = TRUE) + 1.4 * w
# Bandwidth (Silverman / SAS) - distinct n; Kish's effective n when weight only
if (!is.null(wgt_hist) & is.null(freq)) {
n_eff <- (sum(wgt_hist, na.rm = TRUE))^2 / sum(wgt_hist^2, na.rm = TRUE)
} else {
n_eff <- length(rdt)
}
bw <- 1.06 * sdx * n_eff^(-1/5)
# Kernel density estimate
dens <- density(rdt_hist, kernel = "gaussian", bw = bw, n = 512,
from = kernel_lo, to = kernel_hi, na.rm = TRUE)
# Use calculated break to get normal curve
grid <- seq(curve_lo, curve_hi, length.out = 300)
y_norm_percent <- 100 * w * dnorm(grid, mean = mu, sd = sdx)
y_kernel_percent <- 100 * w * dens$y
y_mx <- max(max(y_norm_percent, na.rm = TRUE), max(h$counts, na.rm = TRUE), na.rm = TRUE)
# Initial chart to draw vertical orientation lines underneath bars
hist(Inf, # Putting values here was drawing a single bar, and messing up the plot
breaks = brks,
right = FALSE,
include.lowest = TRUE,
main = paste0("Distribution of", tlbl),
xlab = "",
ylab = "Percent",
xlim = scl,
ylim = c(0, y_mx * 1.025),
plot = TRUE,
axes = FALSE)
# Get axis tick marks and draw orientation lines
aval <- pretty(scl, n = 3)
abline(v = aval, col = "grey90", lwd = 1)
# Create plot using histogram values
plot(h,
col = "#CAD5E5", #"grey85",
border = "grey20",
main = "",
xlab = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.025),
axes = FALSE,
add = TRUE)
# Draw normal curve line
lines(grid, y_norm_percent, col = "steelblue4", lwd = 2)
# Overlay on histogram
lines(dens$x, y_kernel_percent, col = "orangered2", lwd = 2)
# Add x axis label
mtext(xlbl, side = 1, line = par("mar")[1] - 2)
# Add custom axes
axis(side = 1, at = pretty(scl, n = 3), col.ticks = "grey55",
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_histogram2 <- function(dat, var, plt, class, res) {
op <- par("mar")
om <- par("oma")
# 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
# Convert to data frame
dat <- as.data.frame(dat)
dat <- sort(dat, by = class)
# Get analysis variables
cvls <- unique(dat[[class]])
vl1 <- cvls[1]
vl2 <- cvls[2]
# Prepare data
dt1 <- dat[dat[[class]] == vl1, var]
dt2 <- dat[dat[[class]] == vl2, var]
wgt1 <- if (!is.null(plt$weight)) dat[dat[[class]] == vl1, plt$weight] else NULL
wgt2 <- if (!is.null(plt$weight)) dat[dat[[class]] == vl2, plt$weight] else NULL
freq1 <- if (!is.null(plt$freq)) dat[dat[[class]] == vl1, plt$freq] else NULL
freq2 <- if (!is.null(plt$freq)) dat[dat[[class]] == vl2, plt$freq] else NULL
# Apply freq-expansion to mirror SAS hist behavior
if (!is.null(freq1)) {
freq1 <- floor(freq1)
idx1 <- !is.na(dt1) & !is.na(freq1) & freq1 > 0
if (!is.null(wgt1))
idx1 <- idx1 & !is.na(wgt1) & wgt1 > 0
dt1_hist <- rep(dt1[idx1], times = freq1[idx1])
wgt1_hist <- if (!is.null(wgt1)) rep(wgt1[idx1], times = freq1[idx1]) else NULL
} else {
dt1_hist <- dt1
wgt1_hist <- wgt1
}
if (!is.null(freq2)) {
freq2 <- floor(freq2)
idx2 <- !is.na(dt2) & !is.na(freq2) & freq2 > 0
if (!is.null(wgt2))
idx2 <- idx2 & !is.na(wgt2) & wgt2 > 0
dt2_hist <- rep(dt2[idx2], times = freq2[idx2])
wgt2_hist <- if (!is.null(wgt2)) rep(wgt2[idx2], times = freq2[idx2]) else NULL
} else {
dt2_hist <- dt2
wgt2_hist <- wgt2
}
# Output to image file
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Set margins
par(oma = c(5, 1, 2, .75) + 0.1,
mar = c(0, 3, 0, 0) + 0.1,
mfrow = c(2, 1))
#******************************
#* Calculate scales
#******************************
# Overall scale
scl <- range(dat[[var]], na.rm = TRUE)
scl <- c(scl[1] * .875, scl[2] * 1.125)
# Calculate breaks - Histogram 1
brks1 <- pretty(range(dt1, na.rm = TRUE), n = nclass.Sturges(dt1),
min.n = 1, high.u.bias = 3)
# Calculate breaks - Histogram 2
brks2 <- pretty(range(dt2, na.rm = TRUE), n = nclass.Sturges(dt2),
min.n = 1, high.u.bias = 3)
# Calculate breaks - Histogram 1
# brks1 <- get_sas_bins(dt1_hist)
# Calculate breaks - Histogram 2
# brks2 <- get_sas_bins(dt2_hist)
# Scale for Histogram 1 (SAS-style: 1.5 bin-widths of padding)
w1 <- diff(brks1)[1] # bin width
xmin1 <- min(brks1, na.rm = TRUE) - 1.5 * w1
xmax1 <- max(brks1, na.rm = TRUE) + 1.5 * w1
scl1 <- c(xmin1, xmax1)
# Scale for Histogram 2 (SAS-style: 1.5 bin-widths of padding)
w2 <- diff(brks2)[1] # bin width
xmin2 <- min(brks2, na.rm = TRUE) - 1.5 * w2
xmax2 <- max(brks2, na.rm = TRUE) + 1.5 * w2
scl2 <- c(xmin2, xmax2)
# Adjust overall scale if needed by histogram 1
if (scl1[1] < scl[1]) {
scl[1] <- scl1[1]
}
if (scl1[2] > scl[2]) {
scl[2] <- scl1[2]
}
# Adjust overall scale if needed by histogram 2
if (scl2[1] < scl[1]) {
scl[1] <- scl2[1]
}
if (scl2[2] > scl[2]) {
scl[2] <- scl2[2]
}
#******************************
#* Histogram 1
#******************************
# Calculate stats
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
n <- stat_tbl$N[as.character(stat_tbl$CLASS) == as.character(vl1)]
mu <- mean(dt1, na.rm = TRUE)
sdx <- sd(dt1, na.rm = TRUE)
# Histogram (Percent scale) - SAS uses left-closed bins [a, b)
h <- hist(dt1_hist,
breaks = brks1,
right = FALSE,
include.lowest = TRUE,
plot = FALSE)
h$counts <- h$counts / sum(h$counts, na.rm = TRUE) * 100
# Use calculated break to get normal curve
grid <- seq(scl1[1] , scl1[2], length.out = 300)
y_norm_percent <- 100 * w1 * dnorm(grid, mean = mu, sd = sdx)
y_mx <- max(max(y_norm_percent, na.rm = TRUE), max(h$counts, na.rm = TRUE), na.rm = TRUE)
# Initial chart to draw vertical orientation lines underneath bars
hist(Inf,
breaks = brks1,
right = FALSE,
include.lowest = TRUE,
main = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.025),
plot = TRUE,
axes = FALSE)
# Get axis tick marks and draw orientation lines
aval <- pretty(scl, n = 3)
abline(v = aval, col = "grey90", lwd = 1)
# Create plot using histogram values
plot(h,
col = "#CAD5E5", #"grey85",
border = "grey20",
main = "",
xlab = "",
ylab = "Percent",
xlim = scl,
ylim = c(0, y_mx * 1.025), # max(h$counts) * 1.05),
axes = FALSE,
add = TRUE)
# Normal curve overlay (scaled to percent)
lines(grid, y_norm_percent, col = "steelblue4", lwd = 2)
# Bandwidth (Silverman / SAS) - distinct n; Kish's effective n when weight only
if (!is.null(wgt1_hist) & is.null(freq1)) {
n_eff <- (sum(wgt1_hist, na.rm = TRUE))^2 / sum(wgt1_hist^2, na.rm = TRUE)
} else {
n_eff <- length(dt1)
}
bw <- 1.06 * sdx * n_eff^(-1/5)
# Kernel density estimate
dens <- density(dt1_hist, kernel = "gaussian", bw = bw, n = 512, na.rm = TRUE)
# Scale KDE to Percent axis
y_kernel_percent <- 100 * w1 * dens$y
yscl <- seq(50, length(y_kernel_percent) - 25)
# Overlay on histogram
lines(dens$x, y_kernel_percent, col = "orangered2", lwd = 2)
# Add custom axes
axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Add y label
mtext("Percent", side = 2, line = 2.5)
# Add class label
legend("topleft", legend = vl1, bty = "n",
x.intersp = 0, # Spacing between group label and box
y.intersp = 0, # Spacing between content and borders
cex = .9
)
# Frame
box(col = "grey70", lwd = 1)
#******************************
#* Histogram 2
#******************************
# Calculate stats
n <- stat_tbl$N[as.character(stat_tbl$CLASS) == as.character(vl2)]
mu <- mean(dt2, na.rm = TRUE)
sdx <- sd(dt2, na.rm = TRUE)
# Histogram (Percent scale) - SAS uses left-closed bins [a, b)
h <- hist(dt2_hist,
breaks = brks2,
right = FALSE,
include.lowest = TRUE,
plot = FALSE)
h$counts <- h$counts / sum(h$counts, na.rm = TRUE) * 100
# Use calculated break to get normal curve
grid <- seq(scl2[1] , scl2[2], length.out = 300)
y_norm_percent <- 100 * w2 * dnorm(grid, mean = mu, sd = sdx)
y_mx <- max(max(y_norm_percent, na.rm = TRUE),
max(h$counts, na.rm = TRUE), na.rm = TRUE)
# Initial chart to draw vertical orientation lines underneath bars
hist(Inf,
breaks = brks2,
right = FALSE,
include.lowest = TRUE,
main = "",
ylab = "",
xlim = scl,
ylim = c(0, y_mx * 1.025),
plot = TRUE,
axes = FALSE)
# Get axis tick marks and draw orientation lines
aval <- pretty(scl, n = 3)
abline(v = aval, col = "grey90", lwd = 1)
# Create plot using histogram values
plot(h,
col = "#CAD5E5", #"grey85",
border = "grey20",
main = "",
xlab = "",
ylab = "Percent",
xlim = scl,
ylim = c(0, y_mx * 1.025), # max(h$counts) * 1.05),
axes = FALSE,
add = TRUE)
# Normal curve overlay (scaled to percent)
lines(grid, y_norm_percent, col = "steelblue4", lwd = 2)
# Bandwidth (Silverman / SAS) - distinct n; Kish's effective n when weight only
if (!is.null(wgt2_hist) & is.null(freq2)) {
n_eff <- (sum(wgt2_hist, na.rm = TRUE))^2 / sum(wgt2_hist^2, na.rm = TRUE)
} else {
n_eff <- length(dt2)
}
bw <- 1.06 * sdx * n_eff^(-1/5) # 1.06
# Kernel density estimate
dens <- density(dt2_hist, kernel = "gaussian", bw = bw, n = 512, na.rm = TRUE)
# Scale KDE to Percent axis
y_kernel_percent <- 100 * w2 * dens$y
# Overlay on histogram
lines(dens$x, y_kernel_percent, col = "orangered2", lwd = 2)
# Add custom axes (SAS-style: fewer, rounder ticks on x)
axis(side = 1, at = pretty(scl, n = 3), col.ticks = "grey55",
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)
# Add y label
mtext("Percent", side = 2, line = 2.5)
# Add class label
legend("topleft", legend = vl2, bty = "n",
x.intersp = 0, # Spacing between group label and box
y.intersp = 0, # Spacing between content and borders
cex = .9
)
# Frame
box(col = "grey70", lwd = 1)
#******************************
#* Clean up
#******************************
box("outer", col = "grey70", lwd = 1)
# Add titles
mtext(paste0("Distribution of ", var), side = 3, line = par("oma")[3] - 1.75,
outer = TRUE, font = 2, cex = 1.25)
# Add x axis label
mtext(var, side = 1, line = par("oma")[1] - 3.5, outer = TRUE)
# Add legend
mpos <- legend("bottom",
legend = c("Normal ",
"Kernel"),
col = c("steelblue4", "orangered2"),
lwd = 2,
cex = .9,
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, -.38),
xpd = NA)
# Restore margins
par(mar = op, mfrow = c(1, 1), oma = om)
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
# Legend moves to left if mean too far to the right
#' @noRd
render_boxplot1 <- function(dat, var, plt, res) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 2.25 # 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
alph <- (1 - plt$alph) * 100 # Percentage
alpha <- plt$alph # Actual value
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Prepare data
dt <- dat[[var]]
# Set margins
par(mar = c(3, 1, 3, .75) + 0.1)
# Assign labels
xlbl <- var
tlbl <- paste(" ", var)
if (!is.null(plt$varlbl)) {
xlbl <- "Difference"
tlbl <- paste0(" Difference: ", plt$varlbl)
}
# Get scales
if (plt$showh0) {
xscl <- get_scale(dt, .05, plt$h0)
} else {
xscl <- get_scale(dt, .05)
}
# Get boxplot stats
bp <- boxplot_stats1(dt, if (!is.null(plt$freq)) dat[[plt$freq]] else NULL)
if (length(bp$out) > 0) {
if (min(bp$out, na.rm = TRUE) < xscl[1]) {
xscl[1] <- min(bp$out, na.rm = TRUE) * 1.05
}
if (max(bp$out, na.rm = TRUE) > xscl[2]) {
xscl[2] <- max(bp$out, na.rm = TRUE) * 1.05
}
}
# Get stats
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
clm_tbl <- res[[grep("ConfLimits$", names(res))[length(grep("ConfLimits$", names(res)))]]]
n <- stat_tbl$N[1]
mu <- stat_tbl$MEAN[1]
sdx <- stat_tbl$STD[1]
ci <- c(clm_tbl$LCLM[1], clm_tbl$UCLM[1])
## Draw empty plot first (for layering)
plot(dt, rep(1, length(dt)),
type = "n",
xlab = "",
ylab = "",
yaxt = "n",
main = "",
xlim = xscl,
axes = FALSE) #With 95% Confidence Interval for Mean")
# Title
mtext(paste0("Distribution of", tlbl), side = 3,
line = par("mar")[3] - 1.5,
font = 2, cex = 1.25)
# Subtitle
mtext(paste0("With ", alph, "% Confidence Interval for Mean"), side = 3,
line = par("mar")[3] - 2.5,
font = 1)
# X Label
mtext(xlbl, side = 1,
line = par("mar")[1] - 1.5,
font = 1)
## Shaded CI band
usr <- par("usr")
if (is.infinite(ci[1])) ci[1] <- usr[1] - abs(usr[1])*10
if (is.infinite(ci[2])) ci[2] <- usr[2] + abs(usr[2])*10
rect(ci[1], usr[3], ci[2], usr[4],
col = "#B3D2D0",
border = NA)
# Draw axis
aval <- axis(side = 1, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# Orientation lines
abline(v = aval, col = "grey90", lwd = 1)
# Show h0
if (plt$showh0) {
abline(v = plt$h0, col = "grey70", lwd = 2)
}
# Boxplot (horizontal)
bxp(bp,
horizontal = TRUE,
add = TRUE,
boxfill = "#CAD5E5",
border = "grey40",
lty = 1,
whiskcol = "#05379B",
staplecol = "#05379B",
medcol = "#05379B",
medlwd = 1,
boxwex = .65,
outline = FALSE,
axes = FALSE) # Already drawn above
## Mean diamond
points(mu, 1,
pch = 5,
col = "#05379B",
lwd = 1,
cex = 1.25)
# Display outliers
if (length(bp$out) > 0) {
points(bp$out, rep(1, length(bp$out)),
pch = 1,
col = "grey20",
lwd = 1.6,
cex = 1.25)
# Add labels
if (plt$label) {
if (is.null(plt$id)) {
olbls <- bp$onm
} else {
olbls <- dat[[plt$id]][as.integer(bp$onm)]
}
if (length(olbls) > 0) {
text(bp$out, rep(1, length(bp$out)),
labels = olbls,
cex = .9,
pos = 1)
}
}
}
# Legend location - Move to left side if needed
lgnd <- "topright"
pt <- (bp$stats[4] - xscl[1])/diff(xscl)
if (pt > .75) {
lgnd <- "topleft"
}
# Create legend
if (plt$showh0) {
if ( plt$h0 > mu) {
lgnd <- "topleft"
}
legend(lgnd,
legend = c(paste0(alph, "% Confidence"), "Null Value"),
pch = c(15, NA),
fill = c("#B3D2D0", NA),
col = c(NA, "grey70"),
lwd = c(NA, 2),
border = c("grey60", NA),
box.col = "grey80",
inset = c(.01, .02),
x.intersp = .5,
y.intersp = c(0.5, 1.8),
cex = .9,
bty = "o")
} else {
legend(lgnd,
legend = paste0(alph, "% Confidence"),
fill = "#B3D2D0",
border = "grey60",
box.col = "grey80",
inset = c(.01, .02),
x.intersp = .5,
y.intersp = c(0.5, 1.8),
cex = .9,
bty = "o")
}
# 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)
}
# Need to deal with outliers
#' @noRd
render_boxplot2 <- function(dat, var, plt, class, res) {
# browser()
op <- par("mar")
om <- par("oma")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 3.25 # 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
alph <- (1 - plt$alph) * 100 # Percentage
alpha <- plt$alph # Actual value
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Convert to data frame just in case it is a tibble
dat <- as.data.frame(dat)
if (!is.null(plt$id)) {
lbls <- dat[[plt$id]]
} else {
lbls <- NULL
}
# Add rownumbers
rownames(dat) <- seq(1, nrow(dat))
# Get filter variable and values
fvr <- "" # Filter variable
fvls <- c() # Filter values
if (!is.null(class)) {
fvr <- class
dat <- sort(dat, by = fvr)
fvls <- unique(dat[[class]])
fvls <- fvls[!is.na(fvls)]
} else {
stop("Something is wrong")
}
# Prepare data
dt <- dat[[var]]
rwnms <- rownames(dat)
dt1 <- dat[dat[[fvr]] == fvls[1], var]
nms1 <- rwnms[dat[[fvr]] == fvls[1]]
dt2 <- dat[dat[[fvr]] == fvls[2], var]
nms2 <- rwnms[dat[[fvr]] == fvls[2]]
dat1 <- dat[dat[[fvr]] == fvls[1], ]
dat2 <- dat[dat[[fvr]] == fvls[2], ]
# Get number of lines needed to show labels
minlns <- get_line_count(as.character(fvls)) + 1
lmg <- 4
if (minlns > 2) {
lmg <- minlns + 2
}
# Set margins
par(mar = c(4, lmg, 2, .75) + 0.1)
# Get x scale
# xscl <- get_scale(dt, .05)
if (plt$showh0) {
xscl1 <- get_scale(dt1, .05, plt$h0)
} else {
xscl1 <- get_scale(dt1, .05)
}
if (plt$showh0) {
xscl2 <- get_scale(dt2, .05, plt$h0)
} else {
xscl2 <- get_scale(dt2, .05)
}
# Get analysis variables
cvls <- unique(dat[[class]])
cvls <- cvls[!is.na(cvls)]
vl1 <- cvls[1]
vl2 <- cvls[2]
# Get combined scale
xscl <- c()
xscl[1] <- min(xscl1[1], xscl2[1], na.rm = TRUE)
xscl[2] <- max(xscl1[2], xscl2[2], na.rm = TRUE)
# Get boxplot stats for 2 plots
freq1 <- if (!is.null(plt$freq)) dat[dat[[fvr]] == fvls[1], plt$freq] else NULL
freq2 <- if (!is.null(plt$freq)) dat[dat[[fvr]] == fvls[2], plt$freq] else NULL
bp <- boxplot_stats2(dt1, dt2, nms1, nms2, fvls, freq1, freq2)
if (length(bp$out1) > 0) {
if (min(bp$out1, na.rm = TRUE) < xscl[1]) {
xscl[1] <- min(bp$out1, na.rm = TRUE) * 1.05
}
if (max(bp$out1, na.rm = TRUE) > xscl[2]) {
xscl[2] <- max(bp$out1, na.rm = TRUE) * 1.05
}
}
if (length(bp$out2) > 0) {
if (min(bp$out2, na.rm = TRUE) < xscl[1]) {
xscl[1] <- min(bp$out2, na.rm = TRUE) * 1.05
}
if (max(bp$out2, na.rm = TRUE) > xscl[2]) {
xscl[2] <- max(bp$out2, na.rm = TRUE) * 1.05
}
}
# Calculate
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
n1 <- stat_tbl$N[as.character(stat_tbl$CLASS) == as.character(vl1)]
mu1 <- stat_tbl$MEAN[as.character(stat_tbl$CLASS) == as.character(vl1)]
sdx1 <- stat_tbl$STD[as.character(stat_tbl$CLASS) == as.character(vl1)]
n2 <- stat_tbl$N[as.character(stat_tbl$CLASS) == as.character(vl2)]
mu2 <- stat_tbl$MEAN[as.character(stat_tbl$CLASS) == as.character(vl2)]
sdx2 <- stat_tbl$STD[as.character(stat_tbl$CLASS) == as.character(vl2)]
# Draw empty plot first, so we can put ablines() behind the boxes
plot(dt, rep(1, length(dt)),
xlim = xscl,
ylim = c(.5, 2.5),
type = "n",
xlab = "",
ylab = "",
yaxt = "n",
main = "",
axes = FALSE) #With 95% Confidence Interval for Mean")
# Title
mtext(paste0("Distribution of ", var), side = 3,
line = par("mar")[3] - 1.5,
font = 2, cex = 1.25)
# X Label
mtext(var, side = 1,
line = par("mar")[1] - 2,
font = 1)
# Y Label
mtext(fvr, side = 2,
line = par("mar")[2] - 2,
font = 1)
# Draw X axis
aval <- axis(side = 1, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
# # Draw Y axis
axis(side = 2, las = 1, col.ticks = "grey55", labels = fvls,
at = c(2, 1), mgp = c(3, .5, 0), tck = -0.015)
# Orientation lines
abline(v = aval, col = "grey90", lwd = 1)
# Boxplot (horizontal)
bxp(bp,
horizontal = TRUE,
add = TRUE,
boxfill = "#CAD5E5",
border = "grey40",
lty = 1,
whiskcol = "#05379B",
staplecol = "#05379B",
medcol = "#05379B",
medlwd = 1,
show.names = TRUE,
boxwex = .425,
outline = FALSE,
axes = FALSE) # Already drawn above
## Mean diamond
points(c(mu1, mu2), c(2, 1),
pch = 5,
col = "#05379B",
lwd = 1,
cex = 1.25)
# Outliers for first class value
if (length(bp$out1) > 0) {
points(bp$out1, rep(2, length(bp$out1)),
pch = 1,
col = "grey20",
lwd = 1.6,
cex = 1.25)
# Add labels
if (plt$label) {
if (is.null(plt$id)) {
olbls1 <- bp$onm1
} else {
olbls1 <- lbls[as.integer(bp$onm1)]
}
if (length(olbls1) > 0) {
text(bp$out1, rep(2, length(bp$out1)),
labels = olbls1,
cex = .9,
pos = 1)
}
}
}
# Outliers for second class value
if (length(bp$out2) > 0) {
# Plot outlier points
points(bp$out2, rep(1, length(bp$out2)),
pch = 1,
col = "grey20",
lwd = 1.6,
cex = 1.25)
# browser()
# Add labels
if (plt$label) {
if (is.null(plt$id)) {
olbls2 <- bp$onm2
} else {
olbls2 <- lbls[as.integer(bp$onm2)]
}
if (length(olbls2) > 0) {
text(bp$out2, rep(1, length(bp$out2)),
labels = olbls2,
cex = .9,
pos = 1)
}
}
}
# Frame
box(col = "grey70", lwd = 1)
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 <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
# type: pergroup or period
#' @noRd
render_interval1 <- function(dat, var, plt, res) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 2.75 # 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
alph <- (1 - plt$alph) * 100 # Percentage
alpha <- plt$alph # Actual value
# 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")
# Prepare data
dt <- dat[[var]]
# Set margins
par(mar = c(3, 1, 3, .75) + 0.1)
# Calculate basic parameters
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
clm_tbl <- res[[grep("ConfLimits$", names(res))[length(grep("ConfLimits$", names(res)))]]]
n <- stat_tbl$N[1]
mu <- stat_tbl$MEAN[1]
sdx <- stat_tbl$STD[1]
ci <- c(clm_tbl$LCLM[1], clm_tbl$UCLM[1])
# Assign labels
xlbl <- var
tlbl <- paste0(" ", var)
if (!is.null(plt$varlbl)) {
xlbl <- "Difference"
tlbl <- paste0(" Difference: ", plt$varlbl)
}
# Get scales
xscl <- ci
if (is.infinite(xscl[1])) {xscl[1] <- mu-1.5*(xscl[2]-mu)}
if (is.infinite(xscl[2])) {xscl[2] <- mu+1.5*(mu-xscl[1])}
if (plt$showh0) {
xscl <- get_scale(xscl, .001, plt$h0)
} else {
xscl <- get_scale(xscl, .001)
}
# Create plot
plot(ci, rep(1, length(ci)),
type = "n",
yaxt = "n",
ylab = "",
xlab = "",
xlim = xscl,
main = "",
axes = FALSE)
# Title
mtext(paste0("Mean of ", tlbl), side = 3,
line = par("mar")[3] - 1.5,
font = 2, cex = 1.25)
# Subtitle
mtext(paste0("With ", alph, "% Confidence Interval"), side = 3,
line = par("mar")[3] - 2.5,
font = 1)
# X Label
mtext(xlbl, side = 1,
line = par("mar")[1] - 1.5,
font = 1)
# Draw axis
aval <- axis(side = 1, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
## CI line
if (is.infinite(ci[2])) {
# Upper is Inf: Draw arrow pointing right to the edge of the plot
arrows(mu, 1, xscl[2]/1.001, 1, col = "#05379B", lwd = 2, length = 0.1)
segments(ci[1], 1, mu, 1, col = "firebrick3", lwd = 2)
} else if (is.infinite(ci[1])) {
# Lower is -Inf: Draw arrow pointing left to the edge of the plot
arrows(ci[2], 1, xscl[1]/0.999, 1, col = "#05379B", lwd = 2, length = 0.1)
segments(ci[2], 1, mu, 1, col = "firebrick3", lwd = 2)
} else {
# Finite: Draw normal segment
segments(ci[1], 1, ci[2], 1, col = "firebrick3", lwd = 2)
}
## CI end caps
if (is.finite(ci[1])) segments(ci[1], 0.95, ci[1], 1.05, col = "firebrick3", lwd = 2)
if (is.finite(ci[2])) segments(ci[2], 0.95, ci[2], 1.05, col = "firebrick3", lwd = 2)
## Mean diamond
points(mu, 1,
pch = 5,
col = "#05379B",
cex = 1.25)
# Show h0 line
if (plt$showh0) {
abline(v = plt$h0, col = "grey70", lwd = 2)
if (plt$h0 > mu) {
lgnd <- "topleft"
} else {
lgnd <- "topright"
}
legend(lgnd,
horiz = TRUE,
legend = c("Mean", "Null Value"),
pch = c(5, NA),
col = c("#05379B", "grey70"),
lwd = c(NA, 2),
box.col = "grey80",
inset = c(.01, .02),
x.intersp = .5,
y.intersp = c(0.5, 1.8),
pt.cex = 1,
cex = .9,
bty = "o")
} else {
# Create legend
legend("topright",
legend = "Mean",
pch = 5,
col = "#05379B",
box.col = "grey80",
inset = c(.01, .02),
x.intersp = .5,
y.intersp = c(0.5, 1.8),
cex = .9,
pt.cex = 1,
bty = "o")
}
# 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_interval2 <- function(dat, var, plt, res) {
op <- par("mar")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 2.75 # 4.5 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
alph <- (1 - plt$alph) * 100 # Percentage
alpha <- plt$alph # Actual value
# 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")
# Prepare data
dt <- res[[grep("ConfLimits$", names(res))[length(grep("ConfLimits$", names(res)))]]]
# Set margins
par(mar = c(4, 1, 3, .75) + 0.1)
# Get class label
clslbl <- paste0("(", dt$CLASS[1], " - ", dt$CLASS[2], ")")
# Get labels for means
labels <- c("Satterthwaite", "Pooled")
# Pull off data values
mean_diff <- c(dt$MEAN[4], dt$MEAN[3])
lower_ci <- c(dt$LCLM[4], dt$LCLM[3])
upper_ci <- c(dt$UCLM[4], dt$UCLM[3])
# Create y and x values
y <- c(2, 1) # vertical positions
ci <- c(lower_ci, upper_ci)
# Get scales - seed from finite CI bounds plus the means so one-sided tests
# (where all upper_ci or all lower_ci are Inf) produce a readable axis range.
xscl <- c(min(lower_ci, na.rm = TRUE), max(upper_ci, na.rm = TRUE))
mean_mid <- mean(mean_diff, na.rm = TRUE)
if (is.infinite(xscl[1])) { xscl[1] <- mean_mid - 1.5 * (xscl[2] - mean_mid) }
if (is.infinite(xscl[2])) { xscl[2] <- mean_mid + 1.5 * (mean_mid - xscl[1]) }
xscl <- c(min(c(xscl, mean_diff), na.rm = TRUE), max(c(xscl, mean_diff, na.rm = TRUE), na.rm = TRUE))
if (plt$showh0) {
xscl <- get_scale(xscl, .001, plt$h0)
} else {
xscl <- get_scale(xscl, .001)
}
# Create base plot
plot(mean_diff, y,
xlim = xscl,
ylim = c(0.5, 2.5),
xlab = "",
ylab = "",
yaxt = "n",
pch = NA,
main = "",
axes = FALSE)
# Title
mtext(paste("Mean of ", var, " Difference ", clslbl), side = 3,
line = 1.5, cex = 1.25, font = 2)
## Subtitle
sides_lbl <- if (any(is.infinite(upper_ci))) "Upper " else if (any(is.infinite(lower_ci))) "Lower " else ""
mtext(paste0("With ", alph, "% ", sides_lbl, "Confidence Intervals"), side = 3,
line = .5, cex = 1)
# X axis label
mtext("Difference", side = 1, line = 2)
# Draw axis
aval <- axis(side = 1, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
for (k in seq_along(y)) {
if (is.infinite(upper_ci[k])) {
arrows(mean_diff[k], y[k], xscl[2] / 1.001, y[k], col = "#05379B", lwd = 3, length = 0.1)
segments(lower_ci[k], y[k], mean_diff[k], y[k], col = "brown3", lwd = 3)
} else if (is.infinite(lower_ci[k])) {
arrows(mean_diff[k], y[k], xscl[1] / 0.999, y[k], col = "#05379B", lwd = 3, length = 0.1)
segments(upper_ci[k], y[k], mean_diff[k], y[k], col = "brown3", lwd = 3)
} else {
segments(lower_ci[k], y[k], upper_ci[k], y[k], col = "brown3", lwd = 3)
}
if (is.finite(lower_ci[k])) segments(lower_ci[k], y[k] - 0.08, lower_ci[k], y[k] + 0.08, col = "brown3", lwd = 3)
if (is.finite(upper_ci[k])) segments(upper_ci[k], y[k] - 0.08, upper_ci[k], y[k] + 0.08, col = "brown3", lwd = 3)
}
## Mean diamonds
points(mean_diff, y,
pch = 5, # diamond
col = "#05379B",
cex = 1.25)
## Labels above lines
text(mean_diff, y + 0.15, labels,
cex = 0.9)
if (plt$showh0) {
abline(v = plt$h0, col = "grey70", lwd = 2)
if (plt$h0 > mean(ci, na.rm = TRUE)) {
lgnd <- "topleft"
} else {
lgnd <- "topright"
}
legend(lgnd,
horiz = TRUE,
legend = c("Mean", "Null Value"),
pch = c(5, NA),
col = c("#05379B", "grey70"),
lwd = c(NA, 2),
box.col = "grey80",
inset = c(.01, .02),
x.intersp = .5,
y.intersp = c(0.5, 1.8),
pt.cex = 1,
cex = .9,
bty = "o")
} else {
legend("topright",
legend = "Mean",
pch = 5,
col = "#05379B",
box.col = "grey80",
inset = c(.01, .02),
x.intersp = .5,
y.intersp = c(0.5, 1.8),
cex = .9,
bty = "o")
}
# 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)
}
#' @noRd
render_tqqplot1 <- function(dat, var, plt, res) {
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")
# Extract data
rdt <- dat[[var]]
freq <- if (!is.null(plt$freq)) dat[[plt$freq]] else NULL
if (!is.null(freq)) {
freq <- floor(freq)
idx <- !is.na(rdt) & !is.na(freq) & freq > 0
rdt <- rep(rdt[idx], times = freq[idx])
}
# Assign label
ylbl <- var
tlbl <- paste(" ", var)
if (!is.null(plt$varlbl)) {
ylbl <- "Difference"
tlbl <- paste0(": ", plt$varlbl)
}
# 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", tlbl),
xlab = "",
ylab = ylbl,
pch = 1,
col = "#05379B",
cex = 1.3,
axes = FALSE)
# Get stats table
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
# Slope = Standard Deviation, Intercept = Mean
if (!is.null(plt$weight)) {
# For some reason, the qqplot trend line uses the unweighted values.
# So have to calculate them here.
ovar <- dat[[var]]
owgt <- dat[[plt$weight]]
adjmean <- mean(ovar, na.rm = TRUE)
adjstd <- sd(ovar, na.rm = TRUE)
abline(a = adjmean, b = adjstd, col = "grey60")
} else {
abline(a = stat_tbl$MEAN[1], b = stat_tbl$STD[1], 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_tqqplot1_back <- function(dat, var, 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
# rdt <- dat[[var]]
#
# # Assign label
# ylbl <- var
# tlbl <- paste(" ", var)
# if (!is.null(plt$varlbl)) {
# ylbl <- "Difference"
# tlbl <- paste0(": ", plt$varlbl)
# }
#
# # Set margins
# par(mar = c(4, 5, 2, .75) + 0.1)
#
# # Get y scale
# # mx <- max(abs(rdt)) * 1.125 #1.0025
# # scl <- c(-mx, mx)
#
# # Draw plot
# qqnorm(rdt,
# main = paste0("Q-Q Plot of", tlbl),
# xlab = "",
# ylab = ylbl,
# pch = 1, # open circles
# col = "#05379B",
# cex = 1.3,
# # asp = 1,
# # ylim = scl,
# # asp = 1,
# axes = FALSE)
#
# # Add diagonal line
# # This line is a problem. Right now just putting it diagonal
# # between the corners. All qqline() functions are worse.
# usr <- par("usr")
# segments(usr[1], usr[3], usr[2], usr[4], col = "grey60")
# # qqline(rdt, col = "grey60", lwd = 1, qtype = 3)
#
# # 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)
#
#
# }
#' @noRd
render_tqqplot2 <- function(dat, var, plt, class, res) {
op <- par("mar")
om <- par("oma")
# Create temp file path
pth <- tempfile(fileext = ".jpg")
# Standard height and width
hti <- 3.25 # Height in inches
wdi <- 6 # Width in inches
bml <- 6 # bottom margin lines
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Prepare data
dat <- as.data.frame(dat)
dat <- sort(dat, by = class)
# Subset
cvls <- unique(dat[[class]])
dt1 <- dat[dat[[class]] == cvls[1], var]
dt2 <- dat[dat[[class]] == cvls[2], var]
freq1 <- if (!is.null(plt$freq)) dat[dat[[class]] == cvls[1], plt$freq] else NULL
freq2 <- if (!is.null(plt$freq)) dat[dat[[class]] == cvls[2], plt$freq] else NULL
if (!is.null(freq1)) {
freq1 <- floor(freq1)
idx1 <- !is.na(dt1) & !is.na(freq1) & freq1 > 0
dt1 <- rep(dt1[idx1], times = freq1[idx1])
}
if (!is.null(freq2)) {
freq2 <- floor(freq2)
idx2 <- !is.na(dt2) & !is.na(freq2) & freq2 > 0
dt2 <- rep(dt2[idx2], times = freq2[idx2])
}
# Output to image file
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Set margins
par(oma = c(0, 0, 2, .75) + 0.1,
mar = c(4, 4, 0, .0) + 0.1,
mfrow = c(1, 2))
#*******************************
#* QQPlot 1
#*******************************
# Prepare Y data
ydat <- sort(dt1)
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,
axes = FALSE)
# Get stats table
stat_tbl <- res[[grep("Statistics$", names(res))[length(grep("Statistics$", names(res)))]]]
# Slope = Standard Deviation, Intercept = Mean
if (!is.null(plt$weight)) {
# For some reason, the qqplot trend line uses the unweighted values.
# So have to calculate them here.
ovar <- dt1
adjmean <- mean(ovar, na.rm = TRUE)
adjstd <- sd(ovar, na.rm = TRUE)
abline(a = adjmean, b = adjstd, col = "grey60")
} else {
abline(a = stat_tbl$MEAN[as.character(stat_tbl$CLASS) == as.character(cvls[1])],
b = stat_tbl$STD[as.character(stat_tbl$CLASS) == as.character(cvls[1])], 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)
# Y Axis label
mtext(var, side = 2, line = par("mar")[2] - 1.5)
# Add class label
legend("topleft", legend = cvls[1], bty = "n",
x.intersp = 0, # Spacing between group label and box
y.intersp = 0, # Spacing between content and borders
cex = .9
)
# Frame
box(col = "grey70", lwd = 1)
#*******************************
#* QQPlot 2
#*******************************
# Prepare Y data
ydat <- sort(dt2)
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,
axes = FALSE)
# Slope = Standard Deviation, Intercept = Mean
if (!is.null(plt$weight)) {
# For some reason, the qqplot trend line uses the unweighted values.
# So have to calculate them here.
ovar <- dt2
adjmean <- mean(ovar, na.rm = TRUE)
adjstd <- sd(ovar, na.rm = TRUE)
abline(a = adjmean, b = adjstd, col = "grey60")
} else {
abline(a = stat_tbl$MEAN[as.character(stat_tbl$CLASS) == as.character(cvls[2])],
b = stat_tbl$STD[as.character(stat_tbl$CLASS) == as.character(cvls[2])], 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)
# Y Axis label
mtext(var, side = 2, line = par("mar")[2] - 1.5)
# Add class label
legend("topleft", legend = cvls[2], bty = "n",
x.intersp = 0, # Spacing between group label and box
y.intersp = 0, # Spacing between content and borders
cex = .9
)
# Frame
box(col = "grey70", lwd = 1)
#*******************************
#* Clean up
#*******************************
# Title
mtext(paste0("Q-Q Plots of ", var), side = 3, outer = TRUE,
line = par("oma")[3] - 1.75, font = 2, cex = 1.25)
box("outer", col = "grey70", lwd = 1)
par(mar = op, oma = om, mfrow = c(1, 1))
# Close device context
dev.off()
# Put plot in reporter plot object
ret <- create_plot(pth, height = hti, width = wdi)
return(ret)
}
#' @noRd
render_profiles <- function(dat, var, plt, res) {
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
alph <- (1 - plt$alph) * 100 # Percentage
alpha <- plt$alph # Actual value
# 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")
# Prepare data
splt <- trimws(strsplit(plt$varlbl, "-", fixed = TRUE)[[1]])
v1 <- splt[1]
v2 <- splt[2]
vrs <- c(v1, v2)
dt <- as.matrix(dat[ , vrs ])
## X positions (categorical)
x <- c(1, 2)
par(mar = c(3, 3, 3, 3) + 0.1)
## Empty plot frame
plot(x, range(dt, na.rm = TRUE),
type = "n",
xaxt = "n",
xlab = "",
ylab = "",
main = paste0("Paired Profiles for (", v1, ", ", v2, ")"),
axes = FALSE)
# X Axis
axis(1, at = x, labels = c(v1, v2))
# Y Axis
axis(side = 2, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
axis(side = 4, las = 1, col.ticks = "grey55",
mgp = c(3, .5, 0), tck = -0.015)
## Individual subject profiles (thin lines)
for (i in seq_len(nrow(dt))) {
lines(x, dt[i, ],
col = rgb(0.2, 0.4, 0.8, 0.45), # light blue
lwd = 1)
}
## Mean profile (bold line)
mean_profile <- colMeans(dt)
lines(x, mean_profile,
col = "orangered2",
lwd = 3)
## Mean legend
legend("top",
legend = "Mean",
col = "orangered2",
lwd = 3,
box.col = "grey80",
inset = c(.01, .02),
x.intersp = .5,
y.intersp = c(0.5, 1.8),
cex = .9,
bty = "o")
# 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)
}
#' @noRd
render_agreement <- function(dat, var, plt, res) {
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
alph <- (1 - plt$alph) * 100 # Percentage
alpha <- plt$alph # Actual value
# Convert inches to pixels
ht <- hti * 96
wd <- wdi * 96
# Output to image file
jpeg(pth, width = wd, height = ht, quality = 100, units = "px")
# Prepare data
splt <- trimws(strsplit(plt$varlbl, "-", fixed = TRUE)[[1]])
v1 <- splt[1]
v2 <- splt[2]
vrs <- c(v1, v2)
# Get vectors
dt1 <- dat[[v1]]
dt2 <- dat[[v2]]
## Means
m1 <- mean(dt1, na.rm = TRUE)
m2 <- mean(dt2, na.rm = TRUE)
## Axis limits
lims <- range(c(dt1, dt2), na.rm = TRUE)
# Set margins
par(mar = c(4, 10, 2, 6) + 0.1)
## Scatter plot
plot(dt1, dt2,
xlim = lims,
ylim = lims,
pch = 1,
cex = 1.25,
col = "#05379B",
xlab = "",
ylab = v2,
main = paste0("Agreement of ", v2, " and ", v1),
axes = FALSE,
asp = 1) # equal scaling (critical)
# X label
mtext(v1, side = 1, line = par("mar")[1] - 2)
# Axes
axis(side = 1, las = 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)
## Identity (45-degree) line
usr <- par("usr")
segments(usr[1], usr[3], usr[2], usr[4],
col = "grey60", lwd = 1)
## Mean point
points(m1, m2,
pch = 21,
bg = "white",
col = "#05379B",
lwd = 2,
cex = 2.5)
## Legend
legend("right",
legend = "Mean",
pch = 21,
pt.bg = "white",
col = "#05379B",
pt.lwd = 2,
# lwd = 3,
box.col = "grey80",
inset = c(-.2, 0),
x.intersp = 1,
y.intersp = c(0.5, 1.8),
pt.cex = 2.5,
xpd = NA,
bty = "o")
# 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)
}
# Utilities ---------------------------------------------------------------
# Boxplot.stats(data1) can prepare the boxplot data more easily,
# but has no parameter to use SAS-style quantiles. So need to calculate
# everything manually.
#' @noRd
boxplot_stats1 <- function(dt1, freq = NULL) {
# Apply freq expansion if specified, tracking original row indices
if (!is.null(freq)) {
f <- floor(freq)
valid <- !is.na(dt1) & !is.na(f) & f > 0
row_map <- rep(which(valid), times = f[valid])
dt1 <- rep(dt1[valid], times = f[valid])
} else {
row_map <- seq_along(dt1)
}
# Calculate quantile
q1 <- quantile(dt1, probs = c(0, .25, .5, .75, 1), type = 2, na.rm = TRUE)
# Calculate Inter-Quartile Range
iqr1 <- IQR(dt1, type = 2, na.rm = TRUE) * 1.5
# Calculate outlier limits
lm1 <- c(q1[2] - iqr1, q1[4] + iqr1)
# Determine if there are outliers
ol1 <- dt1 <= lm1[1] | dt1 >= lm1[2]
ol1 <- ifelse(is.na(ol1), FALSE, ol1)
out1 <- dt1[ol1]
# Change min and max to remove outliers
onm1 <- c()
if (length(out1) > 0) {
mn1 <- min(dt1[dt1 > lm1[1]], na.rm = TRUE)
mx1 <- max(dt1[dt1 < lm1[2]], na.rm = TRUE)
q1[1] <- mn1
q1[5] <- mx1
onm1 <- row_map[seq(1, length(dt1))[ol1]] # Obs within class: This matches SAS
# onm1 <- nms1[ol1] # Original obs #: This seems more correct
}
# Get length
lndt <- length(dt1)
# Prepare to pass data
bp <- list(
stats = matrix(q1, ncol = 1),
n = lndt,
conf = NULL,
out = out1,
onm = onm1
)
return(bp)
}
# Boxplot.stats(data1) can prepare the boxplot data more easily,
# but has no parameter to use SAS-style quantiles. So need to calculate
# everything manually.
#' @noRd
boxplot_stats2 <- function(dt1, dt2, nms1, nms2, fvls, freq1, freq2) {
# browser()
if (!is.null(freq1)) {
f1 <- floor(freq1)
valid1 <- !is.na(dt1) & !is.na(f1) & f1 > 0
nms1 <- rep(nms1[valid1], times = f1[valid1])
dt1 <- rep(dt1[valid1], times = f1[valid1])
}
if (!is.null(freq2)) {
f2 <- floor(freq2)
valid2 <- !is.na(dt2) & !is.na(f2) & f2 > 0
nms2 <- rep(nms2[valid2], times = f2[valid2])
dt2 <- rep(dt2[valid2], times = f2[valid2])
}
# Calculate quantile
q1 <- quantile(dt1, probs = c(0, .25, .5, .75, 1), type = 2, na.rm = TRUE)
q2 <- quantile(dt2, probs = c(0, .25, .5, .75, 1), type = 2, na.rm = TRUE)
# Calculate Inter-Quartile Range
iqr1 <- IQR(dt1, type = 2, na.rm = TRUE) * 1.5
iqr2 <- IQR(dt2, type = 2, na.rm = TRUE) * 1.5
# Calculate outlier limits
lm1 <- c(q1[2] - iqr1, q1[4] + iqr1)
lm2 <- c(q2[2] - iqr2, q2[4] + iqr2)
# Determine if there are outliers
ol1 <- dt1 <= lm1[1] | dt1 >= lm1[2]
ol1 <- ifelse(is.na(ol1), FALSE, ol1)
out1 <- dt1[ol1]
ol2 <- dt2 <= lm2[1] | dt2 >= lm2[2]
ol2 <- ifelse(is.na(ol2), FALSE, ol2)
out2 <- dt2[ol2]
# Change min and max to remove outliers
onm1 <- c()
if (length(out1) > 0) {
mn1 <- min(dt1[dt1 > lm1[1]], na.rm = TRUE)
mx1 <- max(dt1[dt1 < lm1[2]], na.rm = TRUE)
q1[1] <- mn1
q1[5] <- mx1
onm1 <- nms1[ol1] # Original obs #: This seems more correct
}
# Change min and max to remove outliers
onm2 <- c()
if (length(out2) > 0) {
mn2 <- min(dt1[dt2 > lm2[1]], na.rm = TRUE)
mx2 <- max(dt1[dt2 < lm2[2]], na.rm = TRUE)
q2[1] <- mn2
q2[5] <- mx2
onm2 <- nms2[ol2] # Original obs #: This seems more correct
}
# Get length
lndt <- length(dt1) + length(dt2)
# Prepare to pass data
bp <- list(
stats = cbind(q2, q1),
n = cbind(lndt, lndt),
conf = NULL,
out = numeric(0),
names = fvls,
"out1" = out1,
"out2" = out2,
"onm1" = onm1,
"onm2" = onm2
)
return(bp)
}
#' @noRd
get_sas_bins <- function(x) {
x <- x[!is.na(x)]
rng <- range(x, na.rm = TRUE)
n <- length(x)
n_bins <- ceiling((2 * n)^(1/3))
raw_w <- diff(rng) / n_bins
if (!is.finite(raw_w) || raw_w <= 0) {
return(pretty(rng, n = n_bins, min.n = 1))
}
mag <- 10^floor(log10(raw_w))
cand <- c(1, 2, 2.5, 3, 4, 5) * mag
# SAS rounds DOWN to the nearest nice candidate
le <- cand[cand <= raw_w]
nice_w <- if (length(le) > 0) max(le, na.rm = TRUE) else min(cand, na.rm = TRUE)
# SAS picks edge-aligned vs midpoint-centered bins based on whether
# X_min already lies on the nice_w grid.
m <- floor(rng[1] / nice_w) * nice_w
on_grid <- isTRUE(all.equal(rng[1] - m, 0))
if (on_grid) {
lo <- m - nice_w / 2
} else {
anchor <- 10 * mag
lo <- floor(rng[1] / anchor) * anchor
while (lo + nice_w <= rng[1]) lo <- lo + nice_w
while (lo > rng[1]) lo <- lo - nice_w
}
hi <- lo
while (hi < rng[2]) hi <- hi + nice_w
seq(lo, hi, by = nice_w)
}
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.