# [Picante-commits] r132 - branches/gsoc/R

Wed Jul 9 08:18:39 CEST 2008

```Author: mrhelmus
Date: 2008-07-09 08:18:39 +0200 (Wed, 09 Jul 2008)
New Revision: 132

Modified:
branches/gsoc/R/pblm.r
Log:
cleaned up code

Modified: branches/gsoc/R/pblm.r
===================================================================
--- branches/gsoc/R/pblm.r	2008-07-09 06:17:49 UTC (rev 131)
+++ branches/gsoc/R/pblm.r	2008-07-09 06:18:39 UTC (rev 132)
@@ -1,6 +1,5 @@
pblm<-function(assocs,tree1=NULL,tree2=NULL,covars1=NULL,covars2=NULL,bootstrap=FALSE,nreps=10,maxit=10000,pstart=c(.5,.5)){

-
# Make a vector of associations
A<-as.matrix(as.vector(as.matrix(assocs)))
data.vecs<-A
@@ -21,12 +20,9 @@
pairnames<-c(pairnames,paste(sppnames2[o],sppnames1[u],sep="-"))
}
}
-
-
-
+
#Clean Covariates
#If the covariate applies to both sets, then it should be in the matrix of the longer set
-
covnames<-NULL
C1covs<-NULL
if(is.null(covars1))
@@ -73,7 +69,6 @@
}

-
C2covs<-NULL
if(is.null(covars2))
{
@@ -119,7 +114,7 @@
}

-# Make U
+# Make U, the combined matrix of covariates
U<-NULL
if(is.null(C1) & is.null(C2))
{
@@ -142,7 +137,6 @@
}

# Begin to organize output
-
if(is.null(dim(U)))
{
data.vecs<-data.frame(A)
@@ -152,8 +146,8 @@
}
rownames(data.vecs)<-pairnames

+  ######
# Calculate Star Regression Coefficients
-
#calculate for the star (assuming no phylogenetic correlation)
astar<-solve((t(U)%*%U),(t(U)%*%A))
MSETotal<-cov(A)
@@ -162,11 +156,12 @@
approxCFstar<-rbind(t(astar)-1.96%*%sdaStar, t(astar), t(astar)+1.96%*%sdaStar)
Estar<-A-U%*%astar
MSEStar<-cov(matrix(Estar))
+  #######

if(is.null(tree1) | is.null(tree2))
{
coefs<-approxCFstar
-    rownames(coefs)<-c("upper CI 95%","estimate","lower CI 95%")
+    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))
@@ -205,8 +200,8 @@
invV<-qr.solve(V)

abase<-solve((t(U)%*%invV%*%U),((t(U)%*%invV%*%A)))   #NOTE: Ives in his Matlab code uses a Left matrix division symbol (\)
-    MSEbase<-(t(A-U%*%abase)%*%invV%*%(A-U%*%abase))/(nassocs-1)
-    s2abase<-as.vector(MSEbase)*qr.solve(t(U)%*%invV%*%U)
+    MSEBase<-(t(A-U%*%abase)%*%invV%*%(A-U%*%abase))/(nassocs-1)
+    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
@@ -221,8 +216,7 @@
tau1<-matrix(diag(initV1),nspp1,nspp1) + matrix(diag(initV1),nspp1,nspp1)-2*initV1
tau2<-matrix(diag(initV2),nspp2,nspp2) + matrix(diag(initV2),nspp2,nspp2)-2*initV2

-
-    # The work horse function to estimate ds
+    # The workhorse function to estimate ds
pegls<-function(parameters)
{
d1<-abs(parameters[1])
@@ -263,11 +257,11 @@

#organize output
coefs<-cbind(approxCFfull,approxCFstar,approxCFbase)
-    rownames(coefs)<-c("approx upper CI 95%","estimate","approx lower CI 95%")
+    rownames(coefs)<-c("approx lower CI 95%","estimate","approx upper CI 95%")
colnames(coefs)<-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="-"))
coefs<-t(coefs)
CI.boot<-NULL
-    MSEs<-cbind(data.frame(MSETotal),data.frame(MSEFull), data.frame(MSEStar), data.frame(MSEbase))
+    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))
rownames(residuals)<-pairnames

@@ -322,7 +316,9 @@
}
nr<-dim(bootlist)[1]
nc<-dim(bootlist)[2]
-      alpha<-0.05
+
+      #Calculate bootstrapped CIs
+      alpha<-0.05  # alpha is always 0.05, but could change here
conf<-NULL
for(j in 1:nc)
{
@@ -333,6 +329,7 @@
rownames(signal.strength)<-c("d1","d2")
colnames(signal.strength)<-c("booted lower CI 95%","estimate","booted upper CI 95%")

+      #organize output
CI.boot<-conf
rownames(CI.boot)<-c("d1","d2","intercept",colnames(U)[-1])
colnames(CI.boot)<-c("booted lower CI 95%","booted upper CI 95%")
@@ -340,6 +337,9 @@
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))

} else {
+    ########
+    # If bootstrapping not performed
+
conf<-matrix(NA,2,2)
signal.strength<-data.frame(cbind(conf[1,],dtrue,conf[2,]))
rownames(signal.strength)<-c("d1","d2")
@@ -347,5 +347,4 @@
return(list(MSE=MSEs,signal.strength=signal.strength,coefficients=data.frame(coefs),CI.boot=NULL,variates=data.frame(data.vecs),residuals=residuals,bootvalues=NULL))
}
}
-
-}
\ No newline at end of file
+}
\ No newline at end of file

```