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/
