Dear Pr Harrel,
Thanks for your help at this occasion as well as for previous questions to the list.
I just looked at the example in your intro doc.
However I am dealing with proportions (ie % of patients responding to a given treatment).
In this case I am not sure I can use summarize and then the xYplot.
I am not aware of R graphing tools that can deal directly with proportions adding CI not to mention producing by factor/trellis plots.
This is why I why trying to do it "by hand" (using binconf) with ooplot,
without much success I am afraid.
Best regards,
JL Abitbol, MD
Jean-Louis,
Here is an example.
# Plot proportions and their Wilson confidence limits
set.seed(3)
d <- expand.grid(continent=c('USA','Europe'), year=1999:2001,
reps=1:100)
# Generate binary events from a population probability of 0.2
# of the event, same for all years and continents
d$y <- ifelse(runif(6*100) <= .2, 1, 0)
s <- with(d,
summarize(y, llist(continent,year),
function(y) {
n <- sum(!is.na(y))
s <- sum(y, na.rm=T)
binconf(s, n)
}, type='matrix')
)Dotplot(year ~ Cbind(y) | continent, data=s, ylab='Year',
xlab='Probability')I did have to temporarily override a function in Hmisc to fix a problem. This will be corrected in an upcoming release of Hmisc:
mApply <- function(X, INDEX, FUN=NULL, ..., simplify=TRUE) {
## Matrix tapply
## X: matrix with n rows; INDEX: vector or list of vectors of length n
## FUN: function to operate on submatrices of x by INDEX
## ...: arguments to FUN; simplify: see sapply
## Modification of code by Tony Plate <[EMAIL PROTECTED]> 10Oct02
## If FUN returns more than one number, mApply returns a matrix with
## rows corresponding to unique values of INDEX
nr <- nrow(X)
if(!length(nr)) { ## X not a matrix
r <- tapply(X, INDEX, FUN, ..., simplify=simplify)
if(is.matrix(r)) r <- drop(t(r)) else
if(simplify && is.list(r))
r <- drop(matrix(unlist(r), nrow=length(r),
dimnames=list(names(r),names(r[[1]])), byrow=TRUE))
} else {
idx.list <- tapply(1:nr, INDEX, c)
r <- sapply(idx.list, function(idx,x,fun,...) fun(x[idx,,drop=FALSE],...),
x=X, fun=FUN, ..., simplify=simplify)
if(simplify) r <- drop(t(r))
}
dn <- dimnames(r)
if(length(dn) && !length(dn[[length(dn)]])) {
fx <- FUN(X,...)
dnl <- if(length(names(fx))) names(fx) else dimnames(fx)[[2]]
dn[[length(dn)]] <- dnl
dimnames(r) <- dn
}
if(simplify && is.list(r) && is.array(r)) {ll <- sapply(r, length) maxl <- max(ll) empty <- (1:length(ll))[ll==0] for(i in empty) r[[i]] <- rep(NA, maxl) ## unlist not keep place for NULL entries for nonexistent categories first.not.empty <- ((1:length(ll))[ll > 0])[1] nam <- names(r[[first.not.empty]]) dr <- dim(r)
r <- aperm(array(unlist(r), dim=c(maxl,dr),
dimnames=c(list(nam),dimnames(r))),
c(1+seq(length(dr)), 1))
}
r
}Frank
On Sun, 21 Nov 2004 07:48:58 -0500, "Frank E Harrell Jr" <[EMAIL PROTECTED]> said:
Jean-Louis Abitbol wrote:
Dear All
I am trying to graph a proportion and CI95% by a factor with ooplot (any other better solution ?)
It works well until I try to add the confidence interval.
this is the error message and and a description of the data:
> dat1
PointEst
TT1 1 3.6
TT2 2 5.0
TT3 3 5.8
TT4 4 11.5
TT5 5 7.5
TT5 6 8.7
TT7 7 17.4
dat2
Lower TT1 1 1.0 TT2 2 2.2 TT3 3 2.7 TT4 4 6.7 TT5 5 3.9 TT5 6 4.6 TT7 7 11.5
dat3
Upper TT1 1 12.3 TT2 2 11.2 TT3 3 12.1 TT4 4 19.1 TT5 5 14.2 TT5 6 15.6 TT7 7 25.6
ooplot(dat1,type="barplot",col=rich.colors(7,"temperature"),names.arg=c("X","Y","Z","A","B","C","D"),plot.ci=T,
+ ci.l=dat2,ci.u=dat3, xlab="Treatment", ylab="Percent Normalized
Patients")
Error in ooplot.default(dat1, type = "barplot", col = rich.colors(7,
"temperature"), : 'height' and 'ci.u' must have the same dimensions.
I have tried various ways of supplying ci.l and ci.u (including a vector)
Thanks for the help that anyone can bring,
Regards, JL
One way is to look at the examples for Dotplot in the Hmisc package. Those examples display bootstrap percentile confidence intervals for a
mean.
-- Frank E Harrell Jr Professor and Chair School of Medicine Department of Biostatistics Vanderbilt University
______________________________________________ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
