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 <[email protected]> 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 <[email protected]>wrote:
>
>> 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
>> here<http://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=1000000)
>>
>> 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
>> [email protected]
>> 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 mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
