[Picante-commits] r155 - in branches/gsoc: R man src

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Fri Jul 25 15:17:12 CEST 2008


Author: mrhelmus
Date: 2008-07-25 15:17:12 +0200 (Fri, 25 Jul 2008)
New Revision: 155

Modified:
   branches/gsoc/R/phylostruct.R
   branches/gsoc/R/randomizeSample.R
   branches/gsoc/man/phylostruct.Rd
   branches/gsoc/man/randomizeSample.Rd
   branches/gsoc/src/picante.c
Log:
Edited phylostruct and randomizeSample including the manuals to account for the modified C code. Also cleaned up picante.c.

Modified: branches/gsoc/R/phylostruct.R
===================================================================
--- branches/gsoc/R/phylostruct.R	2008-07-24 23:06:10 UTC (rev 154)
+++ branches/gsoc/R/phylostruct.R	2008-07-25 13:17:12 UTC (rev 155)
@@ -1,40 +1,39 @@
-phylostruct<-function(samp,tree,env=NULL,metric=c("psv","psr","pse","psc","sppregs"),null.model=c("frequency","richness","both"),runs=10,alpha=0.05,fam="binomial"){
+phylostruct<-function(samp,tree,env=NULL,metric=c("psv","psr","pse","psc","sppregs"),null.model=c("frequency","richness","independentswap","trialswap"),runs=100,it=1000,alpha=0.05,fam="binomial"){
 
   metric<-match.arg(metric)
   null.model<-match.arg(null.model)
   if(metric=="sppregs")
   {
-  nulls<-t(replicate(runs,sppregs(randomizeSample(samp,null.model=null.model),env,tree,fam=fam)$correlations))
+  nulls<-t(replicate(runs,sppregs(randomizeSample(samp,null.model=null.model,it=it),env,tree,fam=fam)$correlations))
   obs<-sppregs(samp,env,tree,fam=fam)$correlations
   mean.null<-apply(nulls,2,mean)
   quantiles.null<-t(apply(nulls,2,quantile,probs=c(alpha/2,1-(alpha/2))))
-  
-  return(list(metric=metric,null.model=null.model,runs=runs,obs=obs,mean.null=mean.null
+  if((null.model!="independentswap")&&(null.model!="trialswap")){it=NA}
+  return(list(metric=metric,null.model=null.model,runs=runs,it=it,obs=obs,mean.null=mean.null
                 ,quantiles.null=quantiles.null,phylo.structure=NULL,nulls=nulls))
 
-
   } else {
 
     nulls<-switch(metric,
-                       psv = replicate(runs,mean(psv(as.matrix(randomizeSample(samp,null.model=null.model)),tree,compute.var=FALSE)[,1])),
-                       psr = replicate(runs,mean(psr(as.matrix(randomizeSample(samp,null.model=null.model)),tree,compute.var=FALSE)[,1])),
-                       pse = replicate(runs,mean(pse(as.matrix(randomizeSample(samp,null.model=null.model)),tree)[,1])),
-                       psc = replicate(runs,mean(psc(as.matrix(randomizeSample(samp,null.model=null.model)),tree)[,1])))
+                       psv = replicate(runs,mean(psv(as.matrix(randomizeSample(samp,null.model=null.model,it=it)),tree,compute.var=FALSE)[,1],na.rm=TRUE)),
+                       psr = replicate(runs,mean(psr(as.matrix(randomizeSample(samp,null.model=null.model,it=it)),tree,compute.var=FALSE)[,1],na.rm=TRUE)),
+                       pse = replicate(runs,mean(pse(as.matrix(randomizeSample(samp,null.model=null.model,it=it)),tree)[,1],na.rm=TRUE)),
+                       psc = replicate(runs,mean(psc(as.matrix(randomizeSample(samp,null.model=null.model,it=it)),tree)[,1],na.rm=TRUE)))
     quantiles.null<-quantile(nulls,probs=c(alpha/2,1-(alpha/2)))
     mean.null<-mean(nulls)
     mean.obs<-switch(metric,
-                       psv = mean(psv(samp,tree,compute.var=FALSE)[,1]),
-                       psr = mean(psr(samp,tree,compute.var=FALSE)[,1]),
-                       pse = mean(pse(samp,tree)[,1]),
-                       psc = mean(psc(samp,tree)[,1]))
+                       psv = mean(psv(samp,tree,compute.var=FALSE)[,1],na.rm=TRUE),
+                       psr = mean(psr(samp,tree,compute.var=FALSE)[,1],na.rm=TRUE),
+                       pse = mean(pse(samp,tree)[,1],na.rm=TRUE),
+                       psc = mean(psc(samp,tree)[,1],na.rm=TRUE))
 
     if(mean.obs<=quantiles.null[1])
     {phylo.structure="underdispersed"
     } else {if(mean.obs>=quantiles.null[2]){
     phylo.structure="overdispersed"} else {phylo.structure="random"}
     }
-    
-    return(list(metric=metric,null.model=null.model,runs=runs,mean.obs=mean.obs,mean.null=mean.null
+    if((null.model!="independentswap")&&(null.model!="trialswap")){it=NA}    
+    return(list(metric=metric,null.model=null.model,runs=runs,it=it,mean.obs=mean.obs,mean.null=mean.null
                 ,quantiles.null=quantiles.null,phylo.structure=phylo.structure,null.means=nulls))
   }
 }

Modified: branches/gsoc/R/randomizeSample.R
===================================================================
--- branches/gsoc/R/randomizeSample.R	2008-07-24 23:06:10 UTC (rev 154)
+++ branches/gsoc/R/randomizeSample.R	2008-07-25 13:17:12 UTC (rev 155)
@@ -3,7 +3,7 @@
 } 
 
 `randomizeSample` <-
-function(samp, null.model=c("frequency","richness","richnessindex","both","trialswap"), lang=c("C","R"), it=1000) {
+function(samp, null.model=c("frequency","richness","independentswap","trialswap"), lang=c("C","R"), it=1000) {
 
     samp <- as.matrix(samp)
     lang <- match.arg(lang)
@@ -27,12 +27,7 @@
     	}
 	}
 
-	if (identical(null.model,"richnessindex")) {
-	    ret1 <- .C("richnessindex", m=as.numeric(samp), as.integer(nrow(samp)), as.integer(ncol(samp)), PACKAGE="picante")
-        return(matrix(ret1$randm,nrow=nrow(samp),dimnames=list(rownames(samp),colnames(samp))))
-	}
-
-	if (identical(null.model,"both")) 
+	if (identical(null.model,"independentswap")) 
     {
         if (identical(lang,"R")) {
             #check for presence-absence and warn until abundance implemented

Modified: branches/gsoc/man/phylostruct.Rd
===================================================================
--- branches/gsoc/man/phylostruct.Rd	2008-07-24 23:06:10 UTC (rev 154)
+++ branches/gsoc/man/phylostruct.Rd	2008-07-25 13:17:12 UTC (rev 155)
@@ -4,8 +4,8 @@
 \description{ Randomize sample/community data matrices to create null distributions of given metrics
 }
 \usage{
-phylostruct(samp, tree, env=NULL, metric=c("psv","psr","pse","psc","sppregs"), null.model=c("frequency","richness","both"),
-            runs=10, alpha=0.05, fam="binomial")
+phylostruct(samp, tree, env=NULL, metric=c("psv","psr","pse","psc","sppregs"), null.model=c("frequency",
+            "richness","independentswap","trialswap"), runs=100, it=1000, alpha=0.05, fam="binomial")
 }
 \arguments{
   \item{samp}{ community data matrix, species as columns, communities as rows }
@@ -14,6 +14,7 @@
   \item{metric}{ if \code{metric="psv"}, \code{"psr"}, \code{"pse"}, or \code{"psc"} compares the observed mean of the respective metric to a null distribution at a given alpha; if \code{metric="sppregs"} compares the three correlations produced by \code{\link{sppregs}} to null distributions }
   \item{null.model}{ permutation procedure used to create the null distribution, see \code{\link{randomizeSample}}}
   \item{runs}{ the number of permutations to create the distribution, a rule of thumb is (number of communities)/alpha  }
+  \item{it}{ the number of swaps for the independent and trial-swap null models, see \code{\link{randomizeSample}} }
   \item{alpha}{ probability value to compare the observed mean/correlations to a null distribution }
   \item{fam}{ as in \code{\link{sppregs}} }
 }
@@ -24,6 +25,7 @@
 \item{metric}{ metric used }
 \item{null.model}{ permutation used }
 \item{runs}{ number of permutations }
+\item{it}{ number of swaps if applicable }
 \item{obs}{ observed mean value of a particular metric or the three observed correlations from \code{\link{sppregs}} }
 \item{mean.null}{ mean(s) of the null distribution(s) }
 \item{quantiles.null}{  quantiles of the null distribution(s) to compare to \code{obs}; determined by \code{alpha}}
@@ -42,5 +44,4 @@
 \author{ Matthew Helmus \email{mrhelmus at gmail.com} }
 \seealso{ \code{\link{psd}} ,\code{\link{sppregs}}, \code{\link{randomizeSample}} }
 
-\keyword{univar}
-
+\keyword{univar}
\ No newline at end of file

Modified: branches/gsoc/man/randomizeSample.Rd
===================================================================
--- branches/gsoc/man/randomizeSample.Rd	2008-07-24 23:06:10 UTC (rev 154)
+++ branches/gsoc/man/randomizeSample.Rd	2008-07-25 13:17:12 UTC (rev 155)
@@ -6,7 +6,7 @@
   Various null models for randomizing community data matrices
 }
 \usage{
-randomizeSample(samp, null.model = c("frequency", "richness", "both"))
+randomizeSample(samp, null.model = c("frequency", "richness", "independentswap","trialswap"), lang=c("C","R"), it=1000)
 }
 
 \arguments{
@@ -14,15 +14,20 @@
   \item{null.model}{ Null model
     \item{frequency}{Randomize community data matrix abundances within species (maintains species occurence frequency)}
     \item{richness}{Randomize community data matrix abundances within samples (maintains sample species richness)}
-    \item{both}{Randomize community data matrix by drawing species from sample pool with probability weighted by occurrence frequency (maintains both species occurrence frequency and sample species richness)}
+    \item{independentswap}{Randomize community data matrix with the independent-swap algorithm (Gotelli 2000) by drawing species from sample pool with probability weighted by occurrence frequency (maintains both species occurrence frequency and sample species richness)}
+    \item{trialswap}{Randomize community data matrix with the trial-swap algorithm (Miklos & Podani 2004) by drawing species from sample pool with probability weighted by occurrence frequency (maintains both species occurrence frequency and sample species richness)}
   }
+  \item{lang}{run code in C or R? C is faster.} 
 }
 \value{
   Randomized community data matrix
 }
-\references{Gotelli, N.J. 2000. Null model analysis of species co-occurrence patterns. Ecology 81: 2606-2621.}
+\references{Gotelli, N.J. 2000. Null model analysis of species co-occurrence patterns. Ecology 81: 2606-2621.
+\cr
+\cr
+Miklos I. & Podani J. 2004. Randomization of presence-absence matrices: Comments and new algorithms. Ecology 85: 86-92.
+}
 \author{ Steve Kembel <skembel at berkeley.edu> }
-\section{Warning }{ Null model \code{both} currently only works with presence-absence data. Convert your data to presence-absence before using this null model (e.g. \code{decostand(x,method="pa")})}
 \examples{
 data(phylocom)
 randomizeSample(phylocom$sample, null.model="richness")}

Modified: branches/gsoc/src/picante.c
===================================================================
--- branches/gsoc/src/picante.c	2008-07-24 23:06:10 UTC (rev 154)
+++ branches/gsoc/src/picante.c	2008-07-25 13:17:12 UTC (rev 155)
@@ -84,80 +84,12 @@
     mattovec(v,m,row,column);
     PutRNGstate();
 }
+
 
-void richnessindex(double *v, int *prow, int *pcolumn)
-{
-    int i, j, k, q, p, b, sum=0, we;
-    int row, column;
-    double **m, **randm;
-    row = *prow;
-    column = *pcolumn;
-    q = column;
-    double *randrow[column];
-
-    m = vectomat(v,row,column);//Change the vector v passed from R back to a Matrix
-
-    /*make list of ONES THE LENGTH OF THE COLUMNS*/
-    int s[column];
-    int t;
-    for (t=0;t<column;t++)
-    {
-        s[t]=1;
-    }
-    //////////////////////////////////
-
-    randm = (double **)R_alloc(row,sizeof(double *)); //Initiate the Random Matrix
-    //randrow = (double *)R_alloc(column,sizeof(double)); //Initiate a vector to hold the community that is being randomized
-
-    /*Begin looping to create random matrix*/
-    for (i=0; i<row; i++)
-    {
-        randm[i] = (double *)R_alloc(column,sizeof(double)); //Initiate the Random community at row i
-
-        for (p=0;p<column;p++) //Get obs community i
-        {
-            randrow[i][p]=m[i][p];
-        }
-
-        for (j=0; j<column-1; j++) // randomly sample one species at a time across all rows of community i
-        {
-            for(b=0; b<q; b++)// Sum vector s to give the number to be passed to intrand
-                {
-                sum = sum + s[i];
-                }
-            k=intrand(sum);
-            randm[i][j]=randrow[i][k];//Assign element randm i,j a random element from the community
-            double *randrow[q-1];//Reinitialize randrow THIS WILL PROB BE A BUG IN THE CODE
-
-            /*Fill randrow with all the elements of community i that have not been sampled yet
-            A bug may also come up here if q=k...*/
-            we=0;
-            while(we<q-1)
-                {
-                for (p=0;p<q-k;p++)
-                    {
-                    randrow[i][we]=m[i][p];
-                    we++;
-                    }
-                for (p=k+1;p<q;p++)
-                    {
-                    randrow[i][we]=m[i][p];
-                    we++;
-                    }
-                }
-            q=q-1;
-            }
-        }
-
-    mattovec(v,randm,row,column);
-
-    PutRNGstate();
-}
-
 void richness(double *v, int *prow, int *pcolumn)
 {
 
-    int i,j,k,l;
+    int i,j,k;
     int row, column;
 	double tmp;
     double **m;
@@ -172,7 +104,7 @@
     {
         for (j=0;j<column;j++)
         {
-            while((k=intrand(column))==j);//choose another column (species) at random
+            k=intrand(column);//choose another column (species) at random
             tmp = m[i][j];
             m[i][j] = m[i][k];
             m[i][k] = tmp;
@@ -185,7 +117,7 @@
 void frequency(double *v, int *prow, int *pcolumn)
 {
 
-    int i,j,k,l;
+    int i,j,k;
     int row, column;
 	double tmp;
     double **m;
@@ -200,7 +132,7 @@
     {
         for (j=0;j<row;j++)
         {
-            while((k=intrand(row))==j);//choose another row (sample) at random
+            k=intrand(row);//choose another row (sample) at random
             tmp = m[j][i];
             m[j][i] = m[k][i];
             m[k][i] = tmp;



More information about the Picante-commits mailing list