[Vegan-commits] r1478 - in branches/1.17: R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Feb 5 19:15:08 CET 2011
Author: jarioksa
Date: 2011-02-05 19:15:08 +0100 (Sat, 05 Feb 2011)
New Revision: 1478
Modified:
branches/1.17/R/as.preston.R
branches/1.17/R/kendall.post.R
branches/1.17/R/nestedtemp.R
branches/1.17/R/ordiplot3d.R
branches/1.17/R/orditkplot.R
branches/1.17/R/plot.cascadeKM.R
branches/1.17/R/plot.cca.R
branches/1.17/R/prestondistr.R
branches/1.17/R/prestonfit.R
branches/1.17/R/print.MOStest.R
branches/1.17/R/rgl.renyiaccum.R
branches/1.17/R/showvarparts.R
branches/1.17/inst/ChangeLog
branches/1.17/man/fisherfit.Rd
Log:
merge r1471 to 1475
Modified: branches/1.17/R/as.preston.R
===================================================================
--- branches/1.17/R/as.preston.R 2011-02-02 17:05:13 UTC (rev 1477)
+++ branches/1.17/R/as.preston.R 2011-02-05 18:15:08 UTC (rev 1478)
@@ -1,16 +1,34 @@
-"as.preston" <-
- function (x, ...)
+`as.preston` <-
+ function (x, tiesplit = FALSE, ...)
{
if (inherits(x, "preston"))
return(x)
if (!identical(all.equal(x, round(x)), TRUE))
stop("function accepts only integers (counts)")
- freq <- x[x > 0]
- freq <- ceiling(log2(freq))
- freq <- table(freq)
- nm <- names(freq)
- freq <- as.vector(freq)
- names(freq) <- nm
+ x <- x[x > 0]
+ if (tiesplit) {
+ ## Assume log2(2^k) == k exactly for integer k
+ xlog2 <- log2(x)
+ ties <- xlog2 == ceiling(xlog2)
+ tiefreq <- table(xlog2[ties])
+ notiefreq <- table(ceiling(xlog2[!ties]))
+ itie <- as.numeric(names(tiefreq)) + 1
+ nitie <- as.numeric(names(notiefreq)) + 1
+ freq <- numeric(max(itie, nitie))
+ ## split tied values between two adjacent octaves
+ freq[itie] <- tiefreq/2
+ freq[itie+1] <- freq[itie+1] + tiefreq/2
+ freq[nitie] <- freq[nitie] + notiefreq
+ } else {
+ xlog2 <- ceiling(log2(x))
+ tmp <- table(xlog2)
+ indx <- as.numeric(names(tmp)) + 1
+ freq <- numeric(max(indx))
+ freq[indx] <- tmp
+ }
+ names(freq) <- seq_along(freq) - 1
+ ## remove empty octaves
+ freq <- freq[freq>0]
class(freq) <- "preston"
freq
}
Modified: branches/1.17/R/kendall.post.R
===================================================================
--- branches/1.17/R/kendall.post.R 2011-02-02 17:05:13 UTC (rev 1477)
+++ branches/1.17/R/kendall.post.R 2011-02-05 18:15:08 UTC (rev 1478)
@@ -43,7 +43,7 @@
##CC# Initialize the counters
counter <- as.list(1:ngr)
for(i in 1:ngr){
- counter[[i]] <- vector(l=n.per.gr[i])
+ counter[[i]] <- vector(length = n.per.gr[i])
}
W.gr <- vector("list",ngr)
if(ngr > 1) spear.gr <- vector("list",ngr)
Modified: branches/1.17/R/nestedtemp.R
===================================================================
--- branches/1.17/R/nestedtemp.R 2011-02-02 17:05:13 UTC (rev 1477)
+++ branches/1.17/R/nestedtemp.R 2011-02-05 18:15:08 UTC (rev 1478)
@@ -11,7 +11,7 @@
ind <- matrix(rep(rr, ncol(x)), nrow=nrow(x))
s <- -colSums((x*ind)^2)
t <- -colSums((nrow(x) - (1-x)*ind + 1)^2)
- st <- rank(s+t, ties="random")
+ st <- rank(s+t, ties.method = "random")
st
}
rowpack <- function(x, cr)
@@ -19,15 +19,15 @@
ind <- matrix(rep(cr, each=nrow(x)), nrow=nrow(x))
s <- -rowSums((x*ind)^2)
t <- -rowSums((ncol(x) - (1-x)*ind + 1)^2)
- st <- rank(s+t, ties="random")
+ st <- rank(s+t, ties.method = "random")
st
}
comm <- ifelse(comm > 0, 1, 0)
## Start with columns, expect if nrow > ncol
if (ncol(comm) >= nrow(comm)) {
- i <- rank(-rowSums(comm), ties="average")
+ i <- rank(-rowSums(comm), ties.method = "average")
} else {
- j <- rank(-colSums(comm), ties="average")
+ j <- rank(-colSums(comm), ties.method = "average")
i <- rowpack(comm, j)
}
## Improve eight times
Modified: branches/1.17/R/ordiplot3d.R
===================================================================
--- branches/1.17/R/ordiplot3d.R 2011-02-02 17:05:13 UTC (rev 1477)
+++ branches/1.17/R/ordiplot3d.R 2011-02-05 18:15:08 UTC (rev 1478)
@@ -31,7 +31,7 @@
mul <- ordiArrowMul(cbind(tmp$x, tmp$y), fill=1)
bp.xyz <- pl$xyz.convert(bp * mul)
orig <- pl$xyz.convert(0, 0, 0)
- arrows(orig$x, orig$y, bp.xyz$x, bp.xyz$y, len = arr.len,
+ arrows(orig$x, orig$y, bp.xyz$x, bp.xyz$y, length = arr.len,
col = arr.col)
}
}
Modified: branches/1.17/R/orditkplot.R
===================================================================
--- branches/1.17/R/orditkplot.R 2011-02-02 17:05:13 UTC (rev 1477)
+++ branches/1.17/R/orditkplot.R 2011-02-05 18:15:08 UTC (rev 1478)
@@ -266,11 +266,11 @@
eps = postscript(file=fname, width=xy$dim[1], height=xy$dim[2],
paper="special", horizontal = FALSE),
pdf = pdf(file=fname, width=xy$dim[1], height=xy$dim[2]),
- png = png(file=fname, width=pixdim[1], height=pixdim[2]),
- jpg = jpeg(file=fname, width=pixdim[1], height=pixdim[2],
+ png = png(filename=fname, width=pixdim[1], height=pixdim[2]),
+ jpg = jpeg(filename=fname, width=pixdim[1], height=pixdim[2],
quality = 100),
- tiff = tiff(file=fname, width=pixdim[1], height=pixdim[2]),
- bmp = bmp(file=fname, width=pixdim[1], height=pixdim[2]),
+ tiff = tiff(filename=fname, width=pixdim[1], height=pixdim[2]),
+ bmp = bmp(filename=fname, width=pixdim[1], height=pixdim[2]),
fig = xfig(file=fname, width=xy$dim[1], height=xy$dim[2]))
plot.orditkplot(xy)
dev.off()
Modified: branches/1.17/R/plot.cascadeKM.R
===================================================================
--- branches/1.17/R/plot.cascadeKM.R 2011-02-02 17:05:13 UTC (rev 1477)
+++ branches/1.17/R/plot.cascadeKM.R 2011-02-05 18:15:08 UTC (rev 1478)
@@ -25,8 +25,6 @@
if (sortg) {
x <- orderingKM(x)
}
- else {
- }
main = (paste("K-means partitions comparison"))
xlab = ("Number of groups in each partition")
ylab = ("Objects")
@@ -42,7 +40,7 @@
grid(nx = nrow(x), ny = max.g - min.g + 1, col = gridcol)
box()
axis(2, seq(min.g - min.g + 1, max.g - min.g + 1, by = 1),
- lab = seq(min.g, max.g, by = 1))
+ labels = seq(min.g, max.g, by = 1))
axis(1)
par(mar = c(5, 2, 5, 1))
par(bg = "white", fg = "black", col = "black")
@@ -51,7 +49,7 @@
0.5, max.g + 0.5), yaxs = "i", yaxt = "n", xlab = "Values")
grid(nx = NULL, ny = max.g - min.g + 1, col = gridcol)
box()
- axis(2, seq(min.g, max.g, by = 1), lab = seq(min.g, max.g,
+ axis(2, seq(min.g, max.g, by = 1), labels = seq(min.g, max.g,
by = 1), col.axis = "black")
axis(1)
maxx = which.max(w[])
Modified: branches/1.17/R/plot.cca.R
===================================================================
--- branches/1.17/R/plot.cca.R 2011-02-02 17:05:13 UTC (rev 1477)
+++ branches/1.17/R/plot.cca.R 2011-02-05 18:15:08 UTC (rev 1478)
@@ -88,7 +88,7 @@
else mul <- 1
attr(g$biplot, "arrow.mul") <- mul
arrows(0, 0, mul * g$biplot[, 1], mul * g$biplot[, 2],
- len = 0.05, col = "blue")
+ length = 0.05, col = "blue")
text(1.1 * mul * g$biplot, rownames(g$biplot), col = "blue")
axis(3, at = c(-mul, 0, mul), labels = rep("", 3), col = "blue")
axis(4, at = c(-mul, 0, mul), labels = c(-1, 0, 1), col = "blue")
Modified: branches/1.17/R/prestondistr.R
===================================================================
--- branches/1.17/R/prestondistr.R 2011-02-02 17:05:13 UTC (rev 1477)
+++ branches/1.17/R/prestondistr.R 2011-02-05 18:15:08 UTC (rev 1478)
@@ -1,9 +1,9 @@
-"prestondistr" <-
+`prestondistr` <-
function (x, truncate = -1, ...)
{
fun <- function(par, x, truncate) {
up <- dnorm(x, par[1], par[2], log = TRUE)
- dn <- pnorm(truncate, par[1], par[2], lower = FALSE)
+ dn <- pnorm(truncate, par[1], par[2], lower.tail = FALSE)
-sum(up - log(dn))
}
x <- x[x > 0]
@@ -11,10 +11,10 @@
p <- c(mean(logx), sd(logx))
sol <- optim(p, fun, x = logx, truncate = truncate)
p <- sol$par
- area <- pnorm(truncate, p[1], p[2], lower = FALSE)
+ area <- pnorm(truncate, p[1], p[2], lower.tail = FALSE)
scale <- length(x)/sqrt(2 * pi)/p[2]/area
p <- c(p, scale)
- oct <- as.preston(x)
+ oct <- as.preston(x, ...)
x <- as.numeric(names(oct))
fit <- p[3] * exp(-(x - p[1])^2/2/p[2]^2)
names(p) <- c("mode", "width", "S0")
Modified: branches/1.17/R/prestonfit.R
===================================================================
--- branches/1.17/R/prestonfit.R 2011-02-02 17:05:13 UTC (rev 1477)
+++ branches/1.17/R/prestonfit.R 2011-02-05 18:15:08 UTC (rev 1478)
@@ -1,9 +1,10 @@
-"prestonfit" <-
- function (x, ...)
+`prestonfit` <-
+ function (x, tiesplit = FALSE, ...)
{
- x <- as.preston(x)
+ x <- as.preston(x, tiesplit = tiesplit)
oct <- as.numeric(names(x))
- fit <- glm(x ~ oct + I(oct^2), family = poisson)
+ fit <- glm(x ~ oct + I(oct^2),
+ family = if (tiesplit) quasipoisson else poisson)
fv <- fitted(fit)
p <- coef(fit)
if (!is.na(p[3]) && p[3] < 0) {
@@ -18,6 +19,8 @@
names(p) <- c("mode", "width", "S0")
out <- list(freq = unclass(x), fitted = fv, coefficients = p)
out$method = "Poisson fit to octaves"
+ if(tiesplit)
+ out$method <- paste("Quasi-", out$method, sep="")
class(out) <- c("prestonfit")
out
}
Modified: branches/1.17/R/print.MOStest.R
===================================================================
--- branches/1.17/R/print.MOStest.R 2011-02-02 17:05:13 UTC (rev 1477)
+++ branches/1.17/R/print.MOStest.R 2011-02-05 18:15:08 UTC (rev 1478)
@@ -8,6 +8,6 @@
if (!x$isBracketed)
cat("***** Caution: hump/pit not bracketed by the data ******\n")
cat("\n")
- printCoefmat(coef(x), has.P=TRUE, na.print="")
+ printCoefmat(coef(x), has.Pvalue = TRUE, na.print = "")
invisible(x)
}
Modified: branches/1.17/R/rgl.renyiaccum.R
===================================================================
--- branches/1.17/R/rgl.renyiaccum.R 2011-02-02 17:05:13 UTC (rev 1477)
+++ branches/1.17/R/rgl.renyiaccum.R 2011-02-05 18:15:08 UTC (rev 1478)
@@ -11,7 +11,7 @@
ylen <- ylim[2] - ylim[1] + 1
colorlut <- rainbow(ylen)
col <- colorlut[1000*y-ylim[1]+1]
- rgl.bg(col="white")
+ rgl.bg(color = "white")
rgl.surface(xp, z, y, color=col)
y <- x[,,5] * rgl.height
##rgl.surface(xp,z,y,color="grey", alpha=0.3)
Modified: branches/1.17/R/showvarparts.R
===================================================================
--- branches/1.17/R/showvarparts.R 2011-02-02 17:05:13 UTC (rev 1477)
+++ branches/1.17/R/showvarparts.R 2011-02-05 18:15:08 UTC (rev 1478)
@@ -15,9 +15,9 @@
box()
symbols(cp, circles = rep(rad, min(parts,3)), inches = FALSE, add=TRUE, ...)
if (parts == 4) {
- symbols(0, 0.2, rectangles=cbind(1, 0.5), inch=FALSE, add=TRUE, ...)
+ symbols(0, 0.2, rectangles=cbind(1, 0.5), inches=FALSE, add=TRUE, ...)
symbols(sqrt(1/2), -sqrt(3/4)+0.2, rectangles=cbind(0.5,0.3),
- inch=FALSE, add=TRUE, ...)
+ inches=FALSE, add=TRUE, ...)
}
nlabs <- switch(parts, 2, 4, 8, 16)
if (missing(labels))
Modified: branches/1.17/inst/ChangeLog
===================================================================
--- branches/1.17/inst/ChangeLog 2011-02-02 17:05:13 UTC (rev 1477)
+++ branches/1.17/inst/ChangeLog 2011-02-05 18:15:08 UTC (rev 1478)
@@ -4,6 +4,12 @@
Version 1.17-7 (opened January 10, 2011)
+ * merged r1476: zip data files.
+
+ * merged r1475: partially matched argunments in functions.
+
+ * merged r1471,3,4: split "ties" in prestonfit (and allies).
+
* merged revs 1444 to 1467:
* fixes in anova.cca scope, subset and NA handling:
Modified: branches/1.17/man/fisherfit.Rd
===================================================================
--- branches/1.17/man/fisherfit.Rd 2011-02-02 17:05:13 UTC (rev 1477)
+++ branches/1.17/man/fisherfit.Rd 2011-02-05 18:15:08 UTC (rev 1478)
@@ -27,13 +27,14 @@
\method{confint}{fisherfit}(object, parm, level = 0.95, ...)
\method{profile}{fisherfit}(fitted, alpha = 0.01, maxsteps = 20, del = zmax/5,
...)
-prestonfit(x, ...)
+prestonfit(x, tiesplit = FALSE, ...)
prestondistr(x, truncate = -1, ...)
\method{plot}{prestonfit}(x, xlab = "Frequency", ylab = "Species", bar.col = "skyblue",
line.col = "red", lwd = 2, ...)
\method{lines}{prestonfit}(x, line.col = "red", lwd = 2, ...)
veiledspec(x, ...)
as.fisher(x, ...)
+as.preston(x, tiesplit = FALSE, ...)
}
\arguments{
@@ -45,6 +46,8 @@
\item{alpha}{The extend of profiling as significance.}
\item{maxsteps}{Maximum number of steps in profiling.}
\item{del}{Step length.}
+ \item{tiesplit}{Split frequencies \eqn{1, 2, 4, 8} etc between adjacent
+ octaves.}
\item{truncate}{Truncation point for log-Normal model, in log2
units. Default value \eqn{-1} corresponds to the left border of zero
Octave. The choice strongly influences the fitting results.}
@@ -52,7 +55,9 @@
\item{bar.col}{Colour of data bars.}
\item{line.col}{Colour of fitted line.}
\item{lwd}{Width of fitted line.}
- \item{\dots}{Other parameters passed to functions. }
+ \item{\dots}{Other parameters passed to functions. Ignored in
+ \code{prestonfit} and \code{tiesplit} passed to \code{as.preston} in
+ \code{prestondistr}. }
}
\details{
In Fisher's logarithmic series the expected
@@ -83,52 +88,73 @@
are similar.
Preston (1948) was not satisfied with Fisher's model which seemed to
- imply infinite species richness, and postulated that rare species is a
- diminishing class and most species are in the middle of frequency
+ imply infinite species richness, and postulated that rare species is
+ a diminishing class and most species are in the middle of frequency
scale. This was achieved by collapsing higher frequency classes into
wider and wider ``octaves'' of doubling class limits: 1, 2, 3--4,
- 5--8, 9--16 etc. occurrences. Any logseries data will look like
- lognormal when plotted this way. The expected frequency \eqn{f} at abundance
- octave \eqn{o} is defined by \eqn{f_o = S_0 \exp(-(\log_2(o) -
- \mu)^2/2/\sigma^2)}{f = S0 exp(-(log2(o)-mu)^2/2/sigma^2)}, where
- \eqn{\mu} is the location of the mode and \eqn{\sigma} the width, both
- in \eqn{\log_2}{log2} scale, and \eqn{S_0}{S0} is the expected number
- of species at mode. The lognormal model is usually truncated on the
- left so that some rare species are not observed. Function
+ 5--8, 9--16 etc. occurrences. It seems that Preston regarded
+ frequencies 1, 2, 4, \emph{etc.}. as \dQuote{tied} between octaves
+ (Williamson & Gaston 2005). This means that only half of the species
+ with frequency 1 are shown in the lowest octave, and the rest are
+ transferred to the second octave. Half of the species from the
+ second octave are transferred to the higher one as well, but this is
+ usually not as large a number of species. This practise makes data
+ look more lognormal by reducing the usually high lowest
+ octaves. This can be achieved by setting argument \code{tiesplit =
+ TRUE}. With \code{tiesplit = FALSE} the frequencies are not split,
+ but all ones are in the lowest octave, all twos in the second, etc.
+ Williamson & Gaston (2005) discuss alternative definitions in
+ detail, and they should be consulted for a critical review of
+ log-Normal model.
+
+ Any logseries data will look like lognormal when plotted in
+ Preston's way. The expected frequency \eqn{f} at abundance octave
+ \eqn{o} is defined by \eqn{f_o = S_0 \exp(-(\log_2(o) -
+ \mu)^2/2/\sigma^2)}{f = S0 exp(-(log2(o)-mu)^2/2/sigma^2)}, where
+ \eqn{\mu} is the location of the mode and \eqn{\sigma} the width,
+ both in \eqn{\log_2}{log2} scale, and \eqn{S_0}{S0} is the expected
+ number of species at mode. The lognormal model is usually truncated
+ on the left so that some rare species are not observed. Function
\code{prestonfit} fits the truncated lognormal model as a second
- degree log-polynomial to the octave pooled data using Poisson
- error. Function \code{prestondistr} fits left-truncated Normal distribution to
- \eqn{\log_2}{log2} transformed non-pooled observations with direct
- maximization of log-likelihood. Function \code{prestondistr} is
- modelled after function \code{\link[MASS]{fitdistr}} which can be used
- for alternative distribution models. The functions have common \code{print},
- \code{plot} and \code{lines} methods. The \code{lines} function adds
- the fitted curve to the octave range with line segments showing the
- location of the mode and the width (sd) of the response.
+ degree log-polynomial to the octave pooled data using Poisson (when
+ \eqn{tiesplit = FALSE}) or quasi-Poisson (when \eqn{tiesplit =
+ TRUE}). error. Function \code{prestondistr} fits left-truncated
+ Normal distribution to \eqn{\log_2}{log2} transformed non-pooled
+ observations with direct maximization of log-likelihood. Function
+ \code{prestondistr} is modelled after function
+ \code{\link[MASS]{fitdistr}} which can be used for alternative
+ distribution models.
- The total
- extrapolated richness from a fitted Preston model can be found with
- function \code{veiledspec}. The function accepts results both from
- \code{prestonfit} and from \code{prestondistr}. If \code{veiledspec} is
- called with a species count vector, it will internally use
- \code{prestonfit}. Function \code{\link{specpool}} provides
- alternative ways of estimating the number of unseen species. In fact,
- Preston's lognormal model seems to be truncated at both ends, and this
- may be the main reason why its result differ from lognormal models
- fitted in Rank--Abundance diagrams with functions
- \code{\link{rad.lognormal}}.
+ The functions have common \code{print}, \code{plot} and \code{lines}
+ methods. The \code{lines} function adds the fitted curve to the
+ octave range with line segments showing the location of the mode and
+ the width (sd) of the response. Function \code{as.preston}
+ transforms abundance data to octaves. Argument \code{tiesplit} will
+ not influence the fit in \code{prestondistr}, but it will influence
+ the barplot of the octaves.
+
+ The total extrapolated richness from a fitted Preston model can be
+ found with function \code{veiledspec}. The function accepts results
+ both from \code{prestonfit} and from \code{prestondistr}. If
+ \code{veiledspec} is called with a species count vector, it will
+ internally use \code{prestonfit}. Function \code{\link{specpool}}
+ provides alternative ways of estimating the number of unseen
+ species. In fact, Preston's lognormal model seems to be truncated at
+ both ends, and this may be the main reason why its result differ
+ from lognormal models fitted in Rank--Abundance diagrams with
+ functions \code{\link{rad.lognormal}}.
}
-\value{
- The function \code{prestonfit} returns an object with fitted
+
+\value{ The function \code{prestonfit} returns an object with fitted
\code{coefficients}, and with observed (\code{freq}) and fitted
(\code{fitted}) frequencies, and a string describing the fitting
- \code{method}. Function \code{prestondistr} omits the entry \code{fitted}.
- The function \code{fisherfit} returns the result of \code{\link{nlm}}, where item
- \code{estimate} is \eqn{\alpha}. The result object is amended with the
- following items:
+ \code{method}. Function \code{prestondistr} omits the entry
+ \code{fitted}. The function \code{fisherfit} returns the result of
+ \code{\link{nlm}}, where item \code{estimate} is \eqn{\alpha}. The
+ result object is amended with the following items:
\item{df.residuals}{Residual degrees of freedom.}
- \item{nuisance}{Parameter \eqn{x}.}
- \item{fisher}{Observed data from \code{as.fisher}.}
+ \item{nuisance}{Parameter \eqn{x}.} \item{fisher}{Observed data
+ from \code{as.fisher}.}
}
\references{
@@ -138,8 +164,8 @@
12: 42--58.
Kempton, R.A. & Taylor, L.R. (1974). Log-series and log-normal
- parameters as diversity discriminators for Lepidoptera. \emph{Journal
- of Animal Ecology} 43: 381--399.
+ parameters as diversity discriminators for
+ Lepidoptera. \emph{Journal of Animal Ecology} 43: 381--399.
Preston, F.W. (1948) The commonness and rarity of
species. \emph{Ecology} 29, 254--283.
@@ -148,21 +174,9 @@
not an appropriate null hypothesis for the species--abundance
distribution. \emph{Journal of Animal Ecology} 74, 409--422.
}
-\author{Bob O'Hara \email{bob.ohara at helsinki.fi} (\code{fisherfit}) and Jari Oksanen. }
-\note{It seems that Preston regarded frequencies 1, 2, 4, \emph{etc.}. as ``tied''
- between octaves. This means that only half of the species with
- frequency 1 were shown in the lowest octave, and the rest were
- transferred to the second octave. Half of the species from the second
- octave were transferred to the higher one as well, but this is usually
- not as large number of species. This practise makes data look more
- lognormal by reducing the usually high lowest octaves, but is too
- unfair to be followed. Therefore the octaves used in this function
- include the upper limit. If you do not accept this, you must change
- the function yourself. Williamson & Gaston (2005) discuss alternative
- definitions in detail, and they should be consulted for a
- critical review of log-Normal model.
-}
+\author{Bob O'Hara \email{bob.ohara at helsinki.fi} (\code{fisherfit})
+and Jari Oksanen. }
\seealso{\code{\link{diversity}}, \code{\link{fisher.alpha}},
\code{\link{radfit}}, \code{\link{specpool}}. Function
More information about the Vegan-commits
mailing list