[Vegan-commits] r360 - in branches/1.13: R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon May 19 09:30:24 CEST 2008
Author: jarioksa
Date: 2008-05-19 09:30:24 +0200 (Mon, 19 May 2008)
New Revision: 360
Removed:
branches/1.13/R/MOStest.R
branches/1.13/R/confint.MOStest.R
branches/1.13/R/fieller.MOStest.R
branches/1.13/R/plot.MOStest.R
branches/1.13/R/print.MOStest.R
branches/1.13/R/profile.MOStest.R
branches/1.13/man/MOStest.Rd
Modified:
branches/1.13/inst/ChangeLog
branches/1.13/inst/NEWS
Log:
MOStest removed from release branches/1.13
Deleted: branches/1.13/R/MOStest.R
===================================================================
--- branches/1.13/R/MOStest.R 2008-05-18 20:10:06 UTC (rev 359)
+++ branches/1.13/R/MOStest.R 2008-05-19 07:30:24 UTC (rev 360)
@@ -1,38 +0,0 @@
-`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
-}
Deleted: branches/1.13/R/confint.MOStest.R
===================================================================
--- branches/1.13/R/confint.MOStest.R 2008-05-18 20:10:06 UTC (rev 359)
+++ branches/1.13/R/confint.MOStest.R 2008-05-19 07:30:24 UTC (rev 360)
@@ -1,6 +0,0 @@
-`confint.MOStest` <-
- function (object, parm = 1, level = 0.95, ...)
-{
- require(MASS) || stop("requires packages MASS")
- confint(profile(object), level = level, ...)
-}
Deleted: branches/1.13/R/fieller.MOStest.R
===================================================================
--- branches/1.13/R/fieller.MOStest.R 2008-05-18 20:10:06 UTC (rev 359)
+++ branches/1.13/R/fieller.MOStest.R 2008-05-19 07:30:24 UTC (rev 360)
@@ -1,31 +0,0 @@
-`fieller.MOStest` <-
- function (object, level = 0.95)
-{
- smodel <- summary(object$mod)
- var <- smodel$cov.scaled
- fam <- family(object$mod)
- od <- smodel$dispersion
- k <- coef(object$mod)
- b2 <- -2 * k[3]
- u <- -k[2]/2/k[3]
- alpha <- (1-level)/2
- limits <- numeric(2)
- names(limits) <- paste(round(100*(c(alpha, 1-alpha)), 1), "%")
- wvar <- var[2,2] * od
- uvar <- 4 * var[3,3] * od
- vvar <- -2 * var[2,3] * od
- z <- qnorm(1 - alpha)
- g <- z^2 * uvar/b2^2
- if (g >= 1) {
- limits <- c(NA, NA)
- }
- else {
- x <- u - g * vvar/uvar
- f <- z/b2
- s <- sqrt(wvar - 2 * u * vvar + u^2 * uvar - g * (wvar -
- vvar^2/uvar))
- limits[1] <- (x - f * s)/(1 - g)
- limits[2] <- (x + f * s)/(1 - g)
- }
- limits
-}
Deleted: branches/1.13/R/plot.MOStest.R
===================================================================
--- branches/1.13/R/plot.MOStest.R 2008-05-18 20:10:06 UTC (rev 359)
+++ branches/1.13/R/plot.MOStest.R 2008-05-19 07:30:24 UTC (rev 360)
@@ -1,35 +0,0 @@
-`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()
-}
Deleted: branches/1.13/R/print.MOStest.R
===================================================================
--- branches/1.13/R/print.MOStest.R 2008-05-18 20:10:06 UTC (rev 359)
+++ branches/1.13/R/print.MOStest.R 2008-05-19 07:30:24 UTC (rev 360)
@@ -1,13 +0,0 @@
-`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)
-}
Deleted: branches/1.13/R/profile.MOStest.R
===================================================================
--- branches/1.13/R/profile.MOStest.R 2008-05-18 20:10:06 UTC (rev 359)
+++ branches/1.13/R/profile.MOStest.R 2008-05-19 07:30:24 UTC (rev 360)
@@ -1,55 +0,0 @@
-`profile.MOStest` <-
- function(fitted, alpha = 0.01, maxsteps = 10, del = zmax/5, ...)
-{
- Pnam <- if(fitted$isHump) "hump" else "pit"
- k <- coef(fitted$mod)
- u <- -k[2]/2/k[3]
- n <- length(residuals(fitted$mod))
- std.error <- fieller.MOStest(fitted, level=0.6)
- std.error <- u - std.error[1]
- if (is.na(std.error))
- std.error <- diff(range(model.matrix(fitted$mod)[,2]))
- OrigDev <- deviance(fitted$mod)
- summ <- summary(fitted$mod)
- DispPar <- summ$dispersion
- fam <- family(fitted$mod)
- Y <- fitted$mod$y
- X <- model.matrix(fitted$mod)[,-3]
- Xi <- X
- if (fam$family %in% c("poisson", "binomial", "Negative Binomial")) {
- zmax <- sqrt(qchisq(1 - alpha/2, 1))
- profName <- "z"
- } else {
- zmax <- sqrt(qf(1 - alpha/2, 1, n - 1))
- profName <- "tau"
- }
- zi <- 0
- prof <- vector("list", length=1)
- names(prof) <- Pnam
- uvi <- u
- for (sgn in c(-1, 1)) {
- step <- 0
- z <- 0
- while((step <- step + 1) < maxsteps && abs(z) < zmax) {
- ui <- u + sgn * step * del * std.error
- Xi[,2] <- (X[,2] - ui)^2
- fm <- glm.fit(x = Xi, y = Y, family=fam,
- control = fitted$mod$control)
- uvi <- c(uvi, ui)
- zz <- (fm$deviance - OrigDev)/DispPar
- z <- sgn * sqrt(zz)
- zi <- c(zi, z)
- }
- si <- order(zi)
- prof[[Pnam]] <- structure(data.frame(zi[si]), names=profName)
- uvi <- as.matrix(uvi)
- colnames(uvi) <- Pnam
- prof[[Pnam]]$par.vals <- uvi[si, , drop=FALSE]
- }
- of <- list()
- of$coefficients <- structure(Pnam, names=Pnam)
- val <- structure(prof, original.fit = of, summary = summ)
- class(val) <- c("profile.MOStest", "profile.glm", "profile")
- val
-}
-
Modified: branches/1.13/inst/ChangeLog
===================================================================
--- branches/1.13/inst/ChangeLog 2008-05-18 20:10:06 UTC (rev 359)
+++ branches/1.13/inst/ChangeLog 2008-05-19 07:30:24 UTC (rev 360)
@@ -3,6 +3,11 @@
VEGAN DEVEL VERSIONS at http://r-forge.r-project.org/
Version 1.13-0 (opened May 14, 2008)
+
+ * Branched version 1.12-15 (rev 354) for a new release.
+
+ * MOStest: removed from the release version. It seemed to work
+ correctly, but strangely: needs more research.
Version 1.12-15 (closed May 14, 2008)
Modified: branches/1.13/inst/NEWS
===================================================================
--- branches/1.13/inst/NEWS 2008-05-18 20:10:06 UTC (rev 359)
+++ branches/1.13/inst/NEWS 2008-05-19 07:30:24 UTC (rev 360)
@@ -5,7 +5,8 @@
GENERAL
- - User visible changes in the development version to be 1.13-0.
+ - User visible changes in the development version to be 1.13-0
+ (rev354).
NEW FUNCTIONS
@@ -13,15 +14,6 @@
Koleff et al. (J. Anim. Ecol., 72, 367-382; 2003), with a plot
function to produce triangular plots.
- - MOStest: Mitchell-Olds--Shaw test (quadratic extreme in an
- interval). The test was suggested for analysis humped species
- richness patterns, but function implements a generalized test,
- which also can be used for analysing location the Gaussian
- optimum in gradient analysis. Moreover, it also can be used to
- replace Tokeshi's test for bimodal species frequency
- patterns. The confidence limits for the location of the optimum
- can be estimated with Fieller theorem or from profile deviance.
-
- mso: Helene Wagner's multiscale ordination or spatial
partitioning of cca and rda. This is taken from the Ecological
Archives with minimal edition with the permission of Helene
Deleted: branches/1.13/man/MOStest.Rd
===================================================================
--- branches/1.13/man/MOStest.Rd 2008-05-18 20:10:06 UTC (rev 359)
+++ branches/1.13/man/MOStest.Rd 2008-05-19 07:30:24 UTC (rev 360)
@@ -1,174 +0,0 @@
-\encoding{UTF-8}
-\name{MOStest}
-\alias{MOStest}
-\alias{print.MOStest}
-\alias{plot.MOStest}
-\alias{fieller.MOStest}
-\alias{profile.MOStest}
-\alias{confint.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. The test was popularized
- in ecology for the analysis of humped species richness patterns
- (Mittelbach et al. 2001), but it is more general. With logarithmic
- link function, the quadratic response defines the Gaussian response
- model of ecological gradients (ter Braak & Looman 1986), and the test
- can be used for inspecting the location of Gaussian optimum within a
- given range of the gradient. It can also be used to replace Tokeshi's
- test of \dQuote{bimodal} species frequency distribution.
-}
-\usage{
-MOStest(x, y, interval, ...)
-\method{plot}{MOStest}(x, which = c(1,2,3,6), ...)
-fieller.MOStest(object, level = 0.95)
-\method{profile}{MOStest}(fitted, alpha = 0.01, maxsteps = 10, del = zmax/5, ...)
-\method{confint}{MOStest}(object, parm = 1, level = 0.95, ...)
-}
-
-\arguments{
- \item{x}{The independent variable or plotting object in \code{plot}. }
- \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}} (minus2). }
- \item{object, fitted}{A result object from \code{MOStest}.}
- \item{level}{The confidence level required.}
- \item{alpha}{Maximum significance level allowed.}
- \item{maxsteps}{Maximum number of steps in the profile.}
- \item{del}{A step length parameter for the profile (see code).}
- \item{parm}{Ignored.}
- \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 \eqn{\mu = b_0 + b_1 x + b_2 x^2}
- with given \code{\link{family}} and link function. If \eqn{b_2 < 0},
- this defines a unimodal curve with highest point at
- \eqn{u = -b_2/(2 b_3)} (ter Braak & Looman 1986). If \eqn{b_2 > 0},
- the parabola has a minimum at \eqn{u} and the response is sometimes
- called \dQuote{bimodal}. The null hypothesis is that the extreme
- point \eqn{u} is located within interval given by points \eqn{p_1} and
- \eqn{p_2}. If the extreme point \eqn{u} is exactly at \eqn{p_1}, then
- \eqn{b_1 = 0} on shifted axis \eqn{x - p_1}. In the test, origin of
- \code{x} is shifted to the values \eqn{p_1} and \eqn{p_2}, and the
- test statistic is the value of the first degree coefficient with its
- significance as estimated by the \code{\link{summary.glm}}
- function(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}}.
-
- Function \code{fieller.MOStest} approximates the confidence limits
- of the location of the extreme point (hump or pit) using Fieller's
- theorem following ter Braak & Looman (1986). The test is based on
- quasideviance except if the \code{\link{family}} is \code{poisson}
- or \code{binomial}. Function \code{profile} evaluates the profile
- deviance of the fitted model, and \code{confint} finds the profile
- based confidence limits following Oksanen et al. (2001).
-
- 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. Further, it can be used to test the location of the
- the Gaussian optimum in ecological gradient analysis (ter Braak &
- Looman 1986, Oksanen et al. 2001).
-}
-
-\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.
-
-ter Braak, C.J.F & Looman, C.W.N 1986. Weighted averaging, logistic
-regression and the Gaussian response model. \emph{Vegetatio} 65,
-3--11.
-}
-\author{Jari Oksanen }
-
-\note{
-Function \code{fieller.MOStest} is based on package \pkg{optgrad} in
-the Ecological Archives
-(\url{http://www.esapubs.org/archive/ecol/E082/015/default.htm})
-accompanying Oksanen et al. (2001). The Ecological Archive package
-\pkg{optgrad} also contains profile deviance method for the location
-of the hump or pit, but the current implementation of \code{profile}
-and \code{confint} rather follow the example of
-\code{\link[MASS]{profile.glm}} and \code{\link[MASS]{confint.glm}} in
-the \pkg{MASS} package.
-}
-
-\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)
-## Infinite confidence limits (NA)
-fieller.MOStest(mod)
-## Finite limits with the profile
-confint(mod)
-plot(profile(mod))
-}
-
-\keyword{ models }
-\keyword{ regression }
More information about the Vegan-commits
mailing list