187 lines
7.5 KiB
Python
187 lines
7.5 KiB
Python
"""单因素随机效应 ANOVA + 广义遗传力 H²（纯 Python，numpy + fdist）。
|
||
|
||
模型：y_ij = μ + a_i + e_ij，a_i ~ N(0, σ²g)，e_ij ~ N(0, σ²e)。
|
||
分组（group）= 家系/组合，用于估计家系方差组分：
|
||
n₀ = (N - Σ n_i²/N)/(k-1) 加权平均每组样本量
|
||
σ²e = MSW 组内均方 = 环境/误差方差
|
||
σ²g = (MSB - MSW)/n₀ 组间方差组分（家系遗传方差）
|
||
H² = σ²g/(σ²g + σ²e) 广义遗传力（单株 basis）
|
||
F = MSB/MSW，p = 1-P(F≤f)
|
||
|
||
鲁棒性：组数/样本不足、组内方差为 0（全组同值）、n₀≤0 时均不崩溃，
|
||
按合理默认（H²=0、F=1、p=0.5）返回并给 warning。
|
||
"""
|
||
from __future__ import annotations
|
||
|
||
import numpy as np
|
||
|
||
from . import fdist
|
||
|
||
|
||
def solve(rows: list[dict]) -> dict:
|
||
"""rows: [{"group": str|int, "value": float}]。
|
||
|
||
返回：方差分析表字段 + σ²g/σ²e + H² + F/p + 各组均值 + warning。
|
||
"""
|
||
if not rows:
|
||
raise ValueError("无分组表型数据")
|
||
y = np.array([float(r["value"]) for r in rows], dtype=float)
|
||
g = np.array([r["group"] for r in rows])
|
||
if len(y) < 2 or len(np.unique(g)) < 2:
|
||
raise ValueError("至少需要 2 个组且总样本 ≥ 2 才能做方差分析")
|
||
n = len(y)
|
||
overall = float(y.mean())
|
||
|
||
group_keys = sorted(set(g))
|
||
k = len(group_keys)
|
||
sizes: list[int] = []
|
||
means: list[float] = []
|
||
for gv in group_keys:
|
||
mask = g == gv
|
||
sizes.append(int(mask.sum()))
|
||
means.append(float(y[mask].mean()))
|
||
|
||
ss_total = float(((y - overall) ** 2).sum())
|
||
ss_between = float(sum(sz * (m - overall) ** 2 for sz, m in zip(sizes, means)))
|
||
ss_within = ss_total - ss_between
|
||
ms_between = ss_between / (k - 1) if k > 1 else 0.0
|
||
ms_within = ss_within / (n - k) if n > k else 0.0
|
||
n0 = (n - sum(sz * sz for sz in sizes) / n) / (k - 1) if k > 1 else 0.0
|
||
sigma_e = ms_within
|
||
sigma_g = (ms_between - ms_within) / n0 if n0 > 0 else 0.0
|
||
h2 = sigma_g / (sigma_g + sigma_e) if (sigma_g + sigma_e) > 1e-12 else 0.0
|
||
h2 = float(np.clip(h2, 0.0, 1.0))
|
||
f = ms_between / ms_within if ms_within > 0 else 1.0
|
||
p = fdist.f_pvalue(f, k - 1, n - k) if (k > 1 and n > k and ms_within > 0) else 0.5
|
||
|
||
warning = None
|
||
if ms_within == 0:
|
||
warning = "组内方差为 0（各组观测完全一致），方差组分不可分，H² 无信息。"
|
||
elif sigma_g <= 1e-12:
|
||
warning = "家系方差组分 ≈ 0，组间差异不显著，H² 接近 0。"
|
||
return {
|
||
"n_groups": k,
|
||
"n": n,
|
||
"overall_mean": round(overall, 4),
|
||
"ss_total": round(ss_total, 4),
|
||
"ss_between": round(ss_between, 4),
|
||
"ss_within": round(ss_within, 4),
|
||
"ms_between": round(ms_between, 4),
|
||
"ms_within": round(ms_within, 4),
|
||
"sigma_g": round(sigma_g, 4),
|
||
"sigma_e": round(sigma_e, 4),
|
||
"h2_broad": h2,
|
||
"f": round(float(f), 4),
|
||
"p_value": round(p, 6),
|
||
"groups": [
|
||
{"group": str(gv), "n": sz, "mean": round(m, 4)}
|
||
for gv, sz, m in zip(group_keys, sizes, means)
|
||
],
|
||
"warning": warning,
|
||
}
|
||
|
||
|
||
def solve_rcbd(rows: list[dict]) -> dict:
|
||
"""RCBD 双因素 ANOVA：y_ij = μ + a_i(组合/家系) + b_j(区组) + e_ij。
|
||
|
||
区组效应序贯析出（Type I：grand mean → group → block），遗传差异为区组校正后
|
||
（SS_group 先于区组入模型、SS_block 为其后增量、残差为两者剔除后）。
|
||
输出区组行 SS_block/MS_block/F_block/p_block 与区组均值。区组数不足或无残余
|
||
自由度时按数据现状返回 warning，不崩溃。
|
||
|
||
rows: [{"group": str|int, "block": str|int, "value": float}]
|
||
"""
|
||
if not rows:
|
||
raise ValueError("无 RCBD 表型数据")
|
||
y = np.array([float(r["value"]) for r in rows], dtype=float)
|
||
g = np.array([str(r["group"]) for r in rows])
|
||
b = np.array([str(r["block"]) for r in rows])
|
||
n = len(y)
|
||
if n < 2:
|
||
raise ValueError("至少需要 2 个观测")
|
||
g_keys = sorted(set(g))
|
||
b_keys = sorted(set(b))
|
||
kg, kb = len(g_keys), len(b_keys)
|
||
if kg < 2:
|
||
raise ValueError("至少需要 2 个组合")
|
||
if kb < 2:
|
||
raise ValueError("至少需要 2 个区组才能做 RCBD 分析")
|
||
|
||
overall = float(y.mean())
|
||
|
||
def design(*factors: tuple[list, np.ndarray]) -> np.ndarray:
|
||
# 截距 + 每个因子 drop 第一个水平的哑变量 → 满秩序贯设计阵
|
||
cols = [np.ones(n)]
|
||
for keys, arr in factors:
|
||
for kv in keys[1:]:
|
||
cols.append((arr == kv).astype(float))
|
||
return np.column_stack(cols)
|
||
|
||
X1 = np.ones((n, 1))
|
||
Xg = design((g_keys, g))
|
||
Xgb = design((g_keys, g), (b_keys, b))
|
||
|
||
def rss(X: np.ndarray) -> float:
|
||
beta, *_ = np.linalg.lstsq(X, y, rcond=None)
|
||
resid = y - X @ beta
|
||
return float(resid @ resid)
|
||
|
||
ss_total = float(((y - overall) ** 2).sum())
|
||
ss_group = max(rss(X1) - rss(Xg), 0.0) # 组合效应（区组校正前）
|
||
ss_block = max(rss(Xg) - rss(Xgb), 0.0) # 区组效应（组合剔除后析出）
|
||
ss_error = max(rss(Xgb), 0.0) # 残余误差
|
||
df_group = kg - 1
|
||
df_block = kb - 1
|
||
df_error = n - kg - kb + 1
|
||
ms_group = ss_group / df_group if df_group > 0 else 0.0
|
||
ms_block = ss_block / df_block if df_block > 0 else 0.0
|
||
ms_error = ss_error / df_error if df_error > 0 else 0.0
|
||
sizes = [int((g == gv).sum()) for gv in g_keys]
|
||
means = [float(y[g == gv].mean()) for gv in g_keys]
|
||
n0 = (n - sum(sz * sz for sz in sizes) / n) / (kg - 1) if kg > 1 else 0.0
|
||
sigma_e = ms_error
|
||
sigma_g = (ms_group - ms_error) / n0 if n0 > 0 else 0.0
|
||
h2 = sigma_g / (sigma_g + sigma_e) if (sigma_g + sigma_e) > 1e-12 else 0.0
|
||
h2 = float(np.clip(h2, 0.0, 1.0))
|
||
f_group = ms_group / ms_error if ms_error > 0 else 1.0
|
||
f_block = ms_block / ms_error if ms_error > 0 else 1.0
|
||
p_group = fdist.f_pvalue(f_group, df_group, df_error) if (df_error > 0 and ms_error > 0) else 0.5
|
||
p_block = fdist.f_pvalue(f_block, df_block, df_error) if (df_error > 0 and ms_error > 0) else 0.5
|
||
|
||
warning = None
|
||
if df_error <= 0:
|
||
warning = "区组设计过饱和：无残余误差自由度（每组合每区组仅 1 株），无法估计误差项。"
|
||
elif ms_error == 0:
|
||
warning = "残余误差为 0（区组+组合完全解释变异），F 值无信息。"
|
||
elif sigma_g <= 1e-12:
|
||
warning = "家系方差组分 ≈ 0，组间差异不显著，H² 接近 0。"
|
||
return {
|
||
"n_groups": kg,
|
||
"n_blocks": kb,
|
||
"n": n,
|
||
"overall_mean": round(overall, 4),
|
||
"ss_total": round(ss_total, 4),
|
||
"ss_between": round(ss_group, 4),
|
||
"ss_block": round(ss_block, 4),
|
||
"ss_within": round(ss_error, 4),
|
||
"ms_between": round(ms_group, 4),
|
||
"ms_block": round(ms_block, 4),
|
||
"ms_within": round(ms_error, 4),
|
||
"sigma_g": round(sigma_g, 4),
|
||
"sigma_e": round(sigma_e, 4),
|
||
"h2_broad": h2,
|
||
"f": round(float(f_group), 4),
|
||
"p_value": round(p_group, 6),
|
||
"f_block": round(float(f_block), 4),
|
||
"p_block": round(p_block, 6),
|
||
"groups": [
|
||
{"group": str(gv), "n": sz, "mean": round(m, 4)}
|
||
for gv, sz, m in zip(g_keys, sizes, means)
|
||
],
|
||
"blocks": [
|
||
{"block": str(bv), "n": int((b == bv).sum()), "mean": round(float(y[b == bv].mean()), 4)}
|
||
for bv in b_keys
|
||
],
|
||
"warning": warning,
|
||
}
|