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

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

Author: mrhelmus
Date: 2008-07-01 08:28:41 +0200 (Tue, 01 Jul 2008)
New Revision: 127

Modified:
   branches/gsoc/R/pblm.r
Log:
working function: edit function to account for multiply categories

Modified: branches/gsoc/R/pblm.r
===================================================================
--- branches/gsoc/R/pblm.r	2008-06-29 14:59:36 UTC (rev 126)
+++ branches/gsoc/R/pblm.r	2008-07-01 06:28:41 UTC (rev 127)
@@ -3,6 +3,7 @@
 
   # Make a vector of associations
   A<-as.matrix(as.vector(as.matrix(assocs)))
+  data.vecs<-A
   
   #numbers of species and interactions
   nassocs<-length(A)
@@ -21,9 +22,13 @@
     }
   }
   
+  
+  
   #Clean Covariates
   #If the covariate applies to both sets, then it should be in the matrix of the longer set 
-  covnames<-"intercept"
+  
+  covnames<-NULL
+  C1covs<-NULL
   if(is.null(covars1))
   {
     C1<-NULL
@@ -31,20 +36,45 @@
     if(is.null(dim(covars1)))
     {
       C1<-matrix(covars1,nspp1,nspp2,byrow=FALSE)
-      C1<-as.matrix(as.vector(C1))
+      if(is.factor(covars1))
+      {
+        C1<-as.matrix(as.vector(C1))
+        C1covs<-cbind(C1covs,C1)
+        C1<-as.matrix(model.matrix(~as.factor(C1)-1)[,-1])
+        colnames(C1)<-paste(rep("covar1",length(levels(covars1))-1),levels(covars1)[-1],sep="-")
+      } else {
+        C1<-as.matrix(as.vector(C1))
+        C1covs<-cbind(C1covs,C1)
+      }
       covnames<-c(covnames,"covar1")
     } else {
       C1<-NULL
       for(i in 1:dim(covars1)[2])
       {
         C1hold<-matrix(covars1[,i],nspp1,nspp2,byrow=FALSE)
-        C1hold<-as.matrix(as.vector(C1hold))
-        C1<-cbind(C1,C1hold)
+        if(is.factor(covars1[,i]))
+        {
+          C1hold<-as.matrix(as.vector(C1hold))
+          C1covs<-cbind(C1covs,C1hold)
+          C1hold<-as.matrix(model.matrix(~as.factor(C1hold)-1)[,-1])
+          colnames(C1hold)<-paste(rep(colnames(covars1)[i],length(levels(covars1[,i]))-1),levels(covars1[,i])[-1],sep="-")
+          C1<-cbind(C1,C1hold)
+        } else { 
+          C1hold<-as.matrix(as.vector(C1hold))
+          C1covs<-cbind(C1covs,C1hold)
+          colnames(C1hold)<-colnames(covars1)[i]
+          C1<-cbind(C1,C1hold)
+        }
         covnames<-c(covnames,colnames(covars1)[i])
       }
     }
+  
+  data.vecs<-cbind(data.vecs,C1covs)
   }
 
+
+ 
+  C2covs<-NULL
   if(is.null(covars2))
   {
     C2<-NULL
@@ -52,18 +82,40 @@
     if(is.null(dim(covars2)))
     {
       C2<-matrix(covars2,nspp1,nspp2,byrow=TRUE)
-      C2<-as.matrix(as.vector(C2))
+      if(is.factor(covars2))
+      {
+        C2<-as.matrix(as.vector(C2))
+        C2covs<-cbind(C2covs,C2)
+        C2<-as.matrix(model.matrix(~as.factor(C2)-1)[,-1])
+        colnames(C2)<-paste(rep("covar2",length(levels(covars2))-1),levels(covars2)[-1],sep="-")
+      } else {
+        C2<-as.matrix(as.vector(C2))
+        C2covs<-cbind(C2covs,C2)
+      }
       covnames<-c(covnames,"covar2")
     } else {
       C2<-NULL
       for(i in 1:dim(covars2)[2])
       {
         C2hold<-matrix(covars2[,i],nspp1,nspp2,byrow=TRUE)
-        C2hold<-as.matrix(as.vector(C2hold))
-        C2<-cbind(C2,C2hold)
+        if(is.factor(covars2[,i]))
+        {
+          C2hold<-as.matrix(as.vector(C2hold))
+          C2covs<-cbind(C2covs,C2hold)
+          C2hold<-as.matrix(model.matrix(~as.factor(C2hold)-1)[,-1])
+          colnames(C2hold)<-paste(rep(colnames(covars2)[i],length(levels(covars2[,i]))-1),levels(covars2[,i])[-1],sep="-")
+          C2<-cbind(C2,C2hold)
+        } else { 
+          C2hold<-as.matrix(as.vector(C2hold))
+          C2covs<-cbind(C2covs,C2hold)
+          colnames(C2hold)<-colnames(covars2)[i]
+          C2<-cbind(C2,C2hold)
+        }
         covnames<-c(covnames,colnames(covars2)[i])
       }
     }
+  
+  data.vecs<-cbind(data.vecs,C2covs)
   }
   
   
@@ -90,8 +142,14 @@
   }
 
   # Begin to organize output
-  data.vecs<-cbind(A,U)
-  colnames(data.vecs)<-c("associations", covnames)
+  
+  if(is.null(dim(U)))
+  {
+    data.vecs<-data.frame(A)
+    colnames(data.vecs)<-"associations"
+  } else {    
+    colnames(data.vecs)<-c("associations", covnames)
+  }
   rownames(data.vecs)<-pairnames
   
   # Calculate Star Regression Coefficients
@@ -110,8 +168,7 @@
     rownames(coefficents)<-c("upper CI 95%","estimate","lower CI 95%")
     colnames(coefficents)<-paste("star",covnames,sep="-")
     MSEs<-data.frame(MSEStar)
-    return(list(MSE=MSEs,signal.strength=NULL,coefficents=data.frame(coefficents),variates.residuals=cbind(data.vecs,data.frame(Estar))))
-    
+    return(list(MSE=MSEs,signal.strength=NULL,coefficents=data.frame(coefficents),variates=data.frame(data.vecs),residuals=data.frame(Estar)))
   } else {
   
     #tree1 is the phylogeny for the rows
@@ -206,9 +263,9 @@
     #organize output
     coefficents<-cbind(approxCFfull,approxCFstar,approxCFbase)
     rownames(coefficents)<-c("upper CI 95%","estimate","lower CI 95%")
-    colnames(coefficents)<-c(paste("full",covnames,sep="-"),paste("star",covnames,sep="-"),paste("base",covnames,sep="-"))
+    colnames(coefficents)<-c(paste("full",c("intercept",colnames(U)[-1]),sep="-"),paste("star",c("intercept",colnames(U)[-1]),sep="-"),paste("base",c("intercept",colnames(U)[-1]),sep="-"))
     MSEs<-cbind(MSEFull, data.frame(MSEStar), data.frame(MSEbase))
-    residuals<-cbind(Efull,Estar,Ebase)
+    residuals<-cbind(data.frame(Efull),data.frame(Estar),data.frame(Ebase))
     rownames(residuals)<-pairnames
   
     #bootstrap CIs
@@ -217,7 +274,7 @@
     {
     return(NULL)
     } else {
-    return(list(MSE=MSEs,signal.strength=cbind(d1,d2),coefficents=coefficents,variates.residuals=cbind(data.vecs,data.frame(Efull),data.frame(Estar),data.frame(Ebase))))
+    return(list(MSE=MSEs,signal.strength=data.frame(cbind(d1,d2)),coefficents=data.frame(coefficents),variates=data.frame(data.vecs),residuals=residuals))
     }
   }