Повний BLUP walk-through: 8 трейтів, mating type, G×L, multi-year
Симулятор проходить повний BLUP-пайплайн від сирих даних до присвоєння класу (CD-043). Вводиш матки з родоводом, mating type мам, пасіки, роки, виміри 8 трейтів → отримуєш EBV per queen per trait, reliability, СПЦ (сумарну племінну цінність) і клас.
y_ijkl = μ + L_j + Y_k + GL_ij + a_i + e_ijkl
| Компонент | Природа | Що означає |
|---|---|---|
| μ | константа | Середнє популяції. |
| L_j | фіксований | Ефект j-ї пасіки. |
| Y_k | фіксований | Ефект k-го року (сезону). |
| GL_ij | випадковий | Взаємодія «мати × пасіка» (~N(0, σ²_GL)). |
| a_i | випадковий | Адитивна генетика i-ї матки (~N(0, σ²_A)). Це EBV. |
| e_ijkl | випадковий | Залишок (~N(0, σ²_e)). |
Введи свої матки/родовід/фенотипи. Або обери пресет.
ȳ_i = (y_i1 + y_i2 + ... + y_in) / n_i
a_ii = 1 (сама з собою) a_ij = a_within[мама.mating_type] (напівсестри/супер-сестри тощо) a_ij = 0 (не споріднені)
σ_p обчислюємо з даних (variance всіх y). σ²_L, σ²_GL — з ANOVA-декомпозиції. σ²_A через h² літ. σ²_e — залишок.
σ²_p = variance(усі ȳ_i) (дисперсія фенотипу — загальна) σ²_L = between-apiary variance (дисперсія ефекту пасіки) σ²_GL = interaction variance (дисперсія взаємодії мати × пасіка) σ²_A = h² · σ²_p (адитивна дисперсія — з літ. h²) σ²_e = σ²_p − σ²_A − σ²_L − σ²_GL (залишкова дисперсія / шум)
μ = mean(усі ȳ_i) L_j = mean(ȳ_i на пасіці j) − μ GL_ij = mean(ȳ у клітинці мати_i, пасіка_j) − μ − L_j − mother_effect_i
d_i = ȳ_i − μ − L_j − GL_ij
(diag(n_i) + λ·A⁻¹) · û = n_i · d_i де λ = σ²_e / σ²_a
r²_i = 1 − PEV_i / σ²_a, PEV_i = діагональ (LHS⁻¹) · σ²_e
EBV_std_i = EBV_i / σ_a
Тут ВЖЕ використовуємо ВСІ 8 трейтів, не тільки обраний. Кроки 1-7 для кожного з інших 7 трейтів виконуються так само (не показуємо детально — тільки результати).
СПЦ_i = 100 + 10 · Σ_t (w_t · EBV_std_i,t)
(сума по всіх трейтах t, w_t — ваги з розділу 3.1)За замовчуванням сортовано за СПЦ спадаючи. Клік на заголовок колонки — змінити сортування. Фільтри — уточнити список.
Розділ для захисту: пояснює чому MME виглядає саме так, звідки береться λ = σ²_e/σ²_a, чому PEV = diag(C⁻¹)·σ²_e, і як BLUP співвідноситься з байєсівським оцінюванням. Це не ще одні обчислення — це логіка позаду формул із §9-§10.
y = Xβ + Zu + e u ~ N(0, A · σ²_a) (апріорний розподіл племінних цінностей) e ~ N(0, I · σ²_e) (залишок)
| Символ | Розмір | Що означає |
|---|---|---|
| y | n×1 | Вектор фенотипів (усі спостереження, не по маток). |
| β | p×1 | Фіксовані ефекти (μ, L_j, Y_k). Це НЕ «те, чим ми керуємо» — це систематичні негенетичні ефекти, які ми спостерігаємо у даних і оцінюємо (пасіка, рік). |
| u | q×1 | Випадкові племінні цінності (по одному на матку). Це саме те, що ми хочемо оцінити — вектор EBV. |
| X | n×p | Матриця інциденції фіксованих ефектів: X[i, j] = 1, якщо спостереження i підпадає під фіксований рівень j. |
| Z | n×q | Матриця інциденції випадкових ефектів: Z[i, j] = 1, якщо спостереження i належить матці j. Це НЕ «вулик-матка», а «спостереження-матка». |
| A | q×q | Матриця споріднення (з §5). Задає структуру коваріацій апріорі: Cov(u_i, u_j) = a_ij · σ²_a. |
| e | n×1 | Залишок: усе, що не пояснено моделлю. |
Розглянемо простий випадок: одна матка, n спостережень, немає родичів, немає фіксованих ефектів. Апріорі u ~ N(0, σ²_a). Умовно на u: y_k = u + e_k, e_k ~ N(0, σ²_e).
û = (n/σ²_e) · Σy_k / (n/σ²_e + 1/σ²_a) = n · d / (n + σ²_e/σ²_a) = n · d / (n + λ) де d = ȳ − μ
λ = 0.5 (чисті дані): û = 1·10 / (1 + 0.5) = 10 / 1.5 = 6.67 λ = 5.0 (шумні дані): û = 1·10 / (1 + 5) = 10 / 6 = 1.67Одне й те саме сире відхилення d=10, а BLUP-оцінки відрізняються у 4 рази — бо байєсівський shrinkage «не довіряє» шумним даним і повертає оцінку ближче до 0 (популяційне середнє).
| X'X X'Z | | β̂ | | X'y | | | | | = | | | Z'X Z'Z + A⁻¹ · λ | | û | | Z'y |
Ці рівняння не постульовані, а виводяться з мінімізації штрафної функції:
Q(β, u) = ‖y − Xβ − Zu‖² / σ²_e + u' A⁻¹ u / σ²_a
| Доданок | Що штрафує | Байєсівський сенс |
|---|---|---|
| ‖y − Xβ − Zu‖² / σ²_e | Погана підгонка даних (residual sum of squares, зважений оберненою дисперсією шуму). | −2·log p(y|β,u) з точністю до константи (правдоподібність). |
| u' A⁻¹ u / σ²_a | Великі u, що суперечать апріорній структурі A. | −2·log p(u) з точністю до константи (апріорний розподіл N(0, Aσ²_a)). |
(Z'Z + A⁻¹ · λ) · û = Z'(y − Xβ̂) → (diag(n_i) + λ · A⁻¹) · û = n_i · d_iЦе та сама формула, що у §9. Спрощення дидактичне; повна MME розв'язує β та u одночасно.
Позначимо LHS = C (coefficient matrix). Тоді C⁻¹ — масштабована коваріаційна матриця похибок оцінювання:
Var(û − u) = C⁻¹ · σ²_e (це стандартний результат mixed model theory) PEV_i = Var(û_i − u_i) = C⁻¹[i,i] · σ²_e
r²_i = 1 − PEV_i / σ²_a
= 1 − C⁻¹[i,i] · σ²_e / σ²_a
= 1 − C⁻¹[i,i] · λ
Три числові вправи, які можна порахувати руками за 5-10 хвилин і перевірити результат у симуляторі (§3-§10). Це сита: якщо руками результат сходиться з симулятором — значить, розумієш формули; якщо ні — треба повернутись до §15 і виявити, де саме розуміння відхилилось.
free → a_within = 0.25.Про parent-offspring у бджолиному BLUP: у стандартній NRM (Wright 1922) parent-offspring = 0.5. Але у спрощеному bee-BLUP (наш симулятор + Bienefeld 2007) drone-батько не відстежується, а мати-цариця часто не оцінюється як окрема запис у A — тому у A-матриці ми враховуємо лише sister-relationships (a_within від mating type матері). Ця вправа саме про це.
Крок 1. Діагональ: a_ii = 1 для всіх (кожна матка з собою).
a_11 = a_22 = a_33 = 1
Крок 2. D1↔D2, D1↔D3, D2↔D3 — усі три пари напівсестер через M1.
За mating type «free» a_within = 0.25 для кожної пари.
a_12 = a_13 = a_23 = 0.25
Крок 3. Симетрія: A[i,j] = A[j,i]. Заповнюємо нижній трикутник.
Крок 4. Результат:
D1 D2 D3
D1 [ 1.00 0.25 0.25 ]
D2 [ 0.25 1.00 0.25 ]
D3 [ 0.25 0.25 1.00 ]
Порівняння з іншими mating types (тільки для розуміння):
free → a_within = 0.25 (наш випадок)
МПТ → a_within = 0.375 (материнська селекція трутнів)
II з групою братів → a_within = 0.50 (штучне осіменіння sister-drones)
II з 1 трутнем → a_within = 0.75 (super-sisters, все одне)
Крок 1. Відхилення від середнього: d_1 = y_1 − μ = 50 − 40 = 10 kg Це «наївна» оцінка — якби ми довіряли одному вимірюванню повністю. Крок 2. Обчислюємо λ (Bayes shrinkage factor): λ = σ²_e / σ²_a = 70 / 30 = 2.333 Крок 3. MME для однієї матки без родичів редукується до: (n + λ) · û = n · d де n = 1 (одне спостереження). Крок 4. Підставляємо: (1 + 2.333) · û = 1 · 10 3.333 · û = 10 û = 10 / 3.333 = 3.0 kg Крок 5. Інтерпретація: Сире відхилення 10 kg, а BLUP стягує до 3.0 kg. Частка «справжньої генетики» = 1 / (1 + λ) = 1 / 3.333 = 0.30 = 30%. Це і є h² в спрощеному випадку: 30/(30+70) = 0.30. Тобто BLUP каже: «з 10 kg відхилення тільки 3 kg — генетика, решта 7 kg — шум».
Крок 1. n = 6, d = 10 (не змінилося), λ = 2.333 (не змінилося). Крок 2. MME: (n + λ) · û = n · d (6 + 2.333) · û = 6 · 10 8.333 · û = 60 û = 60 / 8.333 = 7.2 kg Крок 3. PEV (діагональ LHS⁻¹ · σ²_e): Для 1×1 матриці LHS = 8.333, тому LHS⁻¹ = 1/8.333 = 0.120. PEV = 0.120 · σ²_e = 0.120 · 70 = 8.4 kg² Крок 4. Reliability: r² = 1 − PEV / σ²_a = 1 − 8.4 / 30 = 1 − 0.28 = 0.72 Крок 5. Інтерпретація (важливо для захисту): — При n=1: û = 3.0 (30% від d), r² = 1 − 70/(30·3.333) ≈ 0.30. — При n=6: û = 7.2 (72% від d), r² = 0.72. 6 спостережень майже перебороли байєсівське стягування. r² 0.72 — це відмінно (за CD-043 поріг ≥0.30 для присвоєння класу). Правило пам'яті: чим більше n, тим ближче û до d і тим вища r². Асимптотично при n→∞: û → d, r² → 1.
Обчислюється на клієнті. Дата: .