## confint
##
##' Modified version of the \code{\link[stats]{confint}} function,
##' which displays the coefficients in addition to the CIs, and allows
##' for more control on display parameters. A \code{plot} argument and
##' function allow to graph the coefficients and their CIs.
##'
##' @title Confidence Intervals for Model Parameters
##' @param order Logical. If \code{TRUE}, the results are ordered by
##' descending order on the coefficient value.
##' @param groups A \code{factor} in the sense that
##' \code{as.factor(f)} defines the groups, or a list of such factors
##' in which case their interaction is used for the groups. See
##' \code{\link[base]{split}}.
##' @param plot Whether to plot the results.
##' @param ... Further arguments passed to \code{points}.
##' @return A data frame providing the CI and coefficients.
##' @seealso \code{\link[stats]{confint}} for more details on other
##' parameters.
##' @author Mathieu Basille \email{basille@@ufl.edu}
##' @export
##' @examples
##' ## Example of linear model
##' fit <- lm(100/mpg ~ disp + hp + wt + am, data = mtcars)
##' ## Standard 'confint' function
##' stats::confint(fit)
##' ## Same results with modified function
##' confint(fit)
##' ## Argument 'level'
##' stats::confint(fit, level = .9)
##' confint(fit, level = .9)
##' ## Argument 'order'
##' confint(fit, order = TRUE)
##' ## Argument 'groups'
##' confint(fit, groups = c(3, 1, 1, 1, 2))
##' ## Argument 'level', "'order' and 'groups' simultaneously
##' confint(fit, level = .9, order = TRUE, groups = c(3, 1, 1, 1, 2))
##' ## Argument 'parm'
##' stats::confint(fit, "am")
##' confint(fit, "am")
##' ##'
##' ## Plot of the results
##' plot(confint(fit, order = TRUE, groups = c(3, 1, 1, 1, 2)))
##' confint(fit, order = TRUE, groups = c(3, 1, 1, 1, 2), plot = TRUE)
##' confint(fit, order = TRUE, groups = c(3, 1, 1, 1, 2), plot = TRUE,
##' col = c("blue", "red", "green"), pch = 18, cex = 2)
##' confint(fit, order = TRUE, groups = c(3, 1, 1, 1, 2), level = 0.9,
##' plot = TRUE, add.signif = TRUE)
confint <- function(object, parm, level = 0.95, order = FALSE,
groups, plot = FALSE, ...) {
## Computing the CI using 'stats::confint'
ci <- stats::confint(object = object, parm = parm, level = level,
...)
## Retrieving row and colum names
rnci <- rownames(ci)
cnci <- colnames(ci)
## Include the coefficients, and give it the name 'coef'
ci <- data.frame(ci, coef(object)[parm])
names(ci) <- c(cnci, "coef")
## If no 'groups', just check 'order'
if (missing(groups)) {
## And if 'order', simply order the result by 'coef'
if (order)
ci <- ci[order(ci$coef, decreasing = TRUE), ]
}
## If 'groups'...
else {
## ... and if 'order'
if (order) {
## Split the result according to the groups
sci <- split(ci, groups)
## For each element, sort according to 'coef', and adjust
## the rownames; bind all element outputs together
ci <- do.call("rbind", lapply(names(sci),
function(x) {
tmp <- sci[[x]][order(sci[[x]]$coef, decreasing = TRUE), ]
rownames(tmp) <- paste(x, rownames(tmp), sep = ".")
return(tmp)
}))
}
## ... and no 'order'
else {
## Split the result according to the groups, and bind
## all element outputs together
ci <- do.call("rbind", split(ci, groups))
## Adjust the rownames
rownames(ci) <- paste(groups, rnci, sep = ".")[order(groups)]
}
## Set the attribute 'groups'
attr(ci, "groups") <- groups
}
## Set the class
class(ci) <- c("confint", "data.frame")
## Set the attribute 'level'
attr(ci, "level") <- level
## Set the attribute 'model'
attr(ci, "model") <- as.character(substitute(object))
## If plot, call 'plot.confint'
if (plot)
plot(ci, ...)
return(ci)
}
## plot.confint
##
##' @param x A \code{data.frame} of class \code{confint}.
##' @param mar The number of lines of margin, can be useful if the
##' coefficient names do not fit in the left margin. See
##' \code{\link[graphics]{par}} for more details.
##' @param col The color of each coefficient + CI; gray by default. If
##' \code{"groups"}, the color of each (sorted) group; use a hcl
##' palette by default.
##' @param main The title of the plot.
##' @param pch The symbol to be used for the points. See
##' \code{\link[graphics]{par}} for more details.
##' @param add.signif Logical, whether to add an asterik for variables with
##' CIs non-overlapping with zero.
##' @param ... Further arguments passed to \code{points}.
##' @rdname confint
##' @export
plot.confint <- function(x, mar = c(5, 7, 3, 1) + 0.1, col = NULL,
main = attr(x, "model"), pch = 19, add.signif = FALSE, ...)
{
## Check that 'x' is of class "confint"
if (!inherits(x, "confint"))
stop("x should be of class \"confint\"")
## Invert the order
ci <- x[nrow(x):1, ]
## Define the color. If a "groups" attribute exists...
if (!is.null(attr(ci, "groups"))) {
## ... and the color is not set, use a hcl palette
if (is.null(col))
col <- hcl(h = seq(50, 360, length.out = length(unique(attr(ci,
"groups")))), c = 100, l = 50)
## In both cases, repeat the colors according to the groups
col <- rev(rep(col, table(attr(ci, "groups"))))
}
## But if there is no "groups" attribute, and no color set, use
## 'gray(O.3)'
else if (is.null(col))
col <- gray(0.3)
## Set 'par("mar")'
par(mar = mar)
## Call a plot.default on 'ci' (adjust the y-limits to make it look
## more regular)
plot.default(ci[, 3], 1:nrow(ci), xlim = range(ci), axes = FALSE,
xlab = paste("Confidence intervals, level = ", attr(ci,
"level") * 100, "%", sep = ""), ylab = NA, main = main,
type = "n", ylim = c(.7, nrow(ci) + .3))
## Add a grid
grid()
## Add a vertical line at x = 0
abline(v = 0)
## Add the confidence intervals (colored if does not overlap with zero)
for (i in 1:nrow(ci)) lines(x = c(ci[i, 1], ci[i, 2]), y = c(i,
i), col = ifelse(ci[i, 1] < 0 & ci[i, 2] > 0, "black",
col[i]))
## Add the points
points(ci[, 3], 1:nrow(ci), col = col, pch = pch, ...)
## Add significance
if (add.signif)
points(ci[!(ci[, 1] < 0 & ci[, 2] > 0), 3] + (max(ci) -
min(ci))/30, (1:nrow(ci))[!(ci[, 1] < 0 & ci[, 2] >
0)] + .2, pch = 8)
## Add the axes
axis(2, at = 1:nrow(ci), labels = FALSE)
## Add the row.names of the 'confint' call
mtext(side = 2, text = row.names(ci), at = 1:nrow(ci), col = col,
line = 1, las = 1)
## Add horizontal lines to visually differentiate the groups.
if (length(unique(attr(ci, "groups"))) > 1)
abline(h = (which(!duplicated(rev(sort(attr(ci, "groups"))))) -
0.5)[-1], col = gray(0.3), lty = 2)
## Add the x axis
axis(1)
## Add a bounding box
box()
}
## See: http://christophergandrud.blogspot.ca/2012/09/graphically-comparing-confidence.html
## and: https://github.com/christophergandrud/GreenBook/blob/master/Analysis/BasicAnalysisCoefPlots.R
### Example with ggplot2:
### https://diffuseprior.wordpress.com/2012/04/23/probitlogit-marginal-effects-in-r-2/
### http://is-r.tumblr.com/post/38055963968/possibly-slightly-better-text-analysis-with-lme4
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.