I present you with a function that solves a problem that has bugged me for
many years. I think the problem may be general enough to at least consider
adding this function, or a revamped version of it, to the 'stats'
package,
with the other multinomial functions reside.
I'm using R to export data to text files, which are input data for an
external model written in C++. Parts of the data are age distributions, in
the form of relative frequency in each year:
Year Age1 Age2 ... Age10
1980 0.123 0.234 ... 0.001
... ... ... ... ...
Each row should sum to exactly 1. The problem is that when I preprocess
each line in R as p<-a/sum(a), occasionally a line will sum to 0.999,
1.002, or the like. This could either crash the external model or lead to
wrong conclusions.
I believe similar partitioning is commonly used in a wide variety of
models, making this a general problem for many modellers.
In the past, I have checked every line manually, and then arbitrarily
tweaked one or two values up or down to make the row sum to exactly one,
but two people would tweak differently. Another semi-solution is to write
the values to the text file in a very long format, but this would (1) make
it harder to visually check the numbers and (2) the numbers in the article
or report would no longer match the data files exactly, so other
scientists could not repeat the analysis and get the same results.
Once I implemented a quick and dirty solution, simply setting the last
proportion (Age10 above) as 1 minus the sum of ages 1-9. I quickly stopped
using that approach when I started seeing negative values.
After this introduction, the attached round_multinom.html should make
sense. The algorithm I ended up choosing comes from allocating seats in
elections, so I was tempted to provide that application as well, although
it makes the interface and documentation slightly more confusing.
The working title of this function was a short and catchy vote(), but I
changed it to round_multinom(), even though it's not matrix-oriented like
the other *multinom functions. That would probably be straightforward to
do, but I'll keep it as a vector function during the initial discussion.
I'm curious to hear your impressions and ideas. In the worst case, this is
a not-so-great solution to a marginal problem. In the best case, this
might be worth a short note in the Journal of Statistical Software.
Thanks for your time,
Arni
P.S. In case the mailing list doesn't handle attachments, I've placed
the
same files on http://www.hafro.is/~arnima/ for your convenience.
-------------- next part --------------
\name{round_multinom}
\alias{round_multinom}
\encoding{UTF-8}
\title{
Round Multinomial Proportions (or Allocate Seats from Election
Results)
}
\description{
This function can round multinomial proportions to a given number of
decimal places, while making sure the rounded proportions sum to
exactly one. This is achieved using one of three algorithms that were
originally invented to allocate seats from election results.
It is often necessary to round proportions, e.g. to produce legible
percentages for an article or a presentation. Rounding also takes
place when data are written to a text file, to be analyzed by an
external model. The rounded proportions often fail to add to exactly
one,
\preformatted{
a <- c(67630, 116558, 207536, 251555, 356721)
p <- round(a/sum(a), 3) # 0.068 0.117 0.208 0.252 0.357
sum(p) # 1.002}
which would make the rounded proportions illegal input data for many
models. Instead of manually checking and arbitrarily tweaking
proportions so they add to exactly one, this function can guarantee
that condition, using an unbiased algorithm.
}
\usage{
round_multinom(x, digits = NULL, seats = NULL, method="SL", labels =
names(x))
}
\arguments{
\item{x}{vector containing multinomial proportions or counts.}
\item{digits}{
number of decimal places to use when rounding multinomial
proportions.
}
\item{seats}{number of seats to allocate from election results.}
\item{method}{
string specifying the algorithm to use: \code{"DH"},
\code{"MSL"},
or \code{"SL"}.
}
\item{labels}{optional vector of names for the output vector.}
}
\details{
This function should be called \emph{either} with a \code{digits}
argument to round multinomial proportions, \emph{or} with a
\code{seats} argument to allocate seats from election results, not
both.
The algorithms are variations of the \dQuote{highest averages} method
for allocating seats proportionally from multiparty election results:
\describe{
\item{\code{"DH"}}{
d'Hondt method, involves the series 1, 2, 3, \ldots, \eqn{n}.
Favors big parties.
}
\item{\code{"MSL"}}{
Modified Sainte-Lagu??, involves the series 1, 2.4, 3.8, \ldots,
1.4\eqn{n}--0.4. Favors big parties slightly.
}
\item{\code{"SL"}}{
Sainte-Lagu??. Involves the series 1, 3, 5, \ldots, 2\eqn{n}--1.
Does not favor big or small parties. A reasonable default method
for rounding multinomial proportions.
}
}
}
\value{
Vector of same length as \code{x}, with rounded numbers whose sum is
one (or integers whose sum is \code{seats}).
}
\note{
d'Hondt is used to allocate parliamentary seats in most of Europe and
South America, East Timor, Israel, Japan, and Turkey.
Modified Sainte-Lagu?? is used in Norway and Sweden.
Sainte-Lagu?? is used in Bosnia and Herzegovina, Kosovo, Latvia, and
New Zealand.
}
\author{Arni Magnusson.}
\references{
Balinski, M. and V. Ram??rez. 1999. Parametric methods of
apportionment, rounding and production. \emph{Mathematical Social
Sciences} \bold{37}:107--122.
Diaconis, P. and D. Freedman. 1979. On rounding percentages.
\emph{Journal of the American Statistical Association}
\bold{74}:359--364.
Mosteller, F., C. Youtz, and D. Zahn. 1967. The distribution of sums
of rounded percentages. \emph{Demography} \bold{4}:850--858.
\url{http://en.wikipedia.org/wiki/Highest_averages_method} (accessed
10 Feb 2010).
}
\seealso{
\code{\link{round}} is the standard rounding function in \R.
\code{\link{rmultinom}} generates random multinomial vectors.
}
\examples{
## Guarantee that rounded multinomial proportions sum to 1:
a <- c(67630, 116558, 207536, 251555, 356721)
p <- a / sum(a)
p1 <- round(p, 3) # 0.068 0.117 0.208 0.252 0.357
sum(p1) # 1.002, no good
p2 <- round_multinom(p, 3) # 0.068 0.117 0.207 0.251 0.357
sum(p2) # 1
## The multinomial "proportions" can also be raw counts (a, not p):
p3 <- round_multinom(a, 3) # 0.068 0.117 0.207 0.251 0.357
sum(p3) # 1
## Allocate 9 seats from 178 votes using different methods:
votes <- c(red=66, green=80, blue=32)
round_multinom(votes, seats=9, method="DH") # 4 4 1
round_multinom(votes, seats=9, method="MSL") # 3 4 2
round_multinom(votes, seats=9, method="SL") # 3 4 2
}
\keyword{arith}
\keyword{distribution}
-------------- next part --------------
round_multinom <- function(x, digits=NULL, seats=NULL, method="SL",
labels=names(x))
{
method <- match.arg(toupper(method),
c("DH","MSL","SL"))
if(is.null(digits) && is.null(seats) || !is.null(digits) &&
!is.null(seats))
stop("Please pass a value as 'digits' or 'seats', not
both")
if(!is.null(digits))
{
if(digits<0 || digits>6)
stop("Please pass a positive value (0-6) as 'digits'")
n <- as.integer(10^digits)
}
else
n <- seats
party <- seq_along(x)
series <- switch(method,
DH = 1 + 1 *(seq_len(n)-1), # 1, 2, 3, ..., n
MSL = 1 + 1.4*(seq_len(n)-1), # 1, 2.4, 3.8, ..., 1.4n-0.4
SL = 1 + 2 *(seq_len(n)-1)) # 1, 3, 5, ..., 2n-1
output <- data.frame(party=rep(party,each=n), score=as.numeric(sapply(x,
function(votes) votes/series)))
output <- factor(output$party)[order(-output$score)][seq_len(n)]
output <- as.integer(table(output))
if(!is.null(digits))
output <- output / n
names(output) <- labels
return(output)
}