[R] Visualization of coefficients

Achim Zeileis Achim.Zeileis at uibk.ac.at
Wed Jul 7 11:10:42 CEST 2010


On Wed, 7 Jul 2010, Tal Galili wrote:

> Hello David,
> Thanks to your posting I started looking at the function in the arm package.
>  It appears this function is quite mature, and offers (for example) the
> ability to easily overlap coefficients from several models.

Re: more mature. arm's coefplot() is more flexible in certain respects, 
mine is more convenient in others. The overlay functionality is something 
arm's coefplot() is better in and it also as some further options 
(vertical vs. horizontal etc.). My coefplot() has the advantage that it 
does not need any modification as long as coef() and vcov() methods are 
available. Furthermore, "level" can specify the significance level 
(instead of always using one and two standard errors, respectively).
But it shouldn't be too hard to create a superset of all options.

> I updated the post I published on the subject, so at the end of it I give an
> example of comparing the coef of several models:
> http://www.r-statistics.com/2010/07/visualization-of-regression-coefficient
> s-in-r/

As Allan pointed out in his reply, something fully reproducible would be 
nice. Also, if you keep the example with quasi-complete separation, it 
would be worth pointing this out. (Because the maximum likelihood 
estimator is Infinity in this case.)

Finally, the Poisson model in comparison with the binomial models does not 
make much sense, I guess.

Best,
Z

> Thanks again for the pointer.
> 
> Best,
> 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 Wed, Jul 7, 2010 at 12:02 AM, David Atkins <datkins at u.washington.edu>
> wrote:
> 
>
>       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/Conte
> nt/Wissenschaftsforum/Kolloquien/VisualisierungModellierung__Beitrag,proper
> ty=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)
> 
> 
> ______________________________________________
> 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.
> 
> 
> 
>


More information about the R-help mailing list