Finding the standard errors from the likelihood surface in evol.vcv()

Today I've been working on returning the variances of the parameter estimates in evol.vcv(). Just to remind the readers of this blog, this function fits two or more instantaneous variance-covariance matrices for the Brownian evolutionary process to different pre-specified parts of a phylogenetic tree. (For more details, see my 2009 paper with Dave Collar or a previous blog post here.)

The first problem that I've run into is that my function performs optimization of the likelihood as a function of the Cholesky decomposite matrices of R. This is because the evolutionary variance-covariance matrix are constrained to be positive semi-definite - a constraint which will cause problems for optim(). My solution to this problem was to optimize the Cholesky decomposite matrix, which does not have this constraint. My issue, then, is that if I return the Hessian matrix of partial second derivatives these will be the curvature of the likelihood surface for the Cholesky matrix elements - rather than for the variances and covariances of our evolutionary rate matrices.

My current plan to solve this involves first using optim() to find the MLEs based on the Cholesky decomposed matrices; and then feeding these parameter estimates along with a likelihood function for the full rate matrices into hessian() in the package "numDeriv". If this works I will soon report the result.