Hi, Rafael et al. ?OUwie has some information on the standard errors:
The Hessian matrix is used as a means to estimate the approximate standard > errors of the model parameters and to assess whether they are the maximum > likelihood estimates. The variance-covariance matrix of the estimated > values of alpha and sigma^2 are computed as the inverse of the Hessian > matrix and the standard errors are the square roots of the diagonals of > this matrix. The Hessian is a matrix of second-order derivatives and is > approximated in the R package numDeriv. So, if changes in the value of a > parameter results in sharp changes in the slope around the maximum of the > log-likelihood function, the second-order derivative will be large, the > standard error will be small, and the parameter estimate is considered > stable. On the other hand, if the second-order derivative is nearly zero, > then the change in the slope around the maximum is also nearly zero, > indicating that the parameter value can be moved in any direction without > greatly affecting the log-likelihood. In such situations, the standard > error of the parameter will be large. > > For models that allow alpha and sigma^2 to vary (i.e., OUMV, OUMA, and > OUMVA), the complexity of the model can often times be greater than the > information that is contained within the data. As a result one or many > parameters are poorly estimated, which can cause the function to return a > log-likelihood that is suboptimal. This has great potential for poor model > choice and incorrect biological interpretations. An eigendecomposition of > the Hessian can provide an indication of whether the search returned the > maximum likelihood estimates. If all the eigenvalues of the Hessian are > positive, then the Hessian is positive definite, and all parameter > estimates are considered reliable. However, if there are both positive and > negative eigenvalues, then the objective function is at a saddlepoint and > one or several parameters cannot be estimated adequately. One solution is > to just fit a simpler model. Another is to actually identify the offending > parameters. This can be done through the examination of the eigenvectors. > The row order corresponds to the entries in index.matrix, the columns > correspond to the order of values in eigval, and the larger the value of > the row entry the greater the association between the corresponding > parameter and the eigenvalue. Thus, the largest values in the columns > associated with negative eigenvalues are the parameters that are causing > the objective function to be at a saddlepoint. > You can get better (by which I mean more accurate) estimates by parametric bootstrapping (simulation). Doing a quick summary of your OUM results (code at end): Trait: contr nonfor: 0.12 (-0.12, 0.36) inter: 0.12 (-0.12, 0.36) fores: 0.13 (-0.13, 0.39) Phylogenetic half life: 0.06 Tree height: 28.34 Trait: cshd nonfor: 0.24 (-0.23, 0.72) inter: 0.33 (-0.32, 0.98) fores: 0.2 (-0.19, 0.6) Phylogenetic half life: 6.35 Tree height: 28.34 Trait: dors nonfor: 0.04 (-0.04, 0.12) inter: 0.15 (-0.14, 0.44) fores: 0 (0, 0.01) Phylogenetic half life: 3.68 Tree height: 28.34 Trait: vent nonfor: 0.05 (-0.05, 0.14) inter: 0.32 (-0.3, 0.93) fores: -0.03 (0.03, -0.08) Phylogenetic half life: 5.43 Tree height: 28.34 It does seem that you don't have a lot of ability to tell optima apart. The alpha value is pretty low, though not tremendously (the half life is a good way to look at that). But the raw data also suggest it's not three very clear clumps (beeswarm plot, regimes on the x axis, trait 1 value on the y). It's not phylogenetically corrected, but still a good gut check: I agree that for your data, it may be that simpler models work better. It could also explain why you're getting crazy theta values for models with V or A allowed to vary: without strong attraction, it's hard to figure out what the attraction is towards, and especially for the third regime, with only 11 data points, there probably isn't a lot of info to estimate many parameters for just that regime. Best, Brian library(OUwie) library(ape) library(beeswarm) data <- read.csv("~/Downloads/OUM_PARAMS_ant_f.csv", row.names=1) load("~/Downloads/examplemodel.Rdata") traits <- colnames(data) regimes <- c("nonfor", "inter", "fores") observations <- list() for (trait.index in sequence(length(traits))) { cat(paste("\n\nTrait:", traits[trait.index])) for (regime.index in sequence(length(regimes))) { point.estimate <- data[paste0("theta.",regimes[regime.index]), traits[trait.index]] se <- data[paste0("theta.",regimes[regime.index]), traits[trait.index]] result.vector <- round(c(point.estimate,point.estimate-1.96*se, point.estimate+1.96*se),2) cat(paste0("\n",regimes[regime.index], ": ", result.vector[1], " (", result.vector[2], ", ", result.vector[3], ")")) if (trait.index==1) { observations[[regime.index]] <- bla$data[which(bla$data[,1]==regime.index),2] } } cat(paste0("\nPhylogenetic half life: ", round(log(2)/data["alpha", traits[trait.index]], 2))) cat(paste0("\nTree height: ", round(max(ape::branching.times(bla$phy)),2))) } points.to.plot <- bla$data[,1:2] colnames(points.to.plot) <- c("regime", "trait1") beeswarm(trait1~regime, data=points.to.plot, col="black", pch=20) _______________________________________________________________________ Brian O'Meara, http://www.brianomeara.info, especially Calendar <http://brianomeara.info/calendars/omeara/>, CV <http://brianomeara.info/cv/>, and Feedback <http://brianomeara.info/teaching/feedback/> Associate Professor, Dept. of Ecology & Evolutionary Biology, UT Knoxville Associate Head, Dept. of Ecology & Evolutionary Biology, UT Knoxville Associate Director for Postdoctoral Activities, National Institute for Mathematical & Biological Synthesis <http://www.nimbios.org> (NIMBioS) On Wed, Apr 4, 2018 at 4:07 PM, Jacob Berv <jakeberv.r.sig.ph...@gmail.com> wrote: > Even if those models may fit the data better, they may not necessarily > help Rafael determine whether or not the parameter estimates of interest > are different across regimes (though perhaps BMS might be informative). > > Rafael, maybe you could try fixing the ancestral regimes to match most > likely states for each node from your SIMMAP summary? I wonder if that > might decrease your ‘uncertainty’ in parameter estimates. > > I don’t have a great answer for your main question though, which is a good > one. > > Jake > > > On Apr 4, 2018, at 8:59 PM, William Gearty <wgea...@stanford.edu> wrote: > > > > Hi Rafael, > > > > Have you tried running the BM1, BMS, or OU1 models? If so, how do the > AICc > > values for those compare to those of the OUM model? If not, you should > make > > sure you run those. > > If you find that the these models fit your data better, this could > indicate > > that there isn't a large different across regimes and an OUM model isn't > > necessarily suitable. > > > > Best, > > Will > > > > On Wed, Apr 4, 2018 at 12:30 PM, Rafael S Marcondes < > raf.marcon...@gmail.com > >> wrote: > > > >> Dear all, > >> > >> I'm writing (again!) to ask for help interpreting standard errors of > >> parameter estimates in OUwie models. > >> > >> I'm using OUwie to examine how the evolution of bird plumage color > varies > >> across habitat types (my selective regimes) in a tree of 229 tips. I was > >> hoping to be able to make inferences based on OUMV and OUMVA models, > but I > >> was getting nonsensical theta estimates from those. So I've basically > given > >> up on them for now. > >> > >> But even looking at theta estimates from OUM models, I'm getting really > >> large standard errors, often overlapping the estimates from other > selective > >> regimes. So I was wondering what that means exactly. How are these erros > >> calculated? How much do high errors it limit the biological inferences I > >> can make? I'm more interested in the relative thetas across regimes > than on > >> the exact values (testing the prediction that birds in darker habitats > tend > >> to adapt to darker plumages than birds in more illuminated habitats). > >> > >> I have attached a table averaging parameter estimates and errors from > >> models fitted across a posterior distribution of 100 simmaps for four > >> traits; and one exemplar fitted model from one trait in one of those > >> simmaps. > >> > >> Thanks a lot for any help, > >> > >> > >> *--* > >> *Rafael Sobral Marcondes* > >> PhD Candidate (Systematics, Ecology and Evolution/Ornithology) > >> > >> Museum of Natural Science <http://sites01.lsu.edu/wp/mns/> > >> Louisiana State University > >> 119 Foster Hall > >> Baton Rouge, LA 70803, USA > >> > >> Twitter: @brown_birds <https://twitter.com/brown_birds> > >> > >> > >> _______________________________________________ > >> 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-ph...@r-project.org/ > >> > >> > > > > > > -- > > William Gearty > > PhD Candidate, Paleobiology > > Department of Geological Sciences > > Stanford School of Earth, Energy & Environmental Sciences > > williamgearty.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-ph...@r-project.org/ > > _______________________________________________ > 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-ph...@r-project.org/ > [[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/