Ok, I just posted a new version of evol.vcv() - direct link to code here. The main new feature of this version (v0.2) is that it also computes the variances of the parameter estimates. (The standard errors can then just be computed by taking the square roots.) These variances come from the matrix of partial second derivatives of the likelihood surface at the optimum (i.e., the curvature of the surface), also known as the Hessian matrix. The intuitive relationship between the curvature of the likelihood surface and the uncertainty of our parameter estimates is quite straightforward: if the likelihood surface is very flat (low curvature), then the uncertainty in our parameter estimates is large; conversely if the likelihood surface is very peaked at the optimum (i.e., high curvature), then this implies that uncertainty in our parameter estimates is low.
This function (that is, evol.vcv()) implements the method of Revell & Collar (2009) and was also described in a prior blog post.
In R, normally obtaining the Hessian matrix is easy because the built-in optimizer, optim(), optionally returns its value. Unfortunately, as I outlined in a previous post, this would not do for evol.vcv() because in this function I actually optimized the likelihood with respect to the Cholesky matrices rather than the evolutionary variance-covariances themselves. This was done for computational reasons because the evolutionary variance-covariance matrices (called "rate matrices" in our paper) are subject to annoying constraints which mean that for most random values the likelihood cannot be computed.
The solution, which I figured out yesterday and then implemented today, was to optimize with respect to the Cholesky matrices (as before) - but then differentiate the likelihood surface for the full VCV matrices at the optimum using hessian() in the "numDeriv" package (CRAN page here). This seems to work - and this analysis is implemented in v0.2 of evol.vcv() on my R phylogenetics page. This will also be included in the next version of "phytools".
The variances are not computed by default. To get them, load the present version from source:
> source("http://faculty.umb.edu/liam.revell/phytools/evol.vcv/v0.2/evol.vcv.R") # this should work
for SIMMAP style tree and data matrix X. Good luck, and please report if this function seems to work properly or not.