[Picante-commits] r244 - pkg/R pkg/man www
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri May 31 20:48:08 CEST 2013
Author: skembel
Date: 2013-05-31 20:48:07 +0200 (Fri, 31 May 2013)
New Revision: 244
Modified:
pkg/R/expected.pd.R
pkg/man/color.plot.phylo.Rd
pkg/man/comm.phylo.qr.Rd
pkg/man/evol.distinct.Rd
pkg/man/expected.pd.Rd
pkg/man/pblm.Rd
pkg/man/pglmm.Rd
pkg/man/phylostruct.Rd
pkg/man/randomizeSample.Rd
pkg/man/ses.mntd.Rd
pkg/man/sppregs.Rd
www/index.php
Log:
Fix typos in website, add variance.pd function, fix documentation to pass r-devel check for long lines
Modified: pkg/R/expected.pd.R
===================================================================
--- pkg/R/expected.pd.R 2013-05-29 20:33:11 UTC (rev 243)
+++ pkg/R/expected.pd.R 2013-05-31 18:48:07 UTC (rev 244)
@@ -31,28 +31,138 @@
return(byclass)
}
+#expected variance of PD
+variance.pd <- function(phy, upper.bound=TRUE) {
+ Nedges <- length(phy$edge.length)
+ if (upper.bound) {
+ ead <- .edge.children.var.ub(phy)
+ epd <- .binomial.pd(ead,Ntip(phy))
+ varpd <- .var.pd.ub(ead,Ntip(phy))
+ } else {
+ ead <- .edge.children.var.ub(phy)
+ eadvar <- .edge.children.var(phy)
+ epd <- .binomial.pd(ead,Ntip(phy))
+ varpd <- .var.pd(eadvar,Ntip(phy))
+ }
+ return(data.frame(expected.pd=epd, variance.pd=varpd))
+}
+
+# binomial PD
+.binomial.pd <-function(ead,Ntip){
+ epd=vector(mode="numeric",length=Ntip)
+ for(n in 1:Ntip){
+ epd[n]<-sum(ead$byclass$edge.length*(1-dbinom(0,ead$byclass$num.children,n/Ntip)))
+ }
+ return(epd)
+}
+
+#expected variance in pd from R(k),R(k,l)
+.var.pd <- function(test, Ntip) {
+ varpd=vector(mode="numeric",length=Ntip)
+ for(n in 1:Ntip) {
+ varpd[n]=sum(test$byclass.sq$sum.sq.edge.length*((1-n/Ntip)^test$byclass.sq$num.children)*(1-(1-n/Ntip)^test$byclass.sq$num.children))+2*sum(test$byclass.prod$sum.prods.edge.length*((1-n/Ntip)^test$byclass.prod$num.children)*(1-(1-n/Ntip)^test$byclass.prod$num.desc.children))
+ }
+ return(varpd)
+}
+
+# upper bound on variance in PD---quicker by far than var.pd
+.var.pd.ub<-function(test,Ntip) {
+ varpd <- vector(mode="numeric",length=Ntip)
+ for(n in 1:Ntip) {
+ varpd[n]=sum((test$byclass.sq$sum.sq.edge.length+2*test$byclass.prod$edge.length.times.subclade)*((1-n/Ntip)^test$byclass.sq$num.children)*(1-(1-n/Ntip)^test$byclass.sq$num.children))
+ }
+ return(varpd)
+}
+
+# edge children variance
+.edge.children.var <- function(phy) {
+ Nedges <- length(phy$edge.length)
+ edgelen <- vector(mode = "numeric", length=Nedges)
+ children <- vector(mode = "numeric", length=Nedges)
+ rawdata<-data.frame()
+ for (i in 1:Nedges) {
+ edgelen[i] <- phy$edge.length[i]
+ descnode <- phy$edge[i,2]
+ if (descnode <= Ntip(phy)) {
+ children[i] <- 1
+ }
+ else
+ {
+ desc.phy<-.extract.clade.noreord(phy, descnode)
+ children[i] <- Ntip(desc.phy)
+ desc.Nedges <- length(desc.phy$edge.length)
+ desc.edgelen <- vector(mode = "numeric", length=desc.Nedges)
+ desc.children <- vector(mode = "numeric", length=desc.Nedges)
+ for (j in 1:desc.Nedges) {
+ desc.edgelen[j] <- desc.phy$edge.length[j]
+ desc.descnode <- desc.phy$edge[j,2]
+ if (desc.descnode <= Ntip(desc.phy)) {
+ desc.children[j] <- 1
+ }
+ else
+ {
+ desc.children[j] <- Ntip(.extract.clade.noreord(desc.phy, desc.descnode))
+ }
+ }
+ rawdata <- rbind(rawdata,data.frame(num.children=matrix(children[i],nrow=length(desc.children), ncol=1), num.desc.children=desc.children, prod.edge.length=edgelen[i]*desc.edgelen))
+ }
+ }
+ byclass.prod <- aggregate(rawdata$prod.edge.length, by=list(num.children=rawdata$num.children,num.desc.children=rawdata$num.desc.children), sum)
+ rawdata2<- data.frame(num.children=children, node.edge.length=edgelen)
+ byclass.sq <- aggregate(rawdata2$node.edge.length^2, by=list(num.children=rawdata2$num.children), sum)
+ colnames(byclass.prod)[3] <- "sum.prods.edge.length"
+ colnames(byclass.sq)[2] <- "sum.sq.edge.length"
+ return(list(rawdata=rawdata, byclass.prod=byclass.prod,byclass.sq=byclass.sq))
+}
+
+# calculate components needed for upper bound on variance
+.edge.children.var.ub <- function(phy) {
+ phy <- reorder(phy)
+ Nedges <- length(phy$edge.length)
+ edgelen <- vector(mode = "numeric", length=Nedges)
+ children <- vector(mode = "numeric", length=Nedges)
+ product<-vector(mode = "numeric", length=Nedges)
+ for (i in 1:Nedges) {
+ edgelen[i] <- phy$edge.length[i]
+ descnode <- phy$edge[i,2]
+ if (descnode <= Ntip(phy)) {
+ children[i] <- 1
+ }
+ else
+ {
+ desc.tree<-.extract.clade.noreord(phy, descnode)
+ children[i] <- Ntip(desc.tree)
+ product[i]<-edgelen[i]*(sum(desc.tree$edge.length))
+ }
+ }
+ rawdata <- data.frame(num.children=children, edge.length=edgelen,product=product)
+ byclass <- aggregate(rawdata$edge.length, by=list(num.children=rawdata$num.children), sum)
+ byclass.sq <- aggregate(rawdata$edge.length^2, by=list(num.children=rawdata$num.children), sum)
+ byclass.prod<- aggregate(rawdata$product, by=list(num.children=rawdata$num.children), sum)
+ colnames(byclass)[2] <- "edge.length"
+ colnames(byclass.sq)[2] <- "sum.sq.edge.length"
+ colnames(byclass.prod)[2] <- "edge.length.times.subclade"
+ return(list(rawdata=rawdata, byclass=byclass,byclass.sq=byclass.sq,byclass.prod=byclass.prod))
+}
+
# utility function - modified from ape's extract.clade function
-.extract.clade.noreord <- function (phy, node, root.edge = 0, interactive = FALSE)
+.extract.clade.noreord <- function (phy, node, root.edge = 0)
{
Ntip <- length(phy$tip.label)
ROOT <- Ntip + 1
Nedge <- dim(phy$edge)[1]
wbl <- !is.null(phy$edge.length)
- if (interactive)
- node <- identify(phy)$nodes
- else {
- if (length(node) > 1) {
- node <- node[1]
- warning("only the first value of 'node' has been considered")
- }
- if (is.character(node)) {
- if (is.null(phy$node.label))
- stop("the tree has no node labels")
- node <- which(phy$node.label %in% node) + Ntip
- }
- if (node <= Ntip)
- stop("node number must be greater than the number of tips")
- }
+ if (length(node) > 1) {
+ node <- node[1]
+ warning("only the first value of 'node' has been considered")
+ }
+ if (is.character(node)) {
+ if (is.null(phy$node.label))
+ stop("the tree has no node labels")
+ node <- which(phy$node.label %in% node) + Ntip
+ }
+ if (node <= Ntip)
+ stop("node number must be greater than the number of tips")
if (node == ROOT)
return(phy)
# phy <- reorder(phy)
Modified: pkg/man/color.plot.phylo.Rd
===================================================================
--- pkg/man/color.plot.phylo.Rd 2013-05-29 20:33:11 UTC (rev 243)
+++ pkg/man/color.plot.phylo.Rd 2013-05-31 18:48:07 UTC (rev 244)
@@ -7,7 +7,8 @@
color.plot.phylo(phylo, df, trait, taxa.names,
num.breaks = ifelse(is.factor(df[,trait]),
length(levels(df[,trait])), 12),
- col.names = rainbow(ifelse(length(num.breaks) > 1, length(num.breaks) - 1, num.breaks)),
+ col.names = rainbow(ifelse(length(num.breaks) > 1,
+ length(num.breaks) - 1, num.breaks)),
cut.labs = NULL,
leg.title = NULL,
main = trait,
Modified: pkg/man/comm.phylo.qr.Rd
===================================================================
--- pkg/man/comm.phylo.qr.Rd 2013-05-29 20:33:11 UTC (rev 243)
+++ pkg/man/comm.phylo.qr.Rd 2013-05-31 18:48:07 UTC (rev 244)
@@ -51,5 +51,6 @@
\seealso{ \code{\link{randomizeMatrix}} }
\examples{
data(phylocom)
-comm.phylo.qr(phylocom$sample, phylocom$phylo, metric="cij", null.model="sample.taxa.labels", runs=99)}
+comm.phylo.qr(phylocom$sample, phylocom$phylo, metric="cij",
+ null.model="sample.taxa.labels", runs=99)}
\keyword{univar}
Modified: pkg/man/evol.distinct.Rd
===================================================================
--- pkg/man/evol.distinct.Rd 2013-05-29 20:33:11 UTC (rev 243)
+++ pkg/man/evol.distinct.Rd 2013-05-31 18:48:07 UTC (rev 244)
@@ -10,7 +10,8 @@
Calculates evolutionary distinctiveness measures for a suite of species by: a) equal splits (Redding and Mooers 2006) b) fair proportions (Isaac et al., 2007). Returns a datafram with species identifiers and species scores.
}
\usage{
-evol.distinct(tree, type = c("equal.splits", "fair.proportion"), scale = FALSE, use.branch.lengths = TRUE)
+evol.distinct(tree, type = c("equal.splits", "fair.proportion"),
+ scale = FALSE, use.branch.lengths = TRUE)
}
%- maybe also 'usage' for other objects documented here.
\arguments{
Modified: pkg/man/expected.pd.Rd
===================================================================
--- pkg/man/expected.pd.Rd 2013-05-29 20:33:11 UTC (rev 243)
+++ pkg/man/expected.pd.Rd 2013-05-31 18:48:07 UTC (rev 244)
@@ -1,40 +1,53 @@
\name{expected.pd}
\alias{expected.pd}
+\alias{variance.pd}
\alias{ead}
-\title{ Expected PD and Edge Abundance Distribution of a phylogeny }
+\title{ Expected PD, PD Variance, and Edge Abundance Distribution of a phylogeny }
-\description{ Calculates the expected phylogenetic diversity (Faith's PD) under binomial sampling with a fixed probability of each tip being sampled, and the Edge-length Abundance Distribution of a phylogeny. }
+\description{ Calculates the expected phylogenetic diversity (Faith's PD) and variance of PD under binomial sampling with a fixed probability of each tip being sampled, and the Edge-length Abundance Distribution of a phylogeny. }
\usage{
expected.pd(phy)
+variance.pd(phy, upper.bound=TRUE)
ead(phy)
}
\arguments{
\item{phy}{ phylo object }
+ \item{upper.bound}{ Calculate upper bound of PD variance? (default = TRUE) }
}
\value{
\item{n}{ Expected Number of tips sampled }
\item{expected.pd}{ Expected PD for a given n }
+ \item{variance.pd}{ Variance of PD for a given n }
\item{num.children}{ Number of tips descended from an edge }
\item{edge.length}{ Total phylogenetic edge length for a given number of tips descended from an edge }
}
\details{The function \code{expected.pd} calculates the expected phylogenetic diversity (Faith's PD - total branch length) for all subsets of a phylogeny, based on an analytic solution for expected PD.
+The function \code{variance.pd} additionally calculates the variance of expected PD for all subsets of a phylogeny, based on an analytic solution for expected PD. If argument upper.bound=TRUE, a fast solution for the upper bound of the variance is returned. Otherwise, the exact solution for the variance is returned. Note that the exact solution is much slower than the upper bound solution.
+
The function \code{ead} calculates the edge abundance distribution (EAD), the length of edges with different numbers of descendant tips.
}
-\references{ O'Dwyer, Kembel, and Green. 2012. Phylogenetic diversity theory sheds light on the structure of microbial communities. }
+\references{ J.P. O'Dwyer, S.W. Kembel, and J.L. Green. 2012. Phylogenetic Diversity Theory Sheds Light on the Structure of Microbial Communities. PLoS Comput Biol 8(12): e1002832. }
\author{ Steven Kembel <steve.kembel at gmail.com> and James O'Dwyer <jodwyer at santafe.edu> }
\seealso{ \code{\link{pd}} }
\examples{
-randtree <- rcoal(100)
-plot(expected.pd(randtree), xlab="Number of tips", ylab="Expected PD", type="l")
+randtree <- rcoal(300)
+randtree.pd.ub <- variance.pd(randtree, upper.bound=TRUE)
+randtree.pd.exact <- variance.pd(randtree, upper.bound=FALSE)
+plot(expected.pd(randtree), xlab="Number of tips", ylab="Phylogenetic diversity (PD)", type="l", log="xy")
+lines(randtree.pd.exact$expected.pd+1.96*sqrt(randtree.pd.exact$variance.pd), lty=2)
+lines(randtree.pd.exact$expected.pd-1.96*sqrt(randtree.pd.exact$variance.pd), lty=2)
+lines(randtree.pd.ub$expected.pd+1.96*sqrt(randtree.pd.ub$variance.pd), lty=3)
+lines(randtree.pd.ub$expected.pd-1.96*sqrt(randtree.pd.ub$variance.pd), lty=3)
+legend("bottomright", lty=c(1,2,3), legend=c("Expected PD","95 percent CI (exact)","95 percent CI (upper bound)"))
}
\keyword{univar}
Modified: pkg/man/pblm.Rd
===================================================================
--- pkg/man/pblm.Rd 2013-05-29 20:33:11 UTC (rev 243)
+++ pkg/man/pblm.Rd 2013-05-31 18:48:07 UTC (rev 244)
@@ -1,109 +1,110 @@
-\name{pblm}
-\alias{pblm}
-\alias{pblmpredict}
-\title{ Phylogenetic Bipartite Linear Model }
-\description{
- Fits a linear model to the association strengths of a bipartite data set with or without phylogenetic correlation among the interacting species
-}
-\usage{
-pblm(assocs,tree1=NULL,tree2=NULL,covars1=NULL,covars2=NULL,bootstrap=FALSE,nreps=10,maxit=10000,pstart=c(.5,.5))
-pblmpredict(x,tree1.w.novel=NULL,tree2.w.novel=NULL,predict.originals=FALSE)
-
-}
-
-\arguments{
- \item{assocs}{ A matrix of association strengths among two sets of interacting species }
- \item{tree1}{ A phylo tree object or a phylogenetic covariance matrix for the rows of \code{assocs} }
- \item{tree2}{ A phylo tree object or a phylogenetic covariance matrix for the columns of \code{assocs}}
- \item{covars1}{ A matrix of covariates (e.g., traits) for the row species of \code{assocs} }
- \item{covars2}{ A matrix of covariates (e.g., traits) for the column species of \code{assocs} }
- \item{bootstrap}{ logical, bootstrap confidence intervals of the parameter estimates }
- \item{nreps}{ Number of bootstrap replicated data sets to estimate parameter CIs }
- \item{maxit}{ as in \code{\link{optim}} }
- \item{pstart}{ starting values of the two phylogenetic signal strength parameters passed to \code{\link{optim}} }
- \item{x}{ object of class \code{pblm} }
- \item{tree1.w.novel}{ A phylo tree object or a phylogenetic covariance matrix which corresponds to \code{tree1} of \code{x} with species to predict associations }
- \item{tree2.w.novel}{ A phylo tree object or a phylogenetic covariance matrix which corresponds to \code{tree2} of \code{x} with species to predict associations }
- \item{predict.originals}{ if \code{TRUE} then the associations of each original species in the two phylogenies is predicted }
-}
-
-\details{
- Fit a linear model with covariates using estimated generalized least squares to the association strengths between two sets of interacting species.
- Associations can be either binary or continuous. If phylogenies of the two sets of interacting species are supplied,
- two \emph{phyogenetic signal strength} parameters (\emph{d1} and \emph{d2}), one for each species set, based on an Ornstein-Uhlenbeck model of
- evolution with stabilizing selection are estimated. Values of \emph{d=1} indicate no stabilizing selection and correspond to the Brownian motion model of
- evolution; \emph{0<d<1} represents stabilizing selection; \emph{d=0} depicts the absence of phylogenetic correlation (i.e., a star phylogeny); and \emph{d>1} corresponds
- to disruptive selection where phylogenetic signal is amplified. Confidence intervals for these and the other parameters can be estimated with
- bootstrapping.
- \cr
- \cr
- The function \code{pblmpredict} predicts the associations of novel species following the methods given in appendix B of Ives and Godfray (2006).
- }
-
-\value{
- \item{MSE}{ total, full (each \emph{d} estimated), star (\emph{d=0}), and base (\emph{d=1}) mean squared errors }
- \item{signal.strength}{ two estimates of phylogenetic signal strength }
- \item{coefficients}{ estimated intercept and covariate coefficients with approximate 95 percent CIs for the three model types (full, star, base) }
- \item{CI.boot}{ 95 percent CIs for all parameters }
- \item{variates}{ matrix of model variates (can be used for plotting) }
- \item{residuals}{ matrix of residuals from the three models (full, star and base) }
- \item{predicted}{ predicted associations }
- \item{bootvalues}{ matrix of parameters estimated from the \code{nreps} bootstrap replicated data sets used to calculate CIs }
- \item{phylocovs}{ phylogenetic covariance matricies scaled by the estimated \code{d1} and \code{d2} }
- \item{cors.1}{ correlations among predicted and observed associations for species of \code{tree1}, \code{NA} if \code{predict.originals=FALSE} }
- \item{cors.2}{ correlations among predicted and observed associations for species of \code{tree2}, \code{NA} if \code{predict.originals=FALSE} }
- \item{pred.novels1}{ predicted associations for the novel speices of \code{tree1} }
- \item{pred.novels2}{ predicted associations for the novel speices of \code{tree2} }
-}
-
-
-
-\note{Covariates that apply to both species sets (e.g., sampling site) should be supplied in the covariate matrix of the set with the most species.
-\cr
-\cr
-Bootstrapping CIs is slow due to the function \code{\link{optim}} used to estimate the model parameters. See appendix A in Ives and Godfray (2006)
-for a discussion about this boostrapping procedure
-\cr
-\cr
-If \code{pblmpredict=TRUE} the function does not first remove each species in turn when predicting the associations of the original species as
-is done in Ives and Godfray (2006).}
-
-\references{Ives A.R. & Godfray H.C. (2006) Phylogenetic analysis of trophic associations. The American Naturalist, 168, E1-E14 \cr
-\cr
-Blomberg S.P., Garland T.J. & Ives A.R. (2003) Testing for phylogenetic signal in comparative data: Behavioral traits are more labile. Evolution, 57, 717-745
-}
-
-\author{ Matthew Helmus \email{mrhelmus at gmail.com} }
-\examples{
-\dontrun{
-#load example data from Ives & Godfray (2006)
-data(IvesGodfray)
-
-#net attack rate of parasitoid on host eq.4 in Ives and Godfray
-A<-(-1*log(1-IvesGodfray$interactions[,-28]/t(IvesGodfray$interactions[28])))
-
-# Make tips of the phylogenetic trees contemporaneous by extending tips
-p<-dim(IvesGodfray$host)[1]
-q<-dim(IvesGodfray$parasitoid)[1]
-host.cov.scaled<-IvesGodfray$host
-para.cov.scaled<-IvesGodfray$parasitoid
-for (i in 1:p)
-{
- host.cov.scaled[i,i]<-max(host.cov.scaled)
-}
-for (i in 1:q)
-{
- para.cov.scaled[i,i]<-max(para.cov.scaled)
-}
-
-# scale covariance matrices (this reduces numerical problems caused by
-# determinants going to infinity or zero)
- host.cov.scaled<-host.cov.scaled/(det(as.matrix(host.cov.scaled))^(1/p))
- para.cov.scaled<-para.cov.scaled/(det(as.matrix(para.cov.scaled))^(1/q))
-
-pblm.A <- pblm(sqrt(A),tree1=host.cov.scaled,tree2=para.cov.scaled)
-pblm.A$signal.strength #compare to Ives and Godfray (2006) Table 1 Line 1
-pblm.A$MSE
-}
-}
+\name{pblm}
+\alias{pblm}
+\alias{pblmpredict}
+\title{ Phylogenetic Bipartite Linear Model }
+\description{
+ Fits a linear model to the association strengths of a bipartite data set with or without phylogenetic correlation among the interacting species
+}
+\usage{
+pblm(assocs, tree1=NULL, tree2=NULL, covars1=NULL, covars2=NULL, bootstrap=FALSE,
+ nreps=10, maxit=10000, pstart=c(.5,.5))
+pblmpredict(x, tree1.w.novel=NULL, tree2.w.novel=NULL, predict.originals=FALSE)
+
+}
+
+\arguments{
+ \item{assocs}{ A matrix of association strengths among two sets of interacting species }
+ \item{tree1}{ A phylo tree object or a phylogenetic covariance matrix for the rows of \code{assocs} }
+ \item{tree2}{ A phylo tree object or a phylogenetic covariance matrix for the columns of \code{assocs}}
+ \item{covars1}{ A matrix of covariates (e.g., traits) for the row species of \code{assocs} }
+ \item{covars2}{ A matrix of covariates (e.g., traits) for the column species of \code{assocs} }
+ \item{bootstrap}{ logical, bootstrap confidence intervals of the parameter estimates }
+ \item{nreps}{ Number of bootstrap replicated data sets to estimate parameter CIs }
+ \item{maxit}{ as in \code{\link{optim}} }
+ \item{pstart}{ starting values of the two phylogenetic signal strength parameters passed to \code{\link{optim}} }
+ \item{x}{ object of class \code{pblm} }
+ \item{tree1.w.novel}{ A phylo tree object or a phylogenetic covariance matrix which corresponds to \code{tree1} of \code{x} with species to predict associations }
+ \item{tree2.w.novel}{ A phylo tree object or a phylogenetic covariance matrix which corresponds to \code{tree2} of \code{x} with species to predict associations }
+ \item{predict.originals}{ if \code{TRUE} then the associations of each original species in the two phylogenies is predicted }
+}
+
+\details{
+ Fit a linear model with covariates using estimated generalized least squares to the association strengths between two sets of interacting species.
+ Associations can be either binary or continuous. If phylogenies of the two sets of interacting species are supplied,
+ two \emph{phyogenetic signal strength} parameters (\emph{d1} and \emph{d2}), one for each species set, based on an Ornstein-Uhlenbeck model of
+ evolution with stabilizing selection are estimated. Values of \emph{d=1} indicate no stabilizing selection and correspond to the Brownian motion model of
+ evolution; \emph{0<d<1} represents stabilizing selection; \emph{d=0} depicts the absence of phylogenetic correlation (i.e., a star phylogeny); and \emph{d>1} corresponds
+ to disruptive selection where phylogenetic signal is amplified. Confidence intervals for these and the other parameters can be estimated with
+ bootstrapping.
+ \cr
+ \cr
+ The function \code{pblmpredict} predicts the associations of novel species following the methods given in appendix B of Ives and Godfray (2006).
+ }
+
+\value{
+ \item{MSE}{ total, full (each \emph{d} estimated), star (\emph{d=0}), and base (\emph{d=1}) mean squared errors }
+ \item{signal.strength}{ two estimates of phylogenetic signal strength }
+ \item{coefficients}{ estimated intercept and covariate coefficients with approximate 95 percent CIs for the three model types (full, star, base) }
+ \item{CI.boot}{ 95 percent CIs for all parameters }
+ \item{variates}{ matrix of model variates (can be used for plotting) }
+ \item{residuals}{ matrix of residuals from the three models (full, star and base) }
+ \item{predicted}{ predicted associations }
+ \item{bootvalues}{ matrix of parameters estimated from the \code{nreps} bootstrap replicated data sets used to calculate CIs }
+ \item{phylocovs}{ phylogenetic covariance matricies scaled by the estimated \code{d1} and \code{d2} }
+ \item{cors.1}{ correlations among predicted and observed associations for species of \code{tree1}, \code{NA} if \code{predict.originals=FALSE} }
+ \item{cors.2}{ correlations among predicted and observed associations for species of \code{tree2}, \code{NA} if \code{predict.originals=FALSE} }
+ \item{pred.novels1}{ predicted associations for the novel speices of \code{tree1} }
+ \item{pred.novels2}{ predicted associations for the novel speices of \code{tree2} }
+}
+
+
+
+\note{Covariates that apply to both species sets (e.g., sampling site) should be supplied in the covariate matrix of the set with the most species.
+\cr
+\cr
+Bootstrapping CIs is slow due to the function \code{\link{optim}} used to estimate the model parameters. See appendix A in Ives and Godfray (2006)
+for a discussion about this boostrapping procedure
+\cr
+\cr
+If \code{pblmpredict=TRUE} the function does not first remove each species in turn when predicting the associations of the original species as
+is done in Ives and Godfray (2006).}
+
+\references{Ives A.R. & Godfray H.C. (2006) Phylogenetic analysis of trophic associations. The American Naturalist, 168, E1-E14 \cr
+\cr
+Blomberg S.P., Garland T.J. & Ives A.R. (2003) Testing for phylogenetic signal in comparative data: Behavioral traits are more labile. Evolution, 57, 717-745
+}
+
+\author{ Matthew Helmus \email{mrhelmus at gmail.com} }
+\examples{
+\dontrun{
+#load example data from Ives & Godfray (2006)
+data(IvesGodfray)
+
+#net attack rate of parasitoid on host eq.4 in Ives and Godfray
+A<-(-1*log(1-IvesGodfray$interactions[,-28]/t(IvesGodfray$interactions[28])))
+
+# Make tips of the phylogenetic trees contemporaneous by extending tips
+p<-dim(IvesGodfray$host)[1]
+q<-dim(IvesGodfray$parasitoid)[1]
+host.cov.scaled<-IvesGodfray$host
+para.cov.scaled<-IvesGodfray$parasitoid
+for (i in 1:p)
+{
+ host.cov.scaled[i,i]<-max(host.cov.scaled)
+}
+for (i in 1:q)
+{
+ para.cov.scaled[i,i]<-max(para.cov.scaled)
+}
+
+# scale covariance matrices (this reduces numerical problems caused by
+# determinants going to infinity or zero)
+ host.cov.scaled<-host.cov.scaled/(det(as.matrix(host.cov.scaled))^(1/p))
+ para.cov.scaled<-para.cov.scaled/(det(as.matrix(para.cov.scaled))^(1/q))
+
+pblm.A <- pblm(sqrt(A),tree1=host.cov.scaled,tree2=para.cov.scaled)
+pblm.A$signal.strength #compare to Ives and Godfray (2006) Table 1 Line 1
+pblm.A$MSE
+}
+}
\keyword{univar}
\ No newline at end of file
Modified: pkg/man/pglmm.Rd
===================================================================
--- pkg/man/pglmm.Rd 2013-05-29 20:33:11 UTC (rev 243)
+++ pkg/man/pglmm.Rd 2013-05-31 18:48:07 UTC (rev 244)
@@ -87,10 +87,15 @@
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
+#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)
+# 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.
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.
+out$s
}
}
Modified: pkg/man/phylostruct.Rd
===================================================================
--- pkg/man/phylostruct.Rd 2013-05-29 20:33:11 UTC (rev 243)
+++ pkg/man/phylostruct.Rd 2013-05-31 18:48:07 UTC (rev 244)
@@ -1,47 +1,48 @@
-\name{phylostruct}
-\alias{phylostruct}
-\title{ Permutations to Test for Phylogenetic Signal in Community Composition }
-\description{ Randomize sample/community data matrices to create null distributions of given metrics
-}
-\usage{
-phylostruct(samp, tree, env=NULL, metric=c("psv","psr","pse","psc","sppregs"), null.model=c("frequency",
- "richness","independentswap","trialswap"), runs=100, it=1000, alpha=0.05, fam="binomial")
-}
-\arguments{
- \item{samp}{ community data matrix, species as columns, communities as rows }
- \item{tree}{ phylo tree object or a phylogenetic covariance matrix }
- \item{env}{ environmental data matrix }
- \item{metric}{ if \code{metric="psv"}, \code{"psr"}, \code{"pse"}, or \code{"psc"} compares the observed mean of the respective metric to a null distribution at a given alpha; if \code{metric="sppregs"} compares the three correlations produced by \code{\link{sppregs}} to null distributions }
- \item{null.model}{ permutation procedure used to create the null distribution, see \code{\link{randomizeMatrix}}}
- \item{runs}{ the number of permutations to create the distribution, a rule of thumb is (number of communities)/alpha }
- \item{it}{ the number of swaps for the independent and trial-swap null models, see \code{\link{randomizeMatrix}} }
- \item{alpha}{ probability value to compare the observed mean/correlations to a null distribution }
- \item{fam}{ as in \code{\link{sppregs}} }
-}
-
-\details{The function creates null distributions for the \code{\link{psd}} set of metrics and for the correlations of \code{\link{sppregs}} from observed community data sets.}
-
-\value{
-\item{metric}{ metric used }
-\item{null.model}{ permutation used }
-\item{runs}{ number of permutations }
-\item{it}{ number of swaps if applicable }
-\item{obs}{ observed mean value of a particular metric or the three observed correlations from \code{\link{sppregs}} }
-\item{mean.null}{ mean(s) of the null distribution(s) }
-\item{quantiles.null}{ quantiles of the null distribution(s) to compare to \code{obs}; determined by \code{alpha}}
-\item{phylo.structure}{ if \code{obs} less than (alpha/2), \code{phylo.structure="underdispersed"}; if \code{obs} greater than (1-alpha/2), \code{phylo.structure="overdispersed"}; otherwise \code{phylo.structure="random"} and NULL if \code{metric="sppregs"}}
-\item{nulls}{ null values of the distribution(s)}
-}
-
-\references{ Helmus M.R., Bland T.J., Williams C.K. & Ives A.R. (2007a) Phylogenetic measures of biodiversity. American Naturalist, 169, E68-E83
-\cr
-\cr
-Helmus M.R., Savage K., Diebel M.W., Maxted J.T. & Ives A.R. (2007b) Separating the determinants of phylogenetic community structure. Ecology Letters, 10, 917-925
-\cr
-\cr
-Gotelli N.J. (2000) Null model analysis of species co-occurrence patterns. Ecology, 81, 2606-2621}
-
-\author{ Matthew Helmus \email{mrhelmus at gmail.com} }
-\seealso{ \code{\link{psd}} ,\code{\link{sppregs}}, \code{\link{randomizeMatrix}} }
-
+\name{phylostruct}
+\alias{phylostruct}
+\title{ Permutations to Test for Phylogenetic Signal in Community Composition }
+\description{ Randomize sample/community data matrices to create null distributions of given metrics
+}
+\usage{
+phylostruct(samp, tree, env=NULL, metric=c("psv","psr","pse","psc","sppregs"),
+ null.model=c("frequency", "richness","independentswap","trialswap"),
+ runs=100, it=1000, alpha=0.05, fam="binomial")
+}
+\arguments{
+ \item{samp}{ community data matrix, species as columns, communities as rows }
+ \item{tree}{ phylo tree object or a phylogenetic covariance matrix }
+ \item{env}{ environmental data matrix }
+ \item{metric}{ if \code{metric="psv"}, \code{"psr"}, \code{"pse"}, or \code{"psc"} compares the observed mean of the respective metric to a null distribution at a given alpha; if \code{metric="sppregs"} compares the three correlations produced by \code{\link{sppregs}} to null distributions }
+ \item{null.model}{ permutation procedure used to create the null distribution, see \code{\link{randomizeMatrix}}}
+ \item{runs}{ the number of permutations to create the distribution, a rule of thumb is (number of communities)/alpha }
+ \item{it}{ the number of swaps for the independent and trial-swap null models, see \code{\link{randomizeMatrix}} }
+ \item{alpha}{ probability value to compare the observed mean/correlations to a null distribution }
+ \item{fam}{ as in \code{\link{sppregs}} }
+}
+
+\details{The function creates null distributions for the \code{\link{psd}} set of metrics and for the correlations of \code{\link{sppregs}} from observed community data sets.}
+
+\value{
+\item{metric}{ metric used }
+\item{null.model}{ permutation used }
+\item{runs}{ number of permutations }
+\item{it}{ number of swaps if applicable }
+\item{obs}{ observed mean value of a particular metric or the three observed correlations from \code{\link{sppregs}} }
+\item{mean.null}{ mean(s) of the null distribution(s) }
+\item{quantiles.null}{ quantiles of the null distribution(s) to compare to \code{obs}; determined by \code{alpha}}
+\item{phylo.structure}{ if \code{obs} less than (alpha/2), \code{phylo.structure="underdispersed"}; if \code{obs} greater than (1-alpha/2), \code{phylo.structure="overdispersed"}; otherwise \code{phylo.structure="random"} and NULL if \code{metric="sppregs"}}
+\item{nulls}{ null values of the distribution(s)}
+}
+
+\references{ Helmus M.R., Bland T.J., Williams C.K. & Ives A.R. (2007a) Phylogenetic measures of biodiversity. American Naturalist, 169, E68-E83
+\cr
+\cr
+Helmus M.R., Savage K., Diebel M.W., Maxted J.T. & Ives A.R. (2007b) Separating the determinants of phylogenetic community structure. Ecology Letters, 10, 917-925
+\cr
+\cr
+Gotelli N.J. (2000) Null model analysis of species co-occurrence patterns. Ecology, 81, 2606-2621}
+
+\author{ Matthew Helmus \email{mrhelmus at gmail.com} }
+\seealso{ \code{\link{psd}} ,\code{\link{sppregs}}, \code{\link{randomizeMatrix}} }
+
\keyword{univar}
\ No newline at end of file
Modified: pkg/man/randomizeSample.Rd
===================================================================
--- pkg/man/randomizeSample.Rd 2013-05-29 20:33:11 UTC (rev 243)
+++ pkg/man/randomizeSample.Rd 2013-05-31 18:48:07 UTC (rev 244)
@@ -6,7 +6,8 @@
Various null models for randomizing community data matrices
}
\usage{
-randomizeMatrix(samp, null.model = c("frequency", "richness", "independentswap", "trialswap"),iterations = 1000)
+randomizeMatrix(samp, null.model = c("frequency", "richness",
+ "independentswap", "trialswap"), iterations = 1000)
}
\arguments{
Modified: pkg/man/ses.mntd.Rd
===================================================================
--- pkg/man/ses.mntd.Rd 2013-05-29 20:33:11 UTC (rev 243)
+++ pkg/man/ses.mntd.Rd 2013-05-31 18:48:07 UTC (rev 244)
@@ -7,7 +7,8 @@
Standardized effect size of mean nearest taxon distances in communities. When used with a phylogenetic distance matrix, equivalent to -1 times the Nearest Taxon Index (NTI).
}
\usage{
-ses.mntd(samp, dis, null.model = c("taxa.labels", "richness", "frequency", "sample.pool", "phylogeny.pool", "independentswap", "trialswap"),
+ses.mntd(samp, dis, null.model = c("taxa.labels", "richness", "frequency",
+ "sample.pool", "phylogeny.pool", "independentswap", "trialswap"),
abundance.weighted=FALSE, runs = 999, iterations = 1000)
}
Modified: pkg/man/sppregs.Rd
===================================================================
--- pkg/man/sppregs.Rd 2013-05-29 20:33:11 UTC (rev 243)
+++ pkg/man/sppregs.Rd 2013-05-31 18:48:07 UTC (rev 244)
@@ -1,57 +1,58 @@
-\name{sppregs}
-\alias{sppregs}
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/picante -r 244
More information about the Picante-commits
mailing list