[R] Visualization of coefficients

David Atkins datkins at u.washington.edu
Tue Jul 6 23:02:54 CEST 2010



FYI, there is already a function coefplot in the arm package; for 
example, compare:

 > library(arm)
Loading required package: MASS
Loading required package: Matrix
[snip]
Attaching package: 'arm'

The following object(s) are masked from 'package:coda':

     traceplot

 > data("Mroz", package = "car")
 > fm <- glm(lfp ~ ., data = Mroz, family = binomial)
 > coefplot(fm)

with version below.

cheeres, Dave

 >
 > detach("package:arm")
 > coefplot <- function(object, df = NULL, level = 0.95, parm = NULL,
+    labels = TRUE, xlab = "Coefficient confidence intervals", ylab = "",
+    xlim = NULL, ylim = NULL,
+    las = 1, lwd = 1, lty = c(1, 2), pch = 19, col = 1,
+    length = 0, angle = 30, code = 3, ...)
+ {
+    cf <- coef(object)
+    se <- sqrt(diag(vcov(object)))
+    if(is.null(parm)) parm <- seq_along(cf)
+    if(is.numeric(parm) | is.logical(parm)) parm <- names(cf)[parm]
+    if(is.character(parm)) parm <- which(names(cf) %in% parm)
+    cf <- cf[parm]
+    se <- se[parm]
+    k <- length(cf)
+
+    if(is.null(df)) {
+      df <- if(identical(class(object), "lm")) df.residual(object) else 0
+    }
+
+    critval <- if(df > 0 & is.finite(df)) {
+      qt((1 - level)/2, df = df)
+    } else {
+      qnorm((1 - level)/2)
+    }
+    ci1 <- cf + critval * se
+    ci2 <- cf - critval * se
+
+    lwd <- rep(lwd, length.out = 2)
+    lty <- rep(lty, length.out = 2)
+    pch <- rep(pch, length.out = k)
+    col <- rep(col, length.out = k)
+
+    if(is.null(xlim)) xlim <- range(c(0, min(ci1), max(ci2)))
+    if(is.null(ylim)) ylim <- c(1 - 0.05 * k, 1.05 * k)
+
+    if(isTRUE(labels)) labels <- names(cf)
+    if(identical(labels, FALSE)) labels <- ""
+    labels <- rep(labels, length.out = k)
+
+    plot(0, 0, xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab,
+      axes = FALSE, type = "n", las = las, ...)
+    arrows(ci1, 1:k, ci2, 1:k, lty = lty[1], lwd = lwd[1], col = col,
+      length = length, angle = angle, code = code)
+    points(cf, 1:k, pch = pch, col = col)
+    abline(v = 0, lty = lty[2], lwd = lwd[2])
+    axis(1)
+    axis(2, at = 1:k, labels = labels, las = las)
+    box()
+ }
 >
 >
 > coefplot(fm, parm = -1)



Achim Zeileis wrote:

I've thought about adding a plot() method for the coeftest() function in
the "lmtest" package. Essentially, it relies on a coef() and a vcov()
method being available - and that a central limit theorem holds. For
releasing it as a general function in the package the code is still too
raw, but maybe it's useful for someone on the list. Hence, I've included
it below.

An example would be to visualize all coefficients except the intercept for
the Mroz data:

data("Mroz", package = "car")
fm <- glm(lfp ~ ., data = Mroz, family = binomial)
coefplot(fm, parm = -1)

hth,
Z

coefplot <- function(object, df = NULL, level = 0.95, parm = NULL,
    labels = TRUE, xlab = "Coefficient confidence intervals", ylab = "",
    xlim = NULL, ylim = NULL,
    las = 1, lwd = 1, lty = c(1, 2), pch = 19, col = 1,
    length = 0, angle = 30, code = 3, ...)
{
    cf <- coef(object)
    se <- sqrt(diag(vcov(object)))
    if(is.null(parm)) parm <- seq_along(cf)
    if(is.numeric(parm) | is.logical(parm)) parm <- names(cf)[parm]
    if(is.character(parm)) parm <- which(names(cf) %in% parm)
    cf <- cf[parm]
    se <- se[parm]
    k <- length(cf)

    if(is.null(df)) {
      df <- if(identical(class(object), "lm")) df.residual(object) else 0
    }

    critval <- if(df > 0 & is.finite(df)) {
      qt((1 - level)/2, df = df)
    } else {
      qnorm((1 - level)/2)
    }
    ci1 <- cf + critval * se
    ci2 <- cf - critval * se

    lwd <- rep(lwd, length.out = 2)
    lty <- rep(lty, length.out = 2)
    pch <- rep(pch, length.out = k)
    col <- rep(col, length.out = k)

    if(is.null(xlim)) xlim <- range(c(0, min(ci1), max(ci2)))
    if(is.null(ylim)) ylim <- c(1 - 0.05 * k, 1.05 * k)

    if(isTRUE(labels)) labels <- names(cf)
    if(identical(labels, FALSE)) labels <- ""
    labels <- rep(labels, length.out = k)

    plot(0, 0, xlim = xlim, ylim = ylim, xlab = xlab, ylab = ylab,
      axes = FALSE, type = "n", las = las, ...)
    arrows(ci1, 1:k, ci2, 1:k, lty = lty[1], lwd = lwd[1], col = col,
      length = length, angle = angle, code = code)
    points(cf, 1:k, pch = pch, col = col)
    abline(v = 0, lty = lty[2], lwd = lwd[2])
    axis(1)
    axis(2, at = 1:k, labels = labels, las = las)
    box()
}


On Fri, 2 Jul 2010, Tal Galili wrote:

 > Specifically this link:
 > http://tables2graphs.com/doku.php?id=04_regression_coefficients
 >
 > Great reference Bernd, thank you.
 >
 > Tal
 >
 >
 > ----------------Contact
 > Details:-------------------------------------------------------
 > Contact me: Tal.Galili at gmail.com |  972-52-7275845
 > Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il (Hebrew) |
 > www.r-statistics.com (English)
 > 
----------------------------------------------------------------------------------------------
 >
 >
 >
 >
 > On Fri, Jul 2, 2010 at 10:31 AM, Bernd Weiss <bernd.weiss at 
uni-koeln.de>wrote:
 >
 >> Am 02.07.2010 08:10, schrieb Wincent:
 >>> Dear all,
 >>>
 >>> I try to show a subset of coefficients in my presentation. It seems
 >>> that a "standard" table is not a good way to go. I found figure 9
 >>> (page 9) in this file (
 >>>
 >> 
http://www.destatis.de/jetspeed/portal/cms/Sites/destatis/Internet/DE/Content/Wissenschaftsforum/Kolloquien/VisualisierungModellierung__Beitrag,property=file.pdf
 >>>
 >>>
 >> ) looks pretty good. I wonder if there is any function for such plot?
 >>> Or any suggestion on how to present statistical models in a
 >>> presentation?
 >>
 >> Hi Wincent,
 >>
 >> I guess you are looking for "Using Graphs Instead of Tables in Political
 >> Science" by Kastellec/Leoni <http://tables2graphs.com/doku.php>.
 >>
 >> HTH,
 >>
 >> Bernd
 >>
 >> ______________________________________________
 >> R-help at r-project.org mailing list
 >> https://stat.ethz.ch/mailman/listinfo/r-help
 >> PLEASE do read the posting guide
 >> http://www.R-project.org/posting-guide.html
 >> and provide commented, minimal, self-contained, reproducible code.
 >>
 >
 > 	[[alternative HTML version deleted]]
 >
 > ______________________________________________
 > R-help at r-project.org mailing list
 > https://stat.ethz.ch/mailman/listinfo/r-help
 > PLEASE do read the posting guide 
http://www.R-project.org/posting-guide.html
 > and provide commented, minimal, self-contained, reproducible code.
 >

-- 
Dave Atkins, PhD
Research Associate Professor
Department of Psychiatry and Behavioral Science
University of Washington
datkins at u.washington.edu

Center for the Study of Health and Risk Behaviors (CSHRB)		
1100 NE 45th Street, Suite 300 	
Seattle, WA  98105 	
206-616-3879 	
http://depts.washington.edu/cshrb/
(Mon-Wed)	

Center for Healthcare Improvement, for Addictions, Mental Illness,
   Medically Vulnerable Populations (CHAMMP)
325 9th Avenue, 2HH-15
Box 359911
Seattle, WA 98104
http://www.chammp.org
(Thurs)



More information about the R-help mailing list