[Picante-commits] r85 - branches/gsoc/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jun 4 09:56:08 CEST 2008
Author: mrhelmus
Date: 2008-06-04 09:56:08 +0200 (Wed, 04 Jun 2008)
New Revision: 85
Modified:
branches/gsoc/R/phylodiversity.R
Log:
Fixed Bug: PSV, PSR now accept data of one community
Added: PSE, PSC (similar to NTI), and spp.PSVcalc a function to estimate the effect of each species on the mean value of PSV for a set of communities. This completes all of the stats for the Helmus 2007 AmNat paper.
Modified: branches/gsoc/R/phylodiversity.R
===================================================================
--- branches/gsoc/R/phylodiversity.R 2008-06-02 18:58:53 UTC (rev 84)
+++ branches/gsoc/R/phylodiversity.R 2008-06-04 07:56:08 UTC (rev 85)
@@ -114,31 +114,209 @@
PSVcalc<-function(samp,tree,compute.var=TRUE){
+ # Make samp matrix a pa matrix
+ samp[samp>0]<-1
+ flag=0
+ if (is.null(dim(samp))) #if the samp matrix only has one site
+ {
+ samp<-rbind(samp,samp)
+ flag=2
+ }
-# Make samp matrix a pa matrix
-samp<-as.matrix(samp)
-samp[samp>0]<-1
+ if(is(tree)[1]=="phylo")
+ {
+ # Make sure that the species line up
+ samp<-samp[,tree$tip.label]
+ # Make a correlation matrix of the species pool phylogeny
+ Cmatrix<-vcv.phylo(tree,cor=T)
+ } else {
+ Cmatrix<-tree
+ samp<-samp[,colnames(Cmatrix)]
+ }
+
+ # reduce given Cmatrix to the species observed in samp
+ SR<-rowSums(samp)
+ samp<-samp[SR>1,] # prune out locations with <2 species
+ SR<-SR[SR>1]
+
+ #cull the species that are not found in the samp set after all communities with 1 species are removed
+ preval<-colSums(samp)/sum(samp)
+ indexcov<-preval>0
+ Cmatrix<-Cmatrix[indexcov,indexcov]
+ samp<-samp[,indexcov]
+
+ # numbers of locations and species
+ nlocations<-dim(samp)[1]
+ nspecies<-dim(samp)[2]
+
+ ##################################
+ # calculate observed PSVs
+ #
+ PSVs<-NULL
+
+ for(i in 1:nlocations)
+ {
+ index<-seq(1:nrow(Cmatrix))[samp[i,]==1] #species present
+ n<-length(index) #number of species present
+ C<-Cmatrix[index,index] #C for individual locations
+ PSV<-(n*sum(diag(as.matrix(C)))-sum(C))/(n*(n-1))
+ PSVs<-c(PSVs,PSV)
+ }
+ PSVout<-cbind(PSVs,SR)
+
+ if(flag==2)
+ {
+ PSVout<-PSVout[-2,]
+ return(PSVout)
+ } else {
+ if(compute.var==FALSE)
+ {
+ return(data.frame(PSVout))
+ } else {
+
+ PSVvar<-NULL
+ X<-Cmatrix-(sum(sum(Cmatrix-diag(nspecies))))/(nspecies*(nspecies-1));
+ X<-X-diag(diag(X))
+
+ SS1<-sum(X*X)/2
+
+ SS2<-0;
+ for(i in 1:(nspecies-1)){
+ for(j in (i+1):nspecies){
+ SS2<-SS2+X[i,j]*(sum(X[i,])-X[i,j])
+ }
+ }
+ SS3<--SS1-SS2
+
+ S1<-SS1*2/(nspecies*(nspecies-1))
+ S2<-SS2*2/(nspecies*(nspecies-1)*(nspecies-2))
+
+ if(nspecies==3){
+ S3<-0
+ }else{
+ S3<-SS3*2/(nspecies*(nspecies-1)*(nspecies-2)*(nspecies-3))
+ }
+
+ for(n in 2:nspecies){
+ approxi<-2/(n*(n-1))*(S1 + (n-2)*S2 + (n-2)*(n-3)*S3)
+ PSVvar<-rbind(PSVvar, c(n, approxi))
+ }
+
+ vars<-rep(0,nlocations)
+ PSVout<-cbind(PSVout,vars)
+
+ for (g in 1:nlocations)
+ {
+ PSVout[g,3]<-PSVvar[PSVout[g,2]-1,2]
+ }
+ return(data.frame(PSVout))
+ }
+ }
+}
-if(is(tree)[1]=="phylo")
-{
- # Make sure that the species line up
- samp<-samp[,tree$tip.label]
- # Make a correlation matrix of the species pool phylogeny
- Cmatrix<-vcv.phylo(tree,cor=T)
-} else {
- Cmatrix<-tree
- samp<-samp[,colnames(Cmatrix)]
+
+PSRcalc <- function(samp,tree,compute.var=TRUE){
+ PSVout<-PSVcalc(samp,tree,compute.var)
+ if(length(PSVout)==2)
+ {
+ PSRout<-data.frame(cbind(PSVout[1]*PSVout[2],PSVout[2]))
+ names(PSRout)<-c("PSR","SR")
+ rownames(PSRout)<-""
+ return(PSRout)
+ } else {
+ PSRout<-cbind(PSVout[,1]*PSVout[,2],PSVout[,2])
+ if(compute.var!=TRUE) {
+ colnames(PSRout)<-c("PSR","SR")
+ return(data.frame(PSRout))
+ } else {
+ PSRout<-cbind(PSRout,PSVout[,3]*(PSVout[,2])^2)
+ colnames(PSRout)<-c("PSR","SR","vars")
+ return(data.frame(PSRout))
+ }
+ }
}
-if (is.null(dim(samp))) #if the samp matrix only has one site
-{
- n<-sum(samp)
- C<-Cmatrix
- PSV<-(n*sum(diag(as.matrix(C)))-sum(C))/(n*(n-1))
- PSVout<-c(PSV,n)
-} else {
+PSEcalc<-function(samp,tree){
+ flag=0
+ if (is.null(dim(samp))) #if the samp matrix only has one site
+ {
+ samp<-rbind(samp,samp)
+ flag=2
+ }
+ if(is(tree)[1]=="phylo")
+ {
+ # Make sure that the species line up
+ samp<-samp[,tree$tip.label]
+ # Make a correlation matrix of the species pool phylogeny
+ Cmatrix<-vcv.phylo(tree,cor=T)
+ } else {
+ Cmatrix<-tree
+ samp<-samp[,colnames(Cmatrix)]
+ }
+
# reduce given Cmatrix to the species observed in samp
+ SR<-apply(samp>0,1,sum)
+ # prune out locations with <2 species
+ samp<-samp[SR>1,]
+ SR<-SR[SR>1]
+
+ #cull the species that are not found in the samp set after all communities with 1 species are removed
+ preval<-colSums(samp)/sum(samp)
+ indexcov<-preval>0
+ Cmatrix<-Cmatrix[indexcov,indexcov]
+ samp<-samp[,indexcov]
+
+ # numbers of locations and species
+ nlocations<-dim(samp)[1]
+ nspecies<-dim(samp)[2]
+
+ #################################
+ # calculate observed phylogenetic species evenness
+ PSEs<-NULL
+ for(i in 1:nlocations){
+ index<-seq(1,ncol(Cmatrix))[samp[i,]>0] # species present
+ C<-Cmatrix[index,index] #C for individual locations
+ n<-length(index) #location species richness
+ N<-sum(samp[i,]) #location total abundance
+ M<-samp[i,samp[i,]>0] #species abundance column
+ mbar<-mean(M) #mean species abundance
+ PSE<-(N*t(diag(as.matrix(C)))%*%M-t(M)%*%as.matrix(C)%*%M)/(N^2-N*mbar) #phylogenetic species evenness
+ PSEs<-c(PSEs,PSE)
+ }
+ PSEout=cbind(PSEs,SR)
+ if (flag==2)
+ {
+ PSEout<-PSEout[-2,]
+ return(PSEout)
+ } else {
+ return(data.frame(PSEout))
+ }
+}
+
+
+PSCcalc<-function(samp,tree,compute.var=TRUE){
+ # Make samp matrix a pa matrix
+ samp[samp>0]<-1
+ flag=0
+ if (is.null(dim(samp))) #if the samp matrix only has one site
+ {
+ samp<-rbind(samp,samp)
+ flag=2
+ }
+
+ if(is(tree)[1]=="phylo")
+ {
+ # Make sure that the species line up
+ samp<-samp[,tree$tip.label]
+ # Make a correlation matrix of the species pool phylogeny
+ Cmatrix<-vcv.phylo(tree,cor=T)
+ } else {
+ Cmatrix<-tree
+ samp<-samp[,colnames(Cmatrix)]
+ }
+
+ # reduce given Cmatrix to the species observed in samp
SR<-rowSums(samp)
samp<-samp[SR>1,] # prune out locations with <2 species
SR<-SR[SR>1]
@@ -154,73 +332,71 @@
nspecies<-dim(samp)[2]
##################################
- # calculate observed PSVs
+ # calculate observed PSCs
#
- PSVs<-NULL
+ PSCs<-NULL
for(i in 1:nlocations)
{
index<-seq(1:nrow(Cmatrix))[samp[i,]==1] #species present
n<-length(index) #number of species present
C<-Cmatrix[index,index] #C for individual locations
- PSV<-(n*sum(diag(as.matrix(C)))-sum(C))/(n*(n-1))
- PSVs<-c(PSVs,PSV)
+ diag(C)<--1
+ PSC<-sum(apply(C,1,max))/n
+ PSCs<-c(PSCs,PSC)
}
- PSVout<-cbind(PSVs,SR)
+ PSCout<-cbind(PSCs,SR)
+ if (flag==2)
+ {
+ PSCout<-PSCout[-2,]
+ return(PSCout)
+ } else {
+ return(data.frame(PSCout))
+ }
+
}
-if(compute.var==FALSE) {
- return(data.frame(PSVout))
-} else {
-
- PSVvar<-NULL
- X<-Cmatrix-(sum(sum(Cmatrix-diag(nspecies))))/(nspecies*(nspecies-1));
- X<-X-diag(diag(X))
-
- SS1<-sum(X*X)/2
-
- SS2<-0;
- for(i in 1:(nspecies-1)){
- for(j in (i+1):nspecies){
- SS2<-SS2+X[i,j]*(sum(X[i,])-X[i,j])
- }
+spp.PSVcalc<-function(samp,tree){
+ # Make samp matrix a pa matrix
+ samp[samp>0]<-1
+ if (is.null(dim(samp))) #if the samp matrix only has one site
+ {
+ samp<-rbind(samp,samp)
+ }
+ if(is(tree)[1]=="phylo")
+ {
+ # Make sure that the species line up
+ samp<-samp[,tree$tip.label]
+ # Make a correlation matrix of the species pool phylogeny
+ Cmatrix<-vcv.phylo(tree,cor=T)
+ } else {
+ Cmatrix<-tree
+ samp<-samp[,colnames(Cmatrix)]
}
- SS3<--SS1-SS2
+ # reduce given Cmatrix to the species observed in samp
+ samp<-samp[rowSums(samp)>1,] # prune out locations with <2 species
- S1<-SS1*2/(nspecies*(nspecies-1))
- S2<-SS2*2/(nspecies*(nspecies-1)*(nspecies-2))
+ #cull the species that are not found in the samp set after all communities with 1 species are removed
+ preval<-colSums(samp)/sum(samp)
+ indexcov<-preval>0
+ Cmatrix<-Cmatrix[indexcov,indexcov]
+ samp<-samp[,indexcov]
- if(nspecies==3){
- S3<-0
- }else{
- S3<-SS3*2/(nspecies*(nspecies-1)*(nspecies-2)*(nspecies-3))
- }
+ obs.PSV<-mean(PSVcalc(samp,Cmatrix,compute.var=FALSE)[,1])
- for(n in 2:nspecies){
- approxi<-2/(n*(n-1))*(S1 + (n-2)*S2 + (n-2)*(n-3)*S3)
- PSVvar<-rbind(PSVvar, c(n, approxi))
- }
+ # numbers of locations and species
+ nlocations<-dim(samp)[1]
+ nspecies<-dim(samp)[2]
- vars<-rep(0,nlocations)
- PSVout<-cbind(PSVout,vars)
-
- for (g in 1:nlocations)
+ spp.PSVs<-NULL
+ for (j in 1:nspecies)
{
- PSVout[g,3]<-PSVvar[PSVout[g,2]-1,2]
+ spp.samp<-samp[,-j]
+ spp.Cmatrix<-Cmatrix[-j,-j]
+ spp.PSV<-mean(PSVcalc(spp.samp,spp.Cmatrix,compute.var=FALSE)[,1])
+ spp.PSVs<-c(spp.PSVs,spp.PSV)
}
- return(data.frame(PSVout))
- }
-}
-
-PSRcalc <- function(samp,tree,compute.var=TRUE){
- PSVout<-PSVcalc(samp,tree,compute.var)
- PSRout<-cbind(PSVout[,1]*PSVout[,2],PSVout[,2])
- if(compute.var!=TRUE) {
- colnames(PSRout)<-c("PSR","SR")
- return(data.frame(PSRout))
- } else {
- PSRout<-cbind(PSRout,PSVout[,3]*(PSVout[,2])^2)
- colnames(PSRout)<-c("PSR","SR","vars")
- return(data.frame(PSRout))
- }
+ spp.PSVout<-(spp.PSVs-obs.PSV)/sum(abs(spp.PSVs-obs.PSV))
+ names(spp.PSVout)<-colnames(samp)
+ return(spp.PSVout)
}
\ No newline at end of file
More information about the Picante-commits
mailing list