[19199] in s-news-athena
[S] Weighted histogram script
daemon@ATHENA.MIT.EDU (Claus Andersen 4912)
Mon Sep 13 06:21:56 1999
Date: Mon, 13 Sep 1999 12:14:53 +0200
From: Claus Andersen 4912 <claus.andersen@risoe.dk>
To: s-news@wubios.wustl.edu
Message-Id: <01JFX6I9XPKW9KQK2O@risoe.dk>
Mime-Version: 1.0
Content-Type: text/plain; charset=US-ASCII
Content-Transfer-Encoding: 7BIT
On September 9, 1999 Christian Hoffmann wrote:
"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..."
and he had some suggestions about how to make "histogram" do the job. I happened
to be in the need of a similar function. I work with indoor radon measurements from
a national survey. We have conducted about the same number of measurements in
each administrative unit regardless of the number of houses in those units. To obtain
an estimate of the (representative) national distribution of radon measurements, we
need a weighted histogram.
Below, I list a function called histogram.w which solves the above problem. It has
been modified from original "histogram" supplied with S-Plus 2000. The function heading
differs from the original only by the vector called weights.
I hope the function is useful. Please let me know if any mistakes are found.
- Claus
----------------------------------------------------------
Claus E. Andersen
Risoe National Laboratory
Dept. of Nuclear Safety Research, Build. NUK-125
DK-4000 Roskilde, Denmark
Phone: +45-4677 4677 (main Risoe number)
Phone: +45-4677 4912 (direct)
Fax: +45-4677 4959 (fax)
E-mail: claus.andersen@risoe.dk
Internet: http://www.risoe.dk/nuk/nuk-clan.htm
----------------------------------------------------------
S-plusscript:
histogram.w<-function(formula, weights=1,
data= sys.parent(1), panel = "panel.histogram", type = "percent",
nint= "Sturges", endpoints = range(x[!na.x]),
ylim=c(0, max(y)), xlim = range(X, endpoints),
breaks= if(is.factor(x)) seq(0.5, length = length(levels(x)) + 1) else
do.breaks(endpoints,nint), xlab = labs[1], ylab = labs[2],
groups= NULL, ..., scales = NULL, subset = T)
{
#Purpose: Weighted histogram using traditional trellis graphics.
#Description: Modified histogram that allows observations to be
#weighted. Created from the orginal histogram function supplied
#with MathSoft's S-Plus 2000 (Rev. July 1, 1999). Changes
#from the orginal histogram are indicated below. The weight vector
#gives the weight of each observation supplied to the function.
#For example, weights=2 means that all observations are
#counted twice. This can be seen if type="counts". As another
#example consider: histogram.w(c(1,2,2),weights=c(10,1,1)). Here
#the "1" observation is counted 10 times whereas the two
#"2" observations are counted one time each. Observe that
#histogram.w(c(1,2,2)) or histogram.w(c(1,2,2),weights=1)
#give the same results as histogram(c(1,2,2)).
#
#Name: Claus E. Andersen, Risoe National Lab., Denmark
#E-mail: claus.andersen@risoe.dk
#Date: September 10, 1999
#Revised: September 13, 1999
#########BEGIN OF CHANGES (part 1)
tabulate.data <- function(x,index.list,type,breaks,num.y,nint,
num.plot,isub.pop,w){
# isub.pop is an index vector with the subpopulation
# isuper (used below) is the index of all values used
# in a given panel
y <- rep(0, num.y * num.plot)
p <- 0
for(isuper in index.list) {
k <- intersect(isuper,isub.pop)
if(any(k)) {
{
X <- x[k]
divisor <- if(type == "percent") sum(w[isuper])/100 else 1
if (divisor==0) divisor <- 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))
y[c(1:nint) + m] <- tabulate(bin, length(levels(bin)))/divisor
y[m] <- y[m + 1]
}}
p <- p + 1
}
return(y)
} # end of tabulate.data
#########END OF CHANGES (part 1)
do.breaks<- function(x, n)
{
r <- range(x)
e <- diff(r) * 0.01
seq(r[1] - e, r[2] + e, length = n + 1)
}
Summary.factor <- function(x)
{
x <- as.numeric(x)
NextMethod()
}
if(!missing(groups))
stop("groups will not work properly in panel.histogram")
if(!is.null(list(...)$sub))
sub <- list(...)$sub
else sub <- NULL
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))
x <- factor(x)
else if(!is.null(levels(x)) && !is.factor(x))
x <- factor(x, levels(x))
if(is.factor(x)) {
scales$x$at <- seq(along = levels(x))
scales$x$labels <- levels(x)
}
subset <- eval(substitute(subset), data)
######### START OF CHANGES (Part 2)
weights <- eval(substitute(weights), data)
# Extend weights to same length as x:
w <- rep(weights,length=length(x))
w <- w[subset]
######### END OF CHANGES (Part 2)
n <- length(x)
x <- x[subset]
na.x <- is.na(x)
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
}
else {
if(is.factor(x))
nint <- length(levels(x))
else {
if(is.character(nint))
nint <- switch(casefold(nint),
sturges = nclass.sturges(x),
fd = nclass.fd(x),
scott = nclass.scott(x),
stop("Nint method not recognized"))
else if(is.function(nint))
nint <- nint(x)
}
}
num.y <- nint + 1
y <- rep(0, num.y * num.plot)
index.list <- do.index.trellis(glist, na.x | attr(glist, "na.index"))
#########START OF CHANGES (part 3)
# The data are split into subpopulations according to weight:
wfac <- factor(w)
i <- 1:length(x)
wsum <- 0
for (k in levels(wfac)){
# Select observations with the same weight:
isub <- i[wfac==k]
ysub<-tabulate.data(x,index.list=index.list,type=type,
breaks=breaks,num.y=num.y,nint=nint,num.plot=num.plot,
isub=isub,w=w)
wmult <- as.numeric(as.character(k))
# The subpopulations are multiplied by the weight and added to y:
y <- y + ysub*wmult
wsum <- wsum + wmult
} # end for loop
# NB: The above loop should be vectorized
#########END OF CHANGES (part 3)
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, sub = sub,
groups = eval(substitute(groups), data), scales = scales, ...)
attr(data, "call") <- match.call()
data
}# END histogram.w
----------------------------------------------------------
Claus E. Andersen
Risoe National Laboratory
Dept. of Nuclear Safety Research, Build. NUK-125
DK-4000 Roskilde, Denmark
Phone: +45-4677 4677 (main Risoe number)
Phone: +45-4677 4912 (direct)
Fax: +45-4677 4959 (fax)
E-mail: claus.andersen@risoe.dk
Internet: http://www.risoe.dk/nuk/nuk-clan.htm
----------------------------------------------------------
-----------------------------------------------------------------------
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
