[R-sig-phylo] overall p-value from multiple PGLS regression
Hi all,
I am running a number of PGLS regressions, some of which are multiple
regressions. I am using the nlme package with the corBrownian error
structure. If I build a model M via multiple regression, then I can get a
summary of the model that includes p values as well as a number of other
statistics. But how can I get an overall p-value for this model?
In the same vein, as value, std. error, t-value and p-value are all
returned for each of the traits separately, what does one report for
multiple regression models? I have been unable to find a paper that
actually reported this type of information for such a model.
Here is a toy example:
require(nlme)
require(phytools)
tree - pbtree(n=20)
dat - as.data.frame(fastBM(tree, nsim=4))
colnames(dat) - paste('trait', 1:4, sep='')
M - gls(trait1 ~ trait2 + trait3 + trait4, data=dat,
correlation=corBrownian(1, tree))
summary(M)
Any help/advice would be greatly appreciated. Thanks!
-Pascal
--
Pascal Title
PhD candidate, Rabosky Lab
http://www-personal.umich.edu/~drabosky/Home.html
Dept of Ecology and Evolutionary Biology
University of Michigan
pti...@umich.edu
pascaltitle.weebly.com
[[alternative HTML version deleted]]
___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
Re: [R-sig-phylo] obtaining reasonable values from OUCH
Thanks Marguerite for your reply! To answer your questions, as for where I heard about trait variance approaching sigma2/(2*alpha), I had posted to r-sig-phylo about similar problems a little over a year ago, and Carl Boettiger had mentioned this as a way to see if the parameter estimates were in ballpark range or not. The objective here was to calculate evolutionary rate for several climatic niche axes, summarized as principal component axes. Following, Adams et al. 2009 and Rabosky and Adams 2012, I would fit a multivariate model to multiple PC axes, and sum the diagonal elements from the sigma2 matrix, to get an overall rate value. So, yes, I am after specific values so that I can use them to look for correlations, etc. After some testing, I seem to be able to get good estimates for individual traits, so perhaps I can sum these, rather than sum the diagonal elements of a variance-covariance matrix. Thanks for the input! -Pascal On Tue, Jul 2, 2013 at 4:28 AM, Marguerite Butler mbutler...@gmail.comwrote: Good morning Pascal, We (Cressler, King, and myself) have a paper soon to be submitted that shows that model selection is very robust, but parameter estimates, esp. alpha and sigma are very difficult to estimate, because they have high variance. Although the bias is manageable, it is very hard to say much about a single estimate of either parameter because they can vary widely. However, model selection is very robust, so we recommend placing a lot more confidence in the model selection and not so much stake in specific values of parameter estimates. In any case, one should always assess the robustness of the estimates using parametric bootstrap or other resampling methods. Facilities for bootstraps and simulation are already built into OUCH. This is all for the univariate case, but we have no reason to expect the multivariate case to be any better and in fact is probably worse. In my experience, it is frequently the case that one can get very high values, especially for alpha. It just depends on the data and the structure of the tree. It is especially bad if you violate the model assumptions (traits that disappear for example are not particularly gaussian), and even with transformation sometimes you can't do much about it. You never describe the nature of the data. In any case, it is important to remember that it is ultimately a simple model trying to describe stochastic and deterministic aspects of evolution assuming that the stochastic part is Gaussian. However, even so, you can determine how much you can say by doing the parametric bootstraps. It may be that all you can say is alpha is greater than zero, or that alpha is really huge (but you can't say what the value is specifically). That's OK, but you do need to know how much you can conclude. My questions regarding your original post is: why are you expecting V to be proportional to 2alpha over sigma^2? In our original paper (Butler King 2004), we show that V approaches this value when alpha approaches zero (essentially when the model becomes close to BM). But again, this is very hard to estimate. Furthermore, why are you after specific parameter estimates? Does your hypothesis require alpha and sigma to be particular values? Also, since you're in Michigan, you may try asking Aaron. Hope this helps, Marguerite On Jun 29, 2013, at 3:19 PM, Brian O'Meara bome...@utk.edu wrote: No, it's just univariate (I was responding to your example about doing a single character). For multivariate OU, I think OUCH is the only game in town. Best, Brian ___ Brian O'Meara Assistant Professor Dept. of Ecology Evolutionary Biology U. of Tennessee, Knoxville http://www.brianomeara.info Students wanted: Applications due Dec. 15, annually Postdoc collaborators wanted: Check NIMBioS' website Calendar: http://www.brianomeara.info/calendars/omeara On Sat, Jun 29, 2013 at 9:11 PM, Pascal Title pti...@umich.edu wrote: Thanks for the reply, Brian. Does OUwie implement the multivariate model? From the documentation, I don't see any mention of applying it to multiple traits at the same time. cheers, -Pascal On Sat, Jun 29, 2013 at 11:32 PM, Brian O'Meara bome...@utk.edu wrote: You're at best on the border of how many parameters you can estimate given the size of your dataset. One thing you might try is running OUwie, which also implements that model but uses a different optimizer and starting values, though I wouldn't be surprised if you get similar answers. Trying a grid of values for a contour plot can tell you something about the likelihood surface you're trying to optimize on; it could be very flat. We had this function in OUwie until a license conflict with akima made us have to cut this functionality, but it's easy enough to do. Best, Brian
[R-sig-phylo] nonparametric PGLS
Hi everyone I have character data for species on a phylogeny, but simple transformations like log-transformations or square-root transformations are not proving to be sufficient to get the data to be normal. If this were a non-phylogenetic test, I would resort to something like a Spearman correlation, which is non-parametric. But with phylogenetic generalized least squares, is there any method that is nonparametric? Thanks! -Pascal Title [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Re: [R-sig-phylo] convergence issues with hansen model
Thanks for the response, Liam. I just tried out OUwie and it works.
However, it appears to only fit 1 variable at a time (I only get one sigma
squared value).
When running some tests with OUCH, I found that sigma squared is different
when one variable is fit to an OU model, versus when multiple variables are
fit together. Is that expected?
-Pascal
On Mon, Mar 12, 2012 at 11:39 AM, Liam J. Revell liam.rev...@umb.eduwrote:
You might want to try the OUWie package by Beaulieu O'Meara (
http://cran.r-project.org/**web/packages/OUwie/index.htmlhttp://cran.r-project.org/web/packages/OUwie/index.html
)**. I have not tried it yet, but it promises to do multi-optimum OU
model fitting as well.
All the best, Liam
--
Liam J. Revell
University of Massachusetts Boston
web:
http://faculty.umb.edu/liam.**revell/http://faculty.umb.edu/liam.revell/
email: liam.rev...@umb.edu
blog: http://phytools.blogspot.com
On 3/12/2012 2:31 PM, Pascal Title wrote:
Hi all
I have a dataset of 5 PC axes for a phylogeny with 52 tips and I would
like
to get the variance-covariance matrix under a global OU model of
evolution. I find that OU is strongly favored over BM, based on
fitContinuous in GEIGER.
However, I am having issues with convergence when trying to fit a hansen
model to my data. What can I do to get around this problem? Alternatively,
is there another way to get at the multivariate VCV matrix other than with
the OUCH package?
The error is:
*unsuccessful convergence, code 1, see documentation for `optim'*
*Warning message:*
*In hansen(tree = ot, data = otd[varnames], regimes = otd[regimes], :*
* unsuccessful convergence*
I have posted a small file that contains the following R script along with
the data
herehttp://dl.dropbox.com/u/**34644229/OUhansen.ziphttp://dl.dropbox.com/u/34644229/OUhansen.zip
.
Any advice would be greatly appreciated! Thanks!
require(ouch)
require(ape)
tree-read.nexus('tree.nex')
data-read.csv('data.csv')
rownames(data)-data$X
data$X-NULL
tree
head(data)
colnames(data)-varnames #for hansen()
ot-ape2ouch(tree)
otd-as(ot,data.frame)
data$labels-rownames(data)
otd-merge(otd,data,by=**labels,all=TRUE)
rownames(otd)-otd$nodes
ot-with(otd,ouchtree(nodes=**nodes,ancestors=ancestors,**
times=times,labels=labels))
otd$regimes-as.factor(**global)
h1-hansen(tree=ot,data=otd[**varnames],regimes=otd[**
regimes],sqrt.alpha=c(1,0,0,**0,0,1,0,0,0,1,0,0,1,0,1),**
sigma=c(1,0,0,0,0,1,0,0,0,1,0,**0,1,0,1))
--
Pascal Title
Graduate Student
Burns Lab
http://kevinburnslab.com/
Department of Evolutionary Biology
San Diego State University
5500 Campanile Drive
San Diego, CA 92182-4614
[[alternative HTML version deleted]]
___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Re: [R-sig-phylo] convergence issues with hansen model
That helps alot, Carl, thank you. So to bring this back to my original
question, it sounds like the only way to get sigma squared values from a
multivariate OU model is with OUCH. I'm currently not able to get results
because of convergence issues. Does that mean that I simply won't be able
to get such data? Is this the end of the line? Or are there parameters that
I can tweak to get convergence? I know you can increase the number of
iterations, and you can choose different optimization methods, but I don't
know enough about these methods to know whether or not it is appropriate to
change these, or if this could somehow affect the output.
I really appreciate the help!
-Pascal
On Mon, Mar 12, 2012 at 3:01 PM, Carl Boettiger cboet...@gmail.com wrote:
Just wanted to note that there's a fundamental difference between the
way OUCH is handling the multivariate traits and the way geiger and
similar packages are handling it. If you're giving geiger a data
frame with multiple traits, it is fitting each independently, just
saving you the trouble from passing them in one at a time.
OUCH is actually fitting the multivariate trait model, which is the
reason that you get different results fitting a single trait or
fitting multiple traits simultaneously.
In the OUCH model, dX = alpha ( theta - X) dt + sigma dB_t,
X is a vector, alpha is a matrix, theta is a vector, sigma is a
matrix, and dB_t is a vector of Brownian walks.
Only if the off-diagonal terms of alpha and sigma matrices are zero
will the results be the same.
hope that helps,
-Carl
On Mon, Mar 12, 2012 at 4:56 PM, Pascal Title pascalti...@gmail.com
wrote:
Thanks for the response, Liam. I just tried out OUwie and it works.
However, it appears to only fit 1 variable at a time (I only get one
sigma
squared value).
When running some tests with OUCH, I found that sigma squared is
different
when one variable is fit to an OU model, versus when multiple variables
are
fit together. Is that expected?
-Pascal
On Mon, Mar 12, 2012 at 11:39 AM, Liam J. Revell liam.rev...@umb.edu
wrote:
You might want to try the OUWie package by Beaulieu O'Meara (
http://cran.r-project.org/**web/packages/OUwie/index.html
http://cran.r-project.org/web/packages/OUwie/index.html
)**. I have not tried it yet, but it promises to do multi-optimum OU
model fitting as well.
All the best, Liam
--
Liam J. Revell
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.**revell/
http://faculty.umb.edu/liam.revell/
email: liam.rev...@umb.edu
blog: http://phytools.blogspot.com
On 3/12/2012 2:31 PM, Pascal Title wrote:
Hi all
I have a dataset of 5 PC axes for a phylogeny with 52 tips and I would
like
to get the variance-covariance matrix under a global OU model of
evolution. I find that OU is strongly favored over BM, based on
fitContinuous in GEIGER.
However, I am having issues with convergence when trying to fit a
hansen
model to my data. What can I do to get around this problem?
Alternatively,
is there another way to get at the multivariate VCV matrix other than
with
the OUCH package?
The error is:
*unsuccessful convergence, code 1, see documentation for `optim'*
*Warning message:*
*In hansen(tree = ot, data = otd[varnames], regimes = otd[regimes],
:*
* unsuccessful convergence*
I have posted a small file that contains the following R script along
with
the data herehttp://dl.dropbox.com/u/**34644229/OUhansen.zip
http://dl.dropbox.com/u/34644229/OUhansen.zip
.
Any advice would be greatly appreciated! Thanks!
require(ouch)
require(ape)
tree-read.nexus('tree.nex')
data-read.csv('data.csv')
rownames(data)-data$X
data$X-NULL
tree
head(data)
colnames(data)-varnames #for hansen()
ot-ape2ouch(tree)
otd-as(ot,data.frame)
data$labels-rownames(data)
otd-merge(otd,data,by=**labels,all=TRUE)
rownames(otd)-otd$nodes
ot-with(otd,ouchtree(nodes=**nodes,ancestors=ancestors,**
times=times,labels=labels))
otd$regimes-as.factor(**global)
h1-hansen(tree=ot,data=otd[**varnames],regimes=otd[**
regimes],sqrt.alpha=c(1,0,0,**0,0,1,0,0,0,1,0,0,1,0,1),**
sigma=c(1,0,0,0,0,1,0,0,0,1,0,**0,1,0,1))
--
Pascal Title
Graduate Student
Burns Lab
http://kevinburnslab.com/
Department of Evolutionary Biology
San Diego State University
5500 Campanile Drive
San Diego, CA 92182-4614
[[alternative HTML version deleted]]
___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
--
Carl Boettiger
UC Davis
http://www.carlboettiger.info/
--
Pascal Title
Graduate Student
Burns Lab
http://kevinburnslab.com/
Department of Evolutionary Biology
San Diego State University
5500 Campanile Drive
San Diego, CA 92182-4614
[[alternative HTML version deleted
[R-sig-phylo] Freckleton and Harvey randomization test
Hi all, Is the randomization test of Freckleton and Harvey, 2006, Detecting non-Brownian trait evolution in adaptive radiations from PLOS Biol. currently implemented in any R packages? I can't find it anywhere, but thought I would ask before attempting to write a script myself. Cheers, -Pascal -- Pascal Title Graduate Student Burns Lab http://kevinburnslab.com/ Department of Evolutionary Biology San Diego State University 5500 Campanile Drive San Diego, CA 92182-4614 [[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] differing sigma2 results
Hello all
I am trying to calculate sigma2 values for a clade, using 4 PCA axes as my
characters.
I first wanted to figure out if my data better fit a BM model or an OU1
model of evolution. So I compared AICc values using the fitContinuous
function in Geiger.
I found the following:
BM OU
Score1 2.6 0
Score2 9.1 0
Score3 17.2 0
Score4 9.8 0
However, using the pmc package, I ran
pmc(tree,data,modelA='BM',modelB='OU',nboot=100)
for the first PC axis and found that the BM and OU1 distributions are
mostly overlapping, which I believe tells me that you can't really
differentiate BM and OU with the data at hand (see plot of
pmchttp://dl.dropbox.com/u/34644229/Tangarinae.pdf
).
Therefore, I was planning on calculating sigma2 based on a BM model of
evolution.
But I noticed that I get hugely different results using different methods:
(here is my tree http://dl.dropbox.com/u/34644229/testtree.nex, here is
my data http://dl.dropbox.com/u/34644229/spdata.csv)
With OUCH method:
require(geiger)
require(ouch)
tree-read.nexus('testtree.nex')
spdata-read.csv('spdata.csv',stringsAsFactors=FALSE)
rownames(spdata)-spdata$X
spdata$X-NULL
spdata$labels-NULL
spdata-spdata1
ot-ape2ouch(tree)
otd-as(ot,data.frame)
spdata$labels-rownames(spdata)
otd-merge(otd,spdata,by=labels,all=TRUE)
rownames(otd)-otd$nodes
ot-with(otd,ouchtree(nodes=nodes,ancestors=ancestors,times=times,labels=labels))
b1-brown(tree=ot,data=otd[colnames(spdata1)])
sigma2-coef(b1)$sigma.sq.matrix
sigma2
[,1] [,2][,3] [,4]
[1,] 5.885605498 1.17240962 0.7760799 -0.004528573
[2,] 1.172409617 1.52861569 0.1039507 -0.073455526
[3,] 0.776079917 0.10395074 1.0567364 -0.158017246
[4,] -0.004528573 -0.07345553 -0.1580172 1.236342488
With ic.sigma method from GEIGER:
require(geiger)
tree-read.nexus('testtree.nex')
spdata-read.csv('spdata.csv',stringsAsFactors=FALSE)
rownames(spdata)-spdata$X
spdata$X-NULL
spdata$labels-NULL
sig-ic.sigma(tree,spdata)
Score1 Score2 Score3 Score4
Score1 0.1809881293 0.036052743 0.023865217 -0.0001392581
Score2 0.0360527432 0.047006429 0.003196587 -0.0022588293
Score3 0.0238652170 0.003196587 0.032495680 -0.0048591850
Score4 -0.0001392581 -0.002258829 -0.004859185 0.0380187415
Interestingly, the OUCH result is greater than the GEIGER result by a
factor of 32.
I also simulated a tree and data under BM and tried both methods and they
agree, so somehow this has something to do with my data.
Both methods are estimating these parameters under brownian motion, so
where is there a breakdown?
Any advice would be greatly appreciated!
--
Pascal Title
Graduate Student
Burns Lab
http://kevinburnslab.com/
Department of Evolutionary Biology
San Diego State University
5500 Campanile Drive
San Diego, CA 92182-4614
[[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] max iteration error in dismo maxent
Hello
I am constructing a large number of niche models to perform a background
similarity test (as is found in ENMTools), whereby I take a random sample
of points from the distribution of a species, and then build a niche model
using maxent through the dismo package.
I am using R v2.14.1 through RStudio, all packages are up to date as of
today, and I am using maxent v3.3.3k on Windows 7.
My code for that segment:
for (j in 1:100){
print(paste('spA_bgB--iteration ',j),sep='')
bg-spsample(rangeB,n=nrow(as.data.frame(occB)),type='random')
extract(predictorsB[[1]],bg)-x #to remove points with NA values in
predictors
bg[which(!is.na(x)),]-bg
xm-maxent(x=predictorsB,p=bg,args=c('randomseed=TRUE','maximumiterations=10'))
layerNames(predictorsB)-layerNames(fullExtent)
pxBG-predict(xm,fullExtent,progress='text')
This works most of the time, but occasionally I get the following error,
that brings my R script to a halt:
Error in .local(x, n, type, ...) :
iteration did not converge; try enlarging argument iter
I currently have max iterations in the maxent function at 10. I'm not
sure if the error comes from the maxent function or from the predict
function because every time I rerun either of the two functions, it works
and the error does not repeat. So it seems like simply rerunning the
maxent() and/or predict() functions is enough to keep the script moving.
I tried looking in the source code from maxent() and predict() but did not
see this error message, so I don't know how to work around it. What I'm
hoping to do is add a while loop to have maxent() and predict() keep
attempting to build the model until it works.
Can anyone help me out? I would like to know what this error means, and how
to add a while loop to bypass this error.
Thanks!
-Pascal
--
Pascal Title
Graduate Student
Burns Lab
http://kevinburnslab.com/
Department of Evolutionary Biology
San Diego State University
5500 Campanile Drive
San Diego, CA 92182-4614
[[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] very large sigma squared values
Hello I am trying to calculate sigma squared values under the Hansen OU model of evolution, so as to examine rate of evolution of certain characters. I have used PCA to reduce the dimensionality of my data and would like to get sigma squared values for the top 4 PC scores. I am using the OUCH package, and have succeeded in getting results, but the sigma squared values seem incredibly large. The PC data looks like this: head(spdata) Score1 Score2 Score3 Score4 HapUni -0.145379 -0.6747320 -0.7392243 -1.34553078 PhrPle -4.350328 0.9049652 1.0199397 0.24838370 PhrUni -3.586409 0.9624709 0.0682498 0.23511329 XenPar -3.601509 1.6425747 0.2581034 0.07927267 PhrDor -5.977901 0.5893070 1.9651834 0.54601353 PhrEry -5.707650 1.5844235 1.8148729 -0.31778757 and the resulting variance-covariance matrix looks like this: coef(h1)$sigma.sq.matrix [,1] [,2] [,3] [,4] [1,] 142.79238 112.23323 -42.86343 182.90860 [2,] 112.23323 149.15931 -22.40858 111.64635 [3,] -42.86343 -22.40858 29.09922 -26.99329 [4,] 182.90860 111.64635 -26.99329 388.98963 Am I doing this correctly? Or are these large values the result of incorrect methodology? Any advice or suggestions would be greatly appreciated! The tree and data can be downloaded herehttp://dl.dropbox.com/u/34644229/Clade3.zip . require(ouch) require(ape) #Read in tree and data setwd(/Users/Pascal/Desktop/Spec records/clade3) tree-read.nexus(clade3.nex) spdata-read.csv(clade3_pc.csv) rownames(spdata)-spdata$X spdata$X-NULL head(spdata) #Fit OU model ot-ape2ouch(tree) otd-as(ot,data.frame) spdata$labels-rownames(spdata) otd-merge(otd,spdata,by=labels,all=TRUE) rownames(otd)-otd$nodes ot-with(otd,ouchtree(nodes=nodes,ancestors=ancestors,times=times,labels=labels)) otd$regimes-as.factor(global) h3-hansen(tree=ot,data=otd[colnames(spdata)[1:4]],regimes=otd[regimes],sqrt.alpha=c(1,0,0,0,1,0,0,1,0,1),sigma=c(1,0,0,0,1,0,0,1,0,1),maxit=100) summary(h1) -Pascal Title -- Pascal Title Graduate Student Burns Lab http://kevinburnslab.com/ Department of Evolutionary Biology San Diego State University 5500 Campanile Drive San Diego, CA 92182-4614 [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Re: [R-sig-phylo] very large sigma squared values
Thanks for the quick response! At the suggestion of another reader, I made my tree ultrametric and this greatly improved the parameter estimates. Using the ultrametric tree, I did the quick calculation for one variable that Carl Boettiger suggested, and found 6.18, as compared to the variance of that trait data of 5.83. That seems pretty close. And also, after doing some model selection between OU and BM, OU was strongly favored. Now my alpha and sigma^2 values look like this. $alpha [,1] [,2] [,3] [,4] [1,] 5.948908 2.388587 -2.968473 3.285025 [2,] 2.388587 9.827974 -3.891185 -3.769924 [3,] -2.968473 -3.891185 7.708327 5.824082 [4,] 3.285025 -3.769924 5.824082 78.668861 $sigma.squared [,1] [,2] [,3] [,4] [1,] 73.528557 6.771645 -35.121305 51.29463 [2,] 6.771645 17.450129 -9.548777 -11.52932 [3,] -35.121305 -9.548777 31.007077 -13.57953 [4,] 51.294628 -11.529325 -13.579525 58.19881 Some of the values are still somewhat large, though not nearly as much as before. Does this seem more reasonable? -Pascal On Tue, Jan 10, 2012 at 2:20 PM, Carl Boettiger cboet...@gmail.com wrote: Hi Pascal, Looks like your alpha values are also very large (something to check whenever you get large sigma values). To have an idea if your estimates makes sense, you can compare against the equilibrium variance of a (single-peak) OU model as a back-of-the-envelope thing: sigma^2/2*alpha should be somewhat near the variance of your trait. You are right to be suspicious about the large values of sigma. If both sigma and alpha are going to very large values, it may mean that the system is governed by a strong OU process -- you can check with a model choice against the BM model -- in which case you may not be able to estimate alpha and sigma independently. It can also be due to your phylogeny -- a simple example is a star phylogeny, in which alpha and sigma cannot ever be estimated independently but will behave as you describe. Also check the convergence of the algorithm. Convergence doesn't guarantee that you don't have this problem of being on a simga-alpha likelihood ridge, but non-convergence certainly suggests it. Given the large number of iterations you are using for the nelder-mead optimization, this is also a sign of such a likelihood ridge. -Carl On Tue, Jan 10, 2012 at 2:00 PM, Pascal Title pascalti...@gmail.comwrote: Hello I am trying to calculate sigma squared values under the Hansen OU model of evolution, so as to examine rate of evolution of certain characters. I have used PCA to reduce the dimensionality of my data and would like to get sigma squared values for the top 4 PC scores. I am using the OUCH package, and have succeeded in getting results, but the sigma squared values seem incredibly large. The PC data looks like this: head(spdata) Score1 Score2 Score3 Score4 HapUni -0.145379 -0.6747320 -0.7392243 -1.34553078 PhrPle -4.350328 0.9049652 1.0199397 0.24838370 PhrUni -3.586409 0.9624709 0.0682498 0.23511329 XenPar -3.601509 1.6425747 0.2581034 0.07927267 PhrDor -5.977901 0.5893070 1.9651834 0.54601353 PhrEry -5.707650 1.5844235 1.8148729 -0.31778757 and the resulting variance-covariance matrix looks like this: coef(h1)$sigma.sq.matrix [,1] [,2] [,3] [,4] [1,] 142.79238 112.23323 -42.86343 182.90860 [2,] 112.23323 149.15931 -22.40858 111.64635 [3,] -42.86343 -22.40858 29.09922 -26.99329 [4,] 182.90860 111.64635 -26.99329 388.98963 Am I doing this correctly? Or are these large values the result of incorrect methodology? Any advice or suggestions would be greatly appreciated! The tree and data can be downloaded herehttp://dl.dropbox.com/u/34644229/Clade3.zip . require(ouch) require(ape) #Read in tree and data setwd(/Users/Pascal/Desktop/Spec records/clade3) tree-read.nexus(clade3.nex) spdata-read.csv(clade3_pc.csv) rownames(spdata)-spdata$X spdata$X-NULL head(spdata) #Fit OU model ot-ape2ouch(tree) otd-as(ot,data.frame) spdata$labels-rownames(spdata) otd-merge(otd,spdata,by=labels,all=TRUE) rownames(otd)-otd$nodes ot-with(otd,ouchtree(nodes=nodes,ancestors=ancestors,times=times,labels=labels)) otd$regimes-as.factor(global) h3-hansen(tree=ot,data=otd[colnames(spdata)[1:4]],regimes=otd[regimes],sqrt.alpha=c(1,0,0,0,1,0,0,1,0,1),sigma=c(1,0,0,0,1,0,0,1,0,1),maxit=100) summary(h1) -Pascal Title -- Pascal Title Graduate Student Burns Lab http://kevinburnslab.com/ Department of Evolutionary Biology San Diego State University 5500 Campanile Drive San Diego, CA 92182-4614 [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo -- Carl Boettiger UC Davis http://www.carlboettiger.info/ -- Pascal Title Graduate Student
