BOLT-LMM mixed-model association

StructureVariance components of the linear mixed model#✓fully proved

The variance components of the linear mixed model y = X β + g + e with g ~ N(0, σ_g² K) and e ~ N(0, σ_e² I) on n individuals: a positive semidefinite genetic relationship matrix K, a genetic variance σ_g² ≥ 0 and an environmental variance σ_e² > 0.

structure Model (n : ℕ) where
  /-- The genetic relationship (kinship) matrix. -/
  K : Matrix (Fin n) (Fin n) ℝ
  /-- The genetic variance `σ_g²`. -/
  σg2 : ℝ
  /-- The environmental (residual) variance `σ_e²`. -/
  σe2 : ℝ
  posSemidef_K : K.PosSemidef
  σg2_nonneg : 0 ≤ σg2
  σe2_pos : 0 < σe2
Used by 19

DefinitionMarginal phenotypic covariance matrix#✓fully proved

The marginal phenotypic covariance V = σ_g² K + σ_e² I of y = X β + g + e.

noncomputable def V {n : ℕ} (M : Model n) : Matrix (Fin n) (Fin n) ℝ := M.σg2 • M.K + M.σe2 • 1
Uses 1 · Used by 14

TheoremPositive definiteness of the phenotypic covariance#✓fully proved

V = σ_g² K + σ_e² I is positive definite.

theorem posDef_V {n : ℕ} (M : Model n) : M.V.PosDef
Uses 2 · Used by 8

DefinitionMarginal log-likelihood of the mean#✓fully proved

The marginal log-likelihood of the mean μ: the log of the N(μ, V) density at y, -(n log 2π + log det V + (y - μ)ᵀ V⁻¹ (y - μ)) / 2 (det V > 0 by posDef_V).

noncomputable def logLik {n : ℕ} (M : Model n) (y μ : Fin n → ℝ) : ℝ :=
  -(n * Real.log (2 * Real.pi) + Real.log M.V.det + (y - μ) ⬝ᵥ (M.V⁻¹ *ᵥ (y - μ))) / 2
Uses 2 · Used by 5

DefinitionGeneralized least-squares estimator#✓fully proved

The generalized least-squares estimator β̂ = (Xᵀ V⁻¹ X)⁻¹ Xᵀ V⁻¹ y.

noncomputable def gls {n p : ℕ} (M : Model n) (X : Matrix (Fin n) (Fin p) ℝ) (y : Fin n → ℝ) :
    Fin p → ℝ :=
  (Xᵀ * M.V⁻¹ * X)⁻¹ *ᵥ (Xᵀ *ᵥ (M.V⁻¹ *ᵥ y))
Uses 2 · Used by 8

DefinitionResidual-forming matrix#✓fully proved

The residual-forming matrix P = V⁻¹ - V⁻¹ X (Xᵀ V⁻¹ X)⁻¹ Xᵀ V⁻¹.

noncomputable def proj {n p : ℕ} (M : Model n) (X : Matrix (Fin n) (Fin p) ℝ) :
    Matrix (Fin n) (Fin n) ℝ :=
  M.V⁻¹ - M.V⁻¹ * X * (Xᵀ * M.V⁻¹ * X)⁻¹ * Xᵀ * M.V⁻¹
Uses 2 · Used by 8

TheoremGLS maximizes the marginal likelihood#✓fully proved

For X of full column rank, every mean X β has marginal log-likelihood at most that of the GLS fit X β̂.

theorem logLik_le_logLik_gls {n p : ℕ} {M : Model n} {X : Matrix (Fin n) (Fin p) ℝ}
    (hX : Function.Injective X.mulVec) (y : Fin n → ℝ)
    (β : Fin p → ℝ) : M.logLik y (X *ᵥ β) ≤ M.logLik y (X *ᵥ M.gls X y)
Uses 5 · Used by 1

TheoremThe GLS maximizer is unique#✓fully proved

For X of full column rank, X β attains the maximal marginal log-likelihood only at β = β̂.

theorem logLik_eq_logLik_gls_iff {n p : ℕ} {M : Model n} {X : Matrix (Fin n) (Fin p) ℝ}
    (hX : Function.Injective X.mulVec) (y : Fin n → ℝ)
    (β : Fin p → ℝ) :
    M.logLik y (X *ᵥ β) = M.logLik y (X *ᵥ M.gls X y) ↔ β = M.gls X y
Uses 5

TheoremMaximal marginal log-likelihood#✓fully proved

The maximal marginal log-likelihood is -(n log 2π + log det V + yᵀ P y) / 2.

theorem logLik_gls {n p : ℕ} {M : Model n} {X : Matrix (Fin n) (Fin p) ℝ}
    (hX : Function.Injective X.mulVec) (y : Fin n → ℝ) :
    M.logLik y (X *ᵥ M.gls X y) =
      -(n * Real.log (2 * Real.pi) + Real.log M.V.det + y ⬝ᵥ (M.proj X *ᵥ y)) / 2
Uses 6 · Used by 1

TheoremResidual-forming matrix annihilates the fixed effects#✓fully proved

P X = 0: the residual-forming matrix annihilates the fixed effects.

theorem proj_mul_X {n p : ℕ} {M : Model n} {X : Matrix (Fin n) (Fin p) ℝ}
    (hX : Function.Injective X.mulVec) : M.proj X * X = 0
Uses 4 · Used by 1

TheoremIdentity \(PVP = P\) for the residual-forming matrix#✓fully proved

P V P = P. With P X = 0, this says that under the null model y = X β + ε with Cov ε = V the score xᵀ P y = xᵀ P ε has variance xᵀ P V P x = xᵀ P x.

theorem proj_mul_V_mul_proj {n p : ℕ} {M : Model n} {X : Matrix (Fin n) (Fin p) ℝ}
    (hX : Function.Injective X.mulVec) :
    M.proj X * M.V * M.proj X = M.proj X
Uses 5

DefinitionProfile log-likelihood of a SNP effect#✓fully proved

The profile marginal log-likelihood of a SNP effect γ for the SNP genotype vector x: sup_β logLik y (X β + γ x).

noncomputable def profileLogLik {n p : ℕ} (M : Model n)
    (X : Matrix (Fin n) (Fin p) ℝ) (x y : Fin n → ℝ) (γ : ℝ) : ℝ :=
  ⨆ β : Fin p → ℝ, M.logLik y (X *ᵥ β + γ • x)
Uses 2 · Used by 2

DefinitionScore statistic for one SNP#✓fully proved

The score statistic (xᵀ P y)² / (xᵀ P x) for one SNP x. Its denominator is positive when x is not in the column space of X (info_pos); otherwise the statistic is 0 by the convention a / 0 = 0 and carries no information.

noncomputable def scoreStat {n p : ℕ} (M : Model n) (X : Matrix (Fin n) (Fin p) ℝ)
    (x y : Fin n → ℝ) : ℝ :=
  (x ⬝ᵥ (M.proj X *ᵥ y)) ^ 2 / (x ⬝ᵥ (M.proj X *ᵥ x))
Uses 2

TheoremThe profile log-likelihood of a SNP effect is an exact quadratic#✓fully proved

ℓ(γ) = ℓ(0) + γ xᵀ P y - γ² xᵀ P x / 2, so its score at γ = 0 is xᵀ P y and its information is xᵀ P x.

theorem profileLogLik_eq_quadratic {n p : ℕ} {M : Model n} {X : Matrix (Fin n) (Fin p) ℝ}
    (hX : Function.Injective X.mulVec) (x y : Fin n → ℝ)
    (γ : ℝ) :
    M.profileLogLik X x y γ = M.profileLogLik X x y 0 + γ * (x ⬝ᵥ (M.proj X *ᵥ y)) -
      γ ^ 2 / 2 * (x ⬝ᵥ (M.proj X *ᵥ x))
Uses 9 · Used by 1

TheoremThe score of a SNP effect#✓fully proved

The score of a SNP effect at γ = 0 is xᵀ P y.

theorem hasDerivAt_profileLogLik {n p : ℕ} {M : Model n} {X : Matrix (Fin n) (Fin p) ℝ}
    (hX : Function.Injective X.mulVec) (x y : Fin n → ℝ) :
    HasDerivAt (M.profileLogLik X x y) (x ⬝ᵥ (M.proj X *ᵥ y)) 0
Uses 4

TheoremScore at the null generalized least-squares fit#✓fully proved

The score is evaluated at the null GLS fit: xᵀ P y = xᵀ V⁻¹ (y - X β̂).

theorem score_eq {n p : ℕ} (M : Model n) (X : Matrix (Fin n) (Fin p) ℝ) (x y : Fin n → ℝ) :
    x ⬝ᵥ (M.proj X *ᵥ y) = x ⬝ᵥ (M.V⁻¹ *ᵥ (y - X *ᵥ M.gls X y))
Uses 4

TheoremInformation of a SNP outside the column space of \(X\) is positive#✓fully proved

The information xᵀ P x is positive when the SNP x is not in the column space of X.

theorem info_pos {n p : ℕ} {M : Model n} {X : Matrix (Fin n) (Fin p) ℝ}
    (hX : Function.Injective X.mulVec) {x : Fin n → ℝ}
    (hx : x ∉ Set.range X.mulVec) : 0 < x ⬝ᵥ (M.proj X *ᵥ x)
Uses 5

DefinitionBest linear unbiased predictor of genetic effects#✓fully proved

The BLUP ĝ = G V⁻¹ (y - X β̂) of the genetic effects, with G = σ_g² K (Henderson 1975).

noncomputable def blup {n p : ℕ} (M : Model n) (X : Matrix (Fin n) (Fin p) ℝ)
    (y : Fin n → ℝ) : Fin n → ℝ :=
  (M.σg2 • M.K) *ᵥ (M.V⁻¹ *ᵥ (y - X *ᵥ M.gls X y))
Uses 3 · Used by 1

DefinitionSolution of Henderson's mixed-model equations#✓fully proved

(β, g) solves Henderson's mixed-model equations with R = σ_e² I and G = σ_g² K, i.e. the specialization of MixedModel.IsMMESolution to Z = 1 (every individual has its own genetic effect), G = σ_g² K, R = σ_e² I. In unfolded form this is Xᵀ R⁻¹ X β + Xᵀ R⁻¹ g = Xᵀ R⁻¹ y and R⁻¹ X β + (R⁻¹ + G⁻¹) g = R⁻¹ y. The equations are meaningful when G is invertible (isMMESolution_iff); otherwise G⁻¹ is Lean's junk inverse 0.

def IsMMESolution {n p : ℕ} (M : Model n) (X : Matrix (Fin n) (Fin p) ℝ)
    (y : Fin n → ℝ) (β : Fin p → ℝ) (g : Fin n → ℝ) : Prop :=
  MathlibExt.Biology.StatisticalGenetics.MixedModel.IsMMESolution X 1 (M.σg2 • M.K)
    (M.σe2 • 1) y β g
Uses 2 · Used by 1

TheoremHenderson's mixed-model equations give the GLS estimator and the BLUP#✓fully proved

For positive definite K, σ_g² > 0 and X of full column rank, (β, g) solves the mixed-model equations iff β = β̂ and g = σ_g² K V⁻¹ (y - X β̂). This is the specialization of MixedModel.isMMESolution_iff_gls to Z = 1, G = σ_g² K, R = σ_e² I; the generic theorem gives β = glsMatrix X V *ᵥ y and g = G V⁻¹ (y - X β̂), which are Model.gls and Model.blup by the two conversion lemmas below.

theorem isMMESolution_iff {n p : ℕ} {M : Model n} {X : Matrix (Fin n) (Fin p) ℝ}
    (hK : M.K.PosDef) (hg : 0 < M.σg2) (hX : Function.Injective X.mulVec)
    (y : Fin n → ℝ) (β : Fin p → ℝ) (g : Fin n → ℝ) :
    M.IsMMESolution X y β g ↔ β = M.gls X y ∧ g = M.blup X y
Uses 10

DefinitionTwo-Gaussian mixture prior on a SNP effect#✓fully proved

The BOLT-LMM prior on a SNP effect: the mixture p N(0, v₁) + (1 - p) N(0, v₂) of two centered Gaussians with variances v₁, v₂ (v₂ = 0 gives the spike-and-slab prior with a point mass at 0).

noncomputable def mixturePrior (p v₁ v₂ : ℝ≥0) : Measure ℝ :=
  (p : ℝ≥0∞) • gaussianReal 0 v₁ + ((1 - p : ℝ≥0) : ℝ≥0∞) • gaussianReal 0 v₂
Used by 2

TheoremMixture prior is a probability measure#✓fully proved

The mixture prior is a probability measure for a mixing weight p ≤ 1.

theorem isProbabilityMeasure_mixturePrior {p : ℝ≥0} (hp : p ≤ 1) (v₁ v₂ : ℝ≥0) :
    IsProbabilityMeasure (mixturePrior p v₁ v₂)
Uses 1

TheoremMixture prior has finite second moment#✓fully proved

The mixture prior has a finite second moment.

theorem integrable_sq_mixturePrior (p v₁ v₂ : ℝ≥0) :
    Integrable (fun b : ℝ => b ^ 2) (mixturePrior p v₁ v₂)
Uses 1

StructureAdmissible variational factor#✓fully proved

An admissible variational factor for the prior π: a probability measure with a finite second moment and finite KL divergence from π (q ≪ π with integrable log-likelihood ratio).

structure IsAdmissible (π q : Measure ℝ) : Prop where
  isProbabilityMeasure : IsProbabilityMeasure q
  integrable_sq : Integrable (fun b : ℝ => b ^ 2) q
  absolutelyContinuous : q ≪ π
  integrable_llr : Integrable (llr q π) q
Used by 6

DefinitionMean of a variational factor#✓fully proved

The mean ∫ b dq of a factor.

noncomputable def mean (q : Measure ℝ) : ℝ := ∫ b, b ∂q
Used by 4

DefinitionSquared norm of a genotype column#✓fully proved

The squared norm ‖xⱼ‖² = ∑ᵢ X i j ² of the genotype column of SNP j.

def colNormSq {n m : ℕ} (X : Matrix (Fin n) (Fin m) ℝ) (j : Fin m) : ℝ := ∑ i, X i j ^ 2
Used by 7

DefinitionExpected residual sum of squares under the variational factors#✓fully proved

The expected residual sum of squares ‖y - X μ‖² + ∑ⱼ ‖xⱼ‖² vⱼ of the product of the factors q, with means μⱼ and variances vⱼ (integral_rss_pi).

noncomputable def expectedRSS {n m : ℕ} (X : Matrix (Fin n) (Fin m) ℝ) (y : Fin n → ℝ)
    (q : Fin m → Measure ℝ) : ℝ :=
  ∑ i, (y i - ∑ j, X i j * mean (q j)) ^ 2 + ∑ j, colNormSq X j * variance id (q j)
Uses 2 · Used by 3

TheoremThe expected residual sum of squares of a product of factors#✓fully proved

For probability factors with finite second moments, 𝔼_{∏ qⱼ} ‖y - X β‖² = ‖y - X μ‖² + ∑ⱼ ‖xⱼ‖² vⱼ.

theorem integral_rss_pi {n m : ℕ} (X : Matrix (Fin n) (Fin m) ℝ) (y : Fin n → ℝ)
    (q : Fin m → Measure ℝ) [∀ j, IsProbabilityMeasure (q j)]
    (h2 : ∀ j, Integrable (fun b : ℝ => b ^ 2) (q j)) :
    ∫ β, ∑ i, (y i - ∑ j, X i j * β j) ^ 2 ∂(Measure.pi q) = expectedRSS X y q
Uses 3

DefinitionMean-field evidence lower bound#✓fully proved

The mean-field evidence lower bound -(n/2) log (2π σ_e²) - 𝔼_q ‖y - X β‖² / (2 σ_e²) - ∑ⱼ KL(qⱼ ‖ π), meaningful for σ_e² > 0. It is the ELBO for admissible factors (finite KL); for an inadmissible factor toReal would turn an infinite divergence into 0, so all results assume admissibility.

noncomputable def elbo {n m : ℕ} (X : Matrix (Fin n) (Fin m) ℝ) (y : Fin n → ℝ) (σe2 : ℝ)
    (π : Measure ℝ) (q : Fin m → Measure ℝ) : ℝ :=
  -(n * Real.log (2 * Real.pi * σe2)) / 2 - expectedRSS X y q / (2 * σe2) -
    ∑ j, (klDiv (q j) π).toReal
Uses 1 · Used by 5

DefinitionPartial residual of a SNP#✓fully proved

The partial residual rⱼ = y - ∑_{k ≠ j} xₖ μₖ of the other SNPs' posterior means.

noncomputable def partialResidual {n m : ℕ} (X : Matrix (Fin n) (Fin m) ℝ) (y : Fin n → ℝ)
    (q : Fin m → Measure ℝ) (j : Fin m) (i : Fin n) : ℝ :=
  y i - ∑ k ∈ Finset.univ.erase j, X i k * mean (q k)
Uses 1 · Used by 5

DefinitionCoordinate-ascent potential for one SNP effect#✓fully proved

The CAVI potential fⱼ(b) = (b xⱼᵀ rⱼ - ‖xⱼ‖² b² / 2) / σ_e²: the part of the expected log-likelihood 𝔼_{q_{-j}}[log p(y | β)] that depends on βⱼ = b (exists_elbo_update_eq).

noncomputable def potential {n m : ℕ} (X : Matrix (Fin n) (Fin m) ℝ) (y : Fin n → ℝ) (σe2 : ℝ)
    (q : Fin m → Measure ℝ) (j : Fin m) (b : ℝ) : ℝ :=
  (b * ∑ i, X i j * partialResidual X y q j i - colNormSq X j * b ^ 2 / 2) / σe2
Uses 2 · Used by 5

DefinitionCoordinate-ascent update of one variational factor#✓fully proved

The CAVI update of factor j: the prior tilted by the potential, qⱼ ∝ π · exp fⱼ = exp 𝔼_{q_{-j}}[log p(y, β)].

noncomputable def caviUpdate {n m : ℕ} (X : Matrix (Fin n) (Fin m) ℝ) (y : Fin n → ℝ) (σe2 : ℝ)
    (π : Measure ℝ) (q : Fin m → Measure ℝ) (j : Fin m) : Measure ℝ :=
  π.tilted (potential X y σe2 q j)
Uses 1 · Used by 6

DefinitionOne coordinate-ascent step#✓fully proved

One coordinate-ascent step: replace factor j by its CAVI update.

noncomputable def caviStep {n m : ℕ} (X : Matrix (Fin n) (Fin m) ℝ) (y : Fin n → ℝ) (σe2 : ℝ)
    (π : Measure ℝ) (q : Fin m → Measure ℝ) (j : Fin m) : Fin m → Measure ℝ :=
  Function.update q j (caviUpdate X y σe2 π q j)
Uses 1 · Used by 4

DefinitionCoordinate ascent along an update schedule#✓fully proved

Coordinate ascent along the update schedule s: at time t factor s t is updated.

noncomputable def caviRun {n m : ℕ} (X : Matrix (Fin n) (Fin m) ℝ) (y : Fin n → ℝ) (σe2 : ℝ)
    (π : Measure ℝ) (q : Fin m → Measure ℝ) (s : ℕ → Fin m) : ℕ → Fin m → Measure ℝ
  | 0 => q
  | t + 1 => caviStep X y σe2 π (caviRun X y σe2 π q s t) (s t)
Uses 1 · Used by 1

TheoremThe ELBO as a function of one factor#✓fully proved

With the other factors fixed, the ELBO of q with factor j replaced by an admissible q' is ∫ fⱼ dq' - KL(q' ‖ π) + C.

theorem exists_elbo_update_eq {n m : ℕ} {X : Matrix (Fin n) (Fin m) ℝ} {y : Fin n → ℝ}
    {σe2 : ℝ} {π : Measure ℝ}
    (q : Fin m → Measure ℝ) (j : Fin m) :
    ∃ C : ℝ, ∀ q', IsAdmissible π q' →
      elbo X y σe2 π (Function.update q j q') =
        ∫ b, potential X y σe2 q j b ∂q' - (klDiv q' π).toReal + C
Uses 7 · Used by 1

TheoremCoordinate-ascent update is an admissible factor#✓fully proved

The CAVI update is an admissible factor, for a probability prior with a finite second moment.

theorem isAdmissible_caviUpdate {n m : ℕ} {X : Matrix (Fin n) (Fin m) ℝ} {y : Fin n → ℝ}
    {σe2 : ℝ} {π : Measure ℝ}
    [IsProbabilityMeasure π]
    (hπ : Integrable (fun b : ℝ => b ^ 2) π) (hσ : 0 < σe2) (q : Fin m → Measure ℝ)
    (j : Fin m) : IsAdmissible π (caviUpdate X y σe2 π q j)
Uses 5 · Used by 3

TheoremThe CAVI update maximizes the ELBO over its factor, with exact gap#✓fully proved

Replacing factor j by an admissible q' gives ELBO(caviStep q j) - KL(q' ‖ caviUpdate q j).

theorem elbo_update_eq_sub {n m : ℕ} {X : Matrix (Fin n) (Fin m) ℝ} {y : Fin n → ℝ}
    {σe2 : ℝ} {π : Measure ℝ}
    [IsProbabilityMeasure π] (hπ : Integrable (fun b : ℝ => b ^ 2) π)
    (hσ : 0 < σe2) (q : Fin m → Measure ℝ) (j : Fin m) {q' : Measure ℝ}
    (hq' : IsAdmissible π q') :
    elbo X y σe2 π (Function.update q j q') =
      elbo X y σe2 π (caviStep X y σe2 π q j) - (klDiv q' (caviUpdate X y σe2 π q j)).toReal
Uses 9 · Used by 2

TheoremCoordinate ascent does not decrease the ELBO#✓fully proved

Coordinate ascent does not decrease the ELBO.

theorem elbo_le_elbo_caviStep {n m : ℕ} {X : Matrix (Fin n) (Fin m) ℝ} {y : Fin n → ℝ}
    {σe2 : ℝ} {π : Measure ℝ}
    {q : Fin m → Measure ℝ} [IsProbabilityMeasure π]
    (hπ : Integrable (fun b : ℝ => b ^ 2) π) (hσ : 0 < σe2)
    {j : Fin m} (hqj : IsAdmissible π (q j)) :
    elbo X y σe2 π q ≤ elbo X y σe2 π (caviStep X y σe2 π q j)
Uses 5 · Used by 1

TheoremThe ELBO is nondecreasing along coordinate ascent#✓fully proved

The ELBO is nondecreasing along coordinate ascent, for any update schedule.

theorem monotone_elbo_caviRun {n m : ℕ} {X : Matrix (Fin n) (Fin m) ℝ} {y : Fin n → ℝ}
    {σe2 : ℝ} {π : Measure ℝ}
    {q : Fin m → Measure ℝ} [IsProbabilityMeasure π]
    (hπ : Integrable (fun b : ℝ => b ^ 2) π) (hσ : 0 < σe2)
    (hq : ∀ k, IsAdmissible π (q k)) (s : ℕ → Fin m) :
    Monotone fun t => elbo X y σe2 π (caviRun X y σe2 π q s t)
Uses 6

TheoremStationary points of coordinate ascent#✓fully proved

No admissible change of factor j increases the ELBO iff qⱼ is the expected-log-joint update π · exp fⱼ / Z.

theorem isCoordMax_iff {n m : ℕ} {X : Matrix (Fin n) (Fin m) ℝ} {y : Fin n → ℝ}
    {σe2 : ℝ} {π : Measure ℝ}
    {q : Fin m → Measure ℝ} [IsProbabilityMeasure π]
    (hπ : Integrable (fun b : ℝ => b ^ 2) π) (hσ : 0 < σe2)
    {j : Fin m} (hqj : IsAdmissible π (q j)) :
    (∀ q', IsAdmissible π q' → elbo X y σe2 π (Function.update q j q') ≤ elbo X y σe2 π q) ↔
      q j = caviUpdate X y σe2 π q j
Uses 9
40 fully proved

The dependency explorer needs JavaScript and its graph data.

Download the complete graph data.