[Vegan-commits] r2301 - in branches/2.0: . R inst man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Thu Sep 27 18:37:21 CEST 2012


Author: jarioksa
Date: 2012-09-27 18:37:20 +0200 (Thu, 27 Sep 2012)
New Revision: 2301

Added:
   branches/2.0/R/predict.radline.R
Modified:
   branches/2.0/NAMESPACE
   branches/2.0/R/AIC.radfit.R
   branches/2.0/R/coef.radfit.R
   branches/2.0/R/fitted.radfit.R
   branches/2.0/R/lines.radline.R
   branches/2.0/R/plot.radfit.R
   branches/2.0/R/points.radline.R
   branches/2.0/R/rad.null.R
   branches/2.0/R/rad.preempt.R
   branches/2.0/R/radfit.R
   branches/2.0/R/radfit.data.frame.R
   branches/2.0/R/radlattice.R
   branches/2.0/inst/ChangeLog
   branches/2.0/man/radfit.Rd
Log:
merge r2885,2291:96,98,2300 -- add predict.radfit & upgrade radfit functions

Modified: branches/2.0/NAMESPACE
===================================================================
--- branches/2.0/NAMESPACE	2012-09-27 10:53:27 UTC (rev 2300)
+++ branches/2.0/NAMESPACE	2012-09-27 16:37:20 UTC (rev 2301)
@@ -73,6 +73,7 @@
 S3method(adipart, formula)
 # AIC: stats
 S3method(AIC, radfit)
+S3method(AIC, radfit.frame)
 # RsquareAdj: vegan
 S3method(RsquareAdj, cca)
 S3method(RsquareAdj, default)
@@ -123,6 +124,7 @@
 # coef: stats
 S3method(coef, cca)
 S3method(coef, radfit)
+S3method(coef, radfit.frame)
 S3method(coef, rda)
 # confint: stats -- also uses MASS:::confint.glm & MASS:::profile.glm
 # does this work with namespaces??
@@ -144,6 +146,8 @@
 # deviance: stats
 S3method(deviance, cca)
 S3method(deviance, rda)
+S3method(deviance, radfit)
+S3method(deviance, radfit.frame)
 # drop1: stats
 S3method(drop1, cca)
 # eigenvals: vegan
@@ -170,6 +174,7 @@
 S3method(fitted, cca)
 S3method(fitted, procrustes)
 S3method(fitted, radfit)
+S3method(fitted, radfit.frame)
 S3method(fitted, rda)
 # goodness: vegan
 S3method(goodness, cca)
@@ -192,7 +197,11 @@
 S3method(lines, prestonfit)
 S3method(lines, procrustes)
 S3method(lines, radline)
+S3method(lines, radfit)
 S3method(lines, spantree)
+## logLik: stats
+S3method(logLik, radfit)
+S3method(logLik, radfit.frame)
 # model.frame, model.matrix: stats
 S3method(model.frame, cca)
 S3method(model.matrix, cca)
@@ -273,12 +282,16 @@
 S3method(points, orditkplot)
 S3method(points, procrustes)
 S3method(points, radline)
+S3method(points, radfit)
 # predict: stats
 S3method(predict, cca)
 S3method(predict, decorana)
 S3method(predict, fitspecaccum)
 S3method(predict, humpfit)
 S3method(predict, procrustes)
+S3method(predict, radline)
+S3method(predict, radfit)
+S3method(predict, radfit.frame)
 S3method(predict, rda)
 S3method(predict, specaccum)
 # print: base
@@ -346,6 +359,7 @@
 # radfit: vegan
 S3method(radfit, data.frame)
 S3method(radfit, default)
+S3method(radfit, matrix)
 # rda: vegan
 S3method(rda, default)
 S3method(rda, formula)

Modified: branches/2.0/R/AIC.radfit.R
===================================================================
--- branches/2.0/R/AIC.radfit.R	2012-09-27 10:53:27 UTC (rev 2300)
+++ branches/2.0/R/AIC.radfit.R	2012-09-27 16:37:20 UTC (rev 2301)
@@ -1,6 +1,41 @@
-"AIC.radfit" <-
+### these functions are defined _ex machina_ for radline objects which
+### inherit from glm. Here we define them for radfit objects where
+### object$models is a list of radline objects
+
+`AIC.radfit` <-
     function (object, k = 2, ...) 
 {
-    mods <- object$models
-    unlist(lapply(mods, AIC, k = k))
+    sapply(object$models, AIC, k = k, ...)
 }
+
+`deviance.radfit` <-
+    function(object, ...)
+{
+    sapply(object$models, deviance, ...)
+}
+
+`logLik.radfit` <-
+    function(object, ...)
+{
+    sapply(object$models, logLik, ...)
+}
+
+### Define also for radfit.frames which are lists of radfit objects
+
+`AIC.radfit.frame` <-
+    function(object, k = 2, ...)
+{
+    sapply(object, AIC, k = k, ...)
+}
+
+`deviance.radfit.frame` <-
+    function(object, ...)
+{
+    sapply(object, deviance, ...)
+}
+
+`logLik.radfit.frame` <-
+    function(object, ...)
+{
+    sapply(object, logLik, ...)
+}

Modified: branches/2.0/R/coef.radfit.R
===================================================================
--- branches/2.0/R/coef.radfit.R	2012-09-27 10:53:27 UTC (rev 2300)
+++ branches/2.0/R/coef.radfit.R	2012-09-27 16:37:20 UTC (rev 2301)
@@ -1,4 +1,4 @@
-"coef.radfit" <-
+`coef.radfit` <-
     function (object, ...) 
 {
     out <- sapply(object$models, function(x) if (length(coef(x)) < 
@@ -9,3 +9,9 @@
     colnames(out) <- paste("par", 1:3, sep = "")
     out
 }
+
+`coef.radfit.frame` <-
+    function(object, ...)
+{
+    lapply(object, coef, ...)
+}

Modified: branches/2.0/R/fitted.radfit.R
===================================================================
--- branches/2.0/R/fitted.radfit.R	2012-09-27 10:53:27 UTC (rev 2300)
+++ branches/2.0/R/fitted.radfit.R	2012-09-27 16:37:20 UTC (rev 2301)
@@ -1,5 +1,11 @@
 `fitted.radfit` <-
     function(object, ...)
 {
-    matrix(sapply(object$models, fitted), ncol=length(object$models))
+    sapply(object$models, fitted)
 }
+
+`fitted.radfit.frame` <-
+    function(object, ...)
+{
+    lapply(object, fitted, ...)
+}

Modified: branches/2.0/R/lines.radline.R
===================================================================
--- branches/2.0/R/lines.radline.R	2012-09-27 10:53:27 UTC (rev 2300)
+++ branches/2.0/R/lines.radline.R	2012-09-27 16:37:20 UTC (rev 2301)
@@ -1,4 +1,4 @@
-"lines.radline" <-
+`lines.radline` <-
     function (x, ...) 
 {
     lin <- fitted(x)
@@ -6,3 +6,9 @@
     lines(rnk, lin, ...)
     invisible()
 }
+
+`lines.radfit` <-
+    function(x, ...)
+{
+    matlines(fitted(x), ...)
+}

Modified: branches/2.0/R/plot.radfit.R
===================================================================
--- branches/2.0/R/plot.radfit.R	2012-09-27 10:53:27 UTC (rev 2300)
+++ branches/2.0/R/plot.radfit.R	2012-09-27 16:37:20 UTC (rev 2301)
@@ -1,8 +1,13 @@
-"plot.radfit" <-
+`plot.radfit` <-
     function (x, BIC = FALSE, legend = TRUE, ...) 
 {
     if (length(x$y) == 0)
         stop("No species, nothing to plot")
+    ## if 'type = "n"', do not add legend (other types are not
+    ## supported)
+    type <- match.call(expand.dots = FALSE)$...$type
+    if (is.null(type))
+        type <- ""
     out <- plot(x$y, ...)
     if (length(x$y) == 1)
         return(invisible(out))
@@ -14,7 +19,7 @@
     lwd <- rep(1, ncol(fv))
     lwd[emph] <- 3
     matlines(fv, lty = 1, lwd = lwd, ...)
-    if (legend) {
+    if (legend && type != "n") {
         nm <- names(x$models)
         legend("topright", legend = nm, lty = 1, lwd = lwd, col = 1:6)
     }

Modified: branches/2.0/R/points.radline.R
===================================================================
--- branches/2.0/R/points.radline.R	2012-09-27 10:53:27 UTC (rev 2300)
+++ branches/2.0/R/points.radline.R	2012-09-27 16:37:20 UTC (rev 2301)
@@ -1,4 +1,4 @@
-"points.radline" <-
+`points.radline` <-
     function (x, ...) 
 {
     poi <- x$y
@@ -8,3 +8,9 @@
     class(out) <- "ordiplot"
     invisible(out)
 }
+
+`points.radfit` <-
+    function(x, ...)
+{
+    points.radline(x, ...)
+}

Copied: branches/2.0/R/predict.radline.R (from rev 2285, pkg/vegan/R/predict.radline.R)
===================================================================
--- branches/2.0/R/predict.radline.R	                        (rev 0)
+++ branches/2.0/R/predict.radline.R	2012-09-27 16:37:20 UTC (rev 2301)
@@ -0,0 +1,50 @@
+### predict method for radline, radfit & radfit.frame
+
+### All functions take 'newdata' argument which need not be integer:
+### the functions can interpolate, but not necessarily extrapolate, or
+### the extrapolations may be NaN.
+
+`predict.radline`  <-
+    function(object, newdata, total, ...)
+{
+    ## newdata can be ranks
+    if (missing(newdata))
+        x <- seq_along(object$y)
+    else
+        x <- drop(as.matrix(newdata))
+    ## total number of individuals in the community
+    if (missing(total))
+        total <- sum(object$y)
+    ## adjustment for chagned total in call
+    adj <- total/sum(object$y)
+    nobs <- length(object$y)
+    p <- coef(object)
+    switch(object$model,
+           ## linear interpolation, no extrapolation
+           `Brokenstick` = approx(seq_len(nobs),
+           object$fitted.values, x, ...)$y * adj,
+           `Preemption` = exp(log(total) + log(p) + log(1 - p)*(x-1)),
+           ## NaN when rank outside proportional rank 0...1 
+           `Log-Normal` = {
+               slope <- diff(range(ppoints(nobs)))/(nobs-1)
+               intcpt <- 0.5 - slope * (nobs + 1) / 2
+               xnorm <- -qnorm(intcpt + slope * x)
+               exp(p[1] + p[2]*xnorm)*adj
+           },
+           `Zipf` = exp(log(total) + log(p[1]) + p[2]*log(x)),
+           `Zipf-Mandelbrot` = exp(log(total) + log(p[1]) +
+           p[2]*log(x + p[3]))
+           )
+}
+
+`predict.radfit`<-
+    function(object, ...)
+{
+    sapply(object$models, predict, ...)
+}
+
+`predict.radfit.frame`  <-
+    function(object, ...)
+{
+    lapply(object, predict, ...)
+}

Modified: branches/2.0/R/rad.null.R
===================================================================
--- branches/2.0/R/rad.null.R	2012-09-27 10:53:27 UTC (rev 2300)
+++ branches/2.0/R/rad.null.R	2012-09-27 16:37:20 UTC (rev 2301)
@@ -21,6 +21,7 @@
     }
     residuals <- x - fit
     rdf <- nsp
+    names(fit) <- names(x)
     p <- NA
     names(p) <- "S"
     out <- list(model = "Brokenstick", family=fam, y = x, coefficients = p,

Modified: branches/2.0/R/rad.preempt.R
===================================================================
--- branches/2.0/R/rad.preempt.R	2012-09-27 10:53:27 UTC (rev 2300)
+++ branches/2.0/R/rad.preempt.R	2012-09-27 16:37:20 UTC (rev 2301)
@@ -45,6 +45,7 @@
             aic <- NA
         rdf <- length(x) - 1
     }
+    names(fit) <- names(x)
     names(p) <- c("alpha")
     out <- list(model = "Preemption", family = fam, y = x, coefficients = p, 
                 fitted.values = fit, aic = aic, rank = 1, df.residual = rdf, 

Modified: branches/2.0/R/radfit.R
===================================================================
--- branches/2.0/R/radfit.R	2012-09-27 10:53:27 UTC (rev 2300)
+++ branches/2.0/R/radfit.R	2012-09-27 16:37:20 UTC (rev 2301)
@@ -1,5 +1,5 @@
-"radfit" <-
-    function (...)
+`radfit` <-
+    function (x, ...)
 {
     UseMethod("radfit")
 }

Modified: branches/2.0/R/radfit.data.frame.R
===================================================================
--- branches/2.0/R/radfit.data.frame.R	2012-09-27 10:53:27 UTC (rev 2300)
+++ branches/2.0/R/radfit.data.frame.R	2012-09-27 16:37:20 UTC (rev 2301)
@@ -1,9 +1,9 @@
-"radfit.data.frame" <-
-    function(df, ...)
+`radfit.data.frame` <-
+    function(x, ...)
 {
-    ## df *must* have rownames
-    rownames(df) <- rownames(df, do.NULL = TRUE)
-    out <- apply(df, 1, radfit, ...)
+    ## x *must* have rownames
+    rownames(x) <- rownames(x, do.NULL = TRUE)
+    out <- apply(x, 1, radfit, ...)
     if (length(out) == 1)
         out <- out[[1]]
     else {
@@ -12,3 +12,9 @@
     }
     out
 }
+
+`radfit.matrix` <-
+    function(x, ...)
+{
+    radfit(as.data.frame(x), ...)
+}

Modified: branches/2.0/R/radlattice.R
===================================================================
--- branches/2.0/R/radlattice.R	2012-09-27 10:53:27 UTC (rev 2300)
+++ branches/2.0/R/radlattice.R	2012-09-27 16:37:20 UTC (rev 2301)
@@ -1,6 +1,8 @@
 `radlattice` <-
     function(x, BIC = FALSE, ...)
 {
+    if (!inherits(x, "radfit"))
+        stop("function only works with 'radfit' results for single site")
     require(lattice) || stop("requires package 'lattice'")
     y <- x$y
     fv <- unlist(fitted(x))
@@ -10,7 +12,11 @@
     Abundance <- rep(y, p)
     Rank <- rep(1:n, p)
     Model <- factor(rep(mods, each=n), levels = mods)
-    aic <- AIC(x, BIC = BIC)
+    if (BIC)
+        k <- log(length(y))
+    else
+        k <- 2
+    aic <- AIC(x, k = k)
     col <- trellis.par.get("superpose.line")$col
     if (length(col) > 1)
         col <- col[2]

Modified: branches/2.0/inst/ChangeLog
===================================================================
--- branches/2.0/inst/ChangeLog	2012-09-27 10:53:27 UTC (rev 2300)
+++ branches/2.0/inst/ChangeLog	2012-09-27 16:37:20 UTC (rev 2301)
@@ -3,8 +3,10 @@
 VEGAN RELEASE VERSIONS at http://cran.r-project.org/
 
 Version 2.0-5 (opened June 18, 2012)
-
-	* merge r2288,9: scoping in anova.ccabyaxis and anova.ccabyterm.
+	
+	* merge r2291 thru 2296, 2298, 2300: radfit upgrade.
+	* merge r2287,8: scoping in anova.ccabyaxis and anova.ccabyterm.
+	* merge r2285: add predict.radfit.
 	* merge r2262, 2268:2270, and also r1950: mantel and
 	mantel.partial gained argument na.rm = FALSE. This needed hand
 	editing of merges, and also merging old r1950: beware and test.

Modified: branches/2.0/man/radfit.Rd
===================================================================
--- branches/2.0/man/radfit.Rd	2012-09-27 10:53:27 UTC (rev 2300)
+++ branches/2.0/man/radfit.Rd	2012-09-27 16:37:20 UTC (rev 2301)
@@ -3,22 +3,34 @@
 \alias{radfit.default}
 \alias{radfit.data.frame}
 \alias{AIC.radfit}
+\alias{AIC.radfit.frame}
 \alias{as.rad}
 \alias{coef.radfit}
+\alias{coef.radfit.frame}
+\alias{deviance.radfit}
+\alias{deviance.radfit.frame}
+\alias{logLik, radfit}
+\alias{logLik, radfit.frame}
 \alias{fitted.radfit}
+\alias{fitted.radfit.frame}
 \alias{lines.radline}
+\alias{lines.radfit}
 \alias{plot.radfit.frame}
 \alias{plot.radfit}
 \alias{plot.radline}
 \alias{plot.rad}
 \alias{radlattice}
 \alias{points.radline}
+\alias{points.radfit}
 \alias{summary.radfit.frame}
 \alias{rad.preempt}
 \alias{rad.lognormal}
 \alias{rad.zipf}
 \alias{rad.zipfbrot}
 \alias{rad.null}
+\alias{predict.radline}
+\alias{predict.radfit}
+\alias{predict.radfit.frame}
 
 \title{ Rank -- Abundance or Dominance / Diversity Models}
 \description{
@@ -27,152 +39,167 @@
   Zipf and Zipf-Mandelbrot models of species abundance.
 }
 \usage{
-\method{radfit}{data.frame}(df, ...)
-\method{plot}{radfit.frame}(x, order.by, BIC = FALSE, model, legend = TRUE,
-     as.table = TRUE, ...)
 \method{radfit}{default}(x, ...)
-\method{plot}{radfit}(x, BIC = FALSE, legend = TRUE, ...)  
-radlattice(x, BIC = FALSE, ...)
 rad.null(x, family=poisson, ...)
 rad.preempt(x, family = poisson, ...)
 rad.lognormal(x, family = poisson, ...)
 rad.zipf(x, family = poisson, ...)
 rad.zipfbrot(x, family = poisson, ...)
+\method{predict}{radline}(object, newdata, total, ...)
+\method{plot}{radfit}(x, BIC = FALSE, legend = TRUE, ...)  
+\method{plot}{radfit.frame}(x, order.by, BIC = FALSE, model, legend = TRUE,
+     as.table = TRUE, ...)
 \method{plot}{radline}(x, xlab = "Rank", ylab = "Abundance", type = "b", ...)
-\method{lines}{radline}(x, ...)
-\method{points}{radline}(x, ...)
+radlattice(x, BIC = FALSE, ...)
+\method{lines}{radfit}(x, ...)
+\method{points}{radfit}(x, ...)
 as.rad(x)
 \method{plot}{rad}(x, xlab = "Rank", ylab = "Abundance", log = "y", ...)
 }
 
 \arguments{
-  \item{df}{Data frame where sites are rows and species are columns.}
-  \item{x}{A vector giving species abundances in a site, or an object to
+  \item{x}{Data frame, matrix or a vector giving species abundances, or an object to
     be plotted.}
+
+  \item{family}{Error distribution (passed to \code{\link{glm}}). All
+    alternatives accepting \code{link = "log"} in \code{\link{family}}
+    can be used, although not all make sense.}
+
+  \item{object}{A fitted result object.}
+
+  \item{newdata}{Ranks used for ordinations. All models can
+    interpolate to non-integer \dQuote{ranks} (although this may be
+    approximate), but extrapolation may fail}
+
+  \item{total}{The new total used for predicting abundance. Observed
+    total count is used if this is omitted.}
+
   \item{order.by}{A vector used for ordering sites in plots.}
   \item{BIC}{Use Bayesian Information Criterion, BIC, instead of
-    Akaike's AIC. The penalty for a parameter is \eqn{k = \log(S)}{k =
+    Akaike's AIC. The penalty in BIC is \eqn{k = \log(S)}{k =
       log(S)}  where \eqn{S} is the number of species, whereas AIC uses
     \eqn{k = 2}.} 
-  \item{model}{Show only the specified model. If missing, AIC is used to
-    select the model. The model names (which can be abbreviated) are
-    \code{Preemption}, \code{Lognormal}, \code{Veiled.LN},
-    \code{Zipf}, \code{Mandelbrot}. }
+  \item{model}{Show only the specified model. If missing, AIC is used
+    to select the model. The model names (which can be abbreviated)
+    are \code{Null}, \code{Preemption}, \code{Lognormal}, \code{Zipf},
+    \code{Mandelbrot}. }
   \item{legend}{Add legend of line colours.}
   \item{as.table}{Arrange panels starting from upper left corner (passed
     to \code{\link[lattice]{xyplot}}).}
-  \item{family}{Error distribution (passed to \code{\link{glm}}). All
-    alternatives accepting \code{link = "log"} in \code{\link{family}}
-    can be used, although not all make sense.}
   \item{xlab,ylab}{Labels for \code{x} and \code{y} axes.}
   \item{type}{Type of the plot, \code{"b"} for plotting both observed points
     and fitted lines, \code{"p"} for only points, \code{"l"} for only
     fitted lines, and \code{"n"} for only setting the frame. }
   \item{log}{Use logarithmic scale for given axis. The default
-    \code{log =" y"} gives the traditional plot in community ecology
+    \code{log = "y"} gives the traditional plot of community ecology
     where the pre-emption model is a straight line, and with
     \code{log = "xy"} Zipf model is a straight line. With
     \code{log = ""} both axes are in the original arithmetic scale.}
   \item{\dots}{Other parameters to functions. }
 }
 \details{
-  Rank -- Abundance Dominance (RAD) or Dominance/Diversity plots
+
+  Rank--Abundance Dominance (RAD) or Dominance/Diversity plots
   (Whittaker 1965) display logarithmic species abundances against
-  species rank order. These plots are supposed to be
-  effective in analysing types of abundance distributions in
-  communities. These functions fit some of the most popular models mainly
-  following Wilson (1991). Function \code{as.rad} constructs observed
-  RAD data.
-  Functions \code{rad.XXXX} (where \code{XXXX} is a name) fit
-  the individual models, and
-  function \code{radfit} fits all models.  The
-  argument of the function \code{radfit} can be either a vector for a
-  single community or a data frame where each row represents a
-  distinct community.  All these functions have their own \code{plot}
-  functions. When the argument is a data frame, \code{plot} uses
-  \code{\link[lattice]{Lattice}} graphics, and other \code{plot} functions use
-  ordinary graphics. The ordinary graphics functions return invisibly an
-  \code{\link{ordiplot}} object for observed points, and function
-  \code{\link{identify.ordiplot}} can be used to label selected
-  species. The most complete control of graphics can be achieved
-  with \code{rad.XXXX} methods which have \code{points} and \code{lines}
-  functions to add observed values and fitted models into existing
-  graphs.  Alternatively, \code{radlattice} uses
-  \code{\link[lattice]{Lattice}} graphics to display each
-  \code{radfit} model in a separate panel together with their AIC or
-  BIC values.
+  species rank order. These plots are supposed to be effective in
+  analysing types of abundance distributions in communities. These
+  functions fit some of the most popular models mainly following
+  Wilson (1991).
 
+  Functions \code{rad.null}, \code{rad.preempt}, \code{rad.lognormal},
+  \code{rad.zipf} and \code{zipfbrot} fit the individual models
+  (described below) for a single vector (row of data frame), and
+  function \code{radfit} fits all models. The argument of the function
+  \code{radfit} can be either a vector for a single community or a data
+  frame where each row represents a distinct community.
+  
   Function \code{rad.null} fits a brokenstick model where the expected
-  abundance of species at rank \eqn{r} is \eqn{a_r = (J/S) \sum_{x=r}^S
-    (1/x)}{a[r] = J/S sum(from x=r to S) 1/x} (Pielou 1975), where \eqn{J}
-  is the total number of individuals (site total) and \eqn{S} is the
-  total number of species in the community.  This gives a Null model
-  where the individuals are randomly distributed among observed species,
-  and there are no fitted parameters. 
+  abundance of species at rank \eqn{r} is \eqn{a_r = (J/S)
+  \sum_{x=r}^S (1/x)}{a[r] = J/S sum(from x=r to S) 1/x} (Pielou
+  1975), where \eqn{J} is the total number of individuals (site total)
+  and \eqn{S} is the total number of species in the community.  This
+  gives a Null model where the individuals are randomly distributed
+  among observed species, and there are no fitted parameters.
   Function \code{rad.preempt} fits the niche preemption model,
   a.k.a. geometric series or Motomura model, where the expected
-  abundance \eqn{a} of species at rank \eqn{r} is \eqn{a_r = J \alpha (1 -
-    \alpha)^{r-1}}{a[r] = J*alpha*(1-alpha)^(r-1)}. The only estimated
-  parameter is the preemption coefficient \eqn{\alpha} which gives the
-  decay rate of abundance per rank.
-  The niche preemption model is a straight line in a
-  RAD plot. Function \code{rad.lognormal} fits a log-Normal model which
-  assumes that the logarithmic abundances are distributed Normally, or
-  \eqn{a_r =  \exp( \log \mu + \log \sigma N)}{a[r] = exp(log(mu) +
-    log(sigma) * N)}, where \eqn{N} is a Normal deviate. 
-  Function \code{rad.zipf} fits the Zipf model \eqn{a_r = J p_1
-    r^\gamma}{a[r] = J*p1*r^gamma} where \eqn{p_1}{p1} is the fitted
-  proportion of the most abundant species, and \eqn{\gamma} is a
-  decay coefficient. The
-  Zipf -- Mandelbrot 
-  model (\code{rad.zipfbrot}) adds one parameter: \eqn{a_r = J c
-    (r + \beta)^\gamma}{a[r] = J*c*(r+beta)^gamma} after which
-  \eqn{p_1}{p1} of the Zipf model changes into a meaningless scaling
-  constant \eqn{c}. There are grand narratives about ecological
-  mechanisms behind each model (Wilson 1991), but
-  several alternative and contrasting mechanisms can produce
-  similar models and a good fit does not imply a specific mechanism.
+  abundance \eqn{a} of species at rank \eqn{r} is \eqn{a_r = J \alpha
+  (1 - \alpha)^{r-1}}{a[r] = J*alpha*(1-alpha)^(r-1)}. The only
+  estimated parameter is the preemption coefficient \eqn{\alpha} which
+  gives the decay rate of abundance per rank.  The niche preemption
+  model is a straight line in a RAD plot.  Function
+  \code{rad.lognormal} fits a log-Normal model which assumes that the
+  logarithmic abundances are distributed Normally, or \eqn{a_r = \exp(
+  \log \mu + \log \sigma N)}{a[r] = exp(log(mu) + log(sigma) * N)},
+  where \eqn{N} is a Normal deviate.  Function \code{rad.zipf} fits
+  the Zipf model \eqn{a_r = J p_1 r^\gamma}{a[r] = J*p1*r^gamma} where
+  \eqn{p_1}{p1} is the fitted proportion of the most abundant species,
+  and \eqn{\gamma} is a decay coefficient. The Zipf--Mandelbrot model
+  (\code{rad.zipfbrot}) adds one parameter: \eqn{a_r = J c (r +
+  \beta)^\gamma}{a[r] = J*c*(r+beta)^gamma} after which \eqn{p_1}{p1}
+  of the Zipf model changes into a meaningless scaling constant
+  \eqn{c}. 
 
-  Log-Normal and Zipf models are generalized linear
-  models (\code{\link{glm}}) with logarithmic link function.
-  Zipf-Mandelbrot adds one
-  nonlinear parameter to the Zipf model, and is fitted using
+  Log-Normal and Zipf models are generalized linear models
+  (\code{\link{glm}}) with logarithmic link function.  Zipf--Mandelbrot
+  adds one nonlinear parameter to the Zipf model, and is fitted using
   \code{\link{nlm}} for the nonlinear parameter and estimating other
-  parameters and log-Likelihood with \code{\link{glm}}. Pre-emption
-  model is fitted as purely nonlinear model. There are no estimated
-  parameters in the Null model.  The default
-  \code{\link{family}} is \code{poisson} which is appropriate only for
-  genuine counts (integers), but other families that accept \code{link =
-    "log"} can be used. Family \code{\link{Gamma}} may be
-  appropriate for abundance data, such as cover. The ``best''
-  model is selected by \code{\link{AIC}}. Therefore ``quasi'' families
-  such as \code{\link{quasipoisson}} cannot be used: they do not
-  have \code{\link{AIC}} nor log-Likelihood needed in non-linear
-  models.
+  parameters and log-Likelihood with \code{\link{glm}}. Preemption
+  model is fitted as a purely nonlinear model. There are no estimated
+  parameters in the Null model.  
+
+  The default \code{\link{family}} is \code{poisson} which is
+  appropriate only for genuine counts (integers), but other families
+  that accept \code{link = "log"} can be used. Families
+  \code{\link{Gamma}} or \code{\link{gaussian}} may be appropriate for
+  abundance data, such as cover. The ``best'' model is selected by
+  \code{\link{AIC}}. Therefore ``quasi'' families such as
+  \code{\link{quasipoisson}} cannot be used: they do not have
+  \code{\link{AIC}} nor log-Likelihood needed in non-linear models.
+  
+  All these functions have their own \code{plot} functions. When
+  \code{radfit} was applied for a data frame, \code{plot} uses
+  \code{\link[lattice]{Lattice}} graphics, and other \code{plot}
+  functions use ordinary graphics. The ordinary graphics functions
+  return invisibly an \code{\link{ordiplot}} object for observed points,
+  and function \code{\link{identify.ordiplot}} can be used to label
+  selected species.  Alternatively, \code{radlattice} uses
+  \code{\link[lattice]{Lattice}} graphics to display each \code{radfit}
+  model of a single site in a separate panel together with their AIC or
+  BIC values.
+
+  Function \code{as.rad} is a base function to construct ordered RAD
+  data. Its \code{plot} is used by other RAD \code{plot} functions
+  which pass extra arguments (such as \code{xlab} and \code{log}) to
+  this function.
+
 }
 
 \value{
-  Function \code{rad.XXXX} will return an object of class
-  \code{radline}, which is constructed to resemble results of \code{\link{glm}}
-  and has many (but not all) of its components, even when only
-  \code{\link{nlm}} was used in fitting. At least the following
-  \code{\link{glm}} methods can be applied to the result:
-  \code{\link{fitted}}, \code{\link{residuals.glm}}  with alternatives
-  \code{"deviance"} (default), \code{"pearson"}, \code{"response"},
-  function \code{\link{coef}}, \code{\link{AIC}},
-  \code{\link{extractAIC}}, and \code{\link{deviance}}.
-  Function \code{radfit} applied to a vector will return
-  an object of class \code{radfit} with item \code{y} for the
-  constructed RAD, item \code{family} for the error distribution, and
-  item \code{models} containing each \code{radline} object as an
-  item. In addition, there are special \code{AIC}, \code{coef} and
-  \code{fitted} implementations for \code{radfit} results. 
-  When applied to a data frame
-  \code{radfit} will return an object of class \code{radfit.frame} which
-  is a list of \code{radfit} objects; function \code{summary} can be
-  used to display the results for individual \code{radfit} objects.
-  The functions are still
-  preliminary, and the items in the \code{radline} objects may change.
+
+  Functions \code{rad.null}, \code{rad.preempt}, \code{rad.lognormal},
+  \code{zipf} and \code{zipfbrot} fit each a single RAD model to a
+  single site. The result object has class \code{"radline"} and
+  inherits from \code{\link{glm}}, and can be handled by some (but not
+  all) \code{\link{glm}} methods.
+
+  Function \code{radfit} fits all models either to a single site or to
+  all rows of a data frame or a matrix. When fitted to a single site,
+  the function returns an object of class \code{"radfit"} with items
+  \code{y} (observed values), \code{\link{family}}, and \code{models}
+  which is a list of fitted \code{"radline"} models.  When applied for a
+  data frame or matrix, \code{radfit} function returns an object of
+  class \code{"radfit.frame"} which is a list of \code{"radfit"}
+  objects, each item names by the corresponding row name.
+
+  All result objects (\code{"radline"}, \code{"radfit"},
+  \code{"radfit.frame"}) can be accessed with same method functions.
+  The following methods are available: \code{\link{AIC}},
+  \code{\link{coef}}, \code{\link{deviance}}, \code{\link{logLik}}. In
+  addition the fit results can be accessed with \code{\link{fitted}},
+  \code{\link{predict}} and \code{\link{residuals}} (inheriting from
+  \code{\link{residuals.glm}}).The graphical functions were discussed
+  above in Details.
+
 }
 
 \references{



More information about the Vegan-commits mailing list