Hi Pascal, If you'd like to get a better sense of confidence in the estimates, I'd suggest bootstrapping. ouch has a bootstrap function you can apply directly to your model fit (h3), though it will probably require quite a bit of time to run...
Carl On Tue, Jan 10, 2012 at 3:51 PM, Pascal Title <[email protected]> wrote: > 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 > > -- Carl Boettiger UC Davis http://www.carlboettiger.info/ [[alternative HTML version deleted]] _______________________________________________ R-sig-phylo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
