[Picante-commits] r113 - branches/gsoc/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Jun 21 05:44:55 CEST 2008
Author: mrhelmus
Date: 2008-06-21 05:44:55 +0200 (Sat, 21 Jun 2008)
New Revision: 113
Added:
branches/gsoc/R/plot.sppregs.R
branches/gsoc/R/sppregs.R
Removed:
branches/gsoc/R/sppregs.r
Log:
Added plotting function for sppregs
Added: branches/gsoc/R/plot.sppregs.R
===================================================================
--- branches/gsoc/R/plot.sppregs.R (rev 0)
+++ branches/gsoc/R/plot.sppregs.R 2008-06-21 03:44:55 UTC (rev 113)
@@ -0,0 +1,9 @@
+plot.sppregs<-function(sppreg,rows=c(1,3),cex.mag=1,x.label="phylogenetic correlations",y.label=c("occurrence correlations w/ env","occurrence correlations wo/ env","change in correlations")){
+ par(mfrow=rows,las=1,cex=cex.mag)
+ plot(sppreg$cors.phylo,sppreg$cors.pa,xlab=x.label,ylab=y.label[1],main=paste("cor =",round(cor(sppreg$cors.phylo,sppreg$cors.pa,use="pairwise.complete.obs"),4)))
+ abline(0,0,lty=2)
+ plot(sppreg$cors.phylo,sppreg$cors.resid,xlab=x.label,ylab=y.label[2],main=paste("cor =",round(cor(sppreg$cors.phylo,sppreg$cors.resid,use="pairwise.complete.obs"),4)))
+ abline(0,0,lty=2)
+ plot(sppreg$cors.phylo,sppreg$cors.resid-sppreg$cors.pa,xlab=x.label,ylab=y.label[3],main=paste("cor =",round(cor(sppreg$cors.phylo,sppreg$cors.resid-sppreg$cors.pa,use="pairwise.complete.obs"),4)))
+ abline(0,0,lty=2)
+}
\ No newline at end of file
Added: branches/gsoc/R/sppregs.R
===================================================================
--- branches/gsoc/R/sppregs.R (rev 0)
+++ branches/gsoc/R/sppregs.R 2008-06-21 03:44:55 UTC (rev 113)
@@ -0,0 +1,109 @@
+sppregs<-function(samp,env,tree=NULL,fam="binomial"){
+
+ nplots<-dim(samp)[1] # number of units in the occurence data
+ nspp<-dim(samp)[2] # number of species in the occurence data
+ sppnames<-colnames(samp) # vector of species names
+
+
+ if(is.null(dim(env)))
+ {
+ nenvs<-1
+ envnames<-names(env)
+ 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~",envnames[1]) #Make formula
+ for (t in 2:nenvs)
+ {
+ 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)
+ {
+ y<-samp[,i]
+ mod<-glm(formu,data=cbind(y,env))
+ spp.resids<-cbind(spp.resids,mod$y-mod$fitted.values)
+ spp.fits<-cbind(spp.fits,mod$fit)
+ spp.coef<-cbind(spp.coef,summary(mod)$coef[,1])
+ spp.se<-cbind(spp.se,summary(mod)$coef[,2])
+ }
+
+ cors.resid<-cor(spp.resids)[lower.tri(cor(spp.resids))]
+
+ } else {
+ samp[samp>0]<-1
+ for(i in 1:nspp)
+ {
+ y<-samp[,i]
+ 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])
+ spp.se<-cbind(spp.se,summary(mod)$coef[,2])
+
+ }
+ cor.r<-matrix(0,nrow=nspp,ncol=nspp)
+ for (j in 1:dim(spp.resids)[1])
+ {
+ invv<-matrix((spp.fits[j,]*(1-spp.fits[j,]))^(-.5),nrow=1)
+ invv<-t(invv)%*%invv
+ invv[invv>(10^10)]<-(10^10)
+ r<-matrix(spp.resids[j,],nrow=1)
+ cor.r=cor.r+((t(r)%*%r)*invv)
+ }
+
+ RC=cor.r/dim(spp.resids)[1]
+ cors.resid<-RC[upper.tri(RC)] # Observed correlations of residuals among species
+ }
+
+ colnames(spp.coef)<-colnames(samp)
+ colnames(spp.se)<-colnames(samp)
+ colnames(spp.resids)<-colnames(samp)
+ colnames(spp.fits)<-colnames(samp)
+ pairnames=NULL
+ for (o in 1:(nspp-1))
+ {
+ for (u in (o+1):nspp)
+ {
+ pairnames<-c(pairnames,paste(sppnames[o],sppnames[u],sep="-"))
+ }
+ }
+ cors.pa<-cor(samp)[upper.tri(cor(samp))]
+ names(cors.pa)<-pairnames
+ names(cors.resid)<-pairnames
+
+# 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))
+
+ }
+}
\ No newline at end of file
Deleted: branches/gsoc/R/sppregs.r
===================================================================
--- branches/gsoc/R/sppregs.r 2008-06-21 02:52:24 UTC (rev 112)
+++ branches/gsoc/R/sppregs.r 2008-06-21 03:44:55 UTC (rev 113)
@@ -1,116 +0,0 @@
-sppregs<-function(samp,env,tree=NULL,fam="binomial"){
-
- nplots<-dim(samp)[1] # number of units in the occurence data
- nspp<-dim(samp)[2] # number of species in the occurence data
- sppnames<-colnames(samp) # vector of species names
-
-
- if(is.null(dim(env)))
- {
- nenvs<-1
- envnames<-names(env)
- 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~",envnames[1]) #Make formula
- for (t in 2:nenvs)
- {
- 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)
- {
- y<-samp[,i]
- mod<-glm(formu,data=cbind(y,env))
- spp.resids<-cbind(spp.resids,mod$y-mod$fitted.values)
- spp.fits<-cbind(spp.fits,mod$fit)
- spp.coef<-cbind(spp.coef,summary(mod)$coef[,1])
- spp.se<-cbind(spp.se,summary(mod)$coef[,2])
- }
-
- cors.resid<-cor(spp.resids)[lower.tri(cor(spp.resids))]
-
- } else {
- samp[samp>0]<-1
- for(i in 1:nspp)
- {
- y<-samp[,i]
- 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])
- spp.se<-cbind(spp.se,summary(mod)$coef[,2])
-
- }
- cor.r<-matrix(0,nrow=nspp,ncol=nspp)
- for (j in 1:dim(spp.resids)[1])
- {
- invv<-matrix((spp.fits[j,]*(1-spp.fits[j,]))^(-.5),nrow=1)
- invv<-t(invv)%*%invv
- invv[invv>(10^10)]<-(10^10)
- r<-matrix(spp.resids[j,],nrow=1)
- cor.r=cor.r+((t(r)%*%r)*invv)
- }
-
- RC=cor.r/dim(spp.resids)[1]
- cors.resid<-RC[upper.tri(RC)] # Observed correlations of residuals among species
- }
-
- colnames(spp.coef)<-colnames(samp)
- colnames(spp.se)<-colnames(samp)
- colnames(spp.resids)<-colnames(samp)
- colnames(spp.fits)<-colnames(samp)
- pairnames=NULL
- for (o in 1:(nspp-1))
- {
- for (u in (o+1):nspp)
- {
- pairnames<-c(pairnames,paste(sppnames[o],sppnames[u],sep="-"))
- }
- }
- cors.pa<-cor(samp)[upper.tri(cor(samp))]
- names(cors.pa)<-pairnames
- names(cors.resid)<-pairnames
-
-# 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)
More information about the Picante-commits
mailing list