[Vegan-commits] r2631 - in pkg/vegan: R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Oct 3 09:51:01 CEST 2013
Author: jarioksa
Date: 2013-10-03 09:51:01 +0200 (Thu, 03 Oct 2013)
New Revision: 2631
Removed:
pkg/vegan/R/confint.fisherfit.R
pkg/vegan/R/plot.profile.fisherfit.R
pkg/vegan/R/profile.fisherfit.R
Modified:
pkg/vegan/R/fisher.alpha.R
pkg/vegan/R/fisherfit.R
pkg/vegan/R/print.fisherfit.R
pkg/vegan/inst/ChangeLog
pkg/vegan/man/diversity.Rd
pkg/vegan/man/fisherfit.Rd
Log:
re-write fisherfit and remove its standard errors
Deleted: pkg/vegan/R/confint.fisherfit.R
===================================================================
--- pkg/vegan/R/confint.fisherfit.R 2013-09-30 13:33:14 UTC (rev 2630)
+++ pkg/vegan/R/confint.fisherfit.R 2013-10-03 07:51:01 UTC (rev 2631)
@@ -1,6 +0,0 @@
-"confint.fisherfit" <-
- function (object, parm, level=0.95, ...)
-{
- if (!require(MASS)) stop("Needs packages MASS .. not found")
- confint(profile(object), level=level, ...)
-}
Modified: pkg/vegan/R/fisher.alpha.R
===================================================================
--- pkg/vegan/R/fisher.alpha.R 2013-09-30 13:33:14 UTC (rev 2630)
+++ pkg/vegan/R/fisher.alpha.R 2013-10-03 07:51:01 UTC (rev 2631)
@@ -1,12 +1,12 @@
-"fisher.alpha" <-
- function (x, MARGIN = 1, se = FALSE, ...)
+`fisher.alpha` <-
+ function (x, MARGIN = 1, ...)
{
x <- as.matrix(x)
if(ncol(x) == 1)
x <- t(x)
sol <- apply(x, MARGIN, fisherfit)
out <- unlist(lapply(sol, function(x) x$estimate))
- if (se) {
+ if (FALSE) {
out <- list(alpha = out)
out$se <- unlist(lapply(sol, function(x) sqrt(diag(solve(x$hessian)))[1]))
out$df.residual <- unlist(lapply(sol, df.residual))
Modified: pkg/vegan/R/fisherfit.R
===================================================================
--- pkg/vegan/R/fisherfit.R 2013-09-30 13:33:14 UTC (rev 2630)
+++ pkg/vegan/R/fisherfit.R 2013-10-03 07:51:01 UTC (rev 2631)
@@ -1,24 +1,49 @@
+## Fisher alpha is actually based only on the number of species S and
+## number of individuals.
+
`fisherfit` <-
- function (x, ...)
+ function(x, ...)
{
- Dev.logseries <- function(n.r, p, N) {
- r <- as.numeric(names(n.r))
- x <- N/(N + p)
- logmu <- log(p) + log(x) * r - log(r)
- lhood <- -sum(n.r * (logmu - log(n.r)) + 1) - p * log(1-x)
- attr(lhood, "gradient") <- -sum(n.r)/p - log(1-x)
- lhood
+ nr <- as.fisher(x)
+ S <- sum(nr)
+ N <- sum(x)
+ ## Solve 'x' (Fisher alpha).
+ d1fun <- function(x, S, N) x * log(1 + N/x) - S
+ ## We may need to bracket the interval
+ hi <- 50
+ lo <- 1
+ tries <- 0
+ repeat {
+ sol <- try(uniroot(d1fun, c(lo, hi), S = S, N = N, ...), silent = TRUE)
+ if (inherits(sol, "try-error")) {
+ if(d1fun(hi, S, N) < 0)
+ hi <- 2*hi
+ if(d1fun(lo, S, N) > 0)
+ lo <- lo/2
+ tries <- tries + 1
+ }
+ else break
+ ## alpha can tend to +Inf: set root = NA etc.
+ if (tries > 200) {
+ sol <- list(root = NA, f.root = NA, iter = NA, init.it = NA,
+ estim.prec = NA)
+ break
+ }
}
- tmp <- as.rad(x)
- N <- sum(x)
- tmp <- tmp/N
- p <- 1/sum(tmp^2)
- n.r <- as.fisher(x)
- LSeries <- nlm(Dev.logseries, n.r = n.r, p = p, N = N,
- hessian = TRUE, ...)
- LSeries$df.residual <- sum(x > 0) - 1
- LSeries$nuisance <- N/(N + LSeries$estimate)
- LSeries$fisher <- n.r
- class(LSeries) <- "fisherfit"
- LSeries
+ ## 'extendInt' arg was added in R r63162 | maechler | 2013-07-03
+ ## 11:47:22 +0300 (Wed, 03 Jul 2013). Latest release is R 3.0.2 of
+ ## 2013-09-25, but it still does not have the argument. In the
+ ## future we may switch to the following:
+
+ ##sol <- uniroot(d1fun, c(1,50), extendInt = "yes", S = S, N = N, ...)
+
+ nuisance <- N/(N + sol$root)
+ ## we used nlm() earlier, and the following output is compatible
+ out <- list(estimate = sol$root, hessian = NA,
+ iterations = sol$iter, df.residual = NA,
+ nuisance = nuisance, fisher = nr,
+ estim.prec = sol$estim.prec,
+ code = 2*is.na(sol$estim.prec) + 1)
+ class(out) <- "fisherfit"
+ out
}
Deleted: pkg/vegan/R/plot.profile.fisherfit.R
===================================================================
--- pkg/vegan/R/plot.profile.fisherfit.R 2013-09-30 13:33:14 UTC (rev 2630)
+++ pkg/vegan/R/plot.profile.fisherfit.R 2013-10-03 07:51:01 UTC (rev 2631)
@@ -1,16 +0,0 @@
-`plot.profile.fisherfit` <-
- function (x, type = "l", ...)
-{
- tmp <- attr(x, "original.fit")
- est <- tmp$coefficients
- se <- tmp$std.err
- alpha <- x$alpha[, 1]
- tau <- x$alpha[, 2]
- sp <- spline(tau, alpha)
- plot(sp$x, sp$y, type = type, xlab = "alpha", ylab = "tau",
- ...)
- abline(-est/se, 1/se, lty = 2)
- abline(v = est, lty = 3)
- abline(h = 0, lty = 3)
- invisible()
-}
Modified: pkg/vegan/R/print.fisherfit.R
===================================================================
--- pkg/vegan/R/print.fisherfit.R 2013-09-30 13:33:14 UTC (rev 2630)
+++ pkg/vegan/R/print.fisherfit.R 2013-10-03 07:51:01 UTC (rev 2631)
@@ -1,12 +1,8 @@
-"print.fisherfit" <-
+`print.fisherfit` <-
function (x, ...)
{
cat("\nFisher log series model\n")
- cat("No. of species:", sum(x$fisher), "\n\n")
- out <- cbind(x$estimate, sqrt(diag(solve(x$hessian))))
- colnames(out) <- c("Estimate", "Std. Error")
- rownames(out) <- "alpha"
- printCoefmat(out)
- cat("\n")
+ cat("No. of species:", sum(x$fisher), "\n")
+ cat("Fisher alpha: ", x$estimate, "\n\n")
invisible(x)
}
Deleted: pkg/vegan/R/profile.fisherfit.R
===================================================================
--- pkg/vegan/R/profile.fisherfit.R 2013-09-30 13:33:14 UTC (rev 2630)
+++ pkg/vegan/R/profile.fisherfit.R 2013-10-03 07:51:01 UTC (rev 2631)
@@ -1,43 +0,0 @@
-"profile.fisherfit" <-
- function (fitted, alpha = 0.01, maxsteps = 20, del = zmax/5, ...)
-{
- Dev.logseries <- function(n.r, p, N) {
- r <- as.numeric(names(n.r))
- x <- N/(N + p)
- logmu <- log(p) + log(x) * r - log(r)
- lhood <- -sum(n.r * (logmu - log(n.r)) + 1) - p * log(1 -
- x)
- lhood
- }
- par <- fitted$estimate
- names(par) <- "alpha"
- std.err <- sqrt(diag(solve(fitted$hessian)))
- minll <- fitted$minimum
- nr <- fitted$fisher
- N <- sum(as.numeric(names(nr)) * nr)
- zmax <- sqrt(qchisq(1 - alpha/2, 1))
- zi <- 0
- bi <- par
- for (sgn in c(-1, 1)) {
- step <- 0
- z <- 0
- b <- 0
- while ((step <- step + 1) < maxsteps && abs(z) < zmax) {
- b <- par + sgn * step * del * std.err
- fm <- Dev.logseries(nr, b, N)
- zz <- 2 * (fm - minll)
- if (zz > -0.001)
- zz <- max(zz, 0)
- else stop("profiling has found a better solution, so original fit had not converged")
- z <- sgn * sqrt(zz)
- bi <- c(bi, b)
- zi <- c(zi, z)
- }
- }
- si <- order(bi)
- out <- list()
- out$alpha <- data.frame(tau = zi[si], par.vals = bi[si])
- attr(out, "original.fit") <- list(coefficients = par, std.err = std.err)
- class(out) <- c("profile.fisherfit", "profile.glm", "profile")
- out
-}
Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog 2013-09-30 13:33:14 UTC (rev 2630)
+++ pkg/vegan/inst/ChangeLog 2013-10-03 07:51:01 UTC (rev 2631)
@@ -21,7 +21,18 @@
treated differently in factorfit: now they support null
hypothesis, previously they decreased the P-values.
- * fisherfit: use analytic derivatives in non-linear minimizer nlm().
+ * fisherfit: completely rewritten and estimates of standard error
+ removed: I could not find no justification for these. Actually, it
+ seems that the value of Fisher alpha as estimated in the function
+ was independent of the abundance distribution of species, but will
+ be defined by the number of species (S) and number of individuals
+ (N). Now the Fisher alpha is estimated from the relationship S =
+ alpha*(1 + log(N/alpha)) using function uniroot(). Because of
+ this, standard errors cannot be estimated and they were
+ removed. In addition, functions confint.fisherfit,
+ profile.fisherfit and plot.profile.fisherfit were removed. The
+ estimation of standard errors was also removed in function
+ fisher.alpha (that only calls fisherfit).
* nestednodf: matrix fill was wrongly calculated in weighted
analysis. The nominator was length of 'comm', and if input was a
Modified: pkg/vegan/man/diversity.Rd
===================================================================
--- pkg/vegan/man/diversity.Rd 2013-09-30 13:33:14 UTC (rev 2630)
+++ pkg/vegan/man/diversity.Rd 2013-10-03 07:51:01 UTC (rev 2631)
@@ -20,7 +20,7 @@
drarefy(x, sample)
rarecurve(x, step = 1, sample, xlab = "Sample Size", ylab = "Species",
label = TRUE, ...)
-fisher.alpha(x, MARGIN = 1, se = FALSE, ...)
+fisher.alpha(x, MARGIN = 1, ...)
specnumber(x, groups, MARGIN = 1)
}
@@ -32,7 +32,6 @@
\item{base}{ The logarithm \code{base} used in \code{shannon}.}
\item{sample}{Subsample size for rarefying community, either a single
value or a vector.}
- \item{se}{Estimate standard errors.}
\item{step}{Step size for sample sizes in rarefaction curves.}
\item{xlab, ylab}{Axis labels in plots of rarefaction curves.}
\item{label}{Label rarefaction curves by rownames of \code{x} (logical).}
@@ -86,11 +85,7 @@
\code{fisher.alpha} estimates the \eqn{\alpha} parameter of
Fisher's logarithmic series (see \code{\link{fisherfit}}).
The estimation is possible only for genuine
- counts of individuals. The function can optionally return standard
- errors of \eqn{\alpha}. These should be regarded only as rough
- indicators of the accuracy: the confidence limits of \eqn{\alpha} are
- strongly non-symmetric and the standard errors cannot be used in
- Normal inference.
+ counts of individuals.
Function \code{specnumber} finds the number of species. With
\code{MARGIN = 2}, it finds frequencies of species. If \code{groups}
Modified: pkg/vegan/man/fisherfit.Rd
===================================================================
--- pkg/vegan/man/fisherfit.Rd 2013-09-30 13:33:14 UTC (rev 2630)
+++ pkg/vegan/man/fisherfit.Rd 2013-10-03 07:51:01 UTC (rev 2631)
@@ -2,9 +2,6 @@
\alias{fisherfit}
\alias{as.fisher}
\alias{plot.fisherfit}
-\alias{profile.fisherfit}
-\alias{confint.fisherfit}
-\alias{plot.profile.fisherfit}
\alias{prestonfit}
\alias{prestondistr}
\alias{as.preston}
@@ -25,9 +22,6 @@
}
\usage{
fisherfit(x, ...)
-\method{confint}{fisherfit}(object, parm, level = 0.95, ...)
-\method{profile}{fisherfit}(fitted, alpha = 0.01, maxsteps = 20, del = zmax/5,
- ...)
prestonfit(x, tiesplit = TRUE, ...)
prestondistr(x, truncate = -1, ...)
\method{plot}{prestonfit}(x, xlab = "Frequency", ylab = "Species", bar.col = "skyblue",
@@ -45,12 +39,6 @@
\arguments{
\item{x}{Community data vector for fitting functions or their result
object for \code{plot} functions.}
- \item{object, fitted}{Fitted model.}
- \item{parm}{Not used.}
- \item{level}{The confidence level required.}
- \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
@@ -73,9 +61,8 @@
\details{
In Fisher's logarithmic series the expected
number of species \eqn{f} with \eqn{n} observed individuals is
- \eqn{f_n = \alpha x^n / n} (Fisher et al. 1943). The estimation
- follows Kempton & Taylor (1974) and uses function
- \code{\link{nlm}}. The estimation is possible only for genuine
+ \eqn{f_n = \alpha x^n / n} (Fisher et al. 1943).
+ The estimation is possible only for genuine
counts of individuals. The parameter \eqn{\alpha} is used as a
diversity index, and \eqn{\alpha} and its standard error can be
estimated with a separate function \code{\link{fisher.alpha}}. The
@@ -84,20 +71,6 @@
function \code{as.fisher} transforms abundance data into Fisher
frequency table.
- Function \code{fisherfit} estimates the standard error of
- \eqn{\alpha}{alpha}. However, the confidence limits cannot be directly
- estimated from the standard errors, but you should use function
- \code{confint} based on profile likelihood. Function \code{confint}
- uses function \code{\link[MASS]{confint.glm}} of the \pkg{MASS}
- package, using \code{profile.fisherfit} for the profile
- likelihood. Function \code{profile.fisherfit} follows
- \code{\link[MASS]{profile.glm}} and finds the \eqn{\tau}{tau} parameter or
- signed square root of two times log-Likelihood profile. The profile can
- be inspected with a \code{plot} function which shows the \eqn{\tau}{tau}
- and a dotted line corresponding to the Normal assumption: if standard
- errors can be directly used in Normal inference these two lines
- 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
@@ -162,11 +135,8 @@
\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}.}
-
+ result object is amended with the \code{nuisance} parameter and item
+ \code{fisher} for the observed data from \code{as.fisher}
}
\references{
Fisher, R.A., Corbet, A.S. & Williams, C.B. (1943). The relation
@@ -174,10 +144,6 @@
random sample of animal population. \emph{Journal of Animal Ecology}
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.
-
Preston, F.W. (1948) The commonness and rarity of
species. \emph{Ecology} 29, 254--283.
@@ -186,7 +152,7 @@
distribution. \emph{Journal of Animal Ecology} 74, 409--422.
}
-\author{Bob O'Hara (\code{fisherfit}) and Jari Oksanen. }
+\author{Bob O'Hara and Jari Oksanen. }
\seealso{\code{\link{diversity}}, \code{\link{fisher.alpha}},
\code{\link{radfit}}, \code{\link{specpool}}. Function
@@ -200,8 +166,6 @@
data(BCI)
mod <- fisherfit(BCI[5,])
mod
-plot(profile(mod))
-confint(mod)
# prestonfit seems to need large samples
mod.oct <- prestonfit(colSums(BCI))
mod.ll <- prestondistr(colSums(BCI))
More information about the Vegan-commits
mailing list