[Vegan-commits] r2545 - in branches/2.0: R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jul 9 15:22:12 CEST 2013
Author: jarioksa
Date: 2013-07-09 15:22:12 +0200 (Tue, 09 Jul 2013)
New Revision: 2545
Modified:
branches/2.0/R/ordisurf.R
branches/2.0/R/plot.ordisurf.R
branches/2.0/inst/ChangeLog
branches/2.0/man/ordisurf.Rd
Log:
merge ordisurf upgrades r2490:2493, 2534
Modified: branches/2.0/R/ordisurf.R
===================================================================
--- branches/2.0/R/ordisurf.R 2013-07-09 13:07:03 UTC (rev 2544)
+++ branches/2.0/R/ordisurf.R 2013-07-09 13:22:12 UTC (rev 2545)
@@ -17,17 +17,23 @@
`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",
+ fx = FALSE, 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, lwd.cl = par("lwd"), ...)
{
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 +46,76 @@
}
x1 <- X[, 1]
x2 <- X[, 2]
- if (knots <= 0)
- 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, 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, 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)
+ ## handle fx - allow vector of length up to two
+ if(!(missfx <- missing(fx)) && missing(knots))
+ warning("Requested fixed d.f. splines but without specifying 'knots'.\nSwitching to 'fx = FALSE'.")
+ if (length(fx) > 2L)
+ warning("Length of 'fx' supplied exceeds '2'. Using the first two.")
+ ## expand fx robustly, no matter what length supplied
+ fx <- rep(fx, length.out = 2)
+ ## can't have `fx = TRUE` and `select = TRUE`
+ if(!missfx) { ## fx set by user
+ if((miss.select <- missing(select)) && any(fx)) {
+ warning("'fx = TRUE' requested; using 'select = FALSE'")
+ select <- FALSE
+ } else if(!miss.select && isTRUE(select)){
+ stop("Fixed d.f. splines ('fx = TRUE') incompatible with 'select = TRUE'")
+ }
+ }
+ ## handle knots - allow vector of length up to two
+ if (length(knots) > 2L)
+ warning("Length of 'knots' supplied exceeds '2'. 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."))
+ }
+ ## can't use "cr", "cs", "ps" in 2-d smoother with s()
+ if(isTRUE(isotropic) && any(bs %in% c("cr", "cs", "ps"))) {
+ stop("Bases \"cr\", \"cs\", and \"ps\" not allowed in isotropic smooths.")
+ }
+ ## Build formula
+ if (knots[1] <= 0) {
+ f <- formula(y ~ x1 + x2)
+ } else if (knots[1] == 1) { ## why do we treat this differently?
+ f <- formula(y ~ poly(x1, 1) + poly(x2, 1))
+ } else if (knots[1] == 2) {
+ f <- formula(y ~ poly(x1, 2) + poly(x2, 2) + poly(x1, 1):poly(x2, 1))
+ } else if (isotropic) {
+ f <- formula(paste0("y ~ s(x1, x2, k = ", knots[1],
+ ", bs = \"", bs[1], "\", fx = ", fx[1],")"))
+ } else {
+ if (any(bs %in% c("ad"))) {
+ ## only "ad" for now, but "fs" should also not be allowed
+ f <- formula(paste0("y ~ s(x1, k = ", knots[1],
+ ", bs = \"", bs[1],
+ "\", fx = ", fx[1], ") + s(x2, k = ",
+ knots[2], ", bs = \"", bs[2],
+ "\", fx = ", fx[2], ")"))
+ } else {
+ f <- formula(paste0("y ~ te(x1, x2, k = c(",
+ paste0(knots, collapse = ", "),
+ "), bs = c(",
+ paste0("\"", bs, "\"", collapse = ", "),
+ "), fx = c(",
+ paste0(fx, collapse = ", "),
+ "))"))
+ }
+ }
+ ## fit model
+ mod <- gam(f, 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)
@@ -62,14 +123,17 @@
poly <- chull(cbind(x1,x2))
## Move out points of the convex hull to have contour for all data
## points
- xhull1 <- x1[poly] + sign(x1[poly] - mean(x1[poly])) * diff(range(x1))/(GRID - 1)
- xhull2 <- x2[poly] + sign(x2[poly] - mean(x2[poly])) * diff(range(x2))/(GRID - 1)
+ xhull1 <- x1[poly] + sign(x1[poly] - mean(x1[poly])) *
+ diff(range(x1))/(GRID - 1)
+ xhull2 <- x2[poly] + sign(x2[poly] - mean(x2[poly])) *
+ diff(range(x2))/(GRID - 1)
npol <- length(poly)
np <- nrow(newd)
inpoly <- numeric(np)
- inpoly <- .C("pnpoly", as.integer(npol), as.double(xhull1), as.double(xhull2),
- as.integer(np), as.double(newd[,1]), as.double(newd[,2]),
- inpoly = as.integer(inpoly), PACKAGE="vegan")$inpoly
+ inpoly <- .C("pnpoly", as.integer(npol), as.double(xhull1),
+ as.double(xhull2), as.integer(np), as.double(newd[,1]),
+ as.double(newd[,2]), inpoly = as.integer(inpoly),
+ PACKAGE="vegan")$inpoly
is.na(fit) <- inpoly == 0
if(plot) {
if (!add) {
@@ -88,10 +152,12 @@
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))))
+ 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)
+ drawlabels = !is.null(labcex) && labcex > 0,
+ lwd = lwd.cl)
}
mod$grid <- list(x = xn1, y = xn2, z = matrix(fit, nrow = GRID))
class(mod) <- c("ordisurf", class(mod))
Modified: branches/2.0/R/plot.ordisurf.R
===================================================================
--- branches/2.0/R/plot.ordisurf.R 2013-07-09 13:07:03 UTC (rev 2544)
+++ branches/2.0/R/plot.ordisurf.R 2013-07-09 13:22:12 UTC (rev 2545)
@@ -1,6 +1,7 @@
`plot.ordisurf` <- function(x, what = c("contour","persp","gam"),
add = FALSE, bubble = FALSE, col = "red", cex = 1,
- nlevels = 10, levels, labcex = 0.6, ...) {
+ nlevels = 10, levels, labcex = 0.6,
+ lwd.cl = par("lwd"), ...) {
what <- match.arg(what)
y <- x$model$y
x1 <- x$model$x1
@@ -23,7 +24,8 @@
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)
+ drawlabels = !is.null(labcex) && labcex > 0,
+ lwd = lwd.cl)
} else if(isTRUE(all.equal(what, "persp"))) {
persp(X, Y, Z, col = col, cex = cex, ...)
} else {
Modified: branches/2.0/inst/ChangeLog
===================================================================
--- branches/2.0/inst/ChangeLog 2013-07-09 13:07:03 UTC (rev 2544)
+++ branches/2.0/inst/ChangeLog 2013-07-09 13:22:12 UTC (rev 2545)
@@ -8,9 +8,11 @@
* merge r2539: stressplot return data in input order.
* merge r2538: use expression(R^2) in stressplot.
* merge r2535, 2536: better positioning of arrow labels.
+ * merge r2534: lwd in ordisurf.
* merge r2533: ordipointlabel uses ordiArgAbsorber.
* merge r2504: notation in adipart.Rd.
* merge r2497: avoid visible ~ in .Rnw.
+ * merge r2490 thru 2493: ordisurf new options.
* merge r2484 thru 2486: betadisper fixes (centroid with one
group, call it medoid).
Modified: branches/2.0/man/ordisurf.Rd
===================================================================
--- branches/2.0/man/ordisurf.Rd 2013-07-09 13:07:03 UTC (rev 2544)
+++ branches/2.0/man/ordisurf.Rd 2013-07-09 13:22:12 UTC (rev 2545)
@@ -11,11 +11,13 @@
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", fx = FALSE, 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, lwd.cl = par("lwd"), ...)
\method{ordisurf}{formula}(formula, data, ...)
@@ -23,25 +25,46 @@
\method{plot}{ordisurf}(x, what = c("contour","persp","gam"),
add = FALSE, bubble = FALSE, col = "red", cex = 1,
- nlevels = 10, levels, labcex = 0.6, \dots)
+ nlevels = 10, levels, labcex = 0.6, lwd.cl = par("lwd"), \dots)
}
\arguments{
\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
+ \code{plot.ordisurf} an object of class \code{"ordisurf"} as
returned by \code{ordisurf}.}
- \item{y}{ Variable to be plotted. }
+ \item{y}{Variable to be plotted / modelled as a function of the
+ ordination scores.}
\item{choices}{Ordination axes. }
\item{knots}{Number of initial knots in \code{\link[mgcv]{gam}} (one
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 ordination
+ dimensions (as indicated by \code{choices}) respectively.}
+ \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{bs = "tp"}.}
+ \item{fx}{indicates whether the smoothers are fixed degree of freedom
+ regression splines (\code{fx = FALSE}) or penalised regression
+ splines (\code{fx = TRUE}). Can be a vector of length 2 for
+ anisotropic surfaces (\code{isotropic = FALSE}). It doesn't make
+ sense to use \code{fx = TRUE} \strong{and} \code{select = TRUE} and
+ it is an \strong{error} to do so. A warning is issued if you specify
+ \code{fx = TRUE} and forget to use \code{select = FALSE} though
+ fitting continues using \code{select = FALSE}.}
+ \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,11 +72,14 @@
\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
+ \item{bubble}{Use a \dQuote{bubble plot} for points, or vary the point
diameter by the value of the plotted variable. If \code{bubble} is
numeric, its value is used for the maximum symbol size (as in
\code{cex}), or if \code{bubble = TRUE}, the value of \code{cex} gives
@@ -78,12 +104,14 @@
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}.}
+ complex models. An oft-used value is \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{lwd.cl}{numeric; the \code{lwd} (line width) parameter to use
+ when drawing the contour lines.}
\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
@@ -98,44 +126,52 @@
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
+ \item{\dots}{Other parameters passed to \code{\link{scores}}, or
to the graphical functions. See Note below for exceptions.}
}
\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"}).
+
+ \code{ordisurf()} exposes a large number of options from
+ \code{\link[mgcv]{gam}} for specifying the basis functions used for
+ the surface. If you stray from the defaults, do read the
+ \strong{Notes} section below and relevant documentation in
+ \code{\link[mgcv]{s}} and \code{\link[mgcv]{smooth.terms}}.
+
+ 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.
+ fitted response surface during model fitting, unless \code{fx =
+ TRUE}. 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
- \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
@@ -146,25 +182,46 @@
}
\value{
- Function is usually called for its side effect of drawing the
- contour plot. The function returns the result object of class
+ \code{ordisurf} is usually called for its side effect of drawing the
+ contour plot. The function returns a result object of class
\code{"ordisurf"} that inherits from \code{\link[mgcv]{gam}} used
internally to fit the surface, but adds an item \code{grid} that
contains the data for the grid surface. The item \code{grid} has
elements \code{x} and \code{y} which are vectors of axis coordinates,
and element \code{z} that is a matrix of fitted values for
\code{\link{contour}}. The values outside the convex hull of observed
- points are \code{NA} in \code{z}. The \code{\link[mgcv]{gam}}
- component of the result can be used for further analysis like
- predicting new values (see \code{\link[mgcv]{predict.gam}}).
+ points are indicated as \code{NA} in \code{z}. The
+ \code{\link[mgcv]{gam}} component of the result can be used for
+ further analysis like predicting new values (see
+ \code{\link[mgcv]{predict.gam}}).
}
\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 dimensions are required, 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}}).
+ Cubic regression splines and P splines can \strong{only} be used with
+ \code{isotropic = FALSE}.
+
+ 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.
+
+ To get the old behaviour of \code{ordisurf} use \code{select = FALSE},
+ \code{method = "GCV.Cp"}, \code{fx = FALSE}, and \code{bs = "tp"}. The
+ latter two options are the current defaults.
+
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
@@ -177,6 +234,20 @@
illustration of this.
}
+\section{Warning}{
+ The fitted GAM is a regression model and has the usual assumptions of
+ such models. Of particular note is the assumption of independence of
+ residuals. If the observations are not independent (e.g. they are
+ repeat measures on a set of objects, or from an experimental design,
+ \emph{inter alia}) do not trust the \emph{p}-values from the GAM
+ output.
+
+ If you need further control (i.e. to add additional fixed effects to
+ the model, or use more complex smoothers), extract the ordination
+ scores using the \code{scores} function and then generate your own
+ \code{\link[mgcv]{gam}} call.
+}
+
\references{
Marra, G.P & Wood, S.N. (2011) Practical variable selection for
@@ -200,15 +271,28 @@
vare.mds <- monoMDS(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))
+## as above but without the extra penalties on smooth terms,
+## and using GCV smoothness selection (old behaviour of `ordisurf()`):
+with(varechem, ordisurf(vare.mds, Baresoil,col = "blue", add = TRUE,
+ select = FALSE, method = "GCV.Cp"))
## Cover of Cladina arbuscula
fit <- with(varespec, ordisurf(vare.mds, Cla.arb, family=quasipoisson))
## Get fitted values
calibrate(fit)
+## 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` --- the *default*
+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" with `isotropic = FALSE`
+
## Plot method
plot(fit, what = "contour")
@@ -217,6 +301,25 @@
## or via plot.gam directly
plot.gam(fit, cex = 2, pch = 1, col = "blue")
## 'col' effects all objects drawn...
+
+### controlling the basis functions used
+## Use Duchon splines
+with(varechem, ordisurf(vare.mds, Baresoil, bs = "ds"))
+
+## A fixed degrees of freedom smooth, must use 'select = FALSE'
+with(varechem, ordisurf(vare.mds, Baresoil, knots = 4,
+ fx = TRUE, select = FALSE))
+
+## An anisotropic smoother with cubic regression spline bases
+with(varechem, ordisurf(vare.mds, Baresoil, isotropic = FALSE,
+ bs = "cr", knots = 4))
+
+## An anisotropic smoother with cubic regression spline with
+## shrinkage bases & different degrees of freedom in each dimension
+with(varechem, ordisurf(vare.mds, Baresoil, isotropic = FALSE,
+ bs = "cs", knots = c(3,4), fx = TRUE,
+ select = FALSE))
+
}
\keyword{ multivariate }
\keyword{ aplot }
More information about the Vegan-commits
mailing list