Manuela and David,
perhaps you want to try the OU-transformation syntax I wrote to
translate the MATLAB code developed by Blomberg, Ives & Garland. It
works well with small phylogenies, I never tried it in large ones and am
actually curious to see what would be the outcome.
You will find the details below.
Hope this helps.
Cheers,
Enrico
**************************************************************************************************************************
# This function estimates the parameter d for OU-transformation
(Blomberg et al. 2003).
# Results are identical to those obtained from the Matlab module
RegressionV2 (Lavin et al. 2008) when no independent variables are
included in the model.
# To obtain the optimal branch lengths after OU-transformation, see
phylo.OU.matrix.
# Plot command shows the REML curve as a funtion of d (default plot=T).
# IMPORTANT! In some instances d does not converge because some internal
branches may collapse before the optimal transformation is found (this
is very common with trees obtained with the rcoal command). When this
occurs, other initial branch lengths should be employed. Note that the
crashes coincide with those observed in MATLAB, so it is not a problem
with the code.
library(ape)
'phylo.OU' <-
function(trait,phy,plot=T)
{
V1 <- vcv.phylo(phy)
nspp1 <- dim(V1)[1]
V1 <- V1/det(V1)^(1/nspp1)
initV1 <- V1
U <- cbind(rep(1, nspp1))
tau1 <- matrix(diag(initV1), nspp1, nspp1) +
t(matrix(diag(initV1),nspp1, nspp1)) - 2 * initV1
pgls.ML <- function(d1) {
V1 <- (d1^tau1) * (1 - d1^(2 * initV1))/(1 - d1^2)
V1 <- V1/det(V1)^(1/nspp1)
invV <- solve(V1,diag(nspp1))
a <- solve((t(U) %*% invV %*% U), ((t(U) %*% invV
%*% trait)))
E <- (trait - U %*% a)
s2 <- (t(E) %*% invV %*% E)/(nspp1 - 1)
ML<- (-0.5*nspp1*log(2*pi)) -
(0.5*(nspp1*log((nspp1-1)*s2/nspp1) + log(det(V1)) + nspp1))
-ML
}
est <- optimize(pgls.ML,lower=0,upper=3)
d.ML <- est$minimum
ML <- -(est$objective)
pgls.REML <- function(d1) {
V1 <- (d1^tau1) * (1 - d1^(2 * initV1))/(1 - d1^2)
V1 <- V1/det(V1)^(1/nspp1)
invV <- solve(V1,diag(nspp1))
a <- solve((t(U) %*% invV %*% U), ((t(U) %*% invV
%*% trait)))
E <- (trait - U %*% a)
s2 <- (t(E) %*% invV %*% E)/(nspp1 - 1)
REML <- (-0.5*(nspp1-1)*log(2*pi)) +
(0.5*log(det(t(U)%*%U))) - ((0.5*((nspp1 - 1)*log(s2) +
log(det(V1))+log(det((t(U) %*% invV %*% U)))+(nspp1 - 1))))
-REML
}
est1 <- optimize(pgls.REML,lower=0,upper=3)
d.REML <- est1$minimum
REML <- -(est1$objective)
if (plot){
## Plotting all the data for REML
XX <- matrix(1:200/100+0.001,200,2)
for (i in 1:200){
d1 <- XX[i,1]
V1 <- V1/det(V1)^(1/nspp1)
V1 <- (d1^tau1) * (1 - d1^(2 * initV1))/(1 - d1^2)
invV <- solve(V1,diag(nspp1))
a <- solve((t(U) %*% invV %*% U), ((t(U) %*% invV %*% trait)))
E <- (trait - U %*% a)
s2 <- (t(E) %*% invV %*% E)/(nspp1 - 1)
REML.plot <- (-0.5*(nspp1-1)*log(2*pi)) + (0.5*log(det(t(U)%*%U))) -
((0.5*((nspp1 - 1)*log(s2) + log(det(V1))+log(det((t(U) %*% invV %*%
U)))+(nspp1 - 1))))
XX[i,2]<- REML.plot}
plot(XX[,1],XX[,2],xlab='OU-transformation parameter (d)',ylab='REML')}
else{
}
vcv <- (d.REML^tau1) * (1 - d.REML^(2 * initV1))/(1 - d.REML^2)
data.frame(d.REML,REML, d.ML,ML)
}
**************************************************************************************************************************
David Bapst escribió:
Hi Manuela,
I just ran the fitContinuous last night (lambda model) on a dataset of
1766 graptoloid species and their geologic durations. It took it more
than 12 hours to resolve, but it did come up with an answer, finally.
So, maybe your analysis was not permanently stuck.
However, I have seen fitContinuous get permanently 'stuck' on much
smaller trees (103 graptolites), particularly when fitting the OU
model. I once let it sit for two weeks but nothing happened. (I would
be interested in knowing if there is a way to avoid this, as I run my
analyses over large sets of trees...)
-Dave Bapst, UChicago
On Mon, Nov 15, 2010 at 3:40 AM, Manuela Gonzalez
<manuela.gonzalez.sua...@gmail.com> wrote:
Hello,
I am trying to estimate the phylogenetic "inertia" in a variable
reflecting counts for 3975 mammalian species (with an ultrametric and
fully resolved phylogeny).
I can successfully obtain an estimate of K (as in Blomberg et al 2003)
using the command phylosignal (picante package), but I am having trouble
getting an estimate of lambda with fitContinuous from Geiger package. I
am using this code.
/
lambdaFit=fitContinuous(prunedtree, m1$total_bino, model="lambda") /
The command seems to run fine (no error messages) but it just keeps
running and I am not sure if it is really working. I just stopped R
after more than 60 hours (still running but no results). Is it possible
that running this command for nearly 4000 species is just too much? or
should I be (more) patience?
Any help will be greatly appreciated.
Best regards,
Manuela
--
"I have never let my schooling interfere with my education"
Mark Twain
Manuela Gonzalez Suarez, Ph.D.
Department of Conservation Biology
Estación Biológica de Doñana CSIC
Calle Américo Vespucio s/n,
41092 Sevilla, Spain
http://www.ebd.csic.es/carnivoros/personal/gonzalez/
[[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
--
*******************************************************************
Enrico L. Rezende
Departament de Genètica i de Microbiologia
Facultat de Biociències, Edifici Cn
Universitat Autònoma de Barcelona
08193 Bellaterra (Barcelona)
SPAIN
Telephone: +34 93 581 4705
Fax: +34 93 581 2387
E-mail: enrico.reze...@uab.cat
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
