[Picante-commits] r73 - branches/gsoc/R

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

Author: mrhelmus
Date: 2008-05-28 02:14:18 +0200 (Wed, 28 May 2008)
New Revision: 73

Modified:
   branches/gsoc/R/phylodiversity.R
Log:
Added working version of PSV to phylodiversity.R

Modified: branches/gsoc/R/phylodiversity.R
===================================================================
--- branches/gsoc/R/phylodiversity.R	2008-05-25 22:39:48 UTC (rev 72)
+++ branches/gsoc/R/phylodiversity.R	2008-05-28 00:14:18 UTC (rev 73)
@@ -111,3 +111,103 @@
     data.frame(ntaxa=specnumber(samp),mnnd.obs, mnnd.rand.mean, mnnd.rand.sd, mnnd.obs.rank, 
         mnnd.obs.z, mnnd.obs.p=mnnd.obs.rank/(runs+1),runs=runs, row.names = row.names(samp))
 }
+
+
+PSVcalc<-function(samp,tree,compute.var=TRUE){
+
+# 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)]
+}
+
+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 {
+
+  # 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(compute.var==FALSE) {
+  return(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(PSVout)
+  }
+}