[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