# Re: [R-sig-phylo] Interpretation of standard errors of parameter estimates in OUwie models

```Apparently the plot did not come through correctly. I'm attaching a PDF of
it.```
```
Best,
Brian

_______________________________________________________________________
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:24 PM, Brian O'Meara <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)
>
>
> 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
>> >>
>> >>
>> >>
>> >> _______________________________________________
>> >> 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/
>>
>
>
```

swarm.pdf
```_______________________________________________