A good long time ago now some co-authors & I (Lindenfors et al. 2010) published a 'phylogenetic paired t-test.' This is a really simple method in which we ask if a set of paired traits (say, male and female RBC counts; or summer & winter body mass) differ significantly one from the other across species, while taking into account non-independence due to the phylogeny. The method can also account for (known) uncertainty in either or both trait vectors, and can relax the assumption of pure Brownian motion for the correlation structure of the errors (using Pagel's λ model).

The method has been implemented in the function `phyl.pairedttest`

;
however a *phytools* user recently pointed out some issues with
optimizing the model (which uses likelihood). This turned out to be due to
trait scale. For instance, for (numerically) large mean differences between
the two trait vectors, the corresponding element of the Hessian matrix from
the likelihood surface was sometimes numerically indistinguishable from zero.

I fixed this by using an internal rescaling, which also then required
that I back-transform the estimated parameter by (some function of) the
scaling factor. I was also able to make other improvements and I added
an S3 `print`

method for the object class. All of my updates can
be seen
here.

The following is a simple demo using simulated data:

```
library(phytools)
packageVersion("phytools")
```

```
## [1] '0.6.46'
```

```
phylo.heatmap(tree,cbind(x1,x2,x3),fsize=c(0.5,0.9,0.9),
standardize=TRUE)
```

```
phyl.pairedttest(tree,x1,x2)
```

```
##
## Phylogenetic paired t-test:
##
## t = 0.52022, df = 57, p-value = 0.604925
##
## alternative hypothesis:
## true difference in means is not equal to 0
##
## 95 percent confidence interval on the phylogenetic
## difference in mean:
## [-0.999187, 1.72124]
##
## estimates:
## phylogenetic mean difference = 0.361026
## sig^2 of BM model = 0.879318
## lambda (fixed or estimated) = 1
##
## log-likelihood:
## -78.3204
```

```
phyl.pairedttest(tree,x1,x3)
```

```
##
## Phylogenetic paired t-test:
##
## t = -2.44073, df = 57, p-value = 0.0177879
##
## alternative hypothesis:
## true difference in means is not equal to 0
##
## 95 percent confidence interval on the phylogenetic
## difference in mean:
## [-2.95937, -0.323278]
##
## estimates:
## phylogenetic mean difference = -1.64132
## sig^2 of BM model = 0.826258
## lambda (fixed or estimated) = 0.998969
##
## log-likelihood:
## -77.2733
```

Note that the method should yield (nearly approximately) the same values
as a *standard* paired t-test, if we fix λ to zero. The only
difference is that we compute the standard error of the mean difference
from the Hessian matrix instead of analytically:

```
phyl.pairedttest(tree,x1,x2,lambda=0,fixed=TRUE)
```

```
##
## Phylogenetic paired t-test:
##
## t = 1.33211, df = 58, p-value = 0.188035
##
## alternative hypothesis:
## true difference in means is not equal to 0
##
## 95 percent confidence interval on the phylogenetic
## difference in mean:
## [-0.129139, 0.677096]
##
## estimates:
## phylogenetic mean difference = 0.273978
## sig^2 of BM model = 0.720126
## lambda (fixed or estimated) = 0
##
## log-likelihood:
## -113.078
```

```
t.test(x1,x2,paired=TRUE,var.equal=TRUE)
```

```
##
## Paired t-test
##
## data: x1 and x2
## t = 1.321, df = 59, p-value = 0.1916
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.1410439 0.6890009
## sample estimates:
## mean of the differences
## 0.2739785
```

In this example, the data were simulated with a genuine difference between
`x1`

and `x3`

, but not between `x1`

and
`x2`

as follows:

```
tree<-pbtree(n=60)
x1<-fastBM(tree)
x2<-x1+fastBM(tree)
x3<-x1+fastBM(tree,a=2)
```

That's it.

## No comments:

## Post a Comment