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

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

Author: mrhelmus
Date: 2008-06-21 04:52:24 +0200 (Sat, 21 Jun 2008)
New Revision: 112

Modified:
   branches/gsoc/R/sppregs.r
Log:
Added tree acceptance and fixed bug of one env variable

Modified: branches/gsoc/R/sppregs.r
===================================================================
--- branches/gsoc/R/sppregs.r	2008-06-20 12:42:10 UTC (rev 111)
+++ branches/gsoc/R/sppregs.r	2008-06-21 02:52:24 UTC (rev 112)
@@ -9,26 +9,25 @@
   {
     nenvs<-1
     envnames<-names(env)
-    formu<-paste("y~",env.names[1]) #Make formula
+    if (is.null(envnames)){envnames<-"env"}
+    formu<-paste("y~",envnames) #Make formula
   } else {
     nenvs<-dim(env)[2]       # number of environmental variables
     envnames<-colnames(env)  # vecor of env names
 
-    formu<-paste("y~",env.names[1]) #Make formula
+    formu<-paste("y~",envnames[1]) #Make formula
     for (t in 2:nenvs)
     {
-      formu<-paste(formu,env.names[t],sep="+")
+      formu<-paste(formu,envnames[t],sep="+")
     }
   }
 
-
-
   #Fit either the logistic or the standard regressions
   spp.resids<-NULL
   spp.coef<-NULL
   spp.se<-NULL
   spp.fits<-NULL
-
+               
   if(fam=="gaussian")
   {
       for(i in 1:nspp)
@@ -48,7 +47,7 @@
     for(i in 1:nspp)
     {
       y<-samp[,i]
-      mod<-brglm(formu,data=cbind(y,env))
+      mod<-brglm(formu,data=data.frame(cbind(y,env)))
       spp.resids<-cbind(spp.resids,mod$y-mod$fitted.values)
       spp.fits<-cbind(spp.fits,mod$fitted.values)
       spp.coef<-cbind(spp.coef,summary(mod)$coef[,1])
@@ -84,5 +83,34 @@
   cors.pa<-cor(samp)[upper.tri(cor(samp))]
   names(cors.pa)<-pairnames
   names(cors.resid)<-pairnames
-  return(list(family=fam,residuals=spp.resids,coefficients=spp.coef,std.errors=spp.se,cors.pa=cors.pa,cors.resid=cors.resid))
+
+# If a tree is provided
+  if (is.null(tree))
+  {
+    return(list(family=fam,residuals=spp.resids,coefficients=spp.coef,std.errors=spp.se,cors.pa=cors.pa,cors.resid=cors.resid,cors.phylo=NULL))
+  } else {
+  
+  if(is(tree)[1]=="phylo")
+  {
+    tree<-prune.sample(samp,tree)
+    # Make a correlation matrix of the species pool phylogeny
+    Cmatrix<-vcv.phylo(tree,cor=TRUE)
+  } else {
+    Cmatrix<-tree
+    species<-colnames(samp)
+    Cmatrix<-Cmatrix[species,species]
+  }
+  
+  cors.phylo<-Cmatrix[upper.tri(Cmatrix)]
+  names(cors.phylo)<-pairnames
+  return(list(family=fam,residuals=spp.resids,coefficients=spp.coef,std.errors=spp.se,cors.pa=cors.pa,cors.resid=cors.resid,cors.phylo=cors.phylo))
+ 
+  }
 }
+
+samp<-mthood.abund[,4:12]
+env<-mthood.env[,1]
+fam="binomial"
+tree<-mthood.tree
+
+sppregs(samp,env,tree,fam)