Heritability estimation

We provide heritability estimation for Normal, Bernoulli, Probit, Binomial, and Poisson phenotypes. A standard LMM is used for Normal traits:

\[𝐲 = 𝙼𝛂 + 𝐯 + 𝛆,\]

where

\[𝐯 ∼ 𝓝(𝟎, 𝓋₀𝙺) ~~\text{and}~~ 𝛆 ∼ 𝓝(𝟎, 𝓋₁𝙸).\]

A GLMM is used to model the other type of traits:

\[𝐳 = 𝙼𝛂 + 𝐯 + 𝛆, ~~\text{where}~~ yᵢ|𝐳 ∼ 𝙴𝚡𝚙𝙵𝚊𝚖(𝜇ᵢ=g(zᵢ))\]

and 𝐯 and 𝛆 are defined as before.

In both cases, the parameters are the same: 𝛂, 𝓋₀, and 𝓋₁. They are fitted via restricted maximum likelihood for LMM and via maximum likelihood for GLMM. The covariance-matrix 𝙺 given by the user is normalised before the model is fitted as follows:

K = K / K.diagonal().mean()

limix.her.estimate(y, lik, K, M=None, verbose=True)

Estimate the so-called narrow-sense heritability.

It supports Normal, Bernoulli, Probit, Binomial, and Poisson phenotypes.

Parameters

• y (array_like) – Array of trait values of n individuals.

• lik (tuple, "normal", "bernoulli", "probit", "binomial", "poisson") – Sample likelihood describing the residual distribution. Either a tuple or a string specifying the likelihood is required. The Normal, Bernoulli, Probit, and Poisson likelihoods can be selected by providing a string. Binomial likelihood on the other hand requires a tuple because of the number of trials: ("binomial", array_like). Defaults to "normal".

• K (n×n array_like) – Sample covariance, often the so-called kinship matrix. It might be, for example, the estimated kinship relationship between the individuals. The provided matrix will be normalised as K = K / K.diagonal().mean().

• M (n×c array_like, optional) – Covariates matrix. If an array is passed, it will used as is; no normalisation will be performed. If None is passed, an offset will be used as the only covariate. Defaults to None.

• verbose (bool, optional) – True to display progress and summary; False otherwise.

Returns

Estimated heritability.

Return type

float

Examples

>>> from numpy import dot, exp, sqrt
>>> from numpy.random import RandomState
>>> from limix.her import estimate
>>>
>>> random = RandomState(0)
>>>
>>> G = random.randn(150, 200) / sqrt(200)
>>> K = dot(G, G.T)
>>> z = dot(G, random.randn(200)) + random.randn(150)
>>> y = random.poisson(exp(z))
>>>
>>> print(estimate(y, 'poisson', K, verbose=False))  
0.18311439918863426

Notes

It will raise a ValueError exception if non-finite values are passed. Please, refer to the limix.qc.mean_impute() function for missing value imputation.