Standard errors of complex functions of variance components from REML Information matrix

In airemlf90 and other programs, a byproduct of the AIREML is the inverse of the Information Matrix that, asymptotically, describes the sampling covariance of the variance components. It can also be considered as an approximation to the posterior distribution of the parameters.

Getting standard errors of functions of the variance components is rather painful algebraically as one requires the delta method and Taylor series expansions. A much simpler method , implemented in as OPTION se_covar_function in airemlf90, comes through MonteCarlo sampling and was described by Meyer and Houle. Here I will show how this work.
This a model with three variance components, genetic, permanent and residual variance. Estimates from airemlf90 are in vector ahat:

#genetic, permanent, residual

ahat=c(

  0.11478,    

  0.13552,    

  0.25290,

  )

with AI matrix:

# inverse of AI matrix (Sampling Variance)

AI=matrix(c(

  0.16799E-05, -0.96486E-06, -0.82566E-08,

 -0.96486E-06,  0.96167E-06, -0.37113E-07,

 -0.82566E-08, -0.37113E-07,  0.10864E-06)

 ,ncol=3)

asymptotically, the "true" value has a sampling (or posterior) distribution a~N(ahat,AI), and getting samples from this distribution is easy:

require(MASS)

b=mvrnorm(10000,ahat,AI)

> head(b)

          [,1]      [,2]      [,3]

[1,] 0.1146738 0.1357640 0.2529399

[2,] 0.1163889 0.1342926 0.2528479

[3,] 0.1166155 0.1344342 0.2525161

[4,] 0.1142085 0.1358928 0.2534974

[5,] 0.1136835 0.1361108 0.2530133

[6,] 0.1140485 0.1365707 0.2530573

and from here, we can obtain samples of the heritability and its standard deviation:

h2=b[,1]/(b[,1]+b[,2]+b[,3])

sd(h2)

> 0.002318198

This is very useful, i.e. for complex model such as consideration of genomic and dominance variation.