R/denstrip.R

Defines functions panel.sectioned.density panel.bpstrip panel.vwstrip.normal panel.vwstrip panel.cistrip panel.denstrip.legend panel.denstrip.normal panel.denstrip bpstrip vwstrip.normal vwstrip cistrip sectioned.density seqToIntervals densregion.normal densregion.survfit densregion.default densregion denstrip.legend denstrip.normal denstrip

Documented in bpstrip cistrip densregion densregion.default densregion.normal densregion.survfit denstrip denstrip.legend denstrip.normal panel.bpstrip panel.cistrip panel.denstrip panel.denstrip.legend panel.denstrip.normal panel.sectioned.density panel.vwstrip panel.vwstrip.normal sectioned.density seqToIntervals vwstrip vwstrip.normal

denstrip <- function(x, dens, at, width, horiz=TRUE, colmax, colmin="white", scale=1, gamma=1, 
                     ticks=NULL, tlen=1.5, twd, tcol, mticks=NULL, mlen=1.5, mwd, mcol, 
                     lattice=FALSE,
                     ...)
{
    if (!is.numeric(x)) stop("\'x\' must be numeric")
    if (missing(dens)) {
        ## x assumed to be a sample from a distribution, density estimated
        de <- density(x, ...)
        x <- de$x; dens <- de$y
    }
    else {
        if (!is.numeric(dens)) stop("\'dens\' must be numeric")
        if (length(dens) != length(x)) stop("Lengths of \'dens\' and \'x\' must be the same")
        dens <- dens[order(x)]
        x <- sort(x)
    }
    if (lattice) {
        rect.fn <- panel.rect; seg.fn <- panel.segments
        default.width <- diff(current.panel.limits()[[if(horiz) "ylim" else "xlim"]]) / 30
        default.colmax <- trellis.par.get("add.line")$col
        default.twd <- trellis.par.get("add.line")$lwd; default.mwd <- trellis.par.get("add.line")$lwd*2
    }
    else {
        rect.fn <- rect; seg.fn <- segments;
        default.width <- diff(par("usr")[if(horiz) 3:4 else 1:2]) / 30
        default.colmax <- par("fg")
        default.twd <- par("lwd"); default.mwd <- par("lwd")*2
    }
    if (missing(width)) width <- default.width
    if (missing(colmax)) colmax <- default.colmax
    if (missing(twd)) twd <- default.twd
    if (missing(mwd)) mwd <- default.mwd
    if (missing(tcol)) tcol <- colmax
    if (missing(mcol)) mcol <- colmax
    dens <- dens / max(dens) * scale
    n <- length(x)
    rgbmax <- col2rgb(colmax, alpha=TRUE)
    rgbmin <- if (colmin=="transparent") c(col2rgb(colmax, alpha=FALSE), 0) else col2rgb(colmin, alpha=TRUE)
    if (gamma <= 0) stop("gamma must be greater than 0")
    p <- dens[1:(n-1)] ^ gamma
    if (colmin=="transparent") 
        cols <- rgb(p*rgbmax[1] + (1 - p)*rgbmin[1],
                    p*rgbmax[2] + (1 - p)*rgbmin[2],
                    p*rgbmax[3] + (1 - p)*rgbmin[3],
                    alpha=p*rgbmax[4] + (1 - p)*rgbmin[4],
                    maxColorValue=255)
    else 
        cols <- rgb(p*rgbmax[1] + (1 - p)*rgbmin[1],
                    p*rgbmax[2] + (1 - p)*rgbmin[2],
                    p*rgbmax[3] + (1 - p)*rgbmin[3],
                    alpha=rgbmax[4],
                    maxColorValue=255)
    first.col <- c(TRUE, cols[2:(n-1)] != cols[1:(n-2)])
    next.col <- c(first.col,TRUE); next.col[1] <- FALSE    
    if (horiz) {
        xleft <- x[-n][first.col]; xright=x[next.col];
        ybottom <- at-width/2; ytop <- at+width/2
    }
    else {
        xleft <- at-width/2; xright <- at+width/2;
        ybottom <- x[-n][first.col]; ytop <- x[next.col]
    }
    rect.fn(xleft=xleft, ybottom=ybottom, xright=xright, ytop=ytop,
         border=NA, col = cols[first.col])
    if (!is.null(ticks)){
        if (horiz) { tx0 <- tx1 <- ticks; ty0 <- at-width*tlen/2; ty1 <- at+width*tlen/2 }
        else { tx0 <- at-width*tlen/2; tx1 <- at+width*tlen/2; ty0 <- ty1 <- ticks }
        seg.fn(tx0, ty0, tx1, ty1, lwd=twd, col=tcol)
    }
    if (!is.null(mticks)){
        if (horiz) { tmx0 <- tmx1 <- mticks; tmy0 <- at-width*mlen/2; tmy1 <- at+width*mlen/2  }
        else { tmx0 <- at-width*mlen/2; tmx1 <- at+width*mlen/2; tmy0 <- tmy1 <- mticks }
        seg.fn(tmx0, tmy0, tmx1, tmy1, lwd=mwd, col=mcol)
    }
    invisible()
}

denstrip.normal <- function(mean, sd, log=FALSE, nx=1000, ...){    
    x <-
        if (log) qlnorm(ppoints(nx, a=0), mean, sd)
        else qnorm(ppoints(nx, a=0), mean, sd)
    dens <- if (log) dlnorm(x, mean, sd) else dnorm(x, mean, sd)
    denstrip(x=x, dens=dens, ...)
}

denstrip.legend <- function(x, # central x position 
                            y, # central y position
                            width, len, colmax, colmin="white", gamma=1, horiz=FALSE,
                            max=1, nticks=5, ticks, value.adj=0, cex, main="Density", lattice=FALSE)
{
    if (lattice) {
        poly.fn <- panel.polygon; seg.fn <- panel.segments; text.fn <- panel.text
        default.width <- diff(current.panel.limits()[[if(horiz) "ylim" else "xlim"]]) / 30
        default.len <- diff(trellis.last.object()[[if(horiz) "x.limits" else "y.limits"]]) / 4
        default.colmax <- trellis.par.get("add.line")$col
        default.colmax <- trellis.par.get("background")$col
        default.cex <- trellis.par.get("axis.text")$cex * 0.75
      }
    else {
        poly.fn <- polygon; seg.fn <- segments; text.fn <- text
        default.width <- diff(par("usr")[if(horiz) 3:4 else 1:2]) / 30
        default.len <- diff(par("usr")[if(horiz) 1:2 else 3:4]) / 4
        default.colmax <- par("fg")
        default.cex <- par("cex") * 0.75
    }
    if (missing(width)) width <- default.width
    if (missing(len)) len <- default.len
    if (missing(colmax)) colmax <- default.colmax
    if (missing(cex)) cex <- default.cex
    if (horiz) {        
        pt <- x; at <- y
        xdim <- len; ydim <- width
    }
    else {
      pt <- y; at <- x
      xdim <- width; ydim <- len
    }
    npoints <- 1000 # number of distinct colors. no need to make this an argument
    dx <- seq(pt - len/2, pt + len/2, length=npoints)
    ddens <- seq(0, max, length=npoints)
    denstrip(x=dx, dens=ddens, at=at, width=width, colmax=colmax, colmin=colmin, gamma=gamma, horiz=horiz, lattice=lattice)    
    poly.fn(x = c(x - xdim/2, x + xdim/2, x + xdim/2, x - xdim/2), # draw box around strip
            y = c(y - ydim/2, y - ydim/2, y + ydim/2, y + ydim/2))
    if (missing(ticks)) {
        ticks.at <- seq(pt - len/2, pt + len/2, length=nticks)
        tick.nos <- signif(seq(0, max, length=nticks), 3)
    }
    else {
        ticks.at <- pt - len/2 + ticks / max * len
        tick.nos <- ticks
    }
    if (horiz) seg.fn(ticks.at, at - width/2, ticks.at, at - width*0.75)
    else seg.fn(at + width/2, ticks.at, at + width*0.75, ticks.at)
    if (horiz) text.fn(ticks.at, at - width - value.adj, tick.nos, cex = cex, pos=1)
    else text.fn(at + width + value.adj, ticks.at, tick.nos, cex = cex, pos=4)
    text.fn(x, y + ydim/2, main, cex = cex, pos=3)
    invisible()
}

densregion <- function(x, ...) UseMethod("densregion")

densregion.default <- function(x, # times we have estimates for (vector)
                               y, # vector of values on y axis to show distinct densities at.  same number of distinct y values for each t.
                               z, # matrix of densities at t times y
                               pointwise=FALSE,
                               nlevels=100, # number of distinct densities to show
                               colmax=par("fg"),
                               colmin="white",
                               scale=1,
                               gamma=1,
                               contour=FALSE,
                               ...
                               )
{
    if (pointwise)
        z <- z/apply(z, 1, max)
    qz <- unique(quantile(z, seq(0,1,length=nlevels+1), type=1))
    nlevels <- length(qz) - 1
    zq <- cut(z, qz, include.lowest=TRUE, labels=seq(length=nlevels))
    zq <- matrix(as.numeric(zq), nrow=nrow(z), ncol=ncol(z))
    dens <- tapply(z, zq, mean) # unique densities to plot 
    rgbmax <- col2rgb(colmax, alpha=TRUE)
    rgbmin <- if (colmin=="transparent") c(col2rgb(colmax, alpha=FALSE), 0) else col2rgb(colmin, alpha=TRUE)
    if (gamma <= 0) stop("gamma must be greater than 0")
    p <- (dens / max(dens)) ^ gamma
    if (colmin=="transparent") 
        cols <- rgb(p*rgbmax[1] + (1 - p)*rgbmin[1],
                    p*rgbmax[2] + (1 - p)*rgbmin[2],
                    p*rgbmax[3] + (1 - p)*rgbmin[3],
                    alpha=p*rgbmax[4] + (1 - p)*rgbmin[4], maxColorValue=255)
    else 
        cols <- rgb(p*rgbmax[1] + (1 - p)*rgbmin[1],
                    p*rgbmax[2] + (1 - p)*rgbmin[2],
                    p*rgbmax[3] + (1 - p)*rgbmin[3],
                    alpha=rgbmax[4],
                    maxColorValue=255)
    z <- z[order(x),order(y)]
    x <- sort(x)
    y <- sort(y)
    .filled.contour(as.double(x), as.double(y), scale*z/max(z), qz/max(z), col = cols)
    if (contour)
      contour(x, y, scale*z/max(z), add=TRUE, ...)
    invisible()
}

densregion.survfit <- function(x, ny=20, ...) {
    if (!inherits(x, "survfit")) stop(deparse(substitute(x)), " should be a survfit object")
    if (is.null(x$conf.type))
        stop("No confidence intervals available in \"", deparse(substitute(x)), "\"")
    identity <- function(x)x # this is in base R from 2.7 
    tr <- switch(x$conf.type, "log-log"=function(x)log(-log(x)), "log" = log, plain=identity)
    invtr <- switch(x$conf.type, "log-log"=function(x)exp(-exp(x)), "log" = exp, plain=identity)
    drop <- is.na(x$surv) | is.na(x$lower)
    surv <- x$surv[!drop]; lower <- x$lower[!drop]; time <- x$time[!drop]
    lsurv <- tr(surv)
    lse <- abs((tr(surv) - tr(lower)) / qnorm(0.975))
    n <- length(time)
    if (n==0) stop("No observed events")
    y <- matrix(nrow=n, ncol=ny) 
    ## "ny" ordinates, based on normal quantiles, where density must be calculated for each event time
    for (i in 1:n)
        y[i,] <- invtr(qnorm(ppoints(ny, a=0), lsurv[i], lse[i]))
    ## but actually calculate density at all y for every time, defining a grid on the whole plot region
    yy <- sort(unique(y))
    z <- matrix(nrow=n, ncol=length(yy))
    for (i in 1:n)
        z[i,] <- dnorm(tr(yy), lsurv[i], lse[i])
    densregion.default(x=time, y=yy, z=z, ...)
    invisible()
}

densregion.normal <- function(x, mean, sd, ny=20, ...)
{
    n <- length(x) 
    if (n != length(mean))
        stop(deparse(substitute(mean)), " and ", deparse(substitute(x)), " should be the same length")
    if (n != length(sd))
        stop(deparse(substitute(sd)), " and ", deparse(substitute(x)), " should be the same length")
    y <- matrix(nrow=n, ncol=ny)
    for (i in 1:n)
        y[i,] <- qnorm(ppoints(ny, a=0), mean[i], sd[i])
    yy <- sort(unique(y))
    z <- matrix(nrow=n, ncol=length(yy))
    for (i in 1:n)
        z[i,] <- dnorm(yy, mean[i], sd[i])    
    densregion.default(x=x, y=yy, z=z, ...)
    invisible()
}

seqToIntervals <- function(x){
    x <- sort(unique(as.integer(x)))
    breaks <- x[c(1, diff(x)) > 1]
    groups <- cut(x, c(min(x)-1, breaks - 0.5, max(x)+1), labels = FALSE)
    ranges <- tapply(x, groups, range)
    res <- do.call("rbind", ranges)
    colnames(res) <- c("from","to")
    return(res)
}

sectioned.density <- function(x, dens, at, width, offset, ny,
                              method=c("kernel","frequency"), nx, horiz=TRUE, up.left = TRUE,
                              colmax, colmin="white", gamma=1, lattice=FALSE, ...)
{
    if (lattice) {
        rect.fn <- panel.rect
        default.width <- diff(current.panel.limits()[[if(horiz) "ylim" else "xlim"]]) / 20
        default.colmax <- trellis.par.get("add.line")$col
    }
    else {
        rect.fn <- rect
        default.width <- diff(par("usr")[if(horiz) 3:4 else 1:2]) / 20
        default.colmax <- par("fg")
    }
    if (missing(width)) width <- default.width
    if (missing(colmax)) colmax <- default.colmax
    if (!is.numeric(x)) stop("\'x\' must be numeric")
    if (missing(offset)) offset <- width/3
    if (!up.left) offset <- - offset
    if (!missing(dens)) {
        if (!is.numeric(dens)) stop("\'dens\' must be numeric")
        if (length(dens) != length(x)) stop("Lengths of \'dens\' and \'x\' must be the same")
        de <-  list(x=sort(x), y=dens[order(x)])
    }
    else { 
        method <- match.arg(method)
        if (method=="frequency") {
            if (missing(nx)) nx <- nclass.Sturges(x)
            xcuts <- seq(min(x), max(x), length=nx)
            dens <- table(cut(x, xcuts)) / length(x)
            de <- list(x=sort(x), y=dens[findInterval(x, xcuts,rightmost.closed=TRUE)][order(x)])
        }
        else if (method=="kernel") de <- density(x, ...)
    }
    if (missing(ny)) ny <- nclass.Sturges(de$y)
    ycuts <- seq(0, max(de$y), length=ny+1)
    rgbmax <- col2rgb(colmax, alpha=TRUE)
    rgbmin <- if (colmin=="transparent") c(col2rgb(colmax, alpha=FALSE), 0) else col2rgb(colmin, alpha=TRUE)
    if (gamma <= 0) stop("gamma must be greater than 0")
    p <- seq(0,1,length=ny-1) ^ {1 / gamma}
    cols <- rgb(rgbmax[1] + p*(rgbmin[1] - rgbmax[1]),
                rgbmax[2] + p*(rgbmin[2] - rgbmax[2]),
                rgbmax[3] + p*(rgbmin[3] - rgbmax[3]),
                alpha=rgbmax[4] + p*(rgbmin[4] - rgbmax[4]),
                maxColorValue=255)
    for (i in 2:ny){
        ## draw rectangles for each region of x with density greather than cut-off
        ## one for each contiguous block in ind 
        ind <- which(de$y >= ycuts[i])
        ints <- seqToIntervals(ind) # from R.utils
        for (j in seq(length=nrow(ints))) {
            if (horiz) 
                rect.fn(xleft=de$x[ints[j,1]],
                     xright=de$x[ints[j,2]],
                     ybottom=at + (i-1)*offset - (!up.left)*width,
                     ytop=at + (i-1)*offset + up.left*width,
                     col=cols[i-1], border = (if(i==ny) cols[i-2] else NA) )
            else
                rect.fn(ybottom=de$x[ints[j,1]],
                     ytop=de$x[ints[j,2]],
                     xleft=at - (i-1)*offset - up.left*width,
                     xright=at - (i-1)*offset + (!up.left)*width,
                     col=cols[i-1], border = (if(i==ny) cols[i-2] else NA) )
        }
    }
    invisible()
}

cistrip <- function(x, at, d, horiz=TRUE, pch = 16, cex=1, lattice=FALSE, ...)
{
    if (missing(d)) d <- diff(par("usr")[if(horiz) 3:4 else 1:2]) / 60
    if (is.data.frame(x)) x <- as.matrix(x)
    if (!is.numeric(x)) stop("\'x\' must be numeric")
    if (is.vector(x)) x <- matrix(x, ncol=3)
    n <- nrow(x)
    if (length(at) != n) stop("length of \'at\' should equal the number of estimates in \'x\'")
    pch <- rep(pch, length.out=n)
    cex <- rep(cex, length.out=n)
    if (lattice) {
        points.fn <- panel.points; seg.fn <- panel.segments;
    }
    else {
        points.fn <- points; seg.fn <- segments;
    }
    for (i in 1:n) { 
        if (horiz) { 
            points.fn(x[i,1], at[i], pch=pch[i], cex=cex[i], ...)
            seg.fn(x[i,2], at[i], x[i,3], at[i], ...)
            seg.fn(x[i,2], at[i]-d/2, x[i,2], at[i]+d/2, ...)
            seg.fn(x[i,3], at[i]-d/2, x[i,3], at[i]+d/2, ...)
        }
        else {
            points.fn(at[i], x[i,1], pch=pch[i], cex=cex[i], ...)
            seg.fn(at[i], x[i,2], at[i], x[i,3], ...)
            seg.fn(at[i]-d/2, x[i,2], at[i]+d/2, x[i,2], ...)
            seg.fn(at[i]-d/2, x[i,3], at[i]+d/2, x[i,3], ...)
        }
    }
    invisible()
}

vwstrip <- function(x, dens, at, width, horiz=TRUE, scale=1, limits=c(-Inf,Inf),
                    col="gray", border=NULL, lwd, lty,
                    ticks=NULL, tlen=1, twd, tty, lattice=FALSE, ...) 
{
    if (!is.numeric(x)) stop("\'x\' must be numeric")
    if (missing(dens)) {
        ## x assumed to be a sample from a distribution, density estimated
        de <- density(x, ...)
        x <- de$x; dens <- de$y
    }
    else {
        if (!is.numeric(dens)) stop("\'dens\' must be numeric")
        if (length(dens) != length(x)) stop("Lengths of \'dens\' and \'x\' must be the same")
        dens <- dens[order(x)]
        x <- sort(x)
    }
    if (lattice) {
        poly.fn <- panel.polygon; seg.fn <- panel.segments
        default.width <- diff(current.panel.limits()[[if(horiz) "ylim" else "xlim"]]) / 20
        default.lwd <- default.twd <- trellis.par.get("add.line")$lwd
        default.lty <- default.tty <- trellis.par.get("add.line")$lty
        default.border <- trellis.par.get("add.line")$col
    }
    else {
        poly.fn <- polygon; seg.fn <- segments
        default.width <- diff(par("usr")[if(horiz) 3:4 else 1:2]) / 20
        default.lwd <- default.twd <- par("lwd")
        default.lty <- default.tty <- par("lty")
        default.border <- par("fg")
    }
    dens <- dens[x>limits[1] & x<limits[2]]
    x <- x[x>limits[1] & x<limits[2]]
    if (missing(width)) width <- default.width
    if (missing(lwd)) lwd <- default.lwd
    if (missing(lty)) lty <- default.lty
    if (missing(twd)) twd <- default.twd
    if (missing(tty)) tty <- default.tty
    if (is.null(border)) border <- default.border
    dens <- dens / max(dens) * width/2 * scale
    xx <- c(x, rev(x))
    yy <- c(dens, -rev(dens)) + at 
    if (horiz) poly.fn(xx, yy, col=col, border=border, lty=lty, lwd=lwd)
    else poly.fn(yy, xx, col=col, border=border, lty=lty, lwd=lwd)
    if (!is.null(ticks)) {
        dticks <- dens[findInterval(ticks, x)]
        if (horiz) { tx0 <- tx1 <- ticks; ty0 <- at - dticks*tlen; ty1 <- at + dticks*tlen }
        else { ty0 <- ty1 <- ticks; tx0 <- at - dticks*tlen; tx1 <- at + dticks*tlen }
        seg.fn(tx0, ty0, tx1, ty1, lwd=twd, lty=tty)
    }
    invisible()
}

vwstrip.normal <- function(mean, sd, log=FALSE, nx=1000, ...){    
    x <-
        if (log) qlnorm(ppoints(nx, a=0), mean, sd)
        else qnorm(ppoints(nx, a=0), mean, sd)
    dens <- if (log) dlnorm(x, mean, sd) else dnorm(x, mean, sd)
    vwstrip(x=x, dens=dens, ...)
}

bpstrip <- function(x, prob, at, width, horiz=TRUE, scale=1, limits=c(-Inf,Inf),
                    col="gray", border=NULL, lwd, lty,
                    ticks=NULL, tlen=1, twd, tty, lattice=FALSE)
{
    if (!is.numeric(x)) stop("\'x\' must be numeric")
    if (missing(prob)) {
        x <- sort(x)
        prob <- ecdf(x)(x)
    }
    else {
        if (!is.numeric(prob)) stop("\'prob\' must be numeric")
        if (length(prob) != length(x)) stop("Lengths of \'prob\' and \'x\' must be the same")
        prob <- prob[order(x)]
        x <- sort(x)
    }
    if (lattice) {
        poly.fn <- panel.polygon; seg.fn <- panel.segments
        default.width <- diff(current.panel.limits()[[if(horiz) "ylim" else "xlim"]]) / 20
        default.lwd <- default.twd <- trellis.par.get("add.line")$lwd
        default.lty <- default.tty <- trellis.par.get("add.line")$lty
        default.border <- trellis.par.get("add.line")$col
    }
    else {
        poly.fn <- polygon; seg.fn <- segments
        default.width <- diff(par("usr")[if(horiz) 3:4 else 1:2]) / 20
        default.lwd <- default.twd <- par("lwd")
        default.lty <- default.tty <- par("lty")
        default.border <- par("fg")
    }
    prob <- pmin(prob, 1 - prob)
    prob <- prob[x>limits[1] & x<limits[2]]
    x <- x[x>limits[1] & x<limits[2]]
    if (missing(width)) width <- default.width
    if (missing(lwd)) lwd <- default.lwd
    if (missing(lty)) lty <- default.lty
    if (missing(twd)) twd <- default.twd
    if (missing(tty)) tty <- default.tty
    if (is.null(border)) border <- default.border
    prob <- prob / max(prob) * width/2 * scale
    xx <- c(x, rev(x))
    yy <- c(prob, -rev(prob)) + at 
    if (horiz) poly.fn(xx, yy, col=col, border=border, lty=lty, lwd=lwd)
    else poly.fn(yy, xx, col=col, border=border, lty=lty, lwd=lwd)
    if (!is.null(ticks)) {
        pticks <- prob[findInterval(ticks, x)]
        if (horiz) { tx0 <- tx1 <- ticks; ty0 <- at - pticks*tlen; ty1 <- at + pticks*tlen }
        else { ty0 <- ty1 <- ticks; tx0 <- at - pticks*tlen; tx1 <- at + pticks*tlen }
        seg.fn(tx0, ty0, tx1, ty1, lwd=twd, lty=tty)
    }
    invisible()
}

panel.denstrip <- function(...)
{
    denstrip(..., lattice=TRUE)
}

panel.denstrip.normal <- function(...)
{
    denstrip.normal(..., lattice=TRUE)
}

panel.denstrip.legend <- function(...)
{
    denstrip.legend(..., lattice=TRUE)
}

panel.cistrip <- function(...)
{
    cistrip(..., lattice=TRUE)
}

panel.vwstrip <- function(...)
{
    vwstrip(..., lattice=TRUE)
}

panel.vwstrip.normal <- function(...)
{
    vwstrip.normal(..., lattice=TRUE)
}

panel.bpstrip <- function(...)
{
    bpstrip(..., lattice=TRUE)
}

panel.sectioned.density <- function(...)
{
    sectioned.density(..., lattice=TRUE)
}

Try the denstrip package in your browser

Any scripts or data that you put into this service are public.

denstrip documentation built on May 30, 2017, 7:41 a.m.