[Vegan-commits] r1549 - in branches/1.17: R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Mar 19 20:31:33 CET 2011
Author: jarioksa
Date: 2011-03-19 20:31:32 +0100 (Sat, 19 Mar 2011)
New Revision: 1549
Added:
branches/1.17/R/plot.ordisurf.R
Modified:
branches/1.17/R/metaMDSiter.R
branches/1.17/R/metaMDSrotate.R
branches/1.17/R/ordisurf.R
branches/1.17/inst/ChangeLog
branches/1.17/man/metaMDS.Rd
branches/1.17/man/ordisurf.Rd
branches/1.17/man/vegdist.Rd
Log:
merge revs 1530 thru 1548:
- ordisurf got formula method, new plot method, and gained arguments to control
gam fitting
- metaMDSrotate can handle >2D configuration
- metaMDSiter got more flexible 'previous.best'
- typo in vegdist.Rd
Modified: branches/1.17/R/metaMDSiter.R
===================================================================
--- branches/1.17/R/metaMDSiter.R 2011-03-18 13:53:18 UTC (rev 1548)
+++ branches/1.17/R/metaMDSiter.R 2011-03-19 19:31:32 UTC (rev 1549)
@@ -10,11 +10,41 @@
SOL <- FALSE
converged <- FALSE
isotrace <- max(0, trace - 1)
+ ## Previous best or initial configuration
if (!missing(previous.best) && !is.null(previous.best)) {
- s0 <- previous.best
- if (trace)
- cat("Starting from a previous solution\n")
- }
+ ## check if previous.best is from metaMDS or isoMDS
+ if (inherits(previous.best, "metaMDS") ||
+ is.list(previous.best) &&
+ all(c("points", "stress") %in% names(previous.best))) {
+ ## real "previous best"
+ if (NCOL(previous.best$points) == k) {
+ s0 <- previous.best
+ if (trace)
+ cat("Starting from a previous solution\n")
+ } else {
+ init <- previous.best$points
+ nc <- NCOL(init)
+ if (nc > k)
+ init <- init[, 1:k, drop = FALSE]
+ else # nc < k
+ for (i in 1:(k-nc))
+ init <- cbind(init, runif(NROW(init), -0.1, 0.1))
+ # evaluate isoMDS with stress
+ s0 <- isoMDS(dist, init, k = k, maxit = 0)
+ # zero 'tries': this was a new start
+ if (inherits(previous.best, "metaMDS"))
+ previous.best$tries <- 0
+ if (trace)
+ cat(gettextf("Starting from %d-dimensional solution\n", nc))
+ }
+ } else if (is.matrix(previous.best) || is.data.frame(previous.best)) {
+ s0 <- isoMDS(dist, previous.best, k = k, maxit = 0)
+ if (trace)
+ cat("Starting from supplied configuration\n")
+ } else { # an error!
+ stop("'previous.best' of unknown kind")
+ }
+ } # No previous best:
else s0 <- isoMDS(dist, k = k, trace = isotrace)
if (trace)
cat("Run 0 stress", s0$stress, "\n")
@@ -48,10 +78,12 @@
}
flush.console()
}
- if (!missing(previous.best) && !is.null(previous.best$tries))
+ if (!missing(previous.best) && inherits(previous.best, "metaMDS")) {
tries <- tries + previous.best$tries
+ }
out <- list(points = s0$points, dims = k, stress = s0$stress,
- data = attr(dist, "commname"), distance = attr(dist,
- "method"), converged = converged, tries = tries)
+ data = attr(dist, "commname"),
+ distance = attr(dist, "method"), converged = converged,
+ tries = tries)
out
}
Modified: branches/1.17/R/metaMDSrotate.R
===================================================================
--- branches/1.17/R/metaMDSrotate.R 2011-03-18 13:53:18 UTC (rev 1548)
+++ branches/1.17/R/metaMDSrotate.R 2011-03-19 19:31:32 UTC (rev 1549)
@@ -1,27 +1,53 @@
### Rotates metaMDS result so that axis one is parallel to vector 'x'
`metaMDSrotate` <-
- function(object, vec, choices = 1:2, ...)
+ function(object, vec, ...)
{
if (!inherits(object, "metaMDS"))
- stop("function works only with 'metaMDS' results")
- if (length(choices) != 2)
- stop("function can be only used with 2dim plots")
+ stop(gettextf("function works only with 'metaMDS' results"))
+ x <- object$points
+ sp <- object$species
+ N <- NCOL(x)
+ if (N < 2)
+ stop(gettextf("needs at least 2 dimensions"))
vec <- drop(vec)
if (length(dim(vec)) > 1)
- stop("function works only with univariate 'x'")
- ## envfit finds the direction cosine
- rot <- envfit(object, vec, choices = choices, ...)$vectors$arrows
- rot <- drop(rot)
- ## rotation matrix [[sin theta, cos theta] [-cos theta, sin theta]]
- rot <- rbind(rot, rev(rot))
- rot[2,1] <- -rot[2,1]
- ## transpose (inverse) of the rotation matrix
- rot <- t(rot)
- ## Rotation of points and species scores
- object$points[] <- object$points %*% rot
+ stop(gettextf("function works only with univariate 'vec'"))
+ if (!is.numeric(vec))
+ stop(gettextf("'vec' must be numeric"))
+ ## scores must be orthogonal for the next loop to work
+ if (N > 2) {
+ pc <- prcomp(x)
+ x <- pc$x
+ if (!all(is.na(sp)))
+ sp <- sp %*% pc$rotation
+ }
+ ## vectorfit finds the direction cosine. We rotate first axis to
+ ## 'vec' which means that we make other axes orthogonal to 'vec'
+ ## one by one
+ for (k in 2:N) {
+ rot <- vectorfit(x[, c(1,k)], vec, permutations=0)$arrows
+ rot <- drop(rot)
+ ## counterclockwise rotation matrix:
+ ## [cos theta -sin theta]
+ ## [sin theta cos theta]
+ rot <- rbind(rot, rev(rot))
+ rot[1,2] <- -rot[1,2]
+ ## Rotation of points and species scores
+ x[, c(1,k)] <- x[, c(1,k)] %*% rot
+ if (!all(is.na(sp)))
+ sp[, c(1,k)] <- sp[, c(1,k)] %*% rot
+ }
+ ## Rotate 2..N axes to PC
+ if (N > 2 && attr(object$points, "pc")) {
+ pc <- prcomp(x[,-1])
+ x[,-1] <- pc$x
+ if (!all(is.na(sp)))
+ sp[,-1] <- sp[,-1] %*% pc$rotation
+ }
+ ## '[] <-' retains attributes
+ object$points[] <- x
+ object$species[] <- sp
attr(object$points, "pc") <- FALSE
- if (!is.null(object$species))
- object$species[] <- object$species %*% rot
object
}
Modified: branches/1.17/R/ordisurf.R
===================================================================
--- branches/1.17/R/ordisurf.R 2011-03-18 13:53:18 UTC (rev 1548)
+++ branches/1.17/R/ordisurf.R 2011-03-19 19:31:32 UTC (rev 1549)
@@ -1,13 +1,31 @@
`ordisurf` <-
- function (x, y, choices = c(1, 2), knots = 10, family = "gaussian",
- col = "red", thinplate = TRUE, add = FALSE, display = "sites",
+ function(...) UseMethod("ordisurf")
+
+`ordisurf.formula` <-
+ function(formula, data, ...)
+{
+ if (missing(data))
+ data <- parent.frame()
+ x <- formula[[2]]
+ x <- eval.parent(x)
+ formula[[2]] <- NULL
+ y <- drop(as.matrix(model.frame(formula, data, na.action = na.pass)))
+ if (NCOL(y) > 1)
+ stop(gettextf("only one fitted variable allowed in the formula"))
+ ordisurf(x, y, ...)
+}
+
+`ordisurf.default` <-
+ function (x, y, choices = c(1, 2), knots = 10, family = "gaussian",
+ col = "red", thinplate = TRUE, add = FALSE, display = "sites",
w = weights(x), main, nlevels = 10, levels, labcex = 0.6,
- bubble = FALSE, cex = 1, ...)
+ bubble = FALSE, cex = 1, select = FALSE,
+ method = "GCV.Cp", gamma = 1, plot = TRUE, ...)
{
weights.default <- function(object, ...) NULL
GRID = 31
w <- eval(w)
- if (!is.null(w) && length(w) == 1)
+ if (!is.null(w) && length(w) == 1)
w <- NULL
require(mgcv) || stop("Requires package 'mgcv'")
X <- scores(x, choices = choices, display = display, ...)
@@ -26,15 +44,17 @@
mod <- gam(y ~ x1 + x2, family = family, weights = w)
else if (knots == 1)
mod <- gam(y ~ poly(x1, 1) + poly(x2, 1),
- family = family, weights = w)
+ family = family, weights = w, method = method)
else if (knots == 2)
mod <- gam(y ~ poly(x1, 2) + poly(x2, 2) + poly(x1, 1):poly(x2, 1),
- family = family, weights = w)
- else if (thinplate)
- mod <- gam(y ~ s(x1, x2, k = knots), family = family,
- weights = w)
- else mod <- gam(y ~ s(x1, k = knots) + s(x2, k = knots),
- family = family, weights = w)
+ family = family, weights = w, method = method)
+ else if (thinplate)
+ mod <- gam(y ~ s(x1, x2, k = knots), family = family,
+ weights = w, select = select, method = method,
+ gamma = gamma)
+ else mod <- gam(y ~ s(x1, k = knots) + s(x2, k = knots), family = family,
+ weights = w, select = select, method = method,
+ gamma = gamma)
xn1 <- seq(min(x1), max(x1), len=GRID)
xn2 <- seq(min(x2), max(x2), len=GRID)
newd <- expand.grid(x1 = xn1, x2 = xn2)
@@ -51,25 +71,28 @@
as.integer(np), as.double(newd[,1]), as.double(newd[,2]),
inpoly = as.integer(inpoly), PACKAGE="vegan")$inpoly
is.na(fit) <- inpoly == 0
- if (!add) {
- if (bubble) {
- if (is.numeric(bubble))
- cex <- bubble
- cex <- (y - min(y))/diff(range(y)) * (cex-0.4) + 0.4
+ if(plot) {
+ if (!add) {
+ if (bubble) {
+ if (is.numeric(bubble))
+ cex <- bubble
+ cex <- (y - min(y))/diff(range(y)) * (cex-0.4) + 0.4
+ }
+ plot(X, asp = 1, cex = cex, ...)
}
- plot(X, asp = 1, cex = cex, ...)
+ if (!missing(main) || (missing(main) && !add)) {
+ if (missing(main))
+ main <- yname
+ title(main = main)
+ }
+ if (missing(levels))
+ levels <- pretty(range(fit, finite = TRUE), nlevels)
+ ## Only plot surface is select is FALSE or (TRUE and EDF is diff from 0)
+ if(!select || (select && !isTRUE(all.equal(as.numeric(summary(mod)$edf), 0))))
+ contour(xn1, xn2, matrix(fit, nrow=GRID), col = col, add = TRUE,
+ levels = levels, labcex = labcex,
+ drawlabels = !is.null(labcex) && labcex > 0)
}
- if (!missing(main) || (missing(main) && !add)) {
- if (missing(main))
- main <- yname
- title(main = main)
- }
- if (missing(levels))
- levels <- pretty(range(fit, finite = TRUE), nlevels)
-
- contour(xn1, xn2, matrix(fit, nrow=GRID), col = col, add = TRUE,
- levels = levels, labcex = labcex,
- drawlabels = !is.null(labcex) && labcex > 0)
mod$grid <- list(x = xn1, y = xn2, z = matrix(fit, nrow = GRID))
class(mod) <- c("ordisurf", class(mod))
mod
Copied: branches/1.17/R/plot.ordisurf.R (from rev 1548, pkg/vegan/R/plot.ordisurf.R)
===================================================================
--- branches/1.17/R/plot.ordisurf.R (rev 0)
+++ branches/1.17/R/plot.ordisurf.R 2011-03-19 19:31:32 UTC (rev 1549)
@@ -0,0 +1,35 @@
+`plot.ordisurf` <- function(x, what = c("contour","persp","gam"),
+ add = FALSE, bubble = FALSE, col = "red", cex = 1,
+ nlevels = 10, levels, labcex = 0.6, ...) {
+ what <- match.arg(what)
+ y <- x$model$y
+ x1 <- x$model$x1
+ x2 <- x$model$x2
+ X <- x$grid$x
+ Y <- x$grid$y
+ Z <- x$grid$z
+ force(col)
+ force(cex)
+ if(isTRUE(all.equal(what, "contour"))) {
+ if(!add) {
+ if(bubble) {
+ if (is.numeric(bubble))
+ cex <- bubble
+ cex <- (y - min(y))/diff(range(y)) * (cex-0.4) + 0.4
+ }
+ plot(x1, x2, asp = 1, cex = cex, ...)
+ }
+ if (missing(levels))
+ levels <- pretty(range(x$grid$z, finite = TRUE), nlevels)
+ contour(X, Y, Z, col = col, add = TRUE,
+ levels = levels, labcex = labcex,
+ drawlabels = !is.null(labcex) && labcex > 0)
+ } else if(isTRUE(all.equal(what, "persp"))) {
+ persp(X, Y, Z, col = col, cex = cex, ...)
+ } else {
+ class(x) <- class(x)[-1]
+ plot(x, ...) ##col = col, cex = cex, ...)
+ class(x) <- c("ordisurf", class(x))
+ }
+ invisible(x)
+}
Modified: branches/1.17/inst/ChangeLog
===================================================================
--- branches/1.17/inst/ChangeLog 2011-03-18 13:53:18 UTC (rev 1548)
+++ branches/1.17/inst/ChangeLog 2011-03-19 19:31:32 UTC (rev 1549)
@@ -4,8 +4,19 @@
Version 1.17-9 (opened March 9, 2011)
+ * merged 1534: typo in vegdist.Rd.
+
+ * merged 1541,8: metaMDSiter got more flexible 'previous.tries.
+
+ * merged 1538,9, 1542,5: ordisurf gaine new arguments and got a
+ new plot method.
+
+ * merged r1533,7, 1544,7: metaMDSrotate works with > 2 dim.
+
* merged r1532: fix gettextf in d/r/rarefy.
+ * merged r1531,4: ordisurf got formula interface.
+
* merged r1527,8: doc updates in CCorA.Rd & ordihull.Rd.
* merged r1522 thru r1526: implement drarefy; drarefy, rrarefy,
Modified: branches/1.17/man/metaMDS.Rd
===================================================================
--- branches/1.17/man/metaMDS.Rd 2011-03-18 13:53:18 UTC (rev 1548)
+++ branches/1.17/man/metaMDS.Rd 2011-03-19 19:31:32 UTC (rev 1549)
@@ -43,7 +43,7 @@
postMDS(X, dist, pc=TRUE, center=TRUE, halfchange, threshold=0.8,
nthreshold=10, plot=FALSE, ...)
metaMDSredist(object, ...)
-metaMDSrotate(object, vec, choices, ...)
+metaMDSrotate(object, vec, ...)
}
\arguments{
Modified: branches/1.17/man/ordisurf.Rd
===================================================================
--- branches/1.17/man/ordisurf.Rd 2011-03-18 13:53:18 UTC (rev 1548)
+++ branches/1.17/man/ordisurf.Rd 2011-03-19 19:31:32 UTC (rev 1549)
@@ -1,6 +1,9 @@
\name{ordisurf}
\alias{ordisurf}
+\alias{ordisurf.default}
+\alias{ordisurf.formula}
\alias{calibrate.ordisurf}
+\alias{plot.ordisurf}
\title{ Fit and Plot Smooth Surfaces of Variables on Ordination. }
\description{
@@ -8,16 +11,26 @@
plots the result on ordination diagram.
}
\usage{
-ordisurf(x, y, choices=c(1, 2), knots=10, family="gaussian", col="red",
+\method{ordisurf}{default}(x, y, choices=c(1, 2), knots=10, family="gaussian", col="red",
thinplate = TRUE, add = FALSE, display = "sites",
w = weights(x), main, nlevels = 10, levels, labcex = 0.6,
- bubble = FALSE, cex = 1, ...)
+ bubble = FALSE, cex = 1, select = FALSE, method = "GCV.Cp",
+ gamma = 1, plot = TRUE, ...)
+
+\method{ordisurf}{formula}(formula, data, ...)
+
\method{calibrate}{ordisurf}(object, newdata, ...)
+
+\method{plot}{ordisurf}(x, what = c("contour","persp","gam"),
+ add = FALSE, bubble = FALSE, col = "red", cex = 1,
+ nlevels = 10, levels, labcex = 0.6, \dots)
}
\arguments{
- \item{x}{Ordination configuration, either a matrix or a result known
- by \code{\link{scores}}. }
+ \item{x}{For \code{ordisurf} an ordination configuration, either a
+ matrix or a result known by \code{\link{scores}}. For
+ \code{plot.ordisurf} and object of class \code{"ordisurf"} as
+ returned by \code{ordisurf}.}
\item{y}{ Variable to be plotted. }
\item{choices}{Ordination axes. }
\item{knots}{Number of initial knots in \code{\link[mgcv]{gam}} (one
@@ -47,20 +60,72 @@
the maximum. The minimum size will always be \code{cex = 0.4}. The
option only has an effect if \code{add = FALSE}.}
\item{cex}{Character expansion of plotting symbols.}
+ \item{select}{Logical; specify \code{\link[mgcv]{gam}} argument
+ \code{"select"}. If this is \code{TRUE} then \code{\link[mgcv]{gam}} can
+ add an extra penalty to each term so that it can be penalized to
+ zero. This means that the smoothing parameter estimation that is part
+ of fitting can completely remove terms from the model. If the
+ corresponding smoothing parameter is estimated as zero then the extra
+ penalty has no effect.}
+ \item{method}{character; the smoothing parameter estimation
+ method. Options allowed are: \code{"GCV.Cp"} uses GCV for models with
+ unknown scale parameter and Mallows' Cp/UBRE/AIC for models with
+ known scale; \code{"GACV.Cp"} as for \code{"GCV.Cp"} but uses GACV
+ (Generalised Approximate CV) instead of GCV; \code{"REML"} and
+ \code{"ML"} use restricted maximum likelihood or maximum likelihood
+ estimation for both known and unknown scale; and \code{"P-REML"} and
+ \code{"P-ML"} use REML or ML estimation but use a Pearson estimate
+ of the scale.}
+ \item{gamma}{Multiplier to inflate model degrees of freedom in GCV or
+ UBRE/AIC score by. This effectively places an extra penalty on
+ complex models. An oft used value if \code{gamma = 1.4}.}
+ \item{plot}{logical; should any plotting be done by
+ \code{ordisurf}? Useful if all you want is the fitted response
+ surface model.}
+ \item{formula, data}{Alternative definition of the fitted model as
+ \code{x ~ y}, or left-hand side is the ordination \code{x} and
+ right-hand side the single fitted continuous variable
+ \code{y}. The variable \code{y} must be in the working environment
+ or in the data frame or environment given by \code{data}. All
+ other arguments of are passed to the default method.}
\item{object}{An \code{ordisurf} result object.}
- \item{newdata}{Coordinates in two-dimensional ordination for new points.}
- \item{\dots}{ Other graphical parameters. }
+ \item{newdata}{Coordinates in two-dimensional ordination for new
+ points.}
+ \item{what}{character; what type of plot to produce. \code{"contour"}
+ produces a contour plot of the response surface, see
+ \code{\link{contour}} for details. \code{"persp"} produces a
+ perspective plot of the same, see \code{\link{persp}} for
+ details. \code{"gam"} plots the fitted GAM model, an object that
+ inherits from class \code{"gam"} returned by \code{ordisurf}, see
+ \code{\link[mgcv]{plot.gam}}.}
+ \item{\dots}{Other parameters passed to \code{\link[mgcv]{gam}}, or
+ to the graphical functions. See Note below for exceptions.}
}
+
\details{
Function \code{ordisurf} fits a smooth surface using thinplate splines
in \code{\link[mgcv]{gam}}, and uses \code{\link[mgcv]{predict.gam}}
- to find fitted values in a regular grid.
- Function plots the fitted contours with convex hull of data points
- either over an existing ordination diagram or draws a new plot
- The function uses
- \code{\link{scores}} to extract ordination scores, and \code{x} can be
- any result object known by that function.
+ to find fitted values in a regular grid. The smooth surface can be
+ fitted with an extra penalty that allows the entire smoother to be
+ penalized back to 0 degrees of freedom, effectively removing the term
+ from the model. The addition of this extra penalty is invoked by
+ setting argument \code{select} to \code{TRUE}. The function plots the
+ fitted contours with convex hull of data points either over an
+ existing ordination diagram or draws a new plot. If \code{select ==
+ TRUE} and the smooth is effectively penalised out of the model, no
+ contours will be plotted.
+ \code{\link[mgcv]{gam}} determines the degree of smoothness for the
+ fitted response surface during model fitting. Argument \code{method}
+ controls how \code{\link[mgcv]{gam}} performs this smoothness
+ selection. See \code{\link[mgcv]{gam}} for details of the available
+ options. Using \code{"REML"} or \code{"ML"} yields p-values for
+ smooths with the best coverage properties if such things matter to
+ you.
+
+ The function uses \code{\link{scores}} to extract ordination scores,
+ and \code{x} can be any result object known by that function.
+
User can supply a vector of prior weights \code{w}. If the ordination
object has weights, these will be used. In practise this means that
the row totals are used as weights with
@@ -92,11 +157,22 @@
predicting new values (see \code{\link[mgcv]{predict.gam}}).
}
-\author{ Dave Roberts and Jari Oksanen }
+\author{ Dave Roberts, Jari Oksanen and Gavin L. Simpson }
\note{
The default is to use thinplate splines. These make sense in
ordination as they have equal smoothing in all directions and are
- rotation invariant.
+ rotation invariant.
+
+ Graphical arguments supplied to \code{plot.ordisurf} are passed on to
+ the underlying plotting functions, \code{contour}, \code{persp}, and
+ \code{\link[mgcv]{plot.gam}}. The exception to this is that arguments
+ \code{col} and \code{cex} can not currently be passed to
+ \code{\link[mgcv]{plot.gam}} because of a bug in the way that function
+ evaluates arguments when arranging the plot.
+
+ A work-around is to call \code{\link[mgcv]{plot.gam}} directly on the
+ result of a call to \code{ordisurf}. See the Examples for an
+ illustration of this.
}
\seealso{ For basic routines \code{\link[mgcv]{gam}},
@@ -111,10 +187,24 @@
vare.dist <- vegdist(varespec)
vare.mds <- isoMDS(vare.dist)
with(varechem, ordisurf(vare.mds, Baresoil, bubble = 5))
+
+## as above but with extra penalties on smooth terms:
+with(varechem, ordisurf(vare.mds, Baresoil, bubble = 5, col = "blue",
+ add = TRUE, select = TRUE))
+
## Cover of Cladina arbuscula
fit <- with(varespec, ordisurf(vare.mds, Cla.arb, family=quasipoisson))
## Get fitted values
calibrate(fit)
+
+## Plot method
+plot(fit, what = "contour")
+
+## Plotting the "gam" object
+plot(fit, what = "gam") ## 'col' and 'cex' not passed on
+## or via plot.gam directly
+plot.gam(fit, cex = 2, pch = 1, col = "blue")
+## 'col' effects all objects drawn...
}
\keyword{ multivariate }
\keyword{ aplot }
Modified: branches/1.17/man/vegdist.Rd
===================================================================
--- branches/1.17/man/vegdist.Rd 2011-03-18 13:53:18 UTC (rev 1548)
+++ branches/1.17/man/vegdist.Rd 2011-03-19 19:31:32 UTC (rev 1549)
@@ -84,7 +84,7 @@
\cr \tab binary: \eqn{\frac{A+B-2J}{A+B-J}}{(A+B-2*J)/(A+B-J)}
\cr
\code{bray}
- \tab \eqn{d_{jk} = \frac{\sum_i |x_{ij}-x_{ik}|}{\sum_i (x_{ij}+x_{ik})}}{d[jk] = (sum abs(x[ij]-x[ik])/(sum (x[ij]+x[ik]))}
+ \tab \eqn{d_{jk} = \frac{\sum_i |x_{ij}-x_{ik}|}{\sum_i (x_{ij}+x_{ik})}}{d[jk] = (sum abs(x[ij]-x[ik]))/(sum (x[ij]+x[ik]))}
\cr \tab binary: \eqn{\frac{A+B-2J}{A+B}}{(A+B-2*J)/(A+B)}
\cr
\code{kulczynski}
More information about the Vegan-commits
mailing list