Permalink
Browse files

updates to phyl.pairedttest

  • Loading branch information...
liamrevell committed Nov 28, 2017
1 parent 5b0dcd1 commit cb74c348f3a91fe27e5288db08010a3f3c96a932
Showing with 98 additions and 38 deletions.
  1. +4 −4 DESCRIPTION
  2. +1 −0 NAMESPACE
  3. +85 −29 R/phyl.pairedttest.R
  4. +8 −5 man/phyl.pairedttest.Rd
View
@@ -1,6 +1,6 @@
Package: phytools
Version: 0.6-45
Date: 2017-11-17
Version: 0.6-46
Date: 2017-11-28
Title: Phylogenetic Tools for Comparative Biology (and Other Things)
Author: Liam J. Revell
Maintainer: Liam J. Revell <liam.revell@umb.edu>
@@ -13,6 +13,6 @@ ZipData: no
Description: A wide range of functions for phylogenetic analysis. Functionality is concentrated in phylogenetic comparative biology, but also includes a diverse array of methods for visualizing, manipulating, reading or writing, and even inferring phylogenetic trees and data. Included among the functions in phylogenetic comparative biology are various for ancestral state reconstruction, model-fitting, simulation of phylogenies and data, and multivariate analysis. There are a broad range of plotting methods for phylogenies and comparative data which include, but are not restricted to, methods for mapping trait evolution on trees, for projecting trees into phenotypic space or a geographic map, and for visualizing correlated speciation between trees. Finally, there are a number of functions for reading, writing, analyzing, inferring, simulating, and manipulating phylogenetic trees and comparative data not covered by other packages. For instance, there are functions for randomly or non-randomly attaching species or clades to a phylogeny, for estimating supertrees or consensus phylogenies from a set, for simulating trees and phylogenetic data under a range of models, and for a wide variety of other manipulations and analyses that phylogenetic biologists might find useful in their research.
License: GPL (>= 2)
URL: http://github.com/liamrevell/phytools
Packaged: 2017-11-17 12:00:00 EST
Packaged: 2017-11-28 12:00:00 EST
Repository:
Date/Publication: 2017-11-17 12:00:00 EST
Date/Publication: 2017-11-28 12:00:00 EST
View
@@ -130,6 +130,7 @@ S3method(logLik, pgls.Ives)
S3method(print, phyl.cca)
S3method(print, anc.Bayes)
S3method(plot, anc.Bayes)
S3method(print, phyl.pairedttest)
importFrom(animation, ani.options, ani.record, ani.replay, saveVideo)
importFrom(ape, .PlotPhyloEnv, .uncompressTipLabel, ace, all.equal.phylo, as.DNAbin, as.phylo, bind.tree, branching.times, collapse.singles)
View
@@ -1,15 +1,20 @@
# function for phylogenetic paired t-test (Lindenfors et al. 2010)
# written by Liam Revell 2011, 2013, 2015
## function for phylogenetic paired t-test (Lindenfors et al. 2010)
## written by Liam Revell 2011, 2013, 2015
phyl.pairedttest<-function(tree,x1,x2=NULL,se1=NULL,se2=NULL,lambda=1.0,h0=0.0,fixed=FALSE){
# check tree
if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
# convert x1 to a matrix if necessary
phyl.pairedttest<-function(tree,x1,x2=NULL,se1=NULL,se2=NULL,lambda=1.0,
h0=0.0,fixed=FALSE,...){
if(hasArg(tol)) tol<-list(...)$tol
else tol<-1e-12
## check tree
if(!inherits(tree,"phylo"))
stop("tree should be an object of class \"phylo\".")
## convert x1 to a matrix if necessary
if(is.data.frame(x1)) x1<-as.matrix(x1)
# if x2 is NULL
## if x2 is NULL
if(is.null(x2)){
# check to see that x1 has two columns
if(dim(x1)[2]!=2) stop("user must provide two data vectors or matrix with two variables")
## check to see that x1 has two columns
if(dim(x1)[2]!=2)
stop("user must provide two data vectors or matrix with two variables")
else {
x2<-x1[,2]
x1<-x1[,1]
@@ -26,32 +31,83 @@ phyl.pairedttest<-function(tree,x1,x2=NULL,se1=NULL,se2=NULL,lambda=1.0,h0=0.0,f
v2<-rep(0,length(tree$tip))
names(v2)<-tree$tip.label
} else v2<-se2^2
# compute C and sort tree$tip.label
## compute C and sort vectors
C<-vcv.phylo(tree)
x1<-x1[tree$tip.label]
x2<-x2[tree$tip.label]
v1<-v1[tree$tip.label]
v2<-v2[tree$tip.label]
V.diff<-diag(v1+v2); dimnames(V.diff)<-list(tree$tip.label,tree$tip.label)
# compute difference
x1<-x1[rownames(C)]
x2<-x2[rownames(C)]
v1<-v1[rownames(C)]
v2<-v2[rownames(C)]
V.diff<-diag(v1+v2)
dimnames(V.diff)<-list(names(v1),names(v1))
## compute difference
d<-x1-x2
# lambda transformation
## lambda transformation
lambda.transform<-function(C,lambda) lambda*(C-diag(diag(C)))+diag(diag(C))
# likelihood function
likelihood<-function(theta,d,C,V.diff){
## likelihood function
likelihood<-if(fixed) function(theta,d,C,V.diff,fixed.lambda){
theta[1]->sig2
theta[2]->dbar
V<-sig2*lambda.transform(C,fixed.lambda)+V.diff
logL<--t(d-dbar)%*%solve(V,d-dbar)/2-determinant(V)$modulus[1]/2-
length(d)*log(2*pi)/2
-logL[1,1]
} else function(theta,d,C,V.diff){
theta[1]->sig2
theta[2]->lambda
theta[3]->dbar
V<-sig2*lambda.transform(C,lambda)+V.diff
logL<-as.numeric(-t(d-dbar)%*%solve(V,d-dbar)/2-determinant(V)$modulus[1]/2-length(d)*log(2*pi)/2)
return(logL)
logL<--t(d-dbar)%*%solve(V,d-dbar)/2-determinant(V)$modulus[1]/2-
length(d)*log(2*pi)/2
-logL[1,1]
}
## rescale for optimization:
dscale<-1/sqrt(mean(pic(d,multi2di(tree))^2))
d<-d*dscale
V.diff<-V.diff*dscale^2
## maximize the likelihood
if(!fixed) fit<-optim(c(mean(pic(d,multi2di(tree))^2),lambda,mean(d)),likelihood,
d=d,C=C,V.diff=V.diff,method="L-BFGS-B",lower=c(tol,0,-Inf),upper=c(Inf,1,Inf),
hessian=TRUE)
else fit<-optim(c(mean(pic(d,multi2di(tree))^2),mean(d)),likelihood,d=d,C=C,V.diff=V.diff,
fixed.lambda=lambda,method="L-BFGS-B",lower=c(tol,-Inf),upper=c(Inf,Inf),hessian=TRUE)
## run t-test
se.dbar<-if(fixed) sqrt(1/fit$hessian[2,2]) else sqrt(1/fit$hessian[3,3])
t<-if(fixed) (fit$par[2]-h0)/se.dbar else (fit$par[3]-h0)/se.dbar
P<-2*pt(abs(t),df=Ntip(tree)-if(fixed) 2 else 3,lower.tail=FALSE)
obj<-list(dbar=if(fixed) fit$par[2]/dscale else fit$par[3]/dscale,
se=se.dbar/dscale,sig2=fit$par[1]/(dscale^2),
lambda=if(fixed) lambda else fit$par[2],
logL=if(fixed) -likelihood(c(fit$par[1]/(dscale^2),fit$par[2]/dscale),
d/dscale,C,V.diff/(dscale^2),lambda) else -likelihood(c(fit$par[1]/(dscale^2),
fit$par[2],fit$par[3]/dscale),d/dscale,C,V.diff/(dscale^2)),t.dbar=t,
P.dbar=P,df=Ntip(tree)-if(fixed) 2 else 3,h0=h0)
class(obj)<-"phyl.pairedttest"
obj
}
print.phyl.pairedttest<-function(x,...){
if(hasArg(digits)) digits<-list(...)$digits
else digits<-6
rnd<-function(x,digits){
x<-if(abs(x)>1e5) format(x,scientific=TRUE,digits=digits) else
if(abs(x)<=1e5&&abs(x)>1e-5) format(x,digits=digits) else
if(x<=1e-5&&x!=0) format(x,scientific=TRUE,digits=digits) else
format(x,digits=digits)
x
}
# maximize the likelihood
if(!fixed) res=optim(c(mean(pic(d,multi2di(tree))^2),lambda,h0),likelihood,d=d,C=C,V.diff=V.diff,method="L-BFGS-B",lower=c(1e-8,0,-Inf),upper=c(Inf,1,Inf),hessian=TRUE,control=list(fnscale=-1))
else res=optim(c(mean(pic(d,multi2di(tree))^2),lambda,h0),likelihood,d=d,C=C,V.diff=V.diff,method="L-BFGS-B",lower=c(1e-8,lambda-1e-8,-Inf),upper=c(Inf,lambda,Inf),hessian=TRUE,control=list(fnscale=-1))
# test
se.dbar<-sqrt(-1/res$hessian[3,3])
t<-(res$par[3]-h0)/se.dbar
P<-2*pt(abs(t),df=length(tree$tip)-3,lower.tail=F)
return(list(dbar=res$par[3],se=se.dbar,sig2=res$par[1],lambda=round(res$par[2],7),logL=res$value,t.dbar=t,P.dbar=P))
cat("\nPhylogenetic paired t-test:\n\n")
cat(paste(" t = ",rnd(x$t.dbar,digits),", df = ",x$df,
", p-value = ",rnd(x$P.dbar,digits),sep=""))
cat(paste(
"\n\nalternative hypothesis:\n true difference in means is not equal to",
rnd(x$h0,digits)))
cat("\n\n95 percent confidence interval on the phylogenetic\ndifference in mean:\n")
cat(paste(" [",paste(sapply(x$dbar+c(-1.96,1.96)*x$se,rnd,digits=digits),
collapse=", "),"]\n",sep=""))
cat("\nestimates:\n")
cat(" phylogenetic mean difference =",rnd(x$dbar,digits),"\n")
cat(" sig^2 of BM model =",rnd(x$sig2,digits),"\n")
cat(" lambda (fixed or estimated) =",rnd(x$lambda,digits),"\n\n")
cat("log-likelihood:\n")
cat(paste(" ",rnd(x$logL,digits),"\n\n"))
}
View
@@ -3,33 +3,36 @@
\title{Phylogenetic paired t-test}
\usage{
phyl.pairedttest(tree, x1, x2=NULL, se1=NULL, se2=NULL, lambda=1.0, h0=0.0,
fixed=FALSE)
fixed=FALSE, ...)
}
\arguments{
\item{tree}{a phylogeny as an object of class \code{"phylo"}.}
\item{x1}{data vector for first trait, or matrix with two traits in columns.}
\item{x2}{data vector for second trait (or null if \code{x1} is a matrix).}
\item{se1}{standard errors for \code{x1}.}
\item{se2}{standard errors for \code{x2}.}
\item{lambda}{starting value for Pagel's lambda (or fixed value, if \code{fixed=TRUE}).}
\item{lambda}{starting value for Pagel's \lambda (or fixed value, if \code{fixed=TRUE}).}
\item{h0}{null hypothesis (to be tested) for the mean difference between \code{x1} and \code{x2}.}
\item{fixed}{logical value specifying whether or not to optimize lambda.}
\item{...}{optional arguments.}
}
\description{
This function conducts a phylogenetic paired t-test, roughly following Lindenfors et al. (2010; \emph{J. Evol. Biol.}). This is not a phylogenetic ANOVA, in which we want to compare the means of different sets of species on the tree. Instead, we are interested in the difference between two characters, or two measures of a character within a species, and we want to know if this difference is significantly different from zero controlling for the phylogenetic non-independence of species.
}
\details{
Likelihood optimization is performed using \code{\link{optim}} with \code{method="L-BFGS-B"} with box constraints on lambda (0,1).
Likelihood optimization is performed using \code{\link{optim}} with \code{method="L-BFGS-B"} with box constraints on \lambda (0,1).
}
\value{
A list with the following components:
An object of class \code{"phyl.pairedttest"} with the following components:
\item{dbar}{phylogenetic mean difference.}
\item{se}{standard error of \code{dbar}.}
\item{sig2}{estimated evolutionary variance (of the difference).}
\item{lambda}{fitted (or fixed) value of lambda.}
\item{lambda}{fitted (or fixed) value of \lambda.}
\item{logL}{log-likelihood of the fitted model.}
\item{t.dbar}{t-value (\code{(dbar-h0)/se} where \code{se} is computed from the Hessian).}
\item{P.dbar}{P-value.}
\item{df}{the degrees of freedom.}
\item{h0}{the null hypothesis that was tested.}
}
\references{
Lindenfors, P., L. J. Revell, and C. L. Nunn (2010) Sexual dimorphism in primate aerobic capacity: A phylogenetic test. \emph{J. Evol. Biol.}, \bold{23}, 1183-1194.

0 comments on commit cb74c34

Please sign in to comment.