[Vegan-commits] r299 - in pkg: . R inst man

[noreply at r-forge.r-project.org]

Author: gsimpson
Date: 2008-03-30 18:39:52 +0200 (Sun, 30 Mar 2008)
New Revision: 299

Added:
   pkg/R/permute.R
Modified:
   pkg/DESCRIPTION
   pkg/inst/ChangeLog
   pkg/man/permuted.index2.Rd
Log:
New function permute() and example of new API in ?permute

Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION	2008-03-30 15:12:05 UTC (rev 298)
+++ pkg/DESCRIPTION	2008-03-30 16:39:52 UTC (rev 299)
@@ -1,7 +1,7 @@
 Package: vegan
 Title: Community Ecology Package
-Version: 1.12-7
-Date: Mar 25, 2008
+Version: 1.12-8
+Date: Mar 30, 2008
 Author: Jari Oksanen, Roeland Kindt, Pierre Legendre, Bob O'Hara, Gavin L. Simpson, 
   M. Henry H. Stevens  
 Maintainer: Jari Oksanen <jari.oksanen at oulu.fi>

Added: pkg/R/permute.R
===================================================================
--- pkg/R/permute.R	                        (rev 0)
+++ pkg/R/permute.R	2008-03-30 16:39:52 UTC (rev 299)
@@ -0,0 +1,10 @@
+permute <- function(i, n, control) {
+    if(control$complete && !is.null(control$all.perms))
+        perm <- control$all.perms[i,]
+    else {
+        if(control$complete)
+            warning("'$all.perms' is NULL, yet '$complete = TRUE'.\nReturning a random permutation.")
+        perm <- permuted.index2(n, control)
+    }
+    return(perm)
+}

Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog	2008-03-30 15:12:05 UTC (rev 298)
+++ pkg/inst/ChangeLog	2008-03-30 16:39:52 UTC (rev 299)
@@ -2,8 +2,15 @@
 
 VEGAN DEVEL VERSIONS at http://r-forge.r-project.org/
 
-Version 1.12-7 (opened Mar 25, 2008)
+Version 1.12-8 (opened Mar 30, 2008)
 
+	* permute: New high-level untility function for facilitating
+	the production of permutation tests using the new permutation
+	designs allowed by permuted.index2(). An example of the new
+	API is shown in ?permute
+	
+Version 1.12-7 (closed Mar 30, 2008)
+
 	* permuplot: New function, produces a graphical representation
 	of a permutation design given a number of observations and a
 	object returned by permControl(). The function handles all the

Modified: pkg/man/permuted.index2.Rd
===================================================================
--- pkg/man/permuted.index2.Rd	2008-03-30 15:12:05 UTC (rev 298)
+++ pkg/man/permuted.index2.Rd	2008-03-30 16:39:52 UTC (rev 299)
@@ -2,6 +2,7 @@
 \alias{permuted.index2}
 \alias{permControl}
 \alias{print.permControl}
+\alias{permute}
 
 \title{Unrestricted and restricted permutations}
 \description{
@@ -17,6 +18,8 @@
             mirror = FALSE, constant = FALSE,
             ncol = NULL, nrow = NULL,
             all.perms = NULL)
+
+permute(i, n, control)
 }
 
 \arguments{
@@ -45,6 +48,7 @@
     in the spatial grid respectiavly.}
   \item{all.perms}{an object of class \code{allPerms}, the result of a
     call to \code{\link{allPerms}}.}
+  \item{i}{integer; row of \code{control$all.perms} to return.}
 }
 \details{
   \code{permuted.index2} can generate permutations for a wide range of
@@ -71,6 +75,13 @@
   argument \code{constant}. If \code{constant = TRUE} then the same
   permutation will be generated for each level of \code{strata}. The
   default is \code{constant = FALSE}.
+
+  \code{permute} is a higher level utility function for use in a loop
+  within a function implementing a permutation test. The main purpose of
+  \code{permute} is to return the correct permutation in each iteration
+  of the loop, either a random permutation from the current design or
+  the next permutation from \code{control$all.perms} if it is not
+  \code{NULL} and \code{control$complete} is \code{TRUE}.
 }
 \value{
   For \code{permuted.index2} a vector of length \code{n} containing a
@@ -134,6 +145,59 @@
 ## permuting levels of block instead of observations
 permuted.index2(20, permControl(strata = block, type = "strata"))
 
+## Simple function using permute() to assess significance
+## of a t.test  
+pt.test <- function(x, group, control) {
+    ## function to calculate t
+    t.statistic <- function(x, y) {
+        m <- length(x)
+        n <- length(y)
+        xbar <- sum(x) / m
+        ybar <- sum(y) / n
+        xvar <- .Internal(cov(x, NULL, 1, FALSE))
+        yvar <- .Internal(cov(y, NULL, 1, FALSE))
+        pooled <- sqrt(((m-1)*xvar + (n-1)*yvar) / (m+n-2))
+        (xbar - ybar) / (pooled * sqrt(1/m + 1/n))
+    }
+    ## check the control object
+    control <- permCheck(x, control)$control
+    ## number of observations
+    nobs <- getNumObs(x)
+    ## group names
+    lev <- names(table(group))
+    ## vector to hold results, +1 because of observed t
+    t.permu <- numeric(length = control$nperm) + 1
+    ## calculate observed t
+    t.permu[1] <- t.statistic(x[group == lev[1]], x[group == lev[2]])
+    ## generate randomisation distribution of t
+    for(i in seq_along(t.permu)) {
+        ## return a permutation
+        want <- permute(i, nobs, control)
+        ## calculate permuted t
+        t.permu[i+1] <- t.statistic(x[want][group == lev[1]],
+                                    x[want][group == lev[2]])
+    }
+    ## pval from permutation test
+    pval <- sum(abs(t.permu) >= abs(t.permu[1])) / (control$nperm + 1)
+    ## return value
+    return(list(t.stat = t.permu[1], pval = pval))
 }
+
+## generate some data with slightly different means
+set.seed(1234)
+gr1 <- rnorm(20, mean = 9)
+gr2 <- rnorm(20, mean = 10)
+dat <- c(gr1, gr2)
+## grouping variable
+grp <- gl(2, 20, labels = paste("Group", 1:2))
+## create the permutation design
+control <- permControl(type = "free", nperm = 999)
+## perform permutation t test
+perm.val <- pt.test(dat, grp, control)
+perm.val
+
+## compare perm.val with the p-value from t.test()
+t.test(dat ~ grp, var.equal = TRUE)
+}
 \keyword{ htest }
 \keyword{ design }