Hi Luke.
For a fixed value of lambda this is very easy. For this example, you
would just do:
library(ape)
data(bird.orders)
dat<-data.frame(y=rnorm(23),x=rnorm(23))
rownames(dat)<-bird.orders$tip.label
mat<-vcv(bird.orders,cor=TRUE)
# following Blomberg to here
lambda<-0.75 # for instance
fit1<-gls(y~x,correlation=corSymm(lambda*mat[lower.tri(mat)],fixed=TRUE),data=dat)
# compare to corPagel
fit2<-gls(y~x,correlation=corPagel(0.75,bird.orders,fixed=TRUE),data=dat)
summary(fit1)
summary(fit2)
In addition, a trick to get the full correlation/covariance matrix after
lambda transformation (given a set value of lambda) is just:
mat.lambda<-lambda*(mat-diag(diag(mat)))+diag(diag(mat))
I'm not sure whether or not this is exactly what you're looking for, but
I hope it helps.
Best, Liam
--
Liam J. Revell
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.revell/
email: liam.rev...@umb.edu
blog: http://phytools.blogspot.com
On 7/15/2011 4:52 PM, Luke Matthews wrote:
This question picks up on a thread from last week, copied below. I've also
been able to use the approach Simon presented so concisely with gene networks
that are not strictly bifurcating. This only allows for the Brownian model;
that is, no transforms of the entries supplied in the covariance matrix.
I'd like to apply a lambda-type transform to the covariance or correlation
matrix directly. Does anyone know how to co-opt corPagel for use directly with
the covariance matrix, or else an alternative implementation in an nlme
correlation structure?
Thanks
Luke Matthews
[R-sig-phylo] query: how to use an existing covariance matrix directly in a gls
procedure in R?
Simon Blomberg s.blomberg1 at
uq.edu.au<mailto:r-sig-phylo%40r-project.org?Subject=Re%3A%20%5BR-sig-phylo%5D%20query%3A%20how%20to%20use%20an%20existing%20covariance%20matrix%0A%20directly%20in%20a%20gls%20procedure%20in%20R%3F&In-Reply-To=%3C011EDAAAEA824D44AE3EC8D2562AB98402DDFE279F%40UQEXMB06.soe.uq.edu.au%3E>
Fri Jul 8 01:44:48 CEST 2011
* Previous message: [R-sig-phylo] query: how to use an existing covariance matrix
directly in a gls procedure in
R?<https://stat.ethz.ch/pipermail/r-sig-phylo/2011-July/001474.html>
* Next message: [R-sig-phylo] query: how to use an existing covariance matrix
directly in a gls procedure in
R?<https://stat.ethz.ch/pipermail/r-sig-phylo/2011-July/001477.html>
* Messages sorted by: [ date
]<https://stat.ethz.ch/pipermail/r-sig-phylo/2011-July/date.html#1475> [ thread
]<https://stat.ethz.ch/pipermail/r-sig-phylo/2011-July/thread.html#1475> [ subject
]<https://stat.ethz.ch/pipermail/r-sig-phylo/2011-July/subject.html#1475> [ author
]<https://stat.ethz.ch/pipermail/r-sig-phylo/2011-July/author.html#1475>
________________________________
I used to do this before ape existed.
Here is some example code.
library(ape)
data(bird.orders)
# some made up data
dat<- data.frame(y=rnorm(23), x=rnorm(23))
rownames(dat)<- bird.orders$tip.label
mat<- vcv(bird.orders, cor=TRUE)
fit<- gls(y~x, correlation=corSymm(mat[lower.tri(mat)], fixed=TRUE), data=dat)
#assuming ultrametric tree
# compare with corBrownian structure:
fit2<- gls(y~x, correlation=corBrownian(phy=bird.orders), data=dat)
# are the regression coefficients the same?
all.equal(coef(fit), coef(fit2))
[1] TRUE
I hope this helps,
Simon.
________________________________________
From: r-sig-phylo-bounces at
r-project.org<https://stat.ethz.ch/mailman/listinfo/r-sig-phylo> [r-sig-phylo-bounces at
r-project.org<https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>] On Behalf Of Tomlinson,
Kyle [kyle.tomlinson at wur.nl<https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>]
Sent: Friday, July 08, 2011 5:30 AM
To: 'r-sig-phylo at
r-project.org<https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>'
Subject: [R-sig-phylo] query: how to use an existing covariance matrix directly
in a gls procedure in R?
Dear Dr Paradis
I trust that I am allowed to contact you directly like this? If not, my
sincere apologies.
I have constructed the covariance matrix for a phylogenetic tree (using the
vcv() command), which i would like to use directly in a gls() procedure,
instead of using the tree directly as appears to be done in ape (because I only
use a very small part of the tree i have constructed..).
{I am doing this in order to check my own gls procedure and the gls procedure
of PHYLOGr which I dont trust.}
Do you know if this can be done? I have tried to figure it out by considering
the corStruct class options in nlme, but none of these options seems to allow
one to directly input an existing covariance matrix, and I simply dont
understand the object construction methods of R well enough to build a suitable
corStruct object myself.
I hope you can help.
sincerely
Kyle Tomlinson
Resource Ecology Group
Centre for Ecosystem Studies
Wageningen University
Droevendaalsesteeg 3a
6708 PB Wageningen
The Netherlands
Phone: +31 317 485314
Fax: +31 317 484845
email: kyle.tomlinson at
wur.nl<https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>
[[alternative HTML version deleted]]
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo at r-project.org<https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
[[alternative HTML version deleted]]
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
