#' Plot logistic regression
#'
#' Plot combined graphs for logistic regressions
#'
#' @param independ explanatory variable
#' @param depend dependent variable, typically a logical vector
#' @param logi.mod type of fitting, 1 = logistic; 2 = "gaussian" logistic
#' @param type type of representation, "dit" = dit plot; "hist" = histogram
#' @param boxp TRUE = with box plots, FALSE = without
#' @param rug TRUE = with rug plots, FALSE = without
#' @param ylabel y-axis label
#' @param ylabel2 2nd y-axis label
#' @param xlabel x-axix label
#' @param mainlabel overall title for plot
#' @param las.h orientation of axes labels (0 = vertical, 1 = horizontal
#' @param counts add counts above histogram bars
#' @param \dots additional options passed to logi.hist
#'
#' @return A logistic regression plot
#'
#' @references de la Cruz Rot, M. 2005. Improving the Presentation of Results of
#' Logistic Regression with R. ESA Bulletin 86:41-48.
#' \url{https://esapubs.org/bulletin/backissues/086-1/bulletinjan2005.htm}
#'
#' @author M. de la Cruz Rot
#'
#' @examples
#' aq.trans$survived <- aq.trans$fate!="dead"
#' a <- subset(aq.trans, leaf<50 & stage!="recruit", c(leaf,survived))
#' logi.hist.plot(a$leaf, a$survived,
#' type="hist", boxp=FALSE, counts=TRUE, int=10,
#' ylabel="Survival probability", ylabel2="Number of plants",
#' xlab="Number of leaves")
#' b <- glm(survived ~ leaf, binomial, data=a)
#' summary(b)
#'
#' @export
logi.hist.plot <- function(independ, depend, logi.mod = 1,
type = "dit", boxp = TRUE, rug = FALSE,
ylabel = "Probability", ylabel2 = "Frequency", xlabel = "", mainlabel = "",
las.h = 1, counts = FALSE, ...) {
# get the label for the x-axis
# xlabel <- paste(deparse(substitute(independ)))
# define functions:
# set the draw area if no box plots are to be drawn
logi.scater <- function(independ, depend, scater = "n",
x.lab = xlabel, las = las.h) {
plot(independ, depend,
cex = 1, type = scater,
ylab = ylabel, xlab = x.lab, main = mainlabel,
cex.lab = 1.2, las = las
)
}
# add rug plot if desired; you could change pch.rug
# (symbol type) or cex.rug (symbol size)
logi.rug <- function(independ, depend, pch.rug = 16,
cex.rug = 1) {
points(independ, depend, pch = pch.rug, cex = cex.rug)
}
# set the draw area and add box plots; you could change
# cold.box (color of the boxes)
logi.box <- function(independ, depend, col.box = "gray",
x.lab = xlabel, las = las.h) {
plot(independ, depend,
cex = 1, type = "n",
ylim = c(-0.1, 1.1), ylab = ylabel, # change ylabel2 to ylabel Nov 2, 2015
xlab = x.lab, cex.lab = 1.2, las = las
)
indep.1 <- independ[depend == 1]
indep.0 <- independ[depend == 0]
boxplot(indep.1,
horizontal = TRUE, add = TRUE,
at = 1.05, boxwex = 0.1, col = col.box, notch = TRUE
)
boxplot(indep.0,
horizontal = TRUE, add = TRUE,
at = -0.05, boxwex = 0.1, col = col.box, notch = TRUE
)
}
# fit binomial glm and add predicted curve; you could
# change col.cur (color of the curve) or lwd.cur(width
# of the curve)
logi.curve <- function(independ, depend, mod = logi.mod,
col.cur = "red", lwd.cur = 4) {
if (mod == 1) {
mod3 <- glm(depend ~ independ,
family = binomial
)
}
if (mod == 2) {
mod3 <- glm(depend ~ independ +
I(independ^2), family = binomial)
}
x.new <- seq(min(independ), max(independ), len = 100)
y.new <- predict(mod3, data.frame(independ = x.new),
type = "response"
)
lines(x.new, y.new, lwd = lwd.cur, col = col.cur)
}
# add dit plot; you may want to change pch.dit (type of
# points), cex.p (size of points), and incre (space
# between points)
logi.dit <- function(independ, depend, cex.p = 1,
pch.dit = 1, incre = 0.02) {
indep.0 <- independ[depend == 0]
indep.1 <- independ[depend == 1]
uni.plot.0 <- function(x) length(which(indep.0 == x))
uni.plot.1 <- function(x) length(which(indep.1 == x))
# get the number of repeated values of "independ":
cosa.0 <- apply(as.matrix(unique(indep.0)), 1, uni.plot.0)
cosa.1 <- apply(as.matrix(unique(indep.1)), 1, uni.plot.1)
# start ploting:
points(independ, depend, pch = pch.dit, cex = cex.p)
for (i in 1:max(cosa.0)) {
for (j in 1:i) {
points(unique(indep.0)[which(cosa.0 == i + 1)],
rep(0 + incre * j, length(which(cosa.0 == i + 1))),
pch = pch.dit, cex = cex.p
)
}
}
for (i in 1:max(cosa.1)) {
for (j in 1:i) {
points(unique(indep.1)[which(cosa.1 == i + 1)],
rep(1 - incre * j, length(which(cosa.1 == i + 1))),
pch = pch.dit, cex = cex.p
)
}
}
}
# add histograms and frequency axes; you may want to change
# scale.hist (factor to scale histogram height to 0-1
# interval) or col.hist (color of histogram)
logi.hist <- function(independ, depend, scale.hist = 5,
col.hist = "blue", count.hist = TRUE,
intervalo = 0, las.h1 = las.h) {
# get the position of bins
h.br <- hist(independ, plot = FALSE)$br
if (intervalo > 0) {
h.br <- seq(
from = range(h.br)[1],
to = range(h.br)[2], by = intervalo
)
}
h.x <- hist(independ[depend == 0],
breaks = h.br,
plot = FALSE
)$mid
# get counts in each bin
h.y0 <- hist(independ[depend == 0],
breaks = h.br,
plot = FALSE
)$counts
h.y1 <- hist(independ[depend == 1],
breaks = h.br,
plot = FALSE
)$counts
# scale the histogram bars to max desired length:
h.y0n <- h.y0 / (max(c(h.y0, h.y1)) * scale.hist)
h.y1n <- 1 - h.y1 / (max(c(h.y0, h.y1)) * scale.hist)
# draw bottom histogram:
for (i in 1:length(h.y0n)) {
if (h.y0n[i] > 0) {
polygon(c(rep(h.br[i], 2), rep(h.br[i + 1], 2)),
c(0, rep(h.y0n[i], 2), 0),
col = col.hist
)
}
}
# draw top histogram:
for (i in 1:length(h.y1n)) {
if (h.y1n[i] < 1) {
polygon(c(rep(h.br[i], 2), rep(h.br[i + 1], 2)),
c(h.y1n[i], 1, 1, h.y1n[i]),
col = col.hist
)
}
}
# add counts to bins if required:
if (counts == TRUE) {
for (i in 1:length(h.x)) {
text(h.x[i], h.y1n[i], h.y1[i], cex = 1, pos = 1)
text(h.x[i], h.y0n[i], h.y0[i], cex = 1, pos = 3)
}
}
# plot the axes of histograms:
axis.hist <- function(h.y0, h.y1, scale.hist,
las = las.h1) {
tope <- max(c(h.y0, h.y1))
label.down <- c(
0, (ceiling(tope / 10)) * 5,
(ceiling(tope / 10)) * 10
)
label.up <- c(
(ceiling(tope / 10)) * 10,
(ceiling(tope / 10)) * 5, 0
)
at.down <- label.down / (tope * scale.hist)
at.up <- 1 - (label.up / (tope * scale.hist))
at.hist <- c(at.down, at.up)
label.hist <- c(label.down, label.up)
axis(
side = 4, at = at.hist, labels = label.hist,
las = las
)
mtext(ylabel2, side = 4, line = 2, cex = 1.2)
}
axis.hist(h.y0, h.y1, scale.hist)
axis(side = 2, las = las.h1)
}
# set the margins of plot area
old.mar <- par()$mar
par(mar = c(5.1, 4.1, 4.1, 4.1))
# plot the combined graph
if (boxp == TRUE) logi.box(independ, depend)
if (boxp == FALSE) logi.scater(independ, depend)
if (type != "dit") logi.hist(independ, depend, ...)
if (rug == TRUE) logi.rug(independ, depend)
logi.curve(independ, depend)
if (type == "dit") logi.dit(independ, depend)
# reset the margins to old margins
par(mar = old.mar)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.