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

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

Author: mrhelmus
Date: 2008-07-31 13:56:38 +0200 (Thu, 31 Jul 2008)
New Revision: 156

Modified:
   branches/gsoc/R/pblm.R
Log:
Edited pblm function to calculate residuals (correctly) and give predicted values that can be use in the pblm.predict function

Modified: branches/gsoc/R/pblm.R
===================================================================
--- branches/gsoc/R/pblm.R	2008-07-25 13:17:12 UTC (rev 155)
+++ branches/gsoc/R/pblm.R	2008-07-31 11:56:38 UTC (rev 156)
@@ -10,7 +10,6 @@
 	nspp2<-dim(assocs)[2]
   sppnames1<-rownames(assocs)
   sppnames2<-colnames(assocs)
-  
   #make names of species pairs
   pairnames=NULL  # make a vector of pairwise comparison names
   for (o in 1:(nspp2))
@@ -154,17 +153,25 @@
 	s2aStar<-as.vector(MSETotal)*qr.solve((t(U)%*%U))
 	sdaStar<-t(diag(s2aStar)^(.5))
 	approxCFstar<-rbind(t(astar)-1.96%*%sdaStar, t(astar), t(astar)+1.96%*%sdaStar)
-  Estar<-A-U%*%astar
+  Pstar<-U%*%astar
+  Estar<-A-Pstar
   MSEStar<-cov(matrix(Estar))
+  
   #######
-  
   if(is.null(tree1) | is.null(tree2))
   {
     coefs<-approxCFstar
     rownames(coefs)<-c("lower CI 95%","estimate","upper CI 95%")
     colnames(coefs)<-paste("star",c("intercept",colnames(U)[-1]),sep="-")
     MSEs<-cbind(data.frame(MSETotal),data.frame(MSEStar))
-    return(list(MSE=MSEs,signal.strength=NULL,coefficients=data.frame(t(coefs)),CI.boot=NULL,variates=data.frame(data.vecs),residuals=data.frame(Estar),bootvalues=NULL))
+    Pstar<-data.frame(Pstar)
+    colnames(Pstar)<-"star"
+    Estar<-data.frame(Estar)
+    colnames(Estar)<-"star"
+    output<-list(MSE=MSEs,signal.strength=NULL,coefficients=data.frame(t(coefs)),CI.boot=NULL,variates=data.frame(data.vecs),residuals=Estar,predicted=Pstar,bootvalues=NULL,Vfull=NULL)
+    class(output)<-"pblm"
+    return(output)
+    
   } else {
   
     #tree1 is the phylogeny for the rows
@@ -193,7 +200,6 @@
     V1<-as.matrix(V1)
     V2<-as.matrix(V2)
     
-    
     V1<-V1/det(V1)^(1/nspp1)   # scale covariance matrices (this reduces numerical problems caused by
   	V2<-V2/det(V2)^(1/nspp2)   # determinants going to infinity or zero)
   	V<-kronecker(V2,V1)  
@@ -204,8 +210,9 @@
     s2abase<-as.vector(MSEBase)*qr.solve(t(U)%*%invV%*%U)
   	sdabase<-t(diag(s2abase)^(.5))
     approxCFbase<-rbind(t(abase)-1.96%*%sdabase, t(abase), t(abase)+1.96%*%sdabase)
-    Ebase<-A-U%*%abase
-    
+    Pbase<-t(t(U%*%abase)%*%invV)
+    Ebase<-A-Pbase
+  
     ###################
     # Full EGLS estimates of phylogenetic signal
     ##################
@@ -252,7 +259,9 @@
     s2aFull<-as.vector(MSEFull)*qr.solve(t(U)%*%invV%*%U)
   	sdaFull<-t(diag(s2aFull)^(.5))
   	approxCFfull<-rbind(t(aFull)-1.96%*%sdaFull, t(aFull), t(aFull)+1.96%*%sdaFull)
-    Efull<-A-U%*%aFull
+    Pfull<-t(t(U%*%aFull)%*%invV)
+    Efull<-A-Pfull
+
     ########################################
     
     #organize output
@@ -263,10 +272,15 @@
     CI.boot<-NULL
     MSEs<-cbind(data.frame(MSETotal),data.frame(MSEFull), data.frame(MSEStar), data.frame(MSEBase))
     residuals<-cbind(data.frame(Efull),data.frame(Estar),data.frame(Ebase))
+    predicted<-cbind(data.frame(Pfull),data.frame(Pstar),data.frame(Pbase))
     rownames(residuals)<-pairnames
-  
+    rownames(predicted)<-pairnames
+    colnames(predicted)<-c("full","star","base")
+    colnames(residuals)<-c("full","star","base")
+    phylocovs=list(V1=V1,V2=V2)
+    
+    ################
     #bootstrap CIs
-    
     if(bootstrap)
     {
       Vtrue<-V
@@ -334,7 +348,9 @@
       rownames(CI.boot)<-c("d1","d2","intercept",colnames(U)[-1])
       colnames(CI.boot)<-c("booted lower CI 95%","booted upper CI 95%")
       colnames(bootlist)<-c("d1","d2","intercept",colnames(U)[-1])
-      return(list(MSE=MSEs,signal.strength=signal.strength,coefficients=data.frame(coefs),CI.boot=CI.boot,variates=data.frame(data.vecs),residuals=residuals,bootvalues=bootlist))
+      output<-list(MSE=MSEs,signal.strength=signal.strength,coefficients=data.frame(coefs),CI.boot=CI.boot,variates=data.frame(data.vecs),predicted=predicted,residuals=residuals,bootvalues=bootlist,phylocovs=phylocovs)
+      class(output)<-"pblm"
+      return(output)
     
     } else {
     ########
@@ -344,7 +360,9 @@
     signal.strength<-data.frame(cbind(conf[1,],dtrue,conf[2,]))
     rownames(signal.strength)<-c("d1","d2")
     colnames(signal.strength)<-c("booted lower CI 95%","estimate","booted upper CI 95%")
-    return(list(MSE=MSEs,signal.strength=signal.strength,coefficients=data.frame(coefs),CI.boot=NULL,variates=data.frame(data.vecs),residuals=residuals,bootvalues=NULL))
+    output<-list(MSE=MSEs,signal.strength=signal.strength,coefficients=data.frame(coefs),CI.boot=NULL,variates=data.frame(data.vecs),predicted=predicted,residuals=residuals,bootvalues=NULL,phylocovs=phylocovs)
+    class(output)<-"pblm"
+    return(output)
     }
   }                                                                                                       
 }
\ No newline at end of file