Model¶

We also export the underlying inference methods limix uses. Its classes can be accessed via limix module limix.model or via the package glimix_core itself. Both ways should work identically.

Linear Mixed Models¶

class limix.model.lmm.LMM(y, X, QS=None, restricted=False)[source]¶

Fast Linear Mixed Models inference via maximum likelihood.

Examples

>>> from numpy import array
>>> from numpy_sugar.linalg import economic_qs_linear
>>> from glimix_core.lmm import LMM
>>>
>>> X = array([[1, 2], [3, -1]], float)
>>> QS = economic_qs_linear(X)
>>> covariates = array([[1], [1]])
>>> y = array([-1, 2], float)
>>> lmm = LMM(y, covariates, QS)
>>> lmm.fit(verbose=False)
>>> print('%.3f' % lmm.lml())
-3.649

One can also specify which parameters should be fitted:

>>> from numpy import array
>>> from numpy_sugar.linalg import economic_qs_linear
>>> from glimix_core.lmm import LMM
>>>
>>> X = array([[1, 2], [3, -1]], float)
>>> QS = economic_qs_linear(X)
>>> covariates = array([[1], [1]])
>>> y = array([-1, 2], float)
>>> lmm = LMM(y, covariates, QS)
>>> lmm.fix('delta')
>>> lmm.fix('scale')
>>> lmm.delta = 0.5
>>> lmm.scale = 1
>>> lmm.fit(verbose=False)
>>> print('%.3f' % lmm.lml())
-3.832
>>> lmm.unfix('delta')
>>> lmm.fit(verbose=False)
>>> print('%.3f' % lmm.lml())
-3.713

Notes

The LMM model can be equivalently written as

𝐲 ∼ 𝓝(𝚇𝜷, 𝑠((1-𝛿)𝙺 + 𝛿𝙸)),

and we thus have v₀ = s (1 - 𝛿) and v₁ = s 𝛿. Consider the economic eigendecomposition of 𝙺:

𝙺 = [𝚀₀  𝚀₁] [𝚂₀  𝟎] [𝚀₀ᵀ]
             [ 𝟎  𝟎] [𝚀₁ᵀ]

and let

𝙳 = [(1-𝛿)𝚂₀ + 𝛿𝙸₀ 𝟎 ]
[ 𝟎 𝛿𝙸₁].

In order to eliminate the need of 𝚀₁, note that 𝚀𝚀ᵀ = 𝙸 implies that

𝚀₁𝚀₁ᵀ = 𝙸 - 𝚀₀𝚀₀ᵀ.

We will need to solve ((1-𝛿)𝙺 + 𝛿𝙸)𝐱 = 𝐮 for 𝐱. Let 𝙳₀ = ((1-𝛿)𝚂₀ + 𝛿𝙸₀) and let us define

𝙱 = 𝚀₀𝙳₀⁻¹𝚀₀ᵀ                   if 𝛿=0, and
𝙱 = 𝚀₀𝙳₀⁻¹𝚀₀ᵀ + 𝛿⁻¹(𝙸 - 𝚀₀𝚀₀ᵀ)  if 𝛿>0.

We therefore have

𝐱 = 𝙱𝐮.

property X¶

Covariates matrix.

Returns: X – Covariates.
Return type: ndarray

property beta: numpy.ndarray¶

Fixed-effect sizes.

Returns: Optimal fixed-effect sizes.
Return type: effect-sizes

Notes

Setting the derivative of log(p(𝐲)) over effect sizes equal to zero leads to solutions 𝜷 from equation

𝚇ᵀ𝙱𝚇𝜷 = 𝚇ᵀ𝙱𝐲

property beta_covariance¶

Estimates the covariance-matrix of the optimal beta.

Returns: beta-covariance – 𝑠(𝚇ᵀ𝙱𝚇)⁻¹.
Return type: ndarray

References

covariance()[source]¶

Covariance of the prior.

Returns: covariance – v₀𝙺 + v₁𝙸.
Return type: ndarray

property delta: float¶: Variance ratio between K and I.

fit(verbose=True)[source]¶

Maximise the marginal likelihood.

Parameters: verbose (bool, optional) – True for progress output; False otherwise. Defaults to True.

fix(param: str)[source]¶

Disable parameter optimization.

Parameters: param – Possible values are "delta", "beta", and "scale".

get_fast_scanner() → glimix_core.lmm._lmm_scan.FastScanner[source]¶

Return FastScanner for association scan.

Returns: Instance of a class designed to perform very fast association scan.
Return type: fast-scanner

gradient()[source]¶: Not implemented.

lml() → float [source]¶

Log of the marginal likelihood.

Returns: Log of the marginal likelihood.
Return type: lml

Notes

The log of the marginal likelihood is given by

2⋅log(p(𝐲)) = -n⋅log(2π) - n⋅log(s) - log|𝙳|
    -  (𝚀₀ᵀ𝐲)ᵀ(s𝙳₀)⁻¹(𝚀₀ᵀ𝐲)  -  (𝐲)ᵀ(s𝛿)⁻¹(𝐲)  +  (𝚀₀ᵀ𝐲)ᵀ(s𝛿)⁻¹(𝚀₀ᵀ𝐲)
    - (𝚀₀ᵀ𝚇𝜷)ᵀ(s𝙳₀)⁻¹(𝚀₀ᵀ𝚇𝜷) - (𝚇𝜷)ᵀ(s𝛿)⁻¹(𝚇𝜷) + (𝚀₀ᵀ𝚇𝜷)ᵀ(s𝛿)⁻¹(𝚀₀ᵀ𝚇𝜷)
    + 2(𝚀₀ᵀ𝐲)ᵀ(s𝙳₀)⁻¹(𝚇𝜷)    + 2(𝐲)ᵀ(s𝛿)⁻¹(𝚇𝜷) - 2(𝚀₀ᵀ𝐲)ᵀ(s𝛿)⁻¹(𝚀₀ᵀ𝚇𝜷)

By using the optimal 𝜷, the log of the marginal likelihood can be rewritten as:

2⋅log(p(𝐲)) = -n⋅log(2π) - n⋅log(s) - log|𝙳| + (𝚀₀ᵀ𝐲)ᵀ(s𝙳₀)⁻¹𝚀₀ᵀ(𝚇𝜷 - 𝐲)
            + (𝐲)ᵀ(s𝛿)⁻¹(𝚇𝜷 - 𝐲) - (𝚀₀ᵀ𝐲)ᵀ(s𝛿)⁻¹𝚀₀ᵀ(𝚇𝜷 - 𝐲).

In the extreme case where 𝜷 is such that 𝐲 = 𝚇𝜷, the maximum is attained as s→0.

For optimals 𝜷 and s, the log of the marginal likelihood can be further simplified to

2⋅log(p(𝐲; 𝜷, s)) = -n⋅log(2π) - n⋅log s - log|𝙳| - n.

mean()[source]¶

Mean of the prior.

Formally, 𝐦 = 𝚇𝜷.

Returns: mean – Mean of the prior.
Return type: ndarray

property ncovariates: int¶: Number of covariates.

property nsamples: int¶: Number of samples.

property scale: float¶

Scaling factor.

Returns: scale – Scaling factor.
Return type: float

Notes

Setting the derivative of log(p(𝐲; 𝜷)), for which 𝜷 is optimal, over scale equal to zero leads to the maximum

s = n⁻¹(Qᵀ𝐲)ᵀD⁻¹ Qᵀ(𝐲-𝚇𝜷).

In the case of restricted marginal likelihood

s = (n-c)⁻¹(Qᵀ𝐲)ᵀD⁻¹ Qᵀ(𝐲-𝚇𝜷),

where s is the number of covariates.

unfix(param: str)[source]¶

Enable parameter optimization.

Parameters: param – Possible values are "delta", "beta", and "scale".

property v0¶

First variance.

Returns: v0 – s(1 - 𝛿).
Return type: float

property v1: float¶

Second variance.

Returns: s𝛿.
Return type: v1

value() → float [source]¶: Internal use only.

class limix.model.lmm.Kron2Sum(Y, A, X, G, rank=1, restricted=False)[source]¶

LMM for multi-traits fitted via maximum likelihood.

This implementation follows the work published in [CA05]. Let n, c, and p be the number of samples, covariates, and traits, respectively. The outcome variable Y is a n×p matrix distributed according to:

vec(Y) ~ N((A ⊗ X) vec(B), K = C₀ ⊗ GGᵀ + C₁ ⊗ I).

A and X are design matrices of dimensions p×p and n×c provided by the user, where X is the usual matrix of covariates commonly used in single-trait models. B is a c×p matrix of fixed-effect sizes per trait. G is a n×r matrix provided by the user and I is a n×n identity matrices. C₀ and C₁ are both symmetric matrices of dimensions p×p, for which C₁ is guaranteed by our implementation to be of full rank. The parameters of this model are the matrices B, C₀, and C₁.

For implementation purpose, we make use of the following definitions:

𝛃 = vec(B)
M = A ⊗ X
H = MᵀK⁻¹M
Yₓ = LₓY
Yₕ = YₓLₕᵀ
Mₓ = LₓX
Mₕ = (LₕA) ⊗ Mₓ
mₕ = Mₕvec(B)

where Lₓ and Lₕ are defined in glimix_core.cov.Kron2SumCov.

References

CA05: Casale, F. P., Rakitsch, B., Lippert, C., & Stegle, O. (2015). Efficient set tests for the genetic analysis of correlated traits. Nature methods, 12(8), 755.

property A¶

A from the equation 𝐦 = (A ⊗ X) vec(B).

Returns

A –

Return type

ndarray

property B¶

Fixed-effect sizes B from 𝐦 = (A ⊗ X) vec(B).

Returns: fixed-effects – B from 𝐦 = (A ⊗ X) vec(B).
Return type: ndarray

property C0¶

C₀ from equation K = C₀ ⊗ GGᵀ + C₁ ⊗ I.

Returns: C0 – C₀.
Return type: ndarray

property C1¶

C₁ from equation K = C₀ ⊗ GGᵀ + C₁ ⊗ I.

Returns: C1 – C₁.
Return type: ndarray

property M¶

M = (A ⊗ X).

Returns: M – M from M = (A ⊗ X).
Return type: ndarray

property X¶

X from equation M = (A ⊗ X).

Returns: X – X from M = (A ⊗ X).
Return type: ndarray

property beta¶

Fixed-effect sizes 𝛃 = vec(B).

Returns: fixed-effects – 𝛃 from 𝛃 = vec(B).
Return type: ndarray

property beta_covariance¶

Estimates the covariance-matrix of the optimal beta.

Returns: beta-covariance – (MᵀK⁻¹M)⁻¹.
Return type: ndarray

References

covariance()[source]¶

Covariance K = C₀ ⊗ GGᵀ + C₁ ⊗ I.

Returns

covariance –

Return type

ndarray

fit(verbose=True)[source]¶

Maximise the marginal likelihood.

Parameters: verbose (bool, optional) – True for progress output; False otherwise. Defaults to True.

get_fast_scanner()[source]¶

Return FastScanner for association scan.

Returns: Instance of a class designed to perform very fast association scan.
Return type: FastScanner

gradient()[source]¶: Gradient of the log of the marginal likelihood.

lml()[source]¶

Log of the marginal likelihood.

Let 𝐲 = vec(Y), M = A⊗X, and H = MᵀK⁻¹M. The restricted log of the marginal likelihood is given by [R07]:

2⋅log(p(𝐲)) = -(n⋅p - c⋅p) log(2π) + log(｜MᵀM｜) - log(｜K｜) - log(｜H｜)
    - (𝐲-𝐦)ᵀ K⁻¹ (𝐲-𝐦),

where 𝐦 = M𝛃 for 𝛃 = H⁻¹MᵀK⁻¹𝐲.

For implementation purpose, let X = (L₀ ⊗ G) and R = (L₁ ⊗ I)(L₁ ⊗ I)ᵀ. The covariance can be written as:

K = XXᵀ + R.

From the Woodbury matrix identity, we have

𝐲ᵀK⁻¹𝐲 = 𝐲ᵀR⁻¹𝐲 - 𝐲ᵀR⁻¹XZ⁻¹XᵀR⁻¹𝐲,

where Z = I + XᵀR⁻¹X. Note that R⁻¹ = (U₁S₁⁻¹U₁ᵀ) ⊗ I and

XᵀR⁻¹𝐲 = (L₀ᵀW ⊗ Gᵀ)𝐲 = vec(GᵀYWL₀),

where W = U₁S₁⁻¹U₁ᵀ. The term GᵀY can be calculated only once and it will form a r×p matrix. We similarly have

XᵀR⁻¹M = (L₀ᵀWA) ⊗ (GᵀX),

for which GᵀX is pre-computed.

The log-determinant of the covariance matrix is given by

log(｜K｜) = log(｜Z｜) - log(｜R⁻¹｜) = log(｜Z｜) - 2·n·log(｜U₁S₁⁻½｜).

The log of the marginal likelihood can be rewritten as:

2⋅log(p(𝐲)) = -(n⋅p - c⋅p) log(2π) + log(｜MᵀM｜)
- log(｜Z｜) + 2·n·log(｜U₁S₁⁻½｜)
- log(｜MᵀR⁻¹M - MᵀR⁻¹XZ⁻¹XᵀR⁻¹M｜)
- 𝐲ᵀR⁻¹𝐲 + (𝐲ᵀR⁻¹X)Z⁻¹(XᵀR⁻¹𝐲)
- 𝐦ᵀR⁻¹𝐦 + (𝐦ᵀR⁻¹X)Z⁻¹(XᵀR⁻¹𝐦)
+ 2𝐲ᵀR⁻¹𝐦 - 2(𝐲ᵀR⁻¹X)Z⁻¹(XᵀR⁻¹𝐦).

Returns: lml – Log of the marginal likelihood.
Return type: float

References

R07: LaMotte, L. R. (2007). A direct derivation of the REML likelihood function. Statistical Papers, 48(2), 321-327.

mean()[source]¶

Mean 𝐦 = (A ⊗ X) vec(B).

Returns: mean – 𝐦.
Return type: ndarray

property ncovariates¶: Number of covariates, c.

property nsamples¶: Number of samples, n.

property ntraits¶: Number of traits, p.

value()[source]¶: Log of the marginal likelihood.

Generalised Linear Mixed Models¶

class limix.model.glmm.GLMMNormal(eta, tau, X, QS=None)[source]¶

LMM with heterogeneous Normal likelihoods.

We model

\[\tilde{\boldsymbol\mu} \sim \mathcal N(\mathrm X\boldsymbol\beta, v_0 \mathrm K + v_1 \mathrm I + \tilde{\Sigma}),\]

where \(\tilde{\boldsymbol\mu}\) and \(\tilde{\Sigma}\) are typically given by EP approximation to non-Normal likelihood (please, refer to glimix_core.expfam.GLMMExpFam). If that is the case, \(\tilde{\Sigma}\) is a diagonal matrix with non-negative values. Those EP parameters are given to this class via their natural forms: eta and tau.

Parameters

eta (array_like) – \([\tilde{\mu}_0\tilde{\sigma}^{-2}_0 \quad \tilde{\mu}_1\tilde{\sigma}^{-2}_1 \quad\cdots\quad \tilde{\mu}_n\tilde{\sigma}^{-2}_n]\).
tau (array_like) – \([\tilde{\sigma}^{-2}_0\quad\tilde{\sigma}^{-2}_1 \quad\cdots\quad \tilde{\sigma}^{-2}_n]\).
X (array_like) – Covariates.
QS (tuple) – Economic eigendecomposition of \(\mathrm K\).

property beta¶

Fixed-effect sizes.

Returns: \(\boldsymbol\beta\).
Return type: numpy.ndarray

fix(var_name)[source]¶

Prevent a variable to be adjusted.

Parameters: var_name (str) – Variable name.

get_fast_scanner()[source]¶: Return glimix_core.lmm.FastScanner for the current delta.

unfix(var_name)[source]¶

Let a variable be adjusted.

Parameters: var_name (str) – Variable name.

value()[source]¶

Log of the marginal likelihood.

Formally,

\[- \frac{n}{2}\log{2\pi} - \frac{1}{2} \log{\left| v_0 \mathrm K + v_1 \mathrm I + \tilde{\Sigma} \right|} - \frac{1}{2} \left(\tilde{\boldsymbol\mu} - \mathrm X\boldsymbol\beta\right)^{\intercal} \left( v_0 \mathrm K + v_1 \mathrm I + \tilde{\Sigma} \right)^{-1} \left(\tilde{\boldsymbol\mu} - \mathrm X\boldsymbol\beta\right)\]

Returns: \(\log{p(\tilde{\boldsymbol\mu})}\)
Return type: float

class limix.model.glmm.GLMMExpFam(y, lik, X, QS=None, n_int=1000, rtol=1.49e-05, atol=1.49e-08)[source]¶

Generalised Linear Gaussian Processes implementation.

It implements inference over GLMMs via the Expectation Propagation [Min01] algorithm. It currently supports the "Bernoulli", "Probit", "Binomial", and "Poisson" likelihoods. (For heterogeneous Normal likelihood, please refer to glimix_core.glmm.GLMMNormal for a closed-form inference.)

Parameters

y (array_like) – Outcome variable.
lik (tuple) – Likelihood definition. The first item is one of the following likelihood names: "Bernoulli", "Binomial", "Normal", and "Poisson". For Binomial, the second item is an array of outcomes.
X (array_like) – Covariates.
QS (tuple) – Economic eigendecomposition.

Example

>>> from numpy import dot, sqrt, zeros
>>> from numpy.random import RandomState
>>>
>>> from numpy_sugar.linalg import economic_qs
>>>
>>> from glimix_core.glmm import GLMMExpFam
>>>
>>> random = RandomState(0)
>>> nsamples = 10
>>>
>>> X = random.randn(nsamples, 2)
>>> G = random.randn(nsamples, 100)
>>> K = dot(G, G.T)
>>> ntrials = random.randint(1, 100, nsamples)
>>> z = dot(G, random.randn(100)) / sqrt(100)
>>>
>>> successes = zeros(len(ntrials), int)
>>> for i in range(len(ntrials)):
...    successes[i] = sum(z[i] + 0.2 * random.randn(ntrials[i]) > 0)
>>>
>>> QS = economic_qs(K)
>>>
>>> glmm = GLMMExpFam(successes, ('binomial', ntrials), X, QS)
>>> print('Before: %.2f' % glmm.lml())
Before: -16.40
>>> glmm.fit(verbose=False)
>>> print('After: %.2f' % glmm.lml())
After: -13.43

property beta¶

Fixed-effect sizes.

Returns: \(\boldsymbol\beta\).
Return type: numpy.ndarray

covariance()[source]¶

Covariance of the prior.

Returns: \(v_0 \mathrm K + v_1 \mathrm I\).
Return type: numpy.ndarray

fit(verbose=True, factr=100000.0, pgtol=1e-07)[source]¶

Maximise the marginal likelihood.

Parameters

verbose (bool) – True for progress output; False otherwise. Defaults to True.
factr (float, optional) – The iteration stops when (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr * eps, where eps is the machine precision.
pgtol (float, optional) – The iteration will stop when max{|proj g_i | i = 1, ..., n} <= pgtol where pg_i is the i-th component of the projected gradient.

Notes

Please, refer to scipy.optimize.fmin_l_bfgs_b() for further information about factr and pgtol.

fix(var_name)[source]¶

Prevent a variable to be adjusted.

Parameters: var_name (str) – Variable name.

gradient()[source]¶

Gradient of the log of the marginal likelihood.

Returns: Map between variables to their gradient values.
Return type: dict

posteriori_covariance()[source]¶: Covariance of the estimated posteriori.

posteriori_mean()[source]¶

Mean of the estimated posteriori.

This is also the maximum a posteriori estimation of the latent variable.

unfix(var_name)[source]¶

Let a variable be adjusted.

Parameters: var_name (str) – Variable name.

Gaussian Process¶

class limix.model.gp.GP(y, mean, cov)[source]¶

Gaussian Process inference via maximum likelihood.

Parameters

y (array_like) – Outcome variable.
mean (function) – Mean function. (Refer to Mean functions.)
cov (function) – Covariance function. (Refer to Covariance functions.)

Example

>>> from numpy.random import RandomState
>>>
>>> from glimix_core.example import offset_mean
>>> from glimix_core.example import linear_eye_cov
>>> from glimix_core.gp import GP
>>> from glimix_core.random import GPSampler
>>>
>>> random = RandomState(94584)
>>>
>>> mean = offset_mean()
>>> cov = linear_eye_cov()
>>>
>>> y = GPSampler(mean, cov).sample(random)
>>>
>>> gp = GP(y, mean, cov)
>>> print('Before: %.4f' % gp.lml())
Before: -15.5582
>>> gp.fit(verbose=False)
>>> print('After: %.4f' % gp.lml())
After: -13.4791
>>> print(gp)  
GP(...)
  lml: -13.47907874997517
  OffsetMean(): OffsetMean
    offset: 0.7755803668772308
  SumCov(covariances=...): SumCov
    LinearCov(): LinearCov
      scale: 2.061153622438558e-09
    EyeCov(dim=10): EyeCov
      scale: 0.8675680523425126

fit(verbose=True, factr=100000.0, pgtol=1e-07)[source]¶

Maximise the marginal likelihood.

Parameters

verbose (bool) – True for progress output; False otherwise. Defaults to True.
factr (float, optional) – The iteration stops when (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr * eps, where eps is the machine precision.
pgtol (float, optional) – The iteration will stop when max{|proj g_i | i = 1, ..., n} <= pgtol where pg_i is the i-th component of the projected gradient.

Notes

Please, refer to scipy.optimize.fmin_l_bfgs_b() for further information about factr and pgtol.

lml()[source]¶

Log of the marginal likelihood.

Returns: \(\log p(\mathbf y)\)
Return type: float

Generalised Gaussian Process¶

class limix.model.ggp.ExpFamGP(y, lik, mean, cov)[source]¶

Expectation Propagation for Generalised Gaussian Processes.

Parameters

y (array_like) – Outcome variable.
lik_name (str) – Likelihood name.
mean (function) – Mean function. (Refer to Mean functions.)
cov (function) – Covariance function. (Refer to Covariance functions.)

Example

>>> from numpy.random import RandomState
>>>
>>> from glimix_core.example import offset_mean
>>> from glimix_core.example import linear_eye_cov
>>> from glimix_core.ggp import ExpFamGP
>>> from glimix_core.lik import BernoulliProdLik
>>> from glimix_core.link import LogitLink
>>> from glimix_core.random import GGPSampler
>>>
>>> random = RandomState(1)
>>>
>>> lik = BernoulliProdLik(LogitLink())
>>> mean = offset_mean()
>>> cov = linear_eye_cov()
>>>
>>> y = GGPSampler(lik, mean, cov).sample(random)
>>>
>>> ggp = ExpFamGP(y, 'bernoulli', mean, cov)
>>> print('Before: %.4f' % ggp.lml())
Before: -6.9802
>>> ggp.fit(verbose=False)
>>> print('After: %.2f' % ggp.lml())
After: -6.55

fit(verbose=True, factr=100000.0, pgtol=1e-07)[source]¶

Maximise the marginal likelihood.

Parameters

verbose (bool) – True for progress output; False otherwise. Defaults to True.
factr (float, optional) – The iteration stops when (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr * eps, where eps is the machine precision.
pgtol (float, optional) – The iteration will stop when max{|proj g_i | i = 1, ..., n} <= pgtol where pg_i is the i-th component of the projected gradient.

Notes

Please, refer to scipy.optimize.fmin_l_bfgs_b() for further information about factr and pgtol.

lml()[source]¶

Log of the marginal likelihood.

Returns: \(\log p(\mathbf y)\)
Return type: float

Covariance¶

class limix.model.cov.EyeCov(dim)[source]¶

Identity covariance function, K = s·I.

The parameter s is the scale of the matrix.

Example

>>> from glimix_core.cov import EyeCov
>>>
>>> cov = EyeCov(2)
>>> cov.scale = 2.5
>>> print(cov.value())
[[2.5 0. ]
 [0.  2.5]]
>>> g = cov.gradient()
>>> print(g['logscale'])
[[2.5 0. ]
 [0.  2.5]]
>>> cov.name = "I"
>>> print(cov)
EyeCov(dim=2): I
  scale: 2.5

Parameters: dim (int) – Matrix dimension, d.

property dim¶

Dimension of the matrix, d.

It corresponds to the number of rows and to the number of columns.

gradient()[source]¶

Derivative of the covariance matrix over log(s), s⋅I.

Returns: logscale – s⋅I, for scale s and a d×d identity matrix I.
Return type: ndarray

property scale¶: Scale parameter.

value()[source]¶

Covariance matrix.

Returns: K – s⋅I, for scale s and a d×d identity matrix I.
Return type: ndarray

class limix.model.cov.FreeFormCov(dim)[source]¶

General definite positive matrix, K = LLᵀ + ϵI.

A d×d covariance matrix K will have ((d+1)⋅d)/2 parameters defining the lower triangular elements of a Cholesky matrix L such that:

K = LLᵀ + ϵI,

for a very small positive number ϵ. That additional term is necessary to avoid singular and ill conditioned covariance matrices.

Example

>>> from glimix_core.cov import FreeFormCov
>>>
>>> cov = FreeFormCov(2)
>>> cov.L = [[1., .0], [0.5, 2.]]
>>> print(cov.gradient()["Lu"])
[[[0.  2.  0. ]
  [1.  0.5 0. ]]

 [[1.  0.5 0. ]
  [1.  0.  8. ]]]
>>> cov.name = "K"
>>> print(cov)
FreeFormCov(dim=2): K
  L: [[1.  0. ]
      [0.5 2. ]]

property L¶

Lower-triangular matrix L such that K = LLᵀ + ϵI.

Returns: L – Lower-triangular matrix.
Return type: (d, d) ndarray

property Lu¶: Lower-triangular, flat part of L.

eigh()[source]¶

Eigen decomposition of K.

Returns

S (ndarray) – The eigenvalues in ascending order, each repeated according to its multiplicity.
U (ndarray) – Normalized eigenvectors.

fix()[source]¶: Disable parameter optimisation.

gradient()[source]¶

Derivative of the covariance matrix over the parameters of L.

Returns: Lu – Derivative of K over the lower triangular part of L.
Return type: ndarray

listen(func)[source]¶

Listen to parameters change.

Parameters: func (callable) – Function to be called when a parameter changes.

logdet()[source]¶

Log of ｜K｜.

Returns: Log-determinant of K.
Return type: float

property nparams¶: Number of parameters.

property shape¶: Array shape.

unfix()[source]¶: Enable parameter optimisation.

value()[source]¶

Covariance matrix.

Returns: K – Matrix K = LLᵀ + ϵI, for a very small positive number ϵ.
Return type: ndarray

class limix.model.cov.GivenCov(K0)[source]¶

Given covariance function, K = s⋅K₀.

The covariance matrix is the provided matrix K₀ scaled by s: K = s⋅K₀.

Example

>>> from glimix_core.cov import GivenCov
>>> from numpy import dot
>>> from numpy.random import RandomState
>>>
>>> G = RandomState(0).randn(5, 3)
>>> K0 = dot(G, G.T)
>>> cov = GivenCov(K0)
>>> cov.scale = 1.3
>>> cov.name = "K"
>>> print(cov)
GivenCov(K0=...): K
  scale: 1.3

gradient()[source]¶

Derivative of the covariance matrix over log(s).

Returns: logscale – s⋅K₀.
Return type: float

property scale¶: Scale parameter, s.

value()[source]¶

Covariance matrix, s⋅K₀.

Returns: K – s⋅K₀.
Return type: ndarray

class limix.model.cov.LinearCov(X)[source]¶

Linear covariance function, K = s⋅XXᵀ.

The mathematical representation is s⋅XXᵀ, for a n×r matrix X provided by the user and a scaling parameter s.

Example

>>> from glimix_core.cov import LinearCov
>>> from numpy import dot
>>> from numpy.random import RandomState
>>>
>>> X = RandomState(0).randn(2, 3)
>>> cov = LinearCov(X)
>>> cov.scale = 1.3
>>> cov.name = "K"
>>> print(cov)
LinearCov(): K
  scale: 1.3

property X¶: Matrix X from K = s⋅XXᵀ.

fix()[source]¶: Prevent s update during optimization.

gradient()[source]¶

Derivative of the covariance matrix over log(s).

Returns: logscale – s⋅XXᵀ.
Return type: ndarray

property scale¶: Scale parameter.

unfix()[source]¶: Enable s update during optimization.

value()[source]¶

Covariance matrix.

Returns: K – s⋅XXᵀ.
Return type: ndarray

class limix.model.cov.SumCov(covariances)[source]¶

Sum of covariance functions, K = K₀ + K₁ + ⋯.

Example

>>> from glimix_core.cov import LinearCov, SumCov
>>> from numpy.random import RandomState
>>>
>>> random = RandomState(0)
>>> cov_left = LinearCov(random.randn(4, 20))
>>> cov_right = LinearCov(random.randn(4, 15))
>>> cov_left.scale = 0.5
>>>
>>> cov = SumCov([cov_left, cov_right])
>>> cov_left.name = "A"
>>> cov_right.name = "B"
>>> cov.name = "A+B"
>>> print(cov)
SumCov(covariances=...): A+B
  LinearCov(): A
    scale: 0.5
  LinearCov(): B
    scale: 1.0

gradient()[source]¶

Sum of covariance function derivatives.

Returns: ∂K₀ + ∂K₁ + ⋯
Return type: dict

value()[source]¶

Sum of covariance matrices.

Returns: K – K₀ + K₁ + ⋯
Return type: ndarray

class limix.model.cov.LRFreeFormCov(n, m)[source]¶

General semi-definite positive matrix of low rank, K = LLᵀ.

The covariance matrix K is given by LLᵀ, where L is a n×m matrix and n≥m. Therefore, K will have rank(K) ≤ m.

Example

>>> from glimix_core.cov import LRFreeFormCov
>>> cov = LRFreeFormCov(3, 2)
>>> print(cov.L)
[[1. 1.]
 [1. 1.]
 [1. 1.]]
>>> cov.L = [[1, 2], [0, 3], [1, 3]]
>>> print(cov.L)
[[1. 2.]
 [0. 3.]
 [1. 3.]]
>>> cov.name = "F"
>>> print(cov)
LRFreeFormCov(n=3, m=2): F
  L: [[1. 2.]
      [0. 3.]
      [1. 3.]]

property L¶

Matrix L from K = LLᵀ.

Returns: L – Parametric matrix.
Return type: (n, m) ndarray

property Lu¶: Lower-triangular, flat part of L.

fix()[source]¶: Disable parameter optimisation.

gradient()[source]¶

Derivative of the covariance matrix over the lower triangular, flat part of L.

It is equal to

∂K/∂Lᵢⱼ = ALᵀ + LAᵀ,

where Aᵢⱼ is an n×m matrix of zeros except at [Aᵢⱼ]ᵢⱼ=1.

Returns: Lu – Derivative of K over the lower-triangular, flat part of L.
Return type: ndarray

listen(func)[source]¶

Listen to parameters change.

Parameters: func (callable) – Function to be called when a parameter changes.

property nparams¶: Number of parameters.

property shape¶: Array shape.

unfix()[source]¶: Enable parameter optimisation.

value()[source]¶

Covariance matrix.

Returns: K – K = LLᵀ.
Return type: (n, n) ndarray

class limix.model.cov.Kron2SumCov(G, dim, rank)[source]¶

Implements K = C₀ ⊗ GGᵀ + C₁ ⊗ I.

C₀ and C₁ are d×d symmetric matrices. C₀ is a semi-definite positive matrix while C₁ is a positive definite one. G is a n×m matrix and I is a n×n identity matrix. Let M = Uₘ Sₘ Uₘᵀ be the eigen decomposition for any matrix M. The documentation and implementation of this class make use of the following definitions:

X = GGᵀ = Uₓ Sₓ Uₓᵀ
C₁ = U₁ S₁ U₁ᵀ
Cₕ = S₁⁻½ U₁ᵀ C₀ U₁ S₁⁻½
Cₕ = Uₕ Sₕ Uₕᵀ
D = (Sₕ ⊗ Sₓ + Iₕₓ)⁻¹
Lₓ = Uₓᵀ
Lₕ = Uₕᵀ S₁⁻½ U₁ᵀ
L = Lₕ ⊗ Lₓ

The above definitions allows us to write the inverse of the covariance matrix as:

K⁻¹ = LᵀDL,

where D is a diagonal matrix.

Example

>>> from numpy import array
>>> from glimix_core.cov import Kron2SumCov
>>>
>>> G = array([[-1.5, 1.0], [-1.5, 1.0], [-1.5, 1.0]])
>>> Lr = array([[3], [2]], float)
>>> Ln = array([[1, 0], [2, 1]], float)
>>>
>>> cov = Kron2SumCov(G, 2, 1)
>>> cov.C0.L = Lr
>>> cov.C1.L = Ln
>>> print(cov)
Kron2SumCov(G=..., dim=2, rank=1): Kron2SumCov
  LRFreeFormCov(n=2, m=1): C₀
    L: [[3.]
        [2.]]
  FreeFormCov(dim=2): C₁
    L: [[1. 0.]
        [2. 1.]]

property C0¶: Semi-definite positive matrix C₀.

property C1¶: Definite positive matrix C₁.

property D¶: (Sₕ ⊗ Sₓ + Iₕₓ)⁻¹.

property G¶: User-provided matrix G, n×m.

Ge¶: Result of US from the SVD decomposition G = USVᵀ.

LdKL_dot(v)[source]¶

Implements L(∂K)Lᵀv.

The array v can have one or two dimensions and the first dimension has to have size n⋅p.

Let vec(V) = v. We have

L(∂K)Lᵀ⋅v = ((Lₕ∂C₀Lₕᵀ) ⊗ (LₓGGᵀLₓᵀ))vec(V) = vec(LₓGGᵀLₓᵀVLₕ∂C₀Lₕᵀ),

when the derivative is over the parameters of C₀. Similarly,

L(∂K)Lᵀv = ((Lₕ∂C₁Lₕᵀ) ⊗ (LₓLₓᵀ))vec(V) = vec(LₓLₓᵀVLₕ∂C₁Lₕᵀ),

over the parameters of C₁.

property Lh¶: Lₕ.

property Lx¶: Lₓ.

gradient()[source]¶

Gradient of K.

Returns

C0 (ndarray) – Derivative of C₀ over its parameters.
C1 (ndarray) – Derivative of C₁ over its parameters.

gradient_dot(v)[source]¶

Implements ∂K⋅v.

Parameters

v (array_like) – Vector from ∂K⋅v.

Returns

C0.Lu (ndarray) – ∂K⋅v, where the gradient is taken over the C₀ parameters.
C1.Lu (ndarray) – ∂K⋅v, where the gradient is taken over the C₁ parameters.

listen(func)[source]¶

Listen to parameters change.

Parameters: func (callable) – Function to be called when a parameter changes.

logdet()[source]¶

Implements log|K| = - log|D| + n⋅log|C₁|.

Returns: logdet – Log-determinant of K.
Return type: float

logdet_gradient()[source]¶

Implements ∂log|K| = Tr[K⁻¹∂K].

It can be shown that:

∂log|K| = diag(D)ᵀdiag(L(∂K)Lᵀ) = diag(D)ᵀ(diag(Lₕ∂C₀Lₕᵀ)⊗diag(LₓGGᵀLₓᵀ)),

when the derivative is over the parameters of C₀. Similarly,

∂log|K| = diag(D)ᵀdiag(L(∂K)Lᵀ) = diag(D)ᵀ(diag(Lₕ∂C₁Lₕᵀ)⊗diag(I)),

over the parameters of C₁.

Returns

C0 (ndarray) – Derivative of C₀ over its parameters.
C1 (ndarray) – Derivative of C₁ over its parameters.

property nparams¶: Number of parameters.

solve(v)[source]¶

Implements the product K⁻¹⋅v.

Parameters: v (array_like) – Array to be multiplied.
Returns: x – Solution x to the equation K⋅x = y.
Return type: ndarray

value()[source]¶

Covariance matrix K = C₀ ⊗ GGᵀ + C₁ ⊗ I.

Returns: K – C₀ ⊗ GGᵀ + C₁ ⊗ I.
Return type: ndarray

Mean¶

class limix.model.mean.OffsetMean(n)[source]¶

Offset mean function, θ⋅𝟏.

It represents a mean vector 𝐦 = θ⋅𝟏 of size n. The offset is given by the parameter θ.

Example

>>> from glimix_core.mean import OffsetMean
>>>
>>> mean = OffsetMean(3)
>>> mean.offset = 2.0
>>> print(mean.value())
[2. 2. 2.]
>>> print(mean.gradient())
{'offset': array([1., 1., 1.])}
>>> mean.name = "𝐦"
>>> print(mean)
OffsetMean(): 𝐦
  offset: 2.0

fix_offset()[source]¶: Prevent θ update during optimization.

gradient()[source]¶

Gradient of the offset function.

Returns: offset – Vector 𝟏.
Return type: (n,) ndarray

property offset¶: Offset parameter.

unfix_offset()[source]¶: Enable θ update during optimization.

value()[source]¶

Offset mean.

Returns: 𝐦 – θ⋅𝟏.
Return type: (n,) ndarray

class limix.model.mean.LinearMean(X)[source]¶

Linear mean function, X𝜶.

It defines X𝜶, for which X is a n×m matrix provided by the user and 𝜶 is a vector of size m.

Example

>>> from glimix_core.mean import LinearMean
>>>
>>> mean = LinearMean([[1.5, 0.2], [0.5, 0.4]])
>>> mean.effsizes = [1.0, -1.0]
>>> print(mean.value())
[1.3 0.1]
>>> print(mean.gradient()["effsizes"])
[[1.5 0.2]
 [0.5 0.4]]
>>> mean.name = "𝐦"
>>> print(mean)
LinearMean(m=2): 𝐦
  effsizes: [ 1. -1.]

property X¶: An n×m matrix of covariates.

property effsizes¶: Effect-sizes parameter, 𝜶, of size m.

gradient()[source]¶

Gradient of the linear mean function over the effect sizes.

Returns

effsizes –

Return type

(n, m) ndarray

value()[source]¶

Linear mean function.

Returns: 𝐦 – X𝜶.
Return type: (n,) ndarray

class limix.model.mean.SumMean(means)[source]¶

Sum mean function, 𝐟₀ + 𝐟₁ + ….

The mathematical representation is

𝐦 = 𝐟₀ + 𝐟₁ + …

In other words, it is a sum of mean vectors.

Example

>>> from glimix_core.mean import OffsetMean, LinearMean, SumMean
>>>
>>> X = [[5.1, 1.0],
...      [2.1, -0.2]]
>>>
>>> mean0 = LinearMean(X)
>>> mean0.effsizes = [-1.0, 0.5]
>>>
>>> mean1 = OffsetMean(2)
>>> mean1.offset = 2.0
>>>
>>> mean = SumMean([mean0, mean1])
>>>
>>> print(mean.value())
[-2.6 -0.2]
>>> g = mean.gradient()
>>> print(g["SumMean[0].effsizes"])
[[ 5.1  1. ]
 [ 2.1 -0.2]]
>>> print(g["SumMean[1].offset"])
[1. 1.]
>>> mean0.name = "A"
>>> mean1.name = "B"
>>> mean.name = "A+B"
>>> print(mean)
SumMean(means=...): A+B
  LinearMean(m=2): A
    effsizes: [-1.   0.5]
  OffsetMean(): B
    offset: 2.0

gradient()[source]¶

Sum of mean function derivatives.

Returns: ∂𝐦 – ∂𝐟₀ + ∂𝐟₁ + ….
Return type: dict

value()[source]¶

Sum of mean vectors, 𝐟₀ + 𝐟₁ + ….

Returns: 𝐦 – 𝐟₀ + 𝐟₁ + ….
Return type: ndarray

class limix.model.mean.KronMean(A, X)[source]¶

Kronecker mean function, (A⊗X)vec(B).

Let

n be the number of samples;
p the number of traits; and
c the number of covariates.

The mathematical representation is

𝐦 = (A⊗X)vec(B)

where A is a p×p trait design matrix of fixed effects and X is a n×c sample design matrix of fixed effects. B is a c×p matrix of fixed-effect sizes.

property A¶: Matrix A.

property AX¶: A ⊗ X.

property B¶: Effect-sizes parameter, B.

property X¶: Matrix X.

gradient()[source]¶

Gradient of the linear mean function.

Returns: vecB – Derivative of M over vec(B).
Return type: ndarray

property nparams¶: Number of parameters.

value()[source]¶

Kronecker mean function.

Returns: 𝐦 – (A⊗X)vec(B).
Return type: ndarray