[19218] in s-news-athena
[S] Results: weighted histogram
daemon@ATHENA.MIT.EDU (Christian Hoffmann)
Tue Sep 14 08:40:52 1999
Message-Id: <3.0.5.32.19990914143247.009fbb50@mail.wsl.ch>
Date: Tue, 14 Sep 1999 14:32:47 +0200
To: s-news@wubios.wustl.edu
From: Christian Hoffmann <christian.hoffmann@lagonda.wsl.ch>
Mime-Version: 1.0
Content-Type: text/plain; charset="us-ascii"
Hi there,
Here I summarize the help which I got concerning the following problem. My
own solution is given below.
>For the evaluation of data of a stratified survey I would like to extend
the function "histogram" so >that its bars represent the weighted sum of
occurrences of the data in the bins. This is because the >probability for
an individual to be represented in the survey is a different constant for
each stratum. >( => "weights" <- 1/probability). If all "weights" == 1,
then the resulting plot should reproduce the >plot generated by the
original version of "histogram".
> I am using S+3.4 on Unix (Sun/Solaris)
1)
Claus Andersen 4912 <claus.andersen@risoe.dk> furnished a version which
solves the problem exactly as my program below.
2)
Tim Hesterberg <timh@statsci.com>, MathSoft wrote: "I have code that does
weighted histograms, calling a weighted version of the C function
tabulate", see http://www.statsci.com/Hesterberg/tilting/
3)
Peter Perkins <peter@caliban.ucsd.edu> and "Buttrey, Samuel"
<sebuttre@nps.navy.mil> contributed the idea to use tapply.
4)
Frank E Harrell Jr <fharrell@virginia.edu> wrote:
You might look at the wtd.stats in the Hmisc library. I don't have a
weighted histogram but there is a weighted empirical CDF and weighted
table()-type function, see "hesweb1.med.virginia.edu-/biostat/s/unix/" .
-------------------------------------------------------------------------
weighted.histogram <- function(formula, weights = NULL, data =
sys.parent(1), panel = "panel.histogram", type = "percent", nint =
round(log(length(x), base = 2) + 1), endpoints = range(x[!na.x]), ylim =
c(0, max(y)), xlim = range(X, endpoints), breaks = do.breaks(endpoints,
nint), xlab = labs[1], ylab = labs[2], groups = NULL, ..., subset = T)
{
# Present weighted histograms using trellis graphics.
# Created from the orginal histogram function supplied with S+3.4, Changes
are marked below.
# If the weights are missing, the original function will be executed.
# If the weights are all 1, the result will be that of the original function.
# If the weights are different from 1 (arbitrary real), then the histogram
# bars are proportional not to the number of observations, but to the
# sum of the weights corresponding to the observations.
# NA in the weights are NOT allowed.
# If weights are fewer than observations in formula, then weights will be
# recycled to the length of the observations.
#
#Name: Christian W. Hoffmann
#E-mail: hoffmann@WSL.ch
#Date: 1999.09.14
#
do.breaks <- function(x, n)
{
r <- range(x)
e <- diff(r) * 0.01
seq(r[1] - e, r[2] + e, length = n + 1)
}
if(!missing(groups))
stop("groups will not work properly in panel.histogram")
if(!is.null(list(...)$subscripts))
stop("subscripts will not work properly in panel.histogram")
sub.formula <- substitute(formula)
formula <- eval(sub.formula, data)
Z <- do.formula.trellis(formula)
if(!Z$no.response)
stop("formula should be in the form of ~x|g1*g2*...")
expr <- Z$expr
no.given <- is.null(expr)
labs <- c(if(is.numeric(formula)) deparse(sub.formula) else Z$ylab, if(
type == "percent") "Percent of Total" else "Count", "", "")
x <- eval(if(no.given) formula else expr[[1]], data)
if(is.character(x) || !is.null(levels(x)))
stop("first variable in formula must be numeric")
subset <- eval(substitute(subset), data)
n <- length(x)
x <- x[subset]
na.x <- is.na(x)
##### Start of changes part 1
if (!missing(weights)) {
weights <- eval(substitute(weights), data)
# Extend weights to same length as x:
w <- rep(weights,length=length(x))
weights <- w[subset]
}
##### End of changes part 1
glist <- proc.given.trellis(expr[-1], data, n, na.index = na.x, subset =
subset)
index <- attr(glist, "index")
num.plot <- nrow(index)
num.given <- ncol(index)
if(!missing(breaks)) {
if(!missing(nint))
stop("specify either breaks or nint, but not both")
else nint <- length(breaks) - 1
}
num.y <- nint + 1
y <- rep(0, num.y * num.plot)
index.list <- do.index.trellis(glist, na.x | attr(glist, "na.index"))
p <- 0
for(k in index.list) {
if(any(k)) {
X <- x[k]
divisor <- if(type == "percent") length(X)/100 else 1
m <- p * num.y + 1
bin <- cut(X[!is.na(X)], breaks)
if(any(is.na(bin)))
stop(paste(
"breaks do not span the range of x in plot",
p + 1))
##### Start of changes part 2
indy <- c(1:nint) + m
if (missing(weights)) {
y[indy] <- tabulate(bin, length(levels(bin)))/divisor ## is original
}
else { ## sum the weights corresponding to the x's
W <- weights[k]
if (any(is.na(W))) stop ("weights must not contain NA")
divisor <- if (type == "percent") sum(W)/100 else 1
yy <- tapply(W,bin,sum)/divisor ## produces NA insted of 0
y[indy][!is.na(yy)] <- yy[!is.na(yy)] ## use only not-NA
##### End of changes part 2
}
y[m] <- y[m + 1]
}
p <- p + 1
}
grid <- expand.grid(c(list(breaks), lapply(glist, function(x) x$levels)))
X <- grid[[1]]
for(i in 1:num.given) {
glist[[i]]$given <- c(grid[[i + 1]])
glist[[i]]$is.factor <- T
}
attr(glist, "is.factor") <- rep(T, num.given)
data <- packlist.trellis(X, y, glist, panel = panel, labs = labs, xlab
= xlab, ylab = ylab, xlim = xlim, ylim = ylim, groups = eval(
substitute(groups), data), ...)
attr(data, "call") <- match.call()
data
}
-------------------------------------------------------------------------
--Christian
Christian W. Hoffmann
Swiss Federal Institute for Forest, Snow and Landscape Research
CH-8903 Birmensdorf, Switzerland
phone: ++41-1-739 22 77 fax : ++41-1-737 40 80 e-mail: Hoffmann@WSL.CH
WWW: http://www.wsl.ch/land/dynamics/modeling/hoffmann.htm
FTP: ftp://ftp.wsl.ch/pub/hoffmann/
I prefer WWW pages which can be viewed best with *any* browser, contain
*no* cookies, *no* Java, *no* frames, and *no* required images !
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news
