[Vegan-commits] r2490 - in pkg/vegan: . R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Apr 20 05:28:22 CEST 2013
Author: gsimpson
Date: 2013-04-20 05:28:22 +0200 (Sat, 20 Apr 2013)
New Revision: 2490
Modified:
pkg/vegan/DESCRIPTION
pkg/vegan/R/ordisurf.R
pkg/vegan/inst/ChangeLog
pkg/vegan/man/ordisurf.Rd
Log:
big changes to ordisurf in one nice commit that can be reverted easily. See Changelog for details.
Modified: pkg/vegan/DESCRIPTION
===================================================================
--- pkg/vegan/DESCRIPTION 2013-04-09 12:13:31 UTC (rev 2489)
+++ pkg/vegan/DESCRIPTION 2013-04-20 03:28:22 UTC (rev 2490)
@@ -1,7 +1,7 @@
Package: vegan
Title: Community Ecology Package
-Version: 2.1-28
-Date: March 17, 2013
+Version: 2.1-29
+Date: April 19, 2013
Author: Jari Oksanen, F. Guillaume Blanchet, Roeland Kindt, Pierre Legendre,
Peter R. Minchin, R. B. O'Hara, Gavin L. Simpson, Peter Solymos,
M. Henry H. Stevens, Helene Wagner
Modified: pkg/vegan/R/ordisurf.R
===================================================================
--- pkg/vegan/R/ordisurf.R 2013-04-09 12:13:31 UTC (rev 2489)
+++ pkg/vegan/R/ordisurf.R 2013-04-20 03:28:22 UTC (rev 2490)
@@ -17,17 +17,24 @@
`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, select = FALSE,
- method = "GCV.Cp", gamma = 1, plot = TRUE, ...)
+ col = "red", isotropic = TRUE, thinplate = TRUE, bs = "tp",
+ add = FALSE, display = "sites",
+ w = weights(x), main, nlevels = 10, levels,
+ npoints = 31, labcex = 0.6,
+ bubble = FALSE, cex = 1, select = TRUE,
+ method = "REML", gamma = 1, plot = TRUE, ...)
{
weights.default <- function(object, ...) NULL
- GRID = 31
+ if(!missing(thinplate)) {
+ warning("Use of 'thinplate' is deprecated and will soon be removed;\nuse 'isotropic' instead.")
+ isotropic <- thinplate
+ }
+ ## GRID no user-definable - why 31?
+ GRID <- npoints
w <- eval(w)
if (!is.null(w) && length(w) == 1)
w <- NULL
- require(mgcv) || stop("Requires package 'mgcv'")
+ require(mgcv) || stop("Requires package 'mgcv'")
X <- scores(x, choices = choices, display = display, ...)
## The original name of 'y' may be lost in handling NA: save for
## plots
@@ -40,21 +47,52 @@
}
x1 <- X[, 1]
x2 <- X[, 2]
- if (knots <= 0)
+ ## handle knots - allow vector of length up to two
+ if (length(knots) > 2L)
+ warning("Number of knots supplied exceeds '2'. Only using the first two.")
+ ## expand knots robustly, no matter what length supplied
+ knots <- rep(knots, length.out = 2)
+ ## handle the bs - we only allow some of the possible options
+ if (length(bs) > 2L)
+ warning("Number of basis types supplied exceeds '2'. Only using the first two.")
+ bs <- rep(bs, length.out = 2)
+ ## check allowed types
+ BS <- c("tp","ts","cr","cs","ds","ps","ad")
+ want <- match(bs, BS)
+ user.bs <- bs ## store supplied (well expanded supplied ones)
+ bs <- BS[want]
+ if (any(wrong <- is.na(bs))) {
+ stop(paste("Supplied basis type of",
+ paste(sQuote(unique(user.bs[wrong])), collapse = ", "),
+ "not supported."))
+ }
+ ## linear fit - note we match something close to 0
+ if (knots[1] <= 0) {
mod <- gam(y ~ x1 + x2, family = family, weights = w)
- else if (knots == 1)
+ } else if (knots[1] == 1) { ## why do we treat this differently?
mod <- gam(y ~ poly(x1, 1) + poly(x2, 1),
family = family, weights = w, method = method)
- else if (knots == 2)
+ } else if (knots[1] == 2) {
mod <- gam(y ~ poly(x1, 2) + poly(x2, 2) + poly(x1, 1):poly(x2, 1),
family = family, weights = w, method = method)
- else if (thinplate)
- mod <- gam(y ~ s(x1, x2, k = knots), family = family,
+ } else if (isotropic) {
+ mod <- gam(y ~ s(x1, x2, k = knots[1], bs = bs[1]),
+ 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)
+ } else {
+ if (any(bs %in% c("ad"))) { ## only "ad" for now, but "fs" should also not be allowed
+ mod <- gam(y ~ s(x1, k = knots[1], bs = bs[1]) +
+ s(x2, k = knots[2], bs[2]),
+ family = family, weights = w, select = select,
+ method = method, gamma = gamma)
+ } else {
+ mod <- gam(y ~ te(x1, x2, k = knots, bs = bs),
+ 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)
Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog 2013-04-09 12:13:31 UTC (rev 2489)
+++ pkg/vegan/inst/ChangeLog 2013-04-20 03:28:22 UTC (rev 2490)
@@ -2,8 +2,31 @@
VEGAN DEVEL VERSIONS at http://r-forge.r-project.org/
-Version 2.1-28 (opened March 17, 2013)
+Version 2.1-29 (opened April 19, 2013)
+ * ordisurf: significant changes were made to this function:
+
+ - The default for `method` and `select` were changed to `"REML"`
+ and `TRUE` respectivelt.
+
+ - Argument `thinplate` is deprecated in favour of `isotropic`. A
+ warning is now issued if `thinplate` is used.
+
+ - The spline basis for the smoother can now be specified from a
+ subset of those implemented in the mgcv package. This is achieved
+ via the `bs` argument, which defaults to `"tp"` for thin plate
+ regression splines.
+
+ - Argument `knots` and `bs` can now be a vector of length two, one
+ per ordination dimension considered. This is only of use with
+ anisotropic surfaces with `isotropic = FALSE`.
+
+ - The number of locations in each ordination dimension at which the
+ fitted surface is evaluated can now be specified via new argument
+ `npoints`.
+
+Version 2.1-28 (closed April 19, 2013)
+
* betadisper: failed with type = "centroid" when there was only
one group (i.e., in estimating the overall beta diversity in the
data). Reported by Pierre Legendre.
Modified: pkg/vegan/man/ordisurf.Rd
===================================================================
--- pkg/vegan/man/ordisurf.Rd 2013-04-09 12:13:31 UTC (rev 2489)
+++ pkg/vegan/man/ordisurf.Rd 2013-04-20 03:28:22 UTC (rev 2490)
@@ -11,11 +11,12 @@
plots the result on ordination diagram.
}
\usage{
-\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, select = FALSE, method = "GCV.Cp",
- gamma = 1, plot = TRUE, ...)
+\method{ordisurf}{default}(x, y, choices = c(1, 2), knots = 10,
+ family = "gaussian", col = "red", isotropic = TRUE,
+ thinplate = TRUE, bs = "tp", add = FALSE, display = "sites",
+ w = weights(x), main, nlevels = 10, levels, npoints = 31,
+ labcex = 0.6, bubble = FALSE, cex = 1, select = TRUE,
+ method = "REML", gamma = 1, plot = TRUE, ...)
\method{ordisurf}{formula}(formula, data, ...)
@@ -37,11 +38,23 @@
more than degrees of freedom). If \code{knots = 0} or
\code{knots = 1} the function will fit a linear trend surface, and
if \code{knots = 2} the function will fit a quadratic trend surface
- instead of a smooth surface. }
- \item{family}{ Error distribution in \code{\link[mgcv]{gam}}. }
+ instead of a smooth surface. A vector of length 2 is allowed when
+ \code{isotropic = FALSE}, with the first and second elements of
+ \code{knots} refering to the first and second of the ordination
+ dimensions indicated by \code{choices}.}
+ \item{family}{Error distribution in \code{\link[mgcv]{gam}}.}
\item{col}{ Colour of contours. }
- \item{thinplate}{Use thinplate splines in \code{\link[mgcv]{gam}}.}
- \item{add}{Add contours on an existing diagram or draw a new plot. }
+ \item{isotropic, thinplate}{Fit an isotropic smooth surface (i.e. same
+ smoothness in both ordination dimensions) via
+ \code{\link[mgcv]{gam}}. Use of \code{thinplate} is deprecated and
+ will be removed in a future version of the package.}
+ \item{bs}{a two letter character string indicating the smoothing basis
+ to use. (eg \code{"tp"} for thin plate regression spline,
+ \code{"cr"} for cubic regression spline). One of \code{c("tp", "ts",
+ "cr", "cs", "ds", "ps", "ad")}. See
+ \code{\link[mgcv]{smooth.terms}} for an over view of what these
+ refer to. The default is to use thin plate splines \code{"tp"}.}
+ \item{add}{Add contours to an existing diagram or draw a new plot.}
\item{display}{Type of scores known by \code{\link{scores}}: typically
"sites" for ordinary site scores or "lc" for linear combination scores.}
\item{w}{Prior weights on the data. Concerns mainly \code{\link{cca}}
@@ -49,8 +62,11 @@
\item{main}{The main title for the plot, or as default the name of
plotted variable in a new plot.}
\item{nlevels, levels}{Either a vector of \code{levels} for which contours
- are drawn, or suggested number of contours in
- \code{nlevels} if \code{levels} are not supplied.}
+ are drawn, or suggested number of contours in \code{nlevels} if
+ \code{levels} are not supplied.}
+ \item{npoints}{numeric; the number of locations at which to evaluate
+ the fitted surface. This represents the number of locations in each
+ dimension.}
\item{labcex}{Label size in contours. Setting this zero will suppress
labels.}
\item{bubble}{Use \dQuote{bubble plot} for points, or vary the point
@@ -83,7 +99,7 @@
\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
+ \code{x ~ y}, where 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
@@ -104,17 +120,20 @@
\details{
- Function \code{ordisurf} fits a smooth surface using thinplate
+ Function \code{ordisurf} fits a smooth surface using penalised
splines (Wood 2003) in \code{\link[mgcv]{gam}}, and uses
\code{\link[mgcv]{predict.gam}} 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 (see Marra &
Wood, 2011). 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
+ setting argument \code{select} to \code{TRUE}. An alternative is to
+ use a spline basis that includes shrinkage (\code{bs = "ts"} or
+ \code{bs = "cs"}).
+
+ 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
@@ -128,14 +147,13 @@
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
- \code{\link{cca}} or
- \code{\link{decorana}} results. If you do not like this, but want to give
- equal weights to all sites, you should set \code{w = NULL}. The
- behaviour is consistent with \code{\link{envfit}}. For complete
- accordance with constrained \code{\link{cca}}, you should set
+ The 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 \code{\link{cca}}
+ or \code{\link{decorana}} results. If you do not like this, but want
+ to give equal weights to all sites, you should set \code{w =
+ NULL}. The behaviour is consistent with \code{\link{envfit}}. For
+ complete accordance with constrained \code{\link{cca}}, you should set
\code{display = "lc"} (and possibly \code{scaling = 2}).
Function \code{calibrate} returns the fitted values of the response
@@ -161,10 +179,22 @@
\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.
+ The default is to use an isotropic smoother via
+ \code{\link[mgcv]{s}} employing thin plate regression splines
+ (\code{bs = "tp"}). These make sense in ordination as they have
+ equal smoothing in all directions and are rotation invariant. However,
+ if different degrees of smoothness along one dimension, an anisotropic
+ smooth surface may be more applicable. This can be achieved through
+ the use of \code{isotropic = FALSE}, wherein the surface is fitted via
+ a tensor product smoother via \code{\link[mgcv]{te}} (unless \code{bs
+ = "ad"}, in which case separate splines for each dimension are
+ fitted using \code{\link[mgcv]{s}}).
+ Adaptive smooths (\code{bs = "ad"}), especially in two dimensions,
+ require a large number of observations; without many hundreds of
+ observations, the default complexities for the smoother will exceed
+ the number of observations and fitting will fail.
+
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
@@ -217,6 +247,18 @@
## or via plot.gam directly
plot.gam(fit, cex = 2, pch = 1, col = "blue")
## 'col' effects all objects drawn...
+
+## Variable selection via additional shrinkage penalties
+## This allows non-significant smooths to be selected out
+## of the model not just to a linear surface. There are 2
+## options available:
+## - option 1: `select = TRUE`
+with(varechem,
+ ordisurf(vare.mds, Baresoil, method = "REML", select = TRUE))
+## - option 2: use a basis with shrinkage
+with(varechem,
+ ordisurf(vare.mds, Baresoil, method = "REML", bs = "ts"))
+## or bs = "cs"
}
\keyword{ multivariate }
\keyword{ aplot }
More information about the Vegan-commits
mailing list