Here's an updated version of lineplot.CI that will succeed even for
cases where data are not present in all factor combinations. Also, this
version has the option x.cont to specify that the x axis represents a
continuous variable with proportional spacing. A new version of sciplot
with these changes will be posted soon.
## Examples:
source("lineplot.CI.R")
## Generate data
time=c(rep(c(21:30),3),rep(c(1:10),3))
y <- time+rnorm(60,0,1)
factors <- rep(c(1:2),each=30)
## Proportional spacing
lineplot.CI(resp=y, x.factor=time, group=factors, x.cont=TRUE)
## Factorial spacing
lineplot.CI(resp=y, x.factor=time, group=factors)
Manuel
On Fri, 2008-02-15 at 15:18 +0100, Dieter Vanderelst
wrote:> Hi List,
>
> I have a problem plotting data using the lineplot.CI command in the sciplot
package.
>
> I want to plot the data of 2 experimental cases using different lines
(traces). Time is on the X-axis. The tricky thing is that the data collection in
the second case started later than for the first case. This is to say: the first
n data points for the second case are missing.
>
> So far so good. However, when I plot the data using lineplot.CI, the
standard error bars are not aligned correctly with the markers.
>
> I know that this might be difficult to imagine. Here you can find an
example: http://i254.photobucket.com/albums/hh115/MarkerMe/example.png
>
> So, has anybody experienced this problem and solved it before? I think I
could try padding the data of the second case with zeros to eliminate the
missing data. But I hope there is a better solution.
>
> Regards,
> Dieter
> ----------------------------------------------------
> Dieter Vanderelst
> dieter dot vanderelst at emailengine dot org
>
> ______________________________________________
> 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.
--
http://mutualism.williams.edu
-------------- next part --------------
lineplot.CI <-
function(x.factor, response, group = NULL, type = "b", x.cont=FALSE,
legend = TRUE, trace.label = NULL, leg.lab = NULL,
fixed = FALSE, x.leg = NULL, y.leg = NULL, cex.leg = 1, ncol = 1,
pch = c(16,21,15,22,17,24, c(3:14)),
fun = function(x) mean(x, na.rm=TRUE),
ci.fun = function(x) c(fun(x)-se(x), fun(x)+se(x)),
err.width = if(length(levels(as.factor(x.factor)))>10) 0 else .1,
err.col = col, err.lty = 1, xlim = NULL, ylim = NULL,
cex = NULL, lwd = NULL, col = "black", cex.axis = 1,
xaxt = "s", data = NULL, subset = NULL, ...) {
# Set up environment
subset <- eval(substitute(subset), envir=data)
if(!is.null(data)) {
if(!is.null(subset)) data <- subset(data,subset)
x.factor <- eval(substitute(x.factor), envir=data)
response <- eval(substitute(response), envir=data)
group <- eval(substitute(group), envir=data)
}
subset = NULL
#////////////////////////////////////////////////////////////////////#
# Below modified from the interaction.plot function in stats package #
# Take the legend section out to allow greater flexibility #
#////////////////////////////////////////////////////////////////////#
int.plot <- function (x.factor = x.factor, group = group,
response = response, type =
c("l","p", "b"),
legend = legend,
trace.label = deparse(substitute(group)),
fixed = FALSE, xlab = deparse(substitute(x.factor)),
ylab = ylabel, lty = nc:1, pch = NA, xpd = NULL,
leg.bg = par("bg"), leg.bty = "n",
xtick = FALSE,
xlim=xlim, ylim=ylim, axes = TRUE, ...) {
ylabel <- paste(deparse(substitute(fun)), "of ",
deparse(substitute(response)))
type <- match.arg(type)
cells <- tapply(response, list(x.factor, group), fun)
nr <- nrow(cells)
nc <- ncol(cells)
xvals <- if(x.cont) as.numeric(levels(as.factor(x.factor))) else 1:nr
if (is.ordered(x.factor)) {
wn <- getOption("warn")
options(warn = -1)
xnm <- as.numeric(levels(x.factor))
options(warn = wn)
if (!any(is.na(xnm)))
xvals <- xnm
}
xlabs <- rownames(cells)
ylabs <- colnames(cells)
nch <- max(sapply(ylabs, nchar, type = "width"))
if (is.null(xlabs))
xlabs <- as.character(xvals)
if (is.null(ylabs))
ylabs <- as.character(1:nc)
if(is.null(xlim)) {
xlim <- range(xvals)
xleg <- xlim[2] + 0.05 * diff(xlim)
xlim <- xlim + c(-0.2/nr, if
(legend&is.null(x.leg)&is.null(y.leg))
0.2 + 0.02 * nch else 0.2/nr) * diff(xlim)}
else {
xlim
xleg <- xlim[2] - 0.25 * diff(xlim)}
matplot(xvals, cells, ..., type = type, xlim = xlim, ylim = ylim,
xlab = xlab, ylab = ylab, axes = axes, xaxt = "n", col =
col,
lty = lty, lwd=lwd, pch = pch)
if (axes && xaxt != "n") {
axisInt <- function(x, main, sub, lwd, bg, log, asp,
...) axis(1, x, ...)
mgp. <- par("mgp")
if (!xtick)
mgp.[2] <- 0
axisInt(1, at = xvals, labels = xlabs, tick = xtick,
mgp = mgp., xaxt = xaxt, ...)
}
ord <- sort.list(cells[nr, ], decreasing = TRUE)
if (legend) {
yrng <- diff(ylim)
yleg <- ylim[2] - 0.1 * yrng
if (!is.null(xpd) || {
xpd. <- par("xpd")
!is.na(xpd.) && !xpd. && (xpd <- TRUE)
}) {
op <- par(xpd = xpd)
on.exit(par(op))
}
text(xleg, ylim[2] - 0.05 * yrng, paste(" ", trace.label),
adj = 0)
if (!fixed) {
ylabs <- ylabs[ord]
lty <- lty[1 + (ord - 1)%%length(lty)]
col <- col[1 + (ord - 1)%%length(col)]
pch <- pch[ord]
}
}
invisible()
return.data<-if(legend)
list(pch=pch,ord=ord,xleg=xleg,yleg=yleg,ylabs=ylabs,lty=lty,
leg.bty=leg.bty,leg.bg=leg.bg,ord=ord,xvals=xvals,cells=cells)
else list(pch=pch,ord=ord,xvals=xvals)
return(return.data)
}
#////////////////////////////////////////////////////////////////////#
# End Section #
#////////////////////////////////////////////////////////////////////#
# Figure out if we're dealing with 1 or more total groups (i.e., we always
# have 1 "x.factor" and we may have 1 or more "group".
if(is.null(group)) groups = x.factor else {
# If more than 1 "y-factor", combine for plotting
if(length(group[[1]]) > 1)
group <- do.call("paste", c(group, sep = "."))
# "groups" defines the combination of "x.factor" and
"group"
groups <- list(x.factor, group)
}
# Calculate mean and SE's
mn.data <- tapply(response, groups, fun)
CI.data <- tapply(response, groups, ci.fun)
# Determine y-axis plot region
plot.limits = c(min(unlist(CI.data)), max(unlist(CI.data)))
# Draw lines
if(is.null(group)) {
nlevels.x <-
if(x.cont) as.numeric(levels(as.factor(x.factor)))
else 1:nrow(mn.data)
plot(nlevels.x, mn.data, xaxt="n", type=type, col=col, pch=NA,
cex=cex,
cex.axis=cex.axis,
xlim=if(is.null(xlim)) {
c(min(nlevels.x)-.2,max(nlevels.x)+.2)
}
else xlim,
ylim=if(is.null(ylim)) plot.limits else ylim, ...)
if(xaxt!="n") axis(1,label=names(mn.data),
at=nlevels.x, cex.axis=cex.axis)
}
else leg.vals <-
int.plot(x.factor, group, response, type = type,
xlim = xlim, ylim = if(is.null(ylim)) plot.limits else ylim,
cex.axis = cex.axis, trace.label = trace.label, pch = NA,
legend = legend, ...)
# Draw CI's
if(is.null(group)) {
nlevels.x <-
if(x.cont) as.numeric(levels(as.factor(x.factor)))
else 1:nrow(mn.data)
CI.seln <- !is.na(mn.data)
CI.plot <- matrix(unlist(CI.data[CI.seln]),
nrow=sum(CI.seln), byrow=TRUE)
arrows(nlevels.x[CI.seln],
CI.plot[,1],
nlevels.x[CI.seln],
CI.plot[,2],
angle = 90, col = err.col, length = err.width, code = 3,
lwd = lwd, lty=err.lty)
}
else {
nlevels.y <- ncol(mn.data)
for(i in 1:nlevels.y) {
CI.seln <- !is.na(mn.data)[,i]
CI.plot <- matrix(unlist(CI.data[CI.seln,i]),
nrow=sum(CI.seln), byrow=TRUE)
arrows(leg.vals$xvals[CI.seln],
CI.plot[,1],
leg.vals$xvals[CI.seln],
CI.plot[,2],
angle=90, length=err.width,
col=if(length(err.col)>1) err.col[i] else err.col,
lty=if(length(err.lty)>1) err.lty[i] else err.lty,
code=3,lwd=lwd)
}
}
# Draw points (Note: adding points at this point allows points to be in
# the foreground)
if(type %in% c("p", "b")) {
if(is.null(group)) {
nlevels.x <-
if(x.cont) as.numeric(levels(as.factor(x.factor)))
else 1:nrow(mn.data)
points(nlevels.x, mn.data,pch=pch[1],bg="white",cex=cex,col=col)
}
else {
nlevels.y<-dim(mn.data)[2]
for(i in 1:nlevels.y)
points(leg.vals$xvals, mn.data[,i],pch=pch[i],bg="white",
col=if(length(col)>1) col[i] else col,cex=cex)
}
}
#////////////////////////////////////////////////////////////////////#
# Below modified from the interaction.plot function in stats package #
#////////////////////////////////////////////////////////////////////#
if(legend & !is.null(group)) {
legend(x=if(is.null(x.leg)) leg.vals$xleg else x.leg,
y=if(is.null(y.leg)) leg.vals$yleg else y.leg,
legend = if(!is.null(leg.lab)) leg.lab else {
if(fixed) leg.vals$ylabs[leg.vals$ord] else leg.vals$ylabs
},
pch = if(type %in% c("p", "b") & !fixed)
pch[leg.vals$ord] else pch,
col = if(type %in% c("p", "b") &
length(col)>1 & !fixed)
col[leg.vals$ord] else col,
lty = if (type %in% c("l", "b") & fixed)
leg.vals$lty[leg.vals$ord]
else leg.vals$lty,
ncol=ncol, bty = leg.vals$leg.bty, bg = leg.vals$leg.bg,
cex=cex.leg)
}
#////////////////////////////////////////////////////////////////////#
# End Section #
#////////////////////////////////////////////////////////////////////#
invisible(list(vals=mn.data, CI=CI.data))
}
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 189 bytes
Desc: This is a digitally signed message part
Url :
https://stat.ethz.ch/pipermail/r-help/attachments/20080215/cf795c9c/attachment.bin