We have done a grid search for a genetic correlation, from 0 to 1 in steps of 0.1. The other variance components were fixed and only the covariance was varied. This was done making calls to remlf90 with the variance components and OPTION maxrounds 0 . In this manner -2logL is obtained. Technically this gives a so-called profile likelihood. The resulting profile is as below:
So the REML estimate is 0.7, but what can we say about the accuracy of this estimate? This, we solve today here :-)
REML estimation software (VCE, Wombat, AIREMLF90) usually yields s.e. of estimates of variance components and their function such as heritabilities. Good accounts on the theory can be read in
(statistical basis)
Casella, G., & Berger, R. L. (2021). Statistical inference. Chapter 10
Sorensen, D., Gianola, D., & Gianola, D. (2002). Likelihood, Bayesian and MCMC methods in quantitative genetics. (3.26)
(application in REML software)
Meyer, K., & Smith, S. P. (1996). Restricted maximum likelihood estimation for animal models using derivatives of the likelihood. Genetics Selection Evolution, 28(1), 23-49.
Meyer, K., & Houle, D. (2013, October). Sampling based approximation of confidence intervals for functions of genetic covariance matrices. In Proc. Assoc. Advmt. Anim. Breed. Genet (Vol. 20, pp. 523-526).
However, our problem is quite simple. The profile likelihood looks a bit like a parable. The information matrix is minus the expected 2nd derivative of the log-likelihood. To obtain this, we want a function $latex logL=f(r_g)$latex. So, we do an approximation to this function by fitting a quadratic regression $latex logL = a + b r_g +c r_g^2 $latex to the points above. In R (theta= rg throughout):
cbind(theta,minus2logL)
# theta minus2logL
# [1,] 0.10 7172238
# [2,] 0.20 7172210
# [3,] 0.30 7172186
# [4,] 0.40 7172167
# [5,] 0.50 7172152
# [6,] 0.60 7172141
# [7,] 0.70 7172137
# [8,] 0.80 7172140
# [9,] 0.90 7172158
#[10,] 0.99 7172223
res=lm(logL ~ theta + I(theta^2))
# Coefficients:
# (Intercept) theta I(theta^2)
# -1.249e+07 8.040e+02 -6.481e+02
Then from here, I need the 2nd derivative, but this is very simple: $latex \frac{dlogL}{d\theta^2}=2c $latex
Again, in R
# the 2nd derivative of this quadratic is 2c
# for the information matrix (actually a single coefficient)
# I need to pre-multiply by -1 (by definition)
info=-2*res$coefficients[3]
So, info is the information matrix. Now, the asymptotic variance is the inverse of the Information matrix (for instance, Sorensen and Gianola, 3.66) so that (informally as they say):
$latex \hat{\theta}_n \sim N(\theta_0,I^{-1}(\theta_0)) $latex
where I is the information matrix ("info" above). So we invert and then just take the square root to get the s.d.:
# now we invert this matrix
infoi=solve(info)
# this is the asymptotic variance of the parameter.
# so the s.d. is the sq root:
sdasymp=sqrt(infoi)
cat("s.d. of the parameter= ",sdasymp,"\n")
#s.d. of the parameter= 0.04898895
