Aloha Rafael,

I am very glad that you are thinking about power of your analyses and 
confidence in your parameter estimates.  Folks have chimed in with many useful 
responses, however, I believe your original question was whether uncertainty in 
parameter estimates indicate bad models and limits of interpretation?

It’s important to note that model fit and parameter estimates are two different 
aspects of assessment. It is perfectly possible to have robust model selection 
(decent power) but have parameters that have a lot of uncertainty and are hard 
to identify with any precision (large confidence intervals).  That was pretty 
well demonstrated in Cressler, Butler, King 2015. The specific outcome is data 
and model dependent, but as a generalization, model selection is more robust 
than parameter estimation, and within the parameters, the position of the 
thetas are a little bit easier to narrow down than the alphas and sigmas. 

This is not too surprising when you think about the information content of the 
data. We are really trying to fit some relatively complex models with few data 
(one point per species, with phylogenetic structure to boot), and the intuitive 
interpretation is that there are many ways to the same outcome.  If the outcome 
is the particular pattern of phenotypes along the tree, it might be possible to 
get that outcome with strong selection and little stochasticity, or weak 
selection toward optima very far apart (and the observed phenotypes are 
essentially “on the way” toward the optima), or very weak selection and strong 
stochasticity (they were just knocked about that way by drift), etc.  While 
this last scenario is most different than the first two (the first two are more 
OU-type, whereas the last is more BM-like in being dominated by stochasticity), 
my point is that there may not necessarily be a unique set of parameters that 
explain the data. 

Nevertheless, it may still be possible to say something biologically useful. 
First of all, I noticed that you are doing SIMMAP which is great for exploring 
phylogenetic uncertainty, but this is a bit arms-length with regard to OU 
parameter and model uncertainty. A bit apples and oranges. I would strongly 
recommend that you try some parametric bootstrapping on your OU models for both 
model selection and parameter estimation. In OUCH, there are built-in 
facilities to make this really easy, probably there are in OUwie as well (sorry 
I only stick with OUCH :). The key for the model selection bootstrap is that 
you can simulate the data using the best-fit model, but then fit all of the 
alternative models to see how often your best model performs. I like to point 
folks to my former grad student’s paper (Scales et al 2009).  If you want to be 
even more stastically thorough you can do the Botteinger, Coop, and Ralph 2012 
approach which basically does parametric simulations on both the best and 
alternative models.   The parametric bootstraps on the parameter estimates is 
much easier. and you can just run the simulations and fit only the best model 
and draw your CI’s from that. Maybe you’ve done this already. 

With regard to what can be said, as Liam and Brian and others have said, many 
times there are large CIs because we encounter flat or ridged likelihood 
surfaces.  So that may limit your conclusions somewhat to saying something like 
we don’t know what alpha is, precisely, but we know it’s greater than zero.  
Etc.  In your specific case, what you want to know is whether theta 1 > theta 
2, or something like that, right? So without needing to know the precise 
location of the theta’s I would bet that you can look through your parametric 
simulations and see how often that condition is satisfied.  You may find that 
you don’t know what they are exactly, but one is always greater than the other. 

I love your study question by the way. I too study colorful organisms in bright 
and dark habitats. Best of luck!

Hth,
Marguerite

Cressler C., Butler M.A., and King A. A. (2015) Detecting adaptive evolution in 
phylogenetic comparative analysis using the Ornstein-Uhlenbeck model.  Sys. 
Bio. 64(6):953-968. DOI: 10.1093/sysbio/syv043

Scales J.A., King A.A., and Butler M.A. (2009) Running for your life or running 
for your dinner: What drives fiber-type evolution in lizard locomotor muscles? 
Am. Nat.173:543-553. DOI: 10.1086/597613

Boettinger C., Coop G., and P. Ralph (2012) Is your phylogeny informative? 
Measuring the power of comparative methods. Evolution 66-7: 2240-51. DOI: 
10.1111/j.1558-5646.2011.01574.x

> On Apr 4, 2018, at 10:27 AM, Brian O'Meara <omeara.br...@gmail.com> wrote:
> 
> Apparently the plot did not come through correctly. I'm attaching a PDF of it.
> 
> Best,
> Brian
> 
> _______________________________________________________________________
> Brian O'Meara, http://www.brianomeara.info <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:24 PM, Brian O'Meara <omeara.br...@gmail.com 
> <mailto:omeara.br...@gmail.com>> wrote:
> 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 <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 
> <mailto: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 
> > <mailto: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 <mailto: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/ 
> >> <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 
> >> <https://twitter.com/brown_birds>>
> >>
> >>
> >> _______________________________________________
> >> R-sig-phylo mailing list - R-sig-phylo@r-project.org 
> >> <mailto:R-sig-phylo@r-project.org>
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo 
> >> <https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>
> >> Searchable archive at http://www.mail-archive.com/r- 
> >> <http://www.mail-archive.com/r->
> >> sig-ph...@r-project.org/ <http://sig-ph...@r-project.org/>
> >>
> >>
> >
> >
> > --
> > William Gearty
> > PhD Candidate, Paleobiology
> > Department of Geological Sciences
> > Stanford School of Earth, Energy & Environmental Sciences
> > williamgearty.com <http://williamgearty.com/>
> >
> >       [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-phylo mailing list - R-sig-phylo@r-project.org 
> > <mailto:R-sig-phylo@r-project.org>
> > https://stat.ethz.ch/mailman/listinfo/r-sig-phylo 
> > <https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>
> > Searchable archive at 
> > http://www.mail-archive.com/r-sig-phylo@r-project.org/ 
> > <http://www.mail-archive.com/r-sig-phylo@r-project.org/>
> 
> _______________________________________________
> R-sig-phylo mailing list - R-sig-phylo@r-project.org 
> <mailto:R-sig-phylo@r-project.org>
> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo 
> <https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>
> Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/ 
> <http://www.mail-archive.com/r-sig-phylo@r-project.org/>
> 
> 
> <swarm.pdf>_______________________________________________
> R-sig-phylo mailing list - R-sig-phylo@r-project.org 
> <mailto:R-sig-phylo@r-project.org>
> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo 
> <https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>
> Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/ 
> <http://www.mail-archive.com/r-sig-phylo@r-project.org/>
____________________________________________
Marguerite A. Butler
Associate Professor

Department of Biology 
2538 McCarthy Mall, Edmondson Hall 216 
Honolulu, HI 96822

Office: 808-956-4713
Dept: 808-956-8617
Lab:  808-956-5867
FAX:   808-956-4745
http://butlerlab.org
http://manoa.hawaii.edu/biology/people/marguerite-butler
http://www2.hawaii.edu/~mbutler















        [[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/

Reply via email to