plotting RR, 95% CI as table and figure in same plot

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

plotting RR, 95% CI as table and figure in same plot

David Atkins

Hi all--

I am in the process of helping colleagues write up a ms in which we fit
zero-inflated Poisson models.  I would prefer plotting the rate ratios
and 95% CI (as I've found Gelman and others convincing about plotting
tables...), but our journals usually like the numbers themselves.

Thus, I'm looking at a recent JAMA article in which both numbers and
dotplot of RR and 95% CI are presented and wondering about best way to
do this in R.

Essentially, the plot has 3 columns: variable names, RR and 95% CI, and
dotplot of the same.

Using the bioChemists data in the pscl package and errbar function in
Hmisc package, the code below is in the right direction... but still
pretty ugly.

Wondering if folks would have alternative suggestions about how to go
about this, or pointers on cleaning up the code below (eg, I know there
are many functions for plotting errbars/CI).

[And, obviously, there are somethings that would be straightforward to
clean-up such as supplying better variable names, etc., just wanted to
see if there were better overall suggestions before getting too far on
this route.]

Thanks in advance.

cheers, Dave

library(Hmisc)
library(pscl)

## data
data("bioChemists", package = "pscl")
fm_pois <- glm(art ~ ., data = bioChemists, family = poisson)
summary(fm_pois)

### pull out rate-ratios and 95% CI
rr <- exp(cbind(coef(fm_pois), confint(fm_pois)))
rr
### round to 2 decimal places
rr <- round(rr, 2)

### plot
par(mfrow=c(1,3))
plot(0, type = "n", xlim=c(0,2), ylim=c(1,6),
        axes = FALSE, ylab=NULL, xlab=NULL)
text(row.names(rr), x = 1, y = 1:6)

plot(0, type = "n", xlim=c(0,2), ylim=c(1,6),
        axes = FALSE, ylab=NULL, xlab=NULL)
text(paste(rr[,1], " [", rr[,2], ", ", rr[,3], "]", sep = ""), x = 1, y
= 1:6)

errbar(x = factor(row.names(rr)),
                y = rr[,1], yplus = rr[,3],
                yminus = rr[,2])
abline(v = 1, lty =2)

--
Dave Atkins, PhD
Research Associate Professor
Department of Psychiatry and Behavioral Science
University of Washington
[hidden email]

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?
206-897-4210
http://www.chammp.org
(Thurs)

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: plotting RR, 95% CI as table and figure in same plot

Thomas Lumley


You could try the forestplot() function in rmeta, or the original grid code on which it is based,
   http://www.stat.auckland.ac.nz/~paul/RGraphics/chapter1.html

     -thomas

On Mon, 19 Apr 2010, David Atkins wrote:

>
> Hi all--
>
> I am in the process of helping colleagues write up a ms in which we fit
> zero-inflated Poisson models.  I would prefer plotting the rate ratios and
> 95% CI (as I've found Gelman and others convincing about plotting tables...),
> but our journals usually like the numbers themselves.
>
> Thus, I'm looking at a recent JAMA article in which both numbers and dotplot
> of RR and 95% CI are presented and wondering about best way to do this in R.
>
> Essentially, the plot has 3 columns: variable names, RR and 95% CI, and
> dotplot of the same.
>
> Using the bioChemists data in the pscl package and errbar function in Hmisc
> package, the code below is in the right direction... but still pretty ugly.
>
> Wondering if folks would have alternative suggestions about how to go about
> this, or pointers on cleaning up the code below (eg, I know there are many
> functions for plotting errbars/CI).
>
> [And, obviously, there are somethings that would be straightforward to
> clean-up such as supplying better variable names, etc., just wanted to see if
> there were better overall suggestions before getting too far on this route.]
>
> Thanks in advance.
>
> cheers, Dave
>
> library(Hmisc)
> library(pscl)
>
> ## data
> data("bioChemists", package = "pscl")
> fm_pois <- glm(art ~ ., data = bioChemists, family = poisson)
> summary(fm_pois)
>
> ### pull out rate-ratios and 95% CI
> rr <- exp(cbind(coef(fm_pois), confint(fm_pois)))
> rr
> ### round to 2 decimal places
> rr <- round(rr, 2)
>
> ### plot
> par(mfrow=c(1,3))
> plot(0, type = "n", xlim=c(0,2), ylim=c(1,6),
> axes = FALSE, ylab=NULL, xlab=NULL)
> text(row.names(rr), x = 1, y = 1:6)
>
> plot(0, type = "n", xlim=c(0,2), ylim=c(1,6),
> axes = FALSE, ylab=NULL, xlab=NULL)
> text(paste(rr[,1], " [", rr[,2], ", ", rr[,3], "]", sep = ""), x = 1, y =
> 1:6)
>
> errbar(x = factor(row.names(rr)),
> y = rr[,1], yplus = rr[,3],
> yminus = rr[,2])
> abline(v = 1, lty =2)
>
> --
> Dave Atkins, PhD
> Research Associate Professor
> Department of Psychiatry and Behavioral Science
> University of Washington
> [hidden email]
>
> 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?
> 206-897-4210
> http://www.chammp.org
> (Thurs)
>
> ______________________________________________
> [hidden email] 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.
>

Thomas Lumley Assoc. Professor, Biostatistics
[hidden email] University of Washington, Seattle

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: plotting RR, 95% CI as table and figure in same plot

Christos Argyropoulos
In reply to this post by David Atkins

ggplot2 should work (resize to get the plot to the dimensions you need for the paper)

 

library(Hmisc)
library(pscl)
library(ggplot2)

## data
data("bioChemists", package = "pscl")
fm_pois <- glm(art ~ ., data = bioChemists, family = poisson)
summary(fm_pois)

### pull out rate-ratios and 95% CI
rr <- exp(cbind(coef(fm_pois), confint(fm_pois)))
rr
### round to 2 decimal places
rr <- as.data.frame(round(rr, 2))

colnames(rr)<-c("y","ymin","ymax")
rr$labl<-rownames(rr) ## Change this to meaningful labels
rr$x<-1:length(rownames(rr))
gpl<-ggplot(rr,aes(x,y,ymin=ymin,ymax=ymax))
gpl+geom_point()+geom_linerange()+
 geom_hline(aes(yintercept=1),
 linetype="dashed",size=0.5)+
 geom_text(aes(x,y=0.3,label=y,hjust=0),size=3)+
 geom_text(aes(x,y=0.0,label=labl),size=3)+
 geom_text(aes(x,y=0.5,label=paste("[",ymin,",",ymax,"]",sep="")
 ,hjust=0.0),size=3)+
 ylab("Relative Risk")+xlab("")+
 coord_cartesian(ylim=c(-1,1.7))+
 coord_cartesian(xlim=c(0.85,6.15))+
 scale_x_continuous(breaks=NA)+
 scale_y_continuous(breaks=seq(0.8,1.6,.1))+
 opts(
  panel.grid.major = theme_blank(),
  panel.grid.minor=theme_blank(),
  title="",
  panel.background = theme_rect(fill=NA,colour=NA)
 )+
 coord_flip()+
geom_text(aes(x=6.3,y=0.35,label="RR"),size=4)+
geom_text(aes(x=6.3,y=0.60,label="95% CI"),size=4)

 

 

Christos


 

> Date: Mon, 19 Apr 2010 09:29:48 -0700
> From: [hidden email]
> To: [hidden email]
> Subject: [R] plotting RR, 95% CI as table and figure in same plot
>
>
> Hi all--
>
> I am in the process of helping colleagues write up a ms in which we fit
> zero-inflated Poisson models. I would prefer plotting the rate ratios
> and 95% CI (as I've found Gelman and others convincing about plotting
> tables...), but our journals usually like the numbers themselves.
>
> Thus, I'm looking at a recent JAMA article in which both numbers and
> dotplot of RR and 95% CI are presented and wondering about best way to
> do this in R.
>
> Essentially, the plot has 3 columns: variable names, RR and 95% CI, and
> dotplot of the same.
>
> Using the bioChemists data in the pscl package and errbar function in
> Hmisc package, the code below is in the right direction... but still
> pretty ugly.
>
> Wondering if folks would have alternative suggestions about how to go
> about this, or pointers on cleaning up the code below (eg, I know there
> are many functions for plotting errbars/CI).
>
> [And, obviously, there are somethings that would be straightforward to
> clean-up such as supplying better variable names, etc., just wanted to
> see if there were better overall suggestions before getting too far on
> this route.]
>
> Thanks in advance.
>
> cheers, Dave
>
> library(Hmisc)
> library(pscl)
>
> ## data
> data("bioChemists", package = "pscl")
> fm_pois <- glm(art ~ ., data = bioChemists, family = poisson)
> summary(fm_pois)
>
> ### pull out rate-ratios and 95% CI
> rr <- exp(cbind(coef(fm_pois), confint(fm_pois)))
> rr
> ### round to 2 decimal places
> rr <- round(rr, 2)
>
> ### plot
> par(mfrow=c(1,3))
> plot(0, type = "n", xlim=c(0,2), ylim=c(1,6),
> axes = FALSE, ylab=NULL, xlab=NULL)
> text(row.names(rr), x = 1, y = 1:6)
>
> plot(0, type = "n", xlim=c(0,2), ylim=c(1,6),
> axes = FALSE, ylab=NULL, xlab=NULL)
> text(paste(rr[,1], " [", rr[,2], ", ", rr[,3], "]", sep = ""), x = 1, y
> = 1:6)
>
> errbar(x = factor(row.names(rr)),
> y = rr[,1], yplus = rr[,3],
> yminus = rr[,2])
> abline(v = 1, lty =2)
>
> --
> Dave Atkins, PhD
> Research Associate Professor
> Department of Psychiatry and Behavioral Science
> University of Washington
> [hidden email]
>
> 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?
> 206-897-4210
> http://www.chammp.org
> (Thurs)
>
> ______________________________________________
> [hidden email] 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.
     
_________________________________________________________________
Hotmail: Trusted email with powerful SPAM protection.

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: plotting RR, 95% CI as table and figure in same plot

David Atkins
In reply to this post by Thomas Lumley

Thanks to Thomas and Christos for helpful suggestions.

The forestplot (in package rmeta) suggestion seems to work fairly well
for me, though does require a bit of fiddling (no complaints, obviously
using it for a different purpose than it was written).

Below is an example using a slightly "hacked" version of forestplot (and
also using a ZIP model).

[BTW, my "hacks" were to adjust the code so I could set the line weights
and to use circles as opposed to boxes and set the radii.]

cheers, Dave


## data
data("bioChemists", package = "pscl")
fm_zip <- zeroinfl(art ~ ., data = bioChemists)
summary(fm_zip)

### pull out rate-ratios and 95% CI
rr <- exp(cbind(coef(fm_zip), confint(fm_zip)))
rr

### round to 2 decimal places
rr <- format(rr, digits=2)


####################### Alternative: forestplot() from rmeta package
########################
#
library(rmeta)

preds <- c("Intercept","Women","Married","Kids","PhD","Mentor")
tab.txt <- rbind(c("Predictor","RR [95% CI]"), c("", ""),
                                        cbind(preds, paste(rr[1:6,1], " [", rr[1:6,2], ", ", rr[1:6,3], "]",
                                                                                                  sep = "")),
                                        c("", ""),
                                        c("Predictor","OR [95% CI]"), c("", ""),
                                        cbind(preds, paste(rr[7:12,1], " [", rr[7:12,2], ", ", rr[7:12,3],
"]",
                                                                                                  sep = "")))
tab.txt

dat.txt <- rbind(c(NA,NA,NA), c(NA,NA,NA), rr[1:6,],
                                   c(NA,NA,NA), c(NA,NA,NA), c(NA,NA,NA), rr[7:12,])
dat.txt

### NOTE: slightly hacked version of forestplot from rmeta
forestplot2(labeltext = tab.txt,
                  mean = dat.txt[,1], lower = dat.txt[,2], upper = dat.txt[,3],
                  zero=1,
                  is.summary=c(TRUE,rep(FALSE,8),TRUE,rep(FALSE,8)),
                  xlog=FALSE,
                  graphwidth = unit(3, "inches"), lwd= 3, rad = 0.3)



########################### Functions
########################################

forestplot2 <-
function (labeltext, mean, lower, upper, align = NULL, is.summary = FALSE,
     clip = c(-Inf, Inf), xlab = "", zero = 0, graphwidth = unit(2,
         "inches"), col = meta.colors(), xlog = FALSE, xticks = NULL,
     boxsize = NULL, lwd = 1, rad = 0.1, ...)
{
     require("grid") || stop("`grid' package not found")
     require("rmeta") || stop("`rmeta' package not found")
     drawNormalCI <- function(LL, OR, UL, size) {
         size = 0.75 * size
         clipupper <- convertX(unit(UL, "native"), "npc", valueOnly =
TRUE) >
             1
         cliplower <- convertX(unit(LL, "native"), "npc", valueOnly =
TRUE) <
             0
         box <- convertX(unit(OR, "native"), "npc", valueOnly = TRUE)
         clipbox <- box < 0 || box > 1
         if (clipupper || cliplower) {
             ends <- "both"
             lims <- unit(c(0, 1), c("npc", "npc"))
             if (!clipupper) {
                 ends <- "first"
                 lims <- unit(c(0, UL), c("npc", "native"))
             }
             if (!cliplower) {
                 ends <- "last"
                 lims <- unit(c(LL, 1), c("native", "npc"))
             }
             grid.lines(x = lims, y = 0.5, arrow = arrow(ends = ends,
                 length = unit(0.05, "inches")), gp = gpar(col = col$lines))
             if (!clipbox)
                 grid.rect(x = unit(OR, "native"), width = unit(size,
                   "snpc"), height = unit(size, "snpc"), gp = gpar(fill
= col$box,
                   col = col$box))
         }
         else {
             grid.lines(x = unit(c(LL, UL), "native"), y = 0.5,
                 gp = gpar(col = 1, lwd = lwd))
             grid.circle(x = unit(OR, "native"),
              #width = unit(size, "snpc"),
              #height = unit(size, "snpc"),
              r = rad,
              gp = gpar(fill = col$box,
                  col = col$box))
             if ((convertX(unit(OR, "native") + unit(0.5 * size,
                 "lines"), "native", valueOnly = TRUE) > UL) &&
                 (convertX(unit(OR, "native") - unit(0.5 * size,
                   "lines"), "native", valueOnly = TRUE) < LL))
                 grid.lines(x = unit(c(LL, UL), "native"), y = 0.5,
                   gp = gpar(col = col$lines))
         }
     }
     drawSummaryCI <- function(LL, OR, UL, size) {
         grid.polygon(x = unit(c(LL, OR, UL, OR), "native"), y = unit(0.5 +
             c(0, 0.5 * size, 0, -0.5 * size), "npc"), gp = gpar(fill =
col$summary,
             col = col$summary))
     }
     plot.new()
     widthcolumn <- !apply(is.na(labeltext), 1, any)
     nc <- NCOL(labeltext)
     labels <- vector("list", nc)
     if (is.null(align))
         align <- c("l", rep("r", nc - 1))
     else align <- rep(align, length = nc)
     nr <- NROW(labeltext)
     is.summary <- rep(is.summary, length = nr)
     for (j in 1:nc) {
         labels[[j]] <- vector("list", nr)
         for (i in 1:nr) {
             if (is.na(labeltext[i, j]))
                 next
             x <- switch(align[j], l = 0, r = 1, c = 0.5)
             just <- switch(align[j], l = "left", r = "right",
                 c = "center")
             labels[[j]][[i]] <- textGrob(labeltext[i, j], x = x,
                 just = just, gp = gpar(fontface = if (is.summary[i])
                   "bold"
                 else "plain", col = rep(col$text, length = nr)[i]))
         }
     }
     colgap <- unit(8, "mm")
     colwidths <- unit.c(max(unit(rep(1, sum(widthcolumn)), "grobwidth",
         labels[[1]][widthcolumn])), colgap)
     if (nc > 1) {
         for (i in 2:nc) colwidths <- unit.c(colwidths, max(unit(rep(1,
             sum(widthcolumn)), "grobwidth", labels[[i]][widthcolumn])),
             colgap)
     }
     colwidths <- unit.c(colwidths, graphwidth)
     pushViewport(viewport(layout = grid.layout(nr + 1, nc * 2 +
         1, widths = colwidths, heights = unit(c(rep(1, nr), 0.5),
         "lines"))))
     cwidth <- (upper - lower)
     xrange <- c(max(min(lower, na.rm = TRUE), clip[1]), min(max(upper,
         na.rm = TRUE), clip[2]))
     info <- 1/cwidth
     info <- info/max(info[!is.summary], na.rm = TRUE)
     info[is.summary] <- 1
     if (!is.null(boxsize))
         info <- rep(boxsize, length = length(info))
     for (j in 1:nc) {
         for (i in 1:nr) {
             if (!is.null(labels[[j]][[i]])) {
                 pushViewport(viewport(layout.pos.row = i,
layout.pos.col = 2 *
                   j - 1))
                 grid.draw(labels[[j]][[i]])
                 popViewport()
             }
         }
     }
     pushViewport(viewport(layout.pos.col = 2 * nc + 1, xscale = xrange))
     grid.lines(x = unit(zero, "native"), y = 0:1, gp = gpar(col = 1))
     if (xlog) {
         if (is.null(xticks)) {
             ticks <- pretty(exp(xrange))
             ticks <- ticks[ticks > 0]
         }
         else {
             ticks <- xticks
         }
         if (length(ticks)) {
             if (min(lower, na.rm = TRUE) < clip[1])
                 ticks <- c(exp(clip[1]), ticks)
             if (max(upper, na.rm = TRUE) > clip[2])
                 ticks <- c(ticks, exp(clip[2]))
             xax <- xaxisGrob(gp = gpar(cex = 0.6, col = col$axes),
                 at = log(ticks), name = "xax")
             xax1 <- editGrob(xax, gPath("labels"), label = format(ticks,
                 digits = 2))
             grid.draw(xax1)
         }
     }
     else {
         if (is.null(xticks)) {
             grid.xaxis(gp = gpar(cex = 0.6, col = col$axes))
         }
         else if (length(xticks)) {
             grid.xaxis(at = xticks, gp = gpar(cex = 0.6, col = col$axes))
         }
     }
     grid.text(xlab, y = unit(-2, "lines"), gp = gpar(col = col$axes))
     popViewport()
     for (i in 1:nr) {
         if (is.na(mean[i]))
             next
         pushViewport(viewport(layout.pos.row = i, layout.pos.col = 2 *
             nc + 1, xscale = xrange))
         if (is.summary[i])
             drawSummaryCI(lower[i], mean[i], upper[i], info[i])
         else drawNormalCI(lower[i], mean[i], upper[i], info[i])
         popViewport()
     }
     popViewport()
}

Dave Atkins, PhD
Research Associate Professor
Department of Psychiatry and Behavioral Science
University of Washington
[hidden email]

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?
206-897-4210
http://www.chammp.org
(Thurs)


Thomas Lumley wrote:

>
>
> You could try the forestplot() function in rmeta, or the original grid
> code on which it is based,
>   http://www.stat.auckland.ac.nz/~paul/RGraphics/chapter1.html
>
>     -thomas
>
> On Mon, 19 Apr 2010, David Atkins wrote:
>
>>
>> Hi all--
>>
>> I am in the process of helping colleagues write up a ms in which we
>> fit zero-inflated Poisson models.  I would prefer plotting the rate
>> ratios and 95% CI (as I've found Gelman and others convincing about
>> plotting tables...), but our journals usually like the numbers
>> themselves.
>>
>> Thus, I'm looking at a recent JAMA article in which both numbers and
>> dotplot of RR and 95% CI are presented and wondering about best way to
>> do this in R.
>>
>> Essentially, the plot has 3 columns: variable names, RR and 95% CI,
>> and dotplot of the same.
>>
>> Using the bioChemists data in the pscl package and errbar function in
>> Hmisc package, the code below is in the right direction... but still
>> pretty ugly.
>>
>> Wondering if folks would have alternative suggestions about how to go
>> about this, or pointers on cleaning up the code below (eg, I know
>> there are many functions for plotting errbars/CI).
>>
>> [And, obviously, there are somethings that would be straightforward to
>> clean-up such as supplying better variable names, etc., just wanted to
>> see if there were better overall suggestions before getting too far on
>> this route.]
>>
>> Thanks in advance.
>>
>> cheers, Dave
>>
>> library(Hmisc)
>> library(pscl)  
>>
>> ## data
>> data("bioChemists", package = "pscl")
>> fm_pois <- glm(art ~ ., data = bioChemists, family = poisson)
>> summary(fm_pois)
>>
>> ### pull out rate-ratios and 95% CI
>> rr <- exp(cbind(coef(fm_pois), confint(fm_pois)))
>> rr
>> ### round to 2 decimal places
>> rr <- round(rr, 2)
>>
>> ### plot
>> par(mfrow=c(1,3))
>> plot(0, type = "n", xlim=c(0,2), ylim=c(1,6),
>>     axes = FALSE, ylab=NULL, xlab=NULL)
>> text(row.names(rr), x = 1, y = 1:6)
>>
>> plot(0, type = "n", xlim=c(0,2), ylim=c(1,6),
>>     axes = FALSE, ylab=NULL, xlab=NULL)
>> text(paste(rr[,1], " [", rr[,2], ", ", rr[,3], "]", sep = ""), x = 1,
>> y = 1:6)
>>
>> errbar(x = factor(row.names(rr)),
>>         y = rr[,1], yplus = rr[,3],
>>         yminus = rr[,2])
>> abline(v = 1, lty =2)  
>>
>> --
>> Dave Atkins, PhD
>> Research Associate Professor
>> Department of Psychiatry and Behavioral Science
>> University of Washington
>> [hidden email]
>>
>> 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?
>> 206-897-4210
>> http://www.chammp.org
>> (Thurs)
>>
>> ______________________________________________
>> [hidden email] 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.
>>
>
> Thomas Lumley            Assoc. Professor, Biostatistics
> [hidden email]    University of Washington, Seattle
>

______________________________________________
[hidden email] 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.