[Picante-commits] r234 - in pkg: . R inst/doc man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Sep 4 03:20:36 CEST 2012
Author: skembel
Date: 2012-09-04 03:20:36 +0200 (Tue, 04 Sep 2012)
New Revision: 234
Added:
pkg/NAMESPACE
pkg/R/pglmm.R
pkg/man/pglmm.Rd
Modified:
pkg/DESCRIPTION
pkg/R/pblm.R
pkg/R/pblmpredict.R
pkg/R/pcd.R
pkg/R/phylodiversity.R
pkg/R/phylosignal.R
pkg/R/phylostruct.R
pkg/R/sppregs.R
pkg/inst/doc/picante-intro.Rnw
pkg/inst/doc/picante-intro.pdf
pkg/man/picante-package.Rd
pkg/man/traitgram.Rd
Log:
Adding pglmm functions. Fix typos in vignette and documentation. Add NAMESPACE.
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2011-04-22 23:44:26 UTC (rev 233)
+++ pkg/DESCRIPTION 2012-09-04 01:20:36 UTC (rev 234)
@@ -1,11 +1,11 @@
Package: picante
Type: Package
Title: R tools for integrating phylogenies and ecology
-Version: 1.3-0
-Date: 2011-4-22
+Version: 1.4-0
+Date: 2012-8-28
Author: Steven W. Kembel <skembel at uoregon.edu>, David D. Ackerly <dackerly at berkeley.edu>, Simon P. Blomberg <s.blomberg1 at uq.edu.au>, Will K. Cornwell <cornwell at zoology.ubc.ca>, Peter D. Cowan <pdc at berkeley.edu>, Matthew R. Helmus <mrhelmus at wisc.edu>, Helene Morlon <morlon.helene at gmail.com>, Campbell O. Webb <cwebb at oeb.harvard.edu>
Maintainer: Steven W. Kembel <skembel at uoregon.edu>
Depends: ape, vegan, nlme
-Suggests: brglm, circular, quantreg
+Suggests: brglm, circular, corpcor, quantreg, plotrix
Description: Phylocom integration, community analyses, null-models, traits and evolution in R
License: GPL-2
Added: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE (rev 0)
+++ pkg/NAMESPACE 2012-09-04 01:20:36 UTC (rev 234)
@@ -0,0 +1,11 @@
+useDynLib(picante)
+
+# Export all names
+exportPattern(".")
+
+# Import all packages listed as Imports or Depends
+import(
+ ape,
+ vegan,
+ nlme
+)
Modified: pkg/R/pblm.R
===================================================================
--- pkg/R/pblm.R 2011-04-22 23:44:26 UTC (rev 233)
+++ pkg/R/pblm.R 2012-09-04 01:20:36 UTC (rev 234)
@@ -181,7 +181,7 @@
if(is(tree1)[1]=="phylo")
{
if(is.null(tree1$edge.length)){tree1<-compute.brlen(tree1, 1)} #If phylo has no given branch lengths
- V1<-vcv.phylo(tree1,cor=TRUE)
+ V1<-vcv.phylo(tree1,corr=TRUE)
V1<-V1[rownames(assocs),rownames(assocs)]
} else {
V1<-tree1[rownames(assocs),rownames(assocs)]
@@ -190,7 +190,7 @@
if(is(tree2)[1]=="phylo")
{
if(is.null(tree2$edge.length)){tree2<-compute.brlen(tree2, 1)} #If phylo has no given branch lengths
- V2<-vcv.phylo(tree2,cor=TRUE)
+ V2<-vcv.phylo(tree2,corr=TRUE)
V2<-V2[colnames(assocs),colnames(assocs)]
} else {
V2<-tree2[colnames(assocs),colnames(assocs)]
Modified: pkg/R/pblmpredict.R
===================================================================
--- pkg/R/pblmpredict.R 2011-04-22 23:44:26 UTC (rev 233)
+++ pkg/R/pblmpredict.R 2012-09-04 01:20:36 UTC (rev 234)
@@ -72,7 +72,7 @@
if(is(tree1.w.novel)[1]=="phylo")
{
if(is.null(tree1.w.novel$edge.length)){tree1.w.novel<-compute.brlen(tree1.w.novel, 1)} #If phylo has no given branch lengths
- V1<-vcv.phylo(tree1.w.novel,cor=TRUE)
+ V1<-vcv.phylo(tree1.w.novel,corr=TRUE)
} else {
V1<-tree1.w.novel
}
@@ -80,7 +80,7 @@
if(is(tree2.w.novel)[1]=="phylo")
{
if(is.null(tree2.w.novel$edge.length)){tree2.w.novel<-compute.brlen(tree2.w.novel, 1)} #If phylo has no given branch lengths
- V2<-vcv.phylo(tree2.w.novel,cor=TRUE)
+ V2<-vcv.phylo(tree2.w.novel,corr=TRUE)
} else {
V2<-tree2.w.novel
}
Modified: pkg/R/pcd.R
===================================================================
--- pkg/R/pcd.R 2011-04-22 23:44:26 UTC (rev 233)
+++ pkg/R/pcd.R 2012-09-04 01:20:36 UTC (rev 234)
@@ -14,7 +14,7 @@
{
if (is.null(tree$edge.length)) {tree <- compute.brlen(tree, 1)} #If phylo has no given branch lengths
tree <- prune.sample(comm, tree)
- V <- vcv.phylo(tree, cor = TRUE)
+ V <- vcv.phylo(tree, corr = TRUE)
comm <- comm[, tree$tip.label]
} else {
V <- tree
Added: pkg/R/pglmm.R
===================================================================
--- pkg/R/pglmm.R (rev 0)
+++ pkg/R/pglmm.R 2012-09-04 01:20:36 UTC (rev 234)
@@ -0,0 +1,666 @@
+
+###################################################################################################################
+#PGLMM code from Ives A.R. & Helmus M.R. (2011). Generalized linear mixed models for phylogenetic analyses of community structure. Ecological Monographs, 81, 511-525.
+#
+# This code contains the functions pglmm.sim, pglmm.data,
+###################################################################################################################
+
+
+
+#############################################################################################
+# Simulates models used to test the PGLMM models
+# 23 June 2010 MATLAB CODE GLMM_model_simulations.m
+# 6 September 2011 R CODE TRANSLATION - HELMUS
+#############################################################################################
+
+pglmm.sim <- function(tree,nsites=30,modelflag=1,figs=TRUE,second.env=TRUE,compscale = 1)
+{
+
+
+ bspp2 <- NULL
+ Vcomp <- NULL
+ envirogradflag2 <- 0
+
+ if (is(tree)[1] == "phylo")
+ {
+ if (is.null(tree$edge.length))
+ {
+ tree <- compute.brlen(tree, 1)
+ }
+ V <- vcv.phylo(tree, corr = TRUE)
+ } else {
+ V <- tree
+ }
+ nspp<-dim(V)[1]
+ # parameter values for each model
+ if(modelflag==1 | modelflag==2)
+ {
+ Xscale <- 1
+ Mscale <- .5
+ Vscale1 <- 1
+ Vscale2 <- 1
+ b0scale <- .5
+
+ envirogradflag1 <- 1
+ if(second.env){envirogradflag2 <- 1} #envirogradflag2 <- 1
+ elimsitesflag <- 1
+ repulseflag <- 0
+ }
+
+ if(modelflag==3)
+ {
+ Xscale <- 1
+ Mscale <- 0
+ Vscale1 <- 1
+ Vscale2 <- 1
+ compscale <- compscale
+ b0scale <- 0
+
+ # repulsion matrix
+ Vcomp <- solve(V,diag(nspp))
+ Vcomp <- Vcomp/max(Vcomp)
+ Vcomp <- compscale*Vcomp
+ iDcomp <- t(chol(Vcomp))
+ colnames(Vcomp)<-rownames(Vcomp)
+
+ envirogradflag1 <- 1
+ if(second.env){envirogradflag2 <- 1} # envirogradflag2=0;
+ elimsitesflag <- 0
+ repulseflag <- 1
+ }
+
+ if(modelflag==4)
+ {
+ Xscale <- 1
+ Mscale <- .5
+
+ Vscale1 <- .05
+ Vscale2 <- 1
+
+ b0scale <- .5
+
+ envirogradflag1 <- 1
+ if(second.env){envirogradflag2 <- 1} # envirogradflag2=0;
+ elimsitesflag <- 1
+ repulseflag <- 0
+ }
+
+ if(modelflag==5)
+ {
+ Xscale <- 1
+ Mscale <- .5
+
+ Vscale1 <- 1
+ Vscale2 <- 1
+
+ b0scale <- .5
+
+ envirogradflag1 <- 1
+ if(second.env){envirogradflag2 <- 1} #envirogradflag2=1
+
+ elimsitesflag <- 1
+ repulseflag <- 0
+ }
+
+ Vforsim <- V
+ iD <- t(chol(Vforsim))
+
+ # set up environmental gradient among sites
+
+ mx <- t(as.matrix((-(nsites)/2):(nsites/2)))
+ # number of sites
+ m <- length(mx)
+
+
+ #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ #% set up independent variables
+ #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ #% environmental gradient 1
+ if(envirogradflag1 == 1)
+ {
+ e <- iD %*% rnorm(nspp)
+ e <- Vscale1*(e-mean(e))/apply(e,2,sd)
+ bspp1 <- e
+
+ X <- 1/(1+exp(-(b0scale*array(1,c(m,1)) %*% rnorm(nspp) + t(mx) %*% t(e))))
+ X <- Xscale*X
+ Xsmooth <- X
+
+ X1 <- diag(1-Mscale*runif(m)) %*% X
+ }
+
+ # environmental gradient 2
+ if(envirogradflag2 == 1)
+ {
+ e <- iD %*% rnorm(nspp)
+ e <- Vscale2*(e-mean(e))/apply(e,2,sd)
+ bspp2 <- e
+
+ mx2 <- as.matrix(mx[sample(m)])
+ X <- 1/(1+exp(-(b0scale*array(1,c(m,1)) %*% rnorm(nspp) + mx2 %*% t(e))))
+ X <- Xscale*X
+ Xsmooth <- Xsmooth*X
+
+ X2 <- diag(1-Mscale*runif(m)) %*% X
+ X <- X1*X2
+ } else {
+ X <- X1
+ }
+
+ # phylogenetic repulsion
+ if(repulseflag == 1)
+ {
+ bcomp <- NULL
+ for(i in 1:m)
+ {
+ bcomp <- cbind(bcomp, iDcomp %*% rnorm(nspp))
+ }
+ bcomp0 <- 0
+ Xcomp <- exp(bcomp0+bcomp)/(1+exp(bcomp0+bcomp))
+ X <- X*t(Xcomp)
+ }
+
+ #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ #% simulate data
+ #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ #% convert distribution to presence/absence
+ Y <- matrix(0,ncol=nspp,nrow=m)
+ Y[matrix(runif(nspp*m),ncol=nspp) < X] <- 1
+ colnames(Y)<-colnames(X)
+
+ # eliminate sites with 0 spp
+ pick <- (rowSums(Y)>0)
+ Y <- Y[pick,]
+ mx <- mx[pick]
+ m <- length(mx)
+
+ # eliminate sites with 1 spp
+ if(elimsitesflag == 1){
+ pick <- (rowSums(Y)>1)
+ Y <- Y[pick,]
+ mx <- mx[pick]
+ m <- length(mx)
+ }
+
+ if(figs)
+ {
+ if (!require(plotrix))
+ {
+ stop("The 'plotrix' package is required to plot images from this function")
+ }
+
+ if(.Platform$OS.type == "unix") quartz() else windows()
+ par(mfrow=c(5,1),las=1,mar=c(2, 4, 2, 2) - 0.1)
+ matplot(Xsmooth,type="l",ylab="probability",main="Xsmooth")
+ matplot(X,type="l",ylab="probability",main="X")
+ hist(colSums(Y),main="spp per site")
+ hist(rowSums(Y),main="sites per spp")
+ plot(x=mx,y=rowSums(Y),main="SR across gradient",type="l")
+
+ if(.Platform$OS.type == "unix") quartz() else windows()
+ color2D.matplot(1-X/max(X),xlab="species",ylab="sites",main="probabilities")
+
+ if(.Platform$OS.type == "unix") quartz() else windows()
+ color2D.matplot(1-Y,xlab="species",ylab="sites",main="presence-absence")
+ }
+
+ return(list(Vphylo=V,Vcomp=Vcomp,Y=Y,X=X,u=mx,bspp1=bspp1,bspp2=bspp2))
+}
+
+######################################################################################################################
+##Function that organizes the data so that PGLMM can be fit.
+######################################################################################################################
+
+pglmm.data<-function(modelflag=1,sim.dat=NULL,samp=NULL,tree=NULL,traits=NULL,env=NULL,Vcomp=NULL)
+{
+
+ if(!is.null(sim.dat))
+ {
+ tree<-sim.dat$Vphylo
+ Vcomp<-sim.dat$Vcomp
+ samp<-sim.dat$Y
+ traits<-sim.dat$bspp1
+ env<-sim.dat$u
+ }
+
+ is.empty<-function(x){length(x) == 0}
+ if(is.empty(samp))
+ {
+ stop("sample matrix (Y) is empty")
+ }
+
+ if (is(tree)[1] == "phylo")
+ {
+ if (is.null(tree$edge.length))
+ {
+ tree <- compute.brlen(tree, 1)
+ }
+ tree <- prune.sample(samp, tree)
+ samp <- samp[, tree$tip.label]
+ V <- vcv.phylo(tree, corr = TRUE)
+ species <- colnames(samp)
+ preval <- colSums(samp)/sum(sum(samp))
+ species <- species[preval > 0]
+ V <- V[species, species]
+ Vcomp <- Vcomp[species, species]
+ samp <- samp[, colnames(V)]
+ traits <- as.matrix(traits[species,])
+
+
+ } else {
+ V <- tree
+ species <- colnames(samp)
+ preval <- colSums(samp)/sum(sum(samp))
+ species <- species[preval > 0]
+ V <- V[species, species]
+ Vcomp <- Vcomp[species, species]
+ samp <- samp[, colnames(V)]
+ traits <- as.matrix(traits[species,])
+ }
+
+ # X = species independent variables (traits)
+ # U = site independent variables (environment)
+ # Y = site (rows) by species (columns) presence/absence matrix, the dependent variable (binary 0,1)
+ # V = phylogenetic covariance matrix
+ Y<-samp
+ U<-matrix(env,nrow=length(env),ncol=1)
+ X<-traits
+
+ nsites<-dim(Y)[1]
+ nspp<-dim(Y)[2]
+
+ # create dependent variable vector
+ YY<-t(Y)
+ YY<-as.matrix(as.vector(as.matrix(YY)))
+
+ if(modelflag==1)
+ {
+ # set up covariance matrices
+ Vfullspp<-kronecker(diag(nsites),V) #should be nspp*nsites in dimension
+ Vfullsite<-kronecker(diag(nsites),matrix(1,nspp,nspp))
+
+ #VV<-abind(Vfullspp,Vfullsite,along=3) #could potentially also use the abind function in the abind library but I will see if we can get this to work instead
+ VV<-list(Vfullspp=Vfullspp,Vfullsite=Vfullsite)
+
+ # create independent variables
+ XX<-kronecker(matrix(1,nsites,1),diag(nspp))
+ return(list(YY=YY,VV=VV,XX=XX))
+ }
+
+ if(modelflag==2)
+ {
+ # set up covariance matrices
+ u<-scale(U)
+ U<-kronecker(u,matrix(1,nspp,1))
+
+ Vfullspp<-kronecker(matrix(1,nsites,nsites),diag(nspp))
+ VfullsppV<-kronecker(matrix(1,nsites,nsites),V)
+
+ VfullUCU<-diag(as.vector(U)) %*% Vfullspp %*% diag(as.vector(U))
+ VfullUCUV<-diag(as.vector(U)) %*% VfullsppV %*% diag(as.vector(U))
+
+ Vfullsite<-kronecker(diag(nsites),matrix(1,nspp,nspp))
+
+ VV<-list(VfullUCU=VfullUCU,VfullUCUV=VfullUCUV,Vfullsite=Vfullsite)
+
+ # create independent variables
+ XXspp<-kronecker(matrix(1,nsites,1),diag(nspp))
+ XX<-cbind(U,XXspp)
+
+ # create dependent variable vector
+ YY<-as.vector(t(Y))
+
+ return(list(YY=YY,VV=VV,XX=XX))
+ }
+
+ if(modelflag == 3)
+ {
+ # set up covariance matrices
+ u<-scale(U)
+ U<-kronecker(u,matrix(1,nspp,1))
+
+ #Vfullspp<-kronecker(diag(nsites),V)
+
+ #VfullUCU<-diag(as.vector(U)) %*% Vfullspp %*% diag(as.vector(U))
+
+ Vfullsite<-kronecker(diag(nsites),matrix(1,nspp,nspp))
+ if(is.null(Vcomp))
+ {
+ compscale<-1
+ # repulsion matrix
+ Vcomp <- solve(V,diag(nspp))
+ Vcomp <- Vcomp/max(Vcomp)
+ Vcomp <- compscale*Vcomp
+ }
+ Vfullcomp<-kronecker(diag(nsites),Vcomp)
+
+ VV<-list(Vfullcomp=Vfullcomp,Vfullsite=Vfullsite)
+
+ # create independent variables
+ XXspp<-kronecker(matrix(1,nsites,1),diag(nspp))
+ XX<-cbind(XXspp, ((U %*% matrix(1,1,nspp))*XXspp))
+
+ # create dependent variable vector
+ YY<-as.vector(t(Y))
+
+ return(list(YY=YY,VV=VV,XX=XX))
+ }
+
+ if(modelflag == 4)
+ {
+ # set up covariance matrix for traits
+ Vfulltrait <- kronecker(diag(nsites),traits %*% t(traits))
+ traitscale4 <- 100
+ Vfulltrait <- traitscale4*Vfulltrait
+ # set up covariance matrices
+ Vfullsite<-kronecker(diag(nsites),matrix(1,nspp,nspp))
+ VV<-list(Vfulltrait=Vfulltrait,Vfullsite=Vfullsite)
+
+ # create independent variables
+ XXspp<-kronecker(matrix(1,nsites,1),diag(nspp))
+ XX<-XXspp
+
+ # create dependent variable vector
+ YY<-as.vector(t(Y))
+ return(list(YY=YY,VV=VV,XX=XX))
+ }
+
+ if(modelflag == 5)
+ {
+ # set up covariance matrices
+ Vfulltrait <- kronecker(diag(nsites),traits %*% t(traits))
+ traitscale5 <- 10
+ Vfulltrait <- traitscale5*Vfulltrait
+
+ Vfullspp<-kronecker(diag(nsites),V) #should be nspp*nsites in dimension
+
+ Vfullsite<-kronecker(diag(nsites),matrix(1,nspp,nspp))
+
+ VV<-list(Vfulltrait=Vfulltrait,Vfullspp=Vfullspp,Vfullsite=Vfullsite)
+
+ # create independent variables
+ XXspp<-kronecker(matrix(1,nsites,1),diag(nspp))
+ XX<-XXspp
+
+ # create dependent variable vector
+ YY<-as.vector(t(Y))
+ return(list(YY=YY,VV=VV,XX=XX))
+ }
+}
+
+###########################################################################################################
+# Calls a PGLMM for logistic regression, workhorse function
+###########################################################################################################
+#global tH tinvW tVV tX
+tH<<-NULL
+tinvW<<-NULL
+tVV<<-NULL
+tX<<-NULL
+pglmm.fit<-function(dat=NULL,Y=NULL,X=NULL,VV=NULL,sp.init=0.5,maxit=25,exitcountermax=50) # [B,B0,s,S95int,LL,flag]=PGLMM_nosparse_funct(Y,X,VV,s)
+{
+ #####################################################################################################
+ #########REML estimation Function not called directly by the user
+ #####################################################################################################
+ PGLMM.reml<-function(sp)
+ {
+
+
+ sp<-abs(Re(sp))
+ Cd<-matrix(0,dim(tVV[[1]])[1],dim(tVV[[1]])[2])
+ for(i in 1:length(sp))
+ {
+ Cd=Cd + sp[i] * tVV[[i]]
+ }
+ V<- tinvW + Cd
+ invV<- solve(V,diag(x=1,nrow=dim(V)[1],ncol=dim(V)[2]))
+
+ # check to see if V is positive definite
+ if(all( eigen(V)$values >0 ))
+ {
+ cholV<-chol(V)
+ # ML
+ # LL=.5*(2*sum(log(diag(chol(V))))+tH'*invV*tH)
+
+ #REML
+ LL=.5 * (2 * sum(log(diag(cholV))) + t(tH) %*% invV %*% tH + log(det(t(tX) %*% invV %*% tX)))
+ } else {
+ LL<-10^10
+ }
+ LL
+ }
+
+
+ if (!require(corpcor))
+ {
+ stop("The 'corpcor' package is required")
+ }
+
+ # dat = list from PGLMM.data
+ # X = independent variable
+ # Y = dependent variable (binary 0,1)
+ # VV = list containing covariance matrices
+ # s = initial values of s
+
+ # Returns
+ # B_SE = coefficients with SEs from GLMM
+ # B0_SE = coefficients with SEs from logistic regression
+ # s = parameter(s) of the covariance matrix
+ # flag = 1 if convergence achieved, 0 otherwise
+
+ #global tH tinvW tVV tX Lmin S iss
+ if(!is.null(dat))
+ {
+ X<-dat$XX
+ Y<-dat$YY
+ VV<-dat$VV
+ }
+
+ is.empty<-function(x){length(x) == 0}
+ if(any(unlist(lapply(list(X=X,Y=Y,VV=VV),is.empty))))
+ {
+ stop("a data matrix is empty")
+ }
+
+ if(any(unlist(lapply(list(X=X,Y=Y,VV=VV),is.na))))
+ {
+ stop("a data matrix is NA")
+ }
+
+ n<-dim(X)[1]
+ p<-dim(X)[2]
+
+ #initial estimates for the s parameters (scaling parameters of the covariance matrices
+ sp<-matrix(sp.init,length(as.list(VV)),1)
+
+ B0<-matrix(mean(Y),p,1)
+ oldB0<-matrix(10,p,1)
+
+ #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ #% initial values for B
+ #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ counter<-0
+ #sparse<-function(W){as(Matrix(W),"sparseMatrix")}
+ while (matrix((t(oldB0-B0) %*% (oldB0-B0) > 10^-8)) & (counter<100))
+ {
+ oldB0<-B0
+ counter<-counter+1
+
+ mu<-exp(X %*% B0) / (1+exp(X %*% B0))
+ W<-as.vector((mu*(1-mu))^-1)
+ invW<-(diag(W))
+ #invW<-sparse(diag(W))
+
+ Z<-(X %*% B0 + (Y-mu)/(mu*(1-mu)))
+ denom<-(t(X) %*% invW %*% X)
+ #Z<-sparse(X %*% B0 + (Y-mu)/(mu*(1-mu)))
+ #denom<-sparse(t(X) %*% invW %*% X)
+
+ options(warn=-1)
+ if(any(c(is.nan(denom),is.infinite(denom),is.null(denom))))
+ {
+ B0<-solve((t(X)%*% X),(t(X) %*% Y))
+ (counter<-100)
+ } else {
+ #num<-sparse(t(X) %*% invW %*% Z)
+ num<-(t(X) %*% invW %*% Z)
+ B0<-pseudoinverse(denom) %*% num
+ }
+ }
+
+ if (is.nan(B0) | is.infinite(B0))
+ {
+ mu<-matrix(mean(Y),p,1)
+ B0<-log(mu/(1-mu))
+ }
+
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ ## GLMM
+ ##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ # initial estimates
+ B<-B0
+ b<-matrix(0,n,1)
+ #bet<-sparse(rBind(B,b))
+ #mu<-sparse(exp(X %*% B)/(1+exp(X %*% B)))
+ bet<-(rbind(B,b))
+ mu<-(exp(X %*% B)/(1+exp(X %*% B))) # the mean function for the binomial process change this to get at other distributions
+
+ # matrix including fixed covariates and dummies for "random effects"
+ XX<-cbind(X,diag(n))
+ Cdum<-matrix(0,n,n)
+ for(i in 1:length(sp))
+ {
+ Cdum=Cdum+(sp[i]*VV[[i]])
+ }
+ #Cdum<-Matrix(Cdum)
+ est<-t(rbind(sp,B))
+ oldest<-matrix(10^6,dim(est)[1],dim(est)[2])
+
+ exitcounter<-0
+ while (matrix(((est-oldest) %*% t(est-oldest)) > 10^-4) & exitcounter<=exitcountermax)
+ {
+ oldest<-est
+ # Using notation of Breslow & Clayton (1993)
+ # Note: Schall (1991) equation for V incomplete; see McCullagh & Nelder p. 40 and
+ # Breslow & Clayton (1993)
+ W<-as.vector((mu*(1-mu))^-1)
+ #invW<-sparse(diag(W))
+ invW<-(diag(W))
+
+ V<-invW+Cdum #Breslow & Clayton (1993) between eq 10 and 11
+ invV<-solve(V,diag(n)) # needed for eq 10
+
+ #################### I did not deal with this code for when the solve function returns a singular matrix ###################
+ #ww=lastwarn;
+ #lastwarn('noper')
+ #if sum(ww(1:5)=='Matri')==5
+ #'Matrix close to singular: alternative B0 tried'
+ #mu=rand(length(X(1,:)),1);
+ #B0=log(mu./(1-mu));
+ #b=zeros(n,1);
+ #beta=[B0;b];
+
+ #mu=exp(X*B0)./(1+exp(X*B0));
+
+ #invW=diag((mu.*(1-mu)).^-1);
+ #V=invW+C;
+ #invV=V\eye(n);
+
+ #B0'
+ #end
+ #########################################
+
+ Z<-X %*% B + b + (Y-mu)/(mu*(1-mu))
+ denom<-(t(X) %*% invV %*% X) # left side of eq 10
+ #denom<-sparse(t(X) %*% invW %*% X)
+ num<-(t(X) %*% invV %*% Z) #right side of eq 10
+ B<-pseudoinverse(denom) %*% num # solve eq 10 for the fixed effects (alpha in Breslow & Clayton (1993))
+ b<-Cdum %*% invV %*% (Z-X %*% B) #eq 11
+
+ bet<-(rbind(B,b))
+ mu<-exp(XX %*% bet)/(1+exp(XX %*% bet))
+
+ #DEFINE THESE AS GLOBAL tH tinvW tVV tX
+ tH<-Z - X %*% B
+ tinvW<-diag(as.vector((mu*(1-mu))^-1))
+ tX<-X
+ tVV<-VV
+
+ # call to obtain estimates of covariance matrix parameters s
+ #options=optimset('MaxFunEvals',25,'MaxIter',25,'TolX',10^-2,'TolFun',10^-2);
+ #pars<-list(sp=sp,tVV=tVV,tH=tH,tinvW=tinvW,tVV=tVV,tX=tX)
+ sp<-abs(optim(sp,PGLMM.reml,control=list(maxit=maxit,abstol=10^-1))$par)
+
+ if(exitcounter==0){
+ scrn.output<-c(exitcounter,t(sp), t(B[1:4]))
+ names(scrn.output)<-c("exitcounter",paste("sigma",1:length(sp)),paste("B",1:4))
+ print(scrn.output)
+ } else {
+ print(c(exitcounter,t(sp), t(B[1:4])))
+ }
+
+ Cdum<-matrix(0,n,n)
+ for(i in 1:length(sp))
+ {
+ Cdum=Cdum + sp[i] * tVV[[i]]
+ }
+ est<-t(rbind(sp,B))
+ exitcounter<-exitcounter+1
+ }
+
+ # flag cases of non-convergence
+ flag <- "converged"
+ if(exitcounter>=exitcountermax)
+ {
+ flag<-"did not converge, try increasing exitcountermax"
+ }
+
+ ## flag cases of non-convergence
+ #if(is.nan(B)){
+ # return
+ #}
+
+ ##############
+ # compute final estimates and SEs
+ W<-as.vector((mu*(1-mu))^-1)
+ #invW<-sparse(diag(W))
+ invW<-(diag(W))
+ V<-invW+Cdum
+ invV<-solve(V,diag(n))
+
+ Z<-X %*% B + b + (Y-mu)/(mu*(1-mu))
+ B<-solve((t(X) %*% invV %*% X),(t(X) %*% invV %*% Z))
+
+ Chi025<-5.0238
+ Chi05<-3.8414
+
+ #options=optimset('MaxFunEvals',100,'MaxIter',100,'TolX',5*10^-3,'TolFun',10^-4);
+
+ S95int<-NULL
+ S<-sp
+ for(iss in 1:length(sp))
+ {
+ Lmin<-PGLMM.reml(S)
+ #Smin
+ if(sp[iss] > 0.02)
+ {Smin<-optimize(PGLMM.reml,c(0,.9*sp[iss]),tol = 0.0001)$minimum
+ } else {
+ Smin<-0}
+ if(Smin<0.01){Smin<-0}
+ #Smax
+ if(sp[iss] > 0.02){
+ Smax<-optimize(PGLMM.reml,c(1.1*sp[iss],10*sp[iss]),tol = 0.0001)$minimum
+ } else {
+ Smax<-optimize(PGLMM.reml,c(.1,5),tol = 0.0001)$minimum}
+ S95int<-rbind(S95int,c(abs(Smin),abs(Smax)))
+ }
+
+ names(sp)<-names(VV)
+ colnames(S95int)<-c("0.05","0.95")
+ return(list(B=B,B0=B0,s=cbind(sp,S95int),LL=Lmin,flag=flag))
+}
+
+
+###########################END
\ No newline at end of file
Modified: pkg/R/phylodiversity.R
===================================================================
--- pkg/R/phylodiversity.R 2011-04-22 23:44:26 UTC (rev 233)
+++ pkg/R/phylodiversity.R 2012-09-04 01:20:36 UTC (rev 234)
@@ -240,7 +240,7 @@
# Make sure that the species line up
samp<-samp[,tree$tip.label]
# Make a correlation matrix of the species pool phylogeny
- Cmatrix<-vcv.phylo(tree,cor=TRUE)
+ Cmatrix<-vcv.phylo(tree,corr=TRUE)
} else {
Cmatrix<-tree
species<-colnames(samp)
@@ -370,7 +370,7 @@
# Make sure that the species line up
samp<-samp[,tree$tip.label, drop=FALSE]
# Make a correlation matrix of the species pool phylogeny
- Cmatrix<-vcv.phylo(tree,cor=TRUE)
+ Cmatrix<-vcv.phylo(tree,corr=TRUE)
} else {
Cmatrix<-tree
species<-colnames(samp)
@@ -430,7 +430,7 @@
# Make sure that the species line up
samp<-samp[,tree$tip.label]
# Make a correlation matrix of the species pool phylogeny
- Cmatrix<-vcv.phylo(tree,cor=TRUE)
+ Cmatrix<-vcv.phylo(tree,corr=TRUE)
} else {
Cmatrix<-tree
species<-colnames(samp)
@@ -486,7 +486,7 @@
# Make sure that the species line up
samp<-samp[,tree$tip.label]
# Make a correlation matrix of the species pool phylogeny
- Cmatrix<-vcv.phylo(tree,cor=TRUE)
+ Cmatrix<-vcv.phylo(tree,corr=TRUE)
} else {
Cmatrix<-tree
species<-colnames(samp)
Modified: pkg/R/phylosignal.R
===================================================================
--- pkg/R/phylosignal.R 2011-04-22 23:44:26 UTC (rev 233)
+++ pkg/R/phylosignal.R 2012-09-04 01:20:36 UTC (rev 234)
@@ -8,7 +8,7 @@
dat2 <- data.frame(x=x)
mat <- vcv.phylo(phy)
dat2$vars <- diag(mat)
- matc <- vcv.phylo(phy, cor=TRUE) # correlation matrix
+ matc <- vcv.phylo(phy, corr=TRUE) # correlation matrix
ntax <- length(phy$tip.label)
# calculate "phylogenetic" mean via gls
Modified: pkg/R/phylostruct.R
===================================================================
--- pkg/R/phylostruct.R 2011-04-22 23:44:26 UTC (rev 233)
+++ pkg/R/phylostruct.R 2012-09-04 01:20:36 UTC (rev 234)
@@ -4,7 +4,7 @@
null.model<-match.arg(null.model)
if(metric=="sppregs")
{
- nulls<-t(replicate(runs,sppregs(randomizeMatrix(samp,null.model=null.model,it=it),env,tree,fam=fam)$correlations))
+ nulls<-t(replicate(runs,sppregs(randomizeMatrix(samp,null.model=null.model,iterations=it),env,tree,fam=fam)$correlations))
obs<-sppregs(samp,env,tree,fam=fam)$correlations
mean.null<-apply(nulls,2,mean)
quantiles.null<-t(apply(nulls,2,quantile,probs=c(alpha/2,1-(alpha/2))))
@@ -15,10 +15,10 @@
} else {
nulls<-switch(metric,
- psv = replicate(runs,mean(psv(as.matrix(randomizeMatrix(samp,null.model=null.model,it=it)),tree,compute.var=FALSE)[,1],na.rm=TRUE)),
- psr = replicate(runs,mean(psr(as.matrix(randomizeMatrix(samp,null.model=null.model,it=it)),tree,compute.var=FALSE)[,1],na.rm=TRUE)),
- pse = replicate(runs,mean(pse(as.matrix(randomizeMatrix(samp,null.model=null.model,it=it)),tree)[,1],na.rm=TRUE)),
- psc = replicate(runs,mean(psc(as.matrix(randomizeMatrix(samp,null.model=null.model,it=it)),tree)[,1],na.rm=TRUE)))
+ psv = replicate(runs,mean(psv(as.matrix(randomizeMatrix(samp,null.model=null.model,iterations=it)),tree,compute.var=FALSE)[,1],na.rm=TRUE)),
+ psr = replicate(runs,mean(psr(as.matrix(randomizeMatrix(samp,null.model=null.model,iterations=it)),tree,compute.var=FALSE)[,1],na.rm=TRUE)),
+ pse = replicate(runs,mean(pse(as.matrix(randomizeMatrix(samp,null.model=null.model,iterations=it)),tree)[,1],na.rm=TRUE)),
+ psc = replicate(runs,mean(psc(as.matrix(randomizeMatrix(samp,null.model=null.model,iterations=it)),tree)[,1],na.rm=TRUE)))
quantiles.null<-quantile(nulls,probs=c(alpha/2,1-(alpha/2)))
mean.null<-mean(nulls)
mean.obs<-switch(metric,
Modified: pkg/R/sppregs.R
===================================================================
--- pkg/R/sppregs.R 2011-04-22 23:44:26 UTC (rev 233)
+++ pkg/R/sppregs.R 2012-09-04 01:20:36 UTC (rev 234)
@@ -11,7 +11,7 @@
# Make sure that the species line up
samp<-samp[,tree$tip.label]
# Make a correlation matrix of the species pool phylogeny
- Cmatrix<-vcv.phylo(tree,cor=TRUE)
+ Cmatrix<-vcv.phylo(tree,corr=TRUE)
} else {
Cmatrix<-tree
species<-colnames(samp)
Modified: pkg/inst/doc/picante-intro.Rnw
===================================================================
--- pkg/inst/doc/picante-intro.Rnw 2011-04-22 23:44:26 UTC (rev 233)
+++ pkg/inst/doc/picante-intro.Rnw 2012-09-04 01:20:36 UTC (rev 234)
@@ -12,8 +12,8 @@
\date{April 2010}
\begin{document}
-\SweaveOpts{strip.white=true}
-<<echo=false>>=
+\SweaveOpts{strip.white=TRUE}
+<<echo=FALSE>>=
options(width=72)
@
@@ -94,10 +94,10 @@
The following commands set up the layout of the plot to have 2 rows and 3 columns, and then plot a black dot for the species present in each of the six communities:
<<fig=TRUE>>=
-par(mfrow=c(2,3))
+par(mfrow = c(2, 3))
for (i in row.names(comm)) {
- plot(prunedphy, show.tip.label=FALSE, main=i)
- tiplabels(tip=which(comm[i,] > 0), pch=19, cex=2)
+plot(prunedphy, show.tip.label = FALSE, main = i)
+tiplabels(tip = which(prunedphy$tip.label %in% names(which(comm [i, ] > 0))) , pch=19, cex=2)
}
@
Modified: pkg/inst/doc/picante-intro.pdf
===================================================================
(Binary files differ)
Added: pkg/man/pglmm.Rd
===================================================================
--- pkg/man/pglmm.Rd (rev 0)
+++ pkg/man/pglmm.Rd 2012-09-04 01:20:36 UTC (rev 234)
@@ -0,0 +1,96 @@
+\name{pglmm}
+\alias{pglmm}
+\alias{pglmm.sim}
+\alias{pglmm.data}
+\alias{pglmm.fit}
+\title{Phylogenetic Generalized Linear Mixed Model}
+\description{Fit the presence/absence of species across communities to environmental, trait and phylogenetic data with generalized linear mixed models.}
+\usage{
+
+pglmm.sim(tree,nsites=30,modelflag=1,figs=TRUE,second.env=TRUE,compscale = 1)
+pglmm.data(modelflag=1,sim.dat=NULL,samp=NULL,tree=NULL,traits=NULL,env=NULL,Vcomp=NULL)
+pglmm.fit(dat=NULL,Y=NULL,X=NULL,VV=NULL,sp.init=0.5,maxit=25,exitcountermax=50)
+
+}
+
+\arguments{
+ \item{tree}{ Object of class phylo or a phylogenetic covariance matrix }
+ \item{nsites}{ Number of sites to simulate }
+ \item{modelflag}{ A number 1 - 5 indicating data set structure and corresponding to one of the models in Ives and Helmus (2011) }
+ \item{figs}{ Generate figures }
+ \item{compscale}{ compscale }
+ \item{second.env}{ Simulate community data with two environmental covariates }
+ \item{sim.dat}{ Object from \code{pglmm.sim} }
+ \item{samp}{ Species (rows) by site (columns) community data matrix }
+ \item{traits}{ Species X trait data matrix }
+ \item{env}{ Site X environment data matrix }
+ \item{Vcomp}{ Species X species matrix of pairwise repulsion }
+ \item{dat}{ A \code{list} with \code{Y}, \code{X} and \code{VV} from \code{pglmm.sim} or similarly structured in the same way }
+ \item{Y}{ The dependent variable, a (species times site) X 1 matrix of 0 and 1 }
+ \item{X}{ The independent variables, a (species times site) X (trait + environment) matrix }
+ \item{VV}{ A list of the covariance matrices one for each random effect }
+ \item{sp.init}{ Initial values of the variances of the random effects (e.g., phylogenetic signal) }
+ \item{maxit}{ \code{maxit} in \code{\link{optim}} used to estimate the GLMM }
+ \item{exitcountermax}{ Number of iterations to estimate the fixed effect coefficients with penalized quasilikelihood and the variances of the random effects with restricted maximum likelihood }
+}
+
+\details{
+ \emph{Phylogenetic Generalized Linear Mixed Models (PGLMM)} are generalized linear mixed model designed to test for phylogenetic patterns in community structure.
+ Five models are implemented here and are designed to address (1) phylogenetic patterns in community structure,
+ (2) phylogenetic variation in species sensitivities to environmental gradients among communities,
+ (3) phylogenetic repulsion in which closely related species are less likely to co-occur,
+ (4) trait-based variation in species sensitivities to environmental gradients,
+ and (5) the combination traits and phylogeny to explain the variation in species occurrences among communities.
+ Many other models can be designed by differently structuring the independent variable (\code{Y}), dependent variables (\code{X}), and covariance matrices (\code{VV}).
+ This can be done by either editing the \code{pglmm.data} code or designing these objects by hand with custom code.
+}
+
+\value{
+\code{pglmm.sim} returns a list with items:
+
+\item{Vphylo}{ Phylogenetic covariance matrix }
+\item{Vcomp}{ Repulsion matrix }
+\item{Y}{ Community presence/absence matrix }
+\item{X}{ Probabilities of a species being found in a community }
+\item{u}{ Environmental gradient }
+\item{bspp1}{ Species tolerances to environmental gradient 1 (phylogenetic signal) }
+\item{bspp2}{ Species tolerances to environmental gradient 2 (no phylogenetic signal) }
+
+
+\code{pglmm.data} returns a list with items:
+\item{YY}{ Independent variable }
+\item{VV}{ Covariance matrices for the random effects }
+\item{XX}{ Dependent variables }
+
+\code{pglmm.fit} returns a list with items:
+\item{B}{ Coefficient estimates }
+\item{B0}{ Initial estimates of the coefficients }
+\item{s}{ The estimates of the scaling parameters (fitted variances) for the random effects (e.g., estimate of phylogenetic signal in community composition) }
+\item{LL}{ Log likelihood of the final fitted PGLMM }
+\item{flag}{ Did the estimation procedure converge? }
+}
+
+\references{Ives A.R. & Helmus M.R. (2011). Generalized linear mixed models for phylogenetic analyses of community structure. Ecological Monographs, 81, 511-525.}
+\author{ Matthew Helmus <mrhelmus at gmail.com> and Anthony Ives <arives at wisc.edu> }
+\seealso{\code{\link{phylostruct}}}
+\note{
+The estimation procedure can be slow especially for large data sets. Thus, start with a low value for \code{exitcountermax} to test if the data and model seem to be
+correctly structured. Then increase \code{exitcountermax} to get convergence. It behooves the user to fully understand the structure of the covariance matrices
+especially when designing complicated model structures. See Ives and Helmus (2011) for a discussion.
+}
+\examples{
+
+modelflag=1
+sim.dat<-pglmm.sim(stree(16, "balanced"),nsites=30,modelflag=modelflag,second.env=TRUE,compscale=1)
+str(sim.dat)
+
+dat<-pglmm.data(modelflag=modelflag,sim.dat=sim.dat)
+str(dat)
+
+#A low number of iterations, maxit = 25 is probably good for most data sets, but exitcountermax may need to be increased depending on matrix sizes
+out<-pglmm.fit(dat=dat,maxit=25,exitcountermax=30)
+str(out)
+out$s # The first row gives the estimate of phylogenetic signal and the second an estimate of how strongly species richness varies across communities. This later parameter is likely biologically uninformative for most research questions.
+}
+
+\keyword{univar}
Modified: pkg/man/picante-package.Rd
===================================================================
--- pkg/man/picante-package.Rd 2011-04-22 23:44:26 UTC (rev 233)
+++ pkg/man/picante-package.Rd 2012-09-04 01:20:36 UTC (rev 234)
@@ -12,8 +12,8 @@
\tabular{ll}{
Package: \tab picante\cr
Type: \tab Package\cr
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/picante -r 234
More information about the Picante-commits
mailing list