SE of estimated heritability

[kbroman]

There was a question at the R/qtl2 google group about getting an SE on the estimated heritability from est_herit().

Could use equation 2 from Visscher and Goddard (2015), which just depends on eigenvalues of kinship matrix, sample size, and heritability.

[kbroman]

Here's a stab at it:

# estimated SE of estimated heritability
#
# Visscher & Goddard (2015), https://doi.org/10.1534/genetics.114.171017
#
# lambda_i = eigenvalues of K
#
# a = sum[ (lam_i-1)^2 / (1 + h^2(lam_i - 1))^2 ]
#
# b = sum[ (lam_i - 1) / (1 + h^2(lam_i - 1)) ]
#
# var(h_hat^2) = 2 / (a - b^2/N)
#
# N = sample size

se_herit <-
    function(pheno, kinship)
{

    h2 <- as.numeric(est_herit(pheno, kinship))

    d <- decomp_kinship(kinship)

    lam <- d$values
    a <- sapply(h2, function(hsq) sum( (lam-1)^2 / (1 + hsq*(lam-1))^2 ) )
    b <- sapply(h2, function(hsq) sum( (lam-1) / (1 + hsq*(lam-1)) ) )
    n <- length(lam)

    setNames( sqrt(2 / (a - b^2/n)), colnames(pheno) )
}