# 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([, ])
>>> 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([, ])
>>> 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

𝐲 ∼ 𝓝(X𝜷, s((1-𝛿)K + 𝛿I)),


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

$\begin{split}\overbrace{[\mathrm Q₀ \quad \mathrm Q₁]}^{\mathrm Q} \overbrace{\left[\begin{array}{cc} \mathrm S₀ & 𝟎\\ 𝟎 & 𝟎 \end{array}\right]}^{\mathrm S} \left[\begin{array}{c} \mathrm Q₀ᵀ \\ \mathrm Q₁ᵀ \end{array}\right] = \mathrm K\end{split}$

and let

$\begin{split}\mathrm D = \left[ \begin{array}{cc} (1-𝛿)\mathrm S₀ + 𝛿\mathrm I & 𝟎\\ 𝟎 & 𝛿\mathrm I \end{array} \right].\end{split}$

We thus have

((1-𝛿)K + 𝛿I)⁻¹ = QD⁻¹Qᵀ.


A diagonal covariance-matrix can then be used to define an equivalent marginal likelihood:

𝓝(Qᵀ𝐲|QᵀX𝜷, sD).

property X

Covariates matrix.

Returns

X – Covariates.

Return type

ndarray

property beta

Fixed-effect sizes.

Returns

effect-sizes – Optimal fixed-effect sizes.

Return type

numpy.ndarray

Notes

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

(QᵀX)ᵀD⁻¹(QᵀX)𝜷 = (QᵀX)ᵀD⁻¹(Qᵀ𝐲).

property beta_covariance

Estimates the covariance-matrix of the optimal beta.

Returns

beta-covariance – (Xᵀ(s((1-𝛿)K + 𝛿I))⁻¹X)⁻¹.

Return type

ndarray

References

covariance()[source]

Covariance of the prior.

Returns

covariance – v₀K + v₁I.

Return type

ndarray

property delta

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)[source]

Disable parameter optimization.

Parameters

param (str) – Possible values are "delta", "beta", and "scale".

get_fast_scanner()[source]

Return FastScanner for association scan.

Returns

fast-scanner – Instance of a class designed to perform very fast association scan.

Return type

FastScanner

gradient()[source]

Not implemented.

lml()[source]

Log of the marginal likelihood.

Returns

lml – Log of the marginal likelihood.

Return type

float

Notes

The log of the marginal likelihood is given by

2⋅log(p(𝐲)) = -n⋅log(2π) - n⋅log(s) - log|D| - (Qᵀ𝐲)ᵀs⁻¹D⁻¹(Qᵀ𝐲)
+ (Qᵀ𝐲)ᵀs⁻¹D⁻¹(QᵀX𝜷)/2 - (QᵀX𝜷)ᵀs⁻¹D⁻¹(QᵀX𝜷).


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

2⋅log(p(𝐲)) = -n⋅log(2π) - n⋅log(s) - log|D| + (Qᵀ𝐲)ᵀs⁻¹D⁻¹Qᵀ(X𝜷-𝐲).


In the extreme case where 𝜷 is such that 𝐲 = X𝜷, 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|D| - n.

mean()[source]

Mean of the prior.

Formally, 𝐦 = X𝜷.

Returns

mean – Mean of the prior.

Return type

ndarray

property ncovariates

Number of covariates, c.

property nsamples

Number of samples, n.

property scale

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ᵀ(𝐲-X𝜷).


In the case of restricted marginal likelihood

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


where s is the number of covariates.

unfix(param)[source]

Enable parameter optimization.

Parameters

param (str) – Possible values are "delta", "beta", and "scale".

property v0

First variance.

Returns

v0 – s(1 - 𝛿).

Return type

float

property v1

Second variance.

Returns

v1 – s𝛿.

Return type

float

value()[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]

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

unfix(var_name)[source]

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.random import GGPSampler
>>>
>>> random = RandomState(1)
>>>
>>> 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]]
>>> 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.]]
[[[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([, ], 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.

property Ge

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

LdKL_dot(v, v1=None)[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]

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.]
{'offset': array([1., 1., 1.])}
>>> mean.name = "𝐦"
>>> print(mean)
OffsetMean(): 𝐦
offset: 2.0

fix_offset()[source]

Prevent θ update during optimization.

gradient()[source]

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]
[[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]
>>> print(g["SumMean.effsizes"])
[[ 5.1  1. ]
[ 2.1 -0.2]]
>>> print(g["SumMean.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