[Vegan-commits] r348 - in pkg: R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon May 12 20:57:14 CEST 2008
Author: jarioksa
Date: 2008-05-12 20:57:14 +0200 (Mon, 12 May 2008)
New Revision: 348
Added:
pkg/R/confint.MOStest.R
pkg/R/profile.MOStest.R
Modified:
pkg/inst/ChangeLog
pkg/man/MOStest.Rd
Log:
MOStest got profile and confint methods
Added: pkg/R/confint.MOStest.R
===================================================================
--- pkg/R/confint.MOStest.R (rev 0)
+++ pkg/R/confint.MOStest.R 2008-05-12 18:57:14 UTC (rev 348)
@@ -0,0 +1,6 @@
+`confint.MOStest` <-
+ function (object, parm = 1, level = 0.95, ...)
+{
+ require(MASS) || stop("requires packages MASS")
+ confint(profile(object), level = level, ...)
+}
Added: pkg/R/profile.MOStest.R
===================================================================
--- pkg/R/profile.MOStest.R (rev 0)
+++ pkg/R/profile.MOStest.R 2008-05-12 18:57:14 UTC (rev 348)
@@ -0,0 +1,55 @@
+`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: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2008-05-11 13:49:40 UTC (rev 347)
+++ pkg/inst/ChangeLog 2008-05-12 18:57:14 UTC (rev 348)
@@ -21,12 +21,12 @@
* MOStest: added fieller.MOStest for approximate confidence
intervals of the location of the hump or pit. The method is based
on Fieller's theorem following ter Braak & Looman (Vegatatio 65,
- 3-11; 1986) and based on the code from the ESA Ecological Archives
- accompanying Oksanen et al. (Ecology 82, 1191-1197; 2001) and
- published in package optgrad
- (http://www.esapubs.org/archive/ecol/E082/015/default.htm). The
- optgrad package also has profile method, but it is not yet
- implemented.
+ 3-11; 1986) and profile based condifence limits following Oksanen
+ et al. (Ecology 82, 1191-1197; 2001). Both are based on the code
+ from the ESA Ecological Archives accompanying Oksanen et
+ al. (Ecology 82, 1191-1197; 2001) and published in package optgrad
+ (http://www.esapubs.org/archive/ecol/E082/015/default.htm), but
+ profile methods heavily borrow from the MASS package.
Version 1.12-14 (closed May 9, 2008)
Modified: pkg/man/MOStest.Rd
===================================================================
--- pkg/man/MOStest.Rd 2008-05-11 13:49:40 UTC (rev 347)
+++ pkg/man/MOStest.Rd 2008-05-12 18:57:14 UTC (rev 348)
@@ -4,8 +4,9 @@
\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)
@@ -19,6 +20,8 @@
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{
@@ -29,8 +32,12 @@
\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}{A result object from \code{MOStest}.}
+ \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
@@ -68,7 +75,9 @@
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}.
+ 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
@@ -117,7 +126,10 @@
(\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 it is not (yet) implemented in \pkg{vegan}.
+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}}. }
@@ -138,6 +150,9 @@
par(op)
## Infinite confidence limits (NA)
fieller.MOStest(mod)
+## Finite limits with the profile
+confint(mod)
+plot(profile(mod))
}
\keyword{ models }
More information about the Vegan-commits
mailing list