[Vegan-commits] r338 - in pkg: R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed May 7 15:30:37 CEST 2008
Author: jarioksa
Date: 2008-05-07 15:30:37 +0200 (Wed, 07 May 2008)
New Revision: 338
Added:
pkg/R/MOStest.R
pkg/R/plot.MOStest.R
pkg/R/print.MOStest.R
pkg/man/MOStest.Rd
Modified:
pkg/inst/ChangeLog
Log:
Added functions for Mitchell-Olds -- Shaw test
Added: pkg/R/MOStest.R
===================================================================
--- pkg/R/MOStest.R (rev 0)
+++ pkg/R/MOStest.R 2008-05-07 13:30:37 UTC (rev 338)
@@ -0,0 +1,38 @@
+`MOStest` <-
+ function(x, y, interval, ...)
+{
+ if (!missing(interval))
+ interval <- sort(interval)
+ x <- eval(x)
+ m0 <- glm(y ~ x + I(x^2), ...)
+ k <- coef(m0)
+ isHump <- unname(k[3] < 0)
+ hn <- if(isHump) "hump" else "pit"
+ hump <- unname(-k[2]/2/k[3])
+ if (missing(interval))
+ p1 <- min(x)
+ else
+ p1 <- interval[1]
+ if (missing(interval))
+ p2 <- max(x)
+ else
+ p2 <- interval[2]
+ tmp <- glm(y ~ I(x^2 - 2*x*p1) + x, ...)
+ statmin <- coef(summary(tmp))[3, 3:4]
+ tmp <- glm(y ~ I(x^2 - 2*x*p2) + x, ...)
+ statmax <- coef(summary(tmp))[3, 3:4]
+ comb <- 1 - (1-statmin[2])*(1-statmax[2])
+ stats <- rbind(statmin, statmax)
+ rownames(stats) <- paste(hn, c("at min", "at max"))
+ stats <- cbind("min/max" = c(p1,p2), stats)
+ stats <- rbind(stats, "Combined" = c(NA, NA, comb))
+ vec <- c(p1, p2, hump)
+ names(vec) <- c("min", "max", hn)
+ vec <- sort(vec)
+ isBracketed <- names(vec)[2] == hn
+ out <- list(isHump = isHump, isBracketed = isBracketed,
+ hump = vec, family = family(m0), coefficients = stats,
+ mod = m0)
+ class(out) <- "MOStest"
+ out
+}
Added: pkg/R/plot.MOStest.R
===================================================================
--- pkg/R/plot.MOStest.R (rev 0)
+++ pkg/R/plot.MOStest.R 2008-05-07 13:30:37 UTC (rev 338)
@@ -0,0 +1,35 @@
+`plot.MOStest` <-
+ function(x, which = c(1,2,3,6), ...)
+{
+
+ show <- rep(FALSE, 8)
+ show[which] <- TRUE
+ if (show[1]) {
+ X <- x$mod$model$x
+ Y <- x$mod$y
+ xx <- seq(min(X), max(X), len=101)
+ pre <- predict(x$mod, newdata=list(x = xx), se=TRUE)
+ g <- x$mod$family$linkinv
+ fv <- g(pre$fit)
+ hi <- g(pre$fit + 2*pre$se)
+ lo <- g(pre$fit - 2*pre$se)
+ plot(X, Y, ...)
+ matlines(xx, cbind(fv, hi, lo), lty=c(1, 2, 2), lwd=c(2, 1, 1), col=1, ...)
+ }
+ if (show[2]) {
+ require(ellipse) || stop("requires package 'ellipse'")
+ ci <- ellipse(x$mod, which=c(2,3))
+ plot(ci, type="l", lwd=2, xlim=range(ci[,1],0), ylim=range(ci[,2],0), ...)
+ abline(h=0, lty=2, ...)
+ abline(0, -1/2/x$hump["min"], ...)
+ abline(0, -1/2/x$hump["max"], ...)
+ k <- coef(summary(x$mod))[2:3, 1:2]
+ segments(k[1,1] - k[1,2]*2, k[2,1], k[1,1]+k[1,2]*2, k[2,1], lty=3)
+ segments(k[1,1], k[2,1]-k[2,2]*2, k[1,1], k[2,1]+k[2,2]*2, lty=3)
+ }
+ if (any(show[-c(1,2)])) {
+ still <- which(show[-c(1,2)])
+ plot(x$mod, which = still, ...)
+ }
+ invisible()
+}
Added: pkg/R/print.MOStest.R
===================================================================
--- pkg/R/print.MOStest.R (rev 0)
+++ pkg/R/print.MOStest.R 2008-05-07 13:30:37 UTC (rev 338)
@@ -0,0 +1,13 @@
+`print.MOStest` <-
+ function(x, ...)
+{
+ cat("\nMitchell-Olds and Shaw test\n")
+ cat("Null: hump of a quadratic linear predictor is at min or max\n")
+ print(x$family)
+ print(x$hump)
+ if (!x$isBracketed)
+ cat("***** Caution: hump/pit not bracketed by the data ******\n")
+ cat("\n")
+ printCoefmat(coef(x), has.P=TRUE, na.print="")
+ invisible(x)
+}
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2008-05-07 08:05:00 UTC (rev 337)
+++ pkg/inst/ChangeLog 2008-05-07 13:30:37 UTC (rev 338)
@@ -4,6 +4,9 @@
Version 1.12-14 (opened May 7, 2008)
+ * MOStest: new function to implement Mitchell-Olds & Shaw test for
+ the location of quadratic extreme in a defined interval.
+
* capscale: accepts now other dissimilarity function than vegdist,
and optionally uses metaMDSdist to manipulate dissimilarities
similarly as metaMDS. This provides now a one-shot unconstrained
Added: pkg/man/MOStest.Rd
===================================================================
--- pkg/man/MOStest.Rd (rev 0)
+++ pkg/man/MOStest.Rd 2008-05-07 13:30:37 UTC (rev 338)
@@ -0,0 +1,119 @@
+\encoding{UTF-8}
+\name{MOStest}
+\alias{MOStest}
+\alias{print.MOStest}
+\alias{plot.MOStest}
+
+
+\title{ Mitchell-Olds \& Shaw Test for the Location of Quadratic Extreme }
+\description{
+ Mitchell-Olds & Shaw test concerns the location of the highest (hump)
+ or lowest (pit) value of a quadratic curve at given points. Typically,
+ it is used to study whether the quadratic hump or pit is located
+ within a studied interval. The current test is generalized so that it
+ applies generalized linear models (\code{\link{glm}}) with link
+ function instead of simple quadratic curve.
+}
+\usage{
+MOStest(x, y, interval, ...)
+\method{plot}{MOStest}(x, which = c(1,2,3,6), ...)
+}
+
+\arguments{
+ \item{x}{The independent variable. }
+ \item{y}{The dependent variable. }
+ \item{interval}{The two points at which the test statistic is
+ evaluated. If missing, the estremes of \code{x} are used. }
+ \item{which}{Subset of plots produced. Values \code{which=1} and
+ \code{2} define plots specific to \code{MOStest} (see Details), and
+ larger values select a graphs of \code{\link{plot.lm}} (minus 2). }
+ \item{\dots}{ Other variables passed to functions. Function
+ \code{mitchell.olds.test} passes these to \code{\link{glm}} so that
+ these can include \code{\link{family}}. The other functions pass
+ these to underlying graphical functions. }
+}
+\details{
+ The function fits a quadratic curve with given \code{\link{family}}
+ and link function. The origin of \code{x} is shifted to the values
+ given in \code{interval} (defaults to the extremes of \code{x}), and
+ the quadratic model is refitted. If the optimum is is located at the
+ shifted zero, the first degree coefficient of the polynomial will be
+ zero. The test statistic is the value of the first degree coefficient
+ with its significance (Mitchell-Olds & Shaw 1987).
+
+ The test is often presented as a general test for the location of the
+ hump, but it really is dependent on the quadratic fitted curve. If the
+ hump is of different form than quadratic, the test may be
+ insignificant.
+
+ Because of strong assumptions in the test, you should use the support
+ functions to inspect the fit. Function \code{plot(..., which=1)}
+ displays the data points, fitted quadratic model, and its approximate
+ 95\% confidence intervals (2 times SE). Function \code{plot} with
+ \code{which = 2)} (requires \code{\link[ellipse]{ellipse.glm}} in
+ package \pkg{ellipse}) displays the approximate confidence interval of
+ the polynomial coefficients, together with two lines indicating the
+ combinations of the coefficients that produce the evaluated points of
+ \code{x}. Moreover, the cross-hair shows the approximate confidence
+ intervals (2 times SE) for the polynomial coefficients ignoring their
+ correlations. Higher values of \code{which} produce corresponding
+ graphs from \code{\link{plot.lm}}. That is, you must add 2 to the
+ value of \code{which} in \code{\link{plot.lm}}.
+
+ The test is typically used in assessing the significance of diversity
+ hump against productivity gradient (Mittelbach et al. 2001). It also
+ can be used for the location of the pit (deepest points) instead of
+ the Tokeshi test.
+
+ The method is basically similar as the one suggested by Oksanen et
+ al. (2001) for assessing the confidence intervals of the Gaussian
+ response function.
+}
+\value{
+ The function is based on \code{\link{glm}}, and it returns the result
+ of object of \code{glm} amended with the result of the test. The new
+ items in the \code{MOStest} are:
+ \item{isHump }{\code{TRUE} if the response is a
+ hump.}
+ \item{isBracketed}{\code{TRUE} if the hump or the pit is bracketed by
+ the evaluated points.}
+ \item{hump}{Sorted vector of location of the hump or the pit and the
+ points where the test was evaluated.}
+ \item{coefficients}{Table of test statistics and their signficances.}
+}
+\references{
+Mitchell-Olds, T. & Shaw, R.G. 1987. Regression analysis of natural
+selection: statistical inference and biological
+interpretation. \emph{Evolution} 41, 1149--1161.
+
+Mittelbach, G.C. Steiner, C.F., Scheiner, S.M., Gross, K.L., Reynolds,
+H.L., Waide, R.B., Willig, R.M., Dodson, S.I. & Gough, L. 2001. What is
+the observed richness between species richness and productivity?
+\emph{Ecology} 82: 2381--2396.
+
+Oksanen, J., Läärä, E., Tolonen, K. & Warner, B.G. 2001. Confidence
+intervals for the optimum in the Gaussian response
+function. \emph{Ecology} 82, 1191--1197.
+}
+\author{Jari Oksanen }
+
+\seealso{The nointeraction model can be fitted with \code{\link{humpfit}}. }
+\examples{
+## The Al-Mufti data analysed in humpfit():
+mass <- c(140,230,310,310,400,510,610,670,860,900,1050,1160,1900,2480)
+spno <- c(1, 4, 3, 9, 18, 30, 20, 14, 3, 2, 3, 2, 5, 2)
+mod <- MOStest(mass, spno)
+## Insignificant
+mod
+## ... but inadequate shape of the curve
+op <- par(mfrow=c(2,2), mar=c(4,4,1,1)+.1)
+plot(mod)
+## Looks rather like log-link with Poisson error and logarithmic biomass
+mod <- MOStest(log(mass), spno, family=quasipoisson)
+mod
+plot(mod)
+par(op)
+}
+
+\keyword{ models }
+\keyword{ regression }
More information about the Vegan-commits
mailing list