個體育種價估計值使用以下八種統計方法計算：RR-BLUP、BRR、BL、
Bayes A、Bayes B、Bayes C、RKHS、RF 等，各方法均有適合之 R 套件：其中 RR-BLUP 利用 rrBLUP version 4.4，貝氏方法與 RKHS 利用 BGLR version 1.0.5 (Perez and de los Campos 2014)，RF 則利用 randomForest version 4.6-12 (Liaw and Wiener 2002)。
2.5.1 RR-BLUP
GS 預測模型的標準線性模型如式(5)：
𝑦𝑖 = 𝜇 + ∑ 𝑥𝑖𝑗𝛽𝑗+ 𝜀𝑖 (5) 𝑦𝑖為品系 i 的外表型值，𝜇為外表型平均值，𝑥𝑖𝑗為第 i 品系第 j 個分子標誌的基因 型值。𝛽𝑗為第 j 個分子標誌的效應，設為從𝑁(0, 𝜎𝛽2)中抽出的隨機型效應，其中分 子標誌效應變方𝜎𝛽2在 RR-BLUP 的假設下等於𝑉G/𝑁M (𝑉G為遺傳變方；𝑁M為分子 標誌數量)。𝜀𝑖為殘差，取自𝑁(0, 𝜎𝜀2)。𝑥𝑖𝑗以 1、0、-1 構成基因型矩陣，代表雙倍 體基因型 AA、AB、BB (Desta and Ortiz 2014)。利用 Endelman (2011) 發表之 rrBLUP 套件中的 mixed.solve 函數解混合模型估計𝛽𝑗，其先以最大概度法 (maximum likelihood, ML) 或預設的受制最大概度法 (restricted maximum
doi:10.6342/NTU201701806
12
likelihood, REML) 估計兩個變方成分𝜎𝛽2和𝜎𝜀2，再利用頻譜分解 (spectral decomposition) 算法修改𝛽𝑗的常態最小平方估計式：
𝛽̂ = (𝐗′𝐗 + λ𝐈)−1𝐗′𝐲 (6) 𝐗為基因型矩陣，𝐈為單位矩陣，𝐲是外表型向量，脊參數λ = 𝜎𝜀2/𝜎𝛽2，由於分母 較分子小，使得分子標誌效應大幅度壓縮向 0。與一般最小平方法相較多了λ𝐈這 項，使得𝐗′𝐗非奇異 (nonsingular) 反矩陣存在，並減少預測變數間的共線性 (Lorenz et al. 2011)。
2.5.2 BRR
BRR 模型在式(5)的基礎下，另設分子標誌效應變方𝜎𝛽2來自一個具有尺度參 數的反卡方事前分布，效應事前分布為常態分布，分子標誌效應與超參數 (hyper-parameter) 的聯合分布 (join distribution) 如式(7)：
𝑝(𝛽, 𝜎𝛽2) = {∏ 𝑁(𝛽𝑗|0, 𝜎𝛽2)} −2(𝜎𝛽2|𝑑𝑓𝛽, 𝑆𝛽) (7) 其中超參數𝑑𝑓𝛽與𝑆𝛽的值在 BGLR 套件中可自行設定或使用套件內設的規則給出 預設值，本研究參數設定依照預設，自由度為 5，而尺度參數須符合變方分割參 數 R2 表示分子標誌效應變方可解釋外表型變方的比例，預設為 0.5，即預測變數 與殘差各解釋 50%應變數變異 (Perez and de los Campos 2014)。貝氏方法無法以 統計方式取得參數估計值，故而發展出馬可夫鏈蒙地卡羅法 (Markov chain Monte Carlo, MCMC)，此套件藉由吉布斯抽樣法 (Gibbs sampling) 從事後分布中重複抽 樣取得近似參數值的分布，再計算出該分布相應的摘要統計值 (summary
statistics)，估計分子標誌效應 (Lorenz et al. 2011)。參考 Mohammadi et al. (2015) 發表的 PopVar 套件設定，在交叉驗證中抽樣次數 (iteration) 及生樣本數 (burn-in) 分別設為 1500 和 500 以提高計算效率；在基因體選種正式建模時則分別設為 12000 和 3000。此法計算成本高相當耗時，(Perez and de los Campos 2014) 發表 的套件結合 R 內建函數與另外編寫的 C 和 Fortran 程序，以改善電腦運算效率。
doi:10.6342/NTU201701806
13
2.5.3 Bayesian LASSO
BL 模型的分子標誌效應之條件事前分布為常態分布如式(8)：
𝑝(𝛽|𝜏2, 𝜎𝜀2) = ∏ 𝑁(𝛽𝑗|0, 𝜏𝑗2𝜎𝜀2) (8) 𝜎𝜀2取自反卡方分布，𝜏𝑗2為分子標誌變方參數，具分子標誌專一性，使模型針對不 同大小效應估計值壓縮，其來自指數事前分布，以式(9)表示：
𝑝(𝜏2|λ) = ∏ 𝐸𝑥𝑝 (𝜏𝑗2|λ2
2) (9)
在 BGLR 套件中，正則化參數 (regularization parameter) λ 有三種事前分布可選 擇，分別是固定值、λ2~𝐺𝑎𝑚𝑚𝑎(𝑟, 𝑠)或maxλ ~𝐵𝑒𝑡𝑎(𝑝0, 𝜋0)，本研究使用預設的 γ 分布 (gamma distribution)，s 為 1.1，與 BRR 同樣以內定規則求尺度參數解 (Perez and de los Campos 2014; Perez et al. 2010)。綜合以上假設，最後分子標誌效 應的邊際分布為雙指數分布，根據效應估計值大小有不同程度的壓縮，效應越大 壓縮幅度越小 (Perez et al. 2010)。同樣以吉布斯抽樣法估計分子標誌效應，設定 同 BRR 所述。
2.5.4 BayesA、BayesB 與 BayesC
BayesA 為 Meueissen et al. (2001) 提出以貝氏分析放寬 RR-BLUP 遺傳效應均 勻分布在基因體中的假設之方法，每個分子標誌效應 j 取自常態分布𝑁 (0, 𝜎𝛽2𝑗)，
使每個分子標誌以不同程度壓縮至 0 (Lorenz et al. 2011)。其變方參數抽樣自反卡
方分布，超參數之尺度參數來自γ 分布，最後分子標誌效應邊際分布為尺度化 t
分布，聯合分布以式(10)表示 (Perez and de los Campos 2014)：
𝑝(𝛽, 𝜎𝛽2, 𝑆𝛽) = {∏ 𝑁 (𝛽𝑗|0, 𝜎𝛽2𝑗) χ−2(𝜎𝛽2𝑗|𝑑𝑓𝛽, 𝑆𝛽)} 𝐺(𝑆𝛽|𝑟, 𝑠) (10) BayesB 在 BayesA 的基礎下另外設定分子標誌完全沒有效應的機率，其來自β 分布 (beta distribution)，聯合分布如式(11)(Perez and de los Campos 2014)：
doi:10.6342/NTU201701806
14
𝑝(𝛽, 𝜎𝛽2, 𝜋) = {∏[𝜋𝑁(𝛽𝑗|0, 𝜎𝛽2)
+ (1 − 𝜋)1(𝛽𝑗 = 0)]χ−2(𝜎𝛽2𝑗|𝑑𝑓𝛽, 𝑆𝛽)} 𝐵(𝜋|𝑝0, 𝜋0)
× 𝐺(𝑆𝛽|𝑟, 𝑠)
(11)
BayesC 修改自 BayesB，假定全部存在之分子標誌效應事前變方相等，故效應邊 際分布與 BRR 相同為常態分布，聯合分布如式(12)(Lorenz et al. 2011; Perez and de los Campos 2014)：
𝑝(𝛽, 𝜎𝛽2, 𝜋) = {∏[𝜋𝑁(𝛽𝑗|0, 𝜎𝛽2)
+ (1 − 𝜋)1(𝛽𝑗 = 0)]} × χ−2(𝜎𝛽2|𝑑𝑓𝛽, 𝑆𝛽)𝐵(𝜋|𝑝0, 𝜋0)
(12)
本研究使用 BGLR 套件計算，參數皆使用預設值，自由度𝑑𝑓𝛽為 5，γ 分布的形狀 參數 (shape parameter) s 為 1.1，在預設 R2 為 0.5 下以內建算法計算尺度參數 𝑆𝛽，再求得速率參數 (rate parameter) r。β 分布之分子標誌具非 0 效應的事前機率 𝜋0預設為 0.5，事前計數 (prior counts)量𝑝0為 10 (Perez and de los Campos 2014)。
同樣以吉布斯抽樣法估計分子標誌效應，設定同 BRR 所述。
2.5.5 RKHS
RKHS 的模型為非線性迴歸模型：
{𝑦𝑖 = 𝜇 + 𝑢𝑖+ 𝜀𝑖 with
𝑝(𝜇, 𝑢, 𝜀) ∝ 𝑁(𝑢|0, 𝑲𝜎𝑢2)𝑁(𝜀|0, 𝑰𝜎𝜀2) (13) 核函數 K 將分子標誌資料轉換成個體間的平方歐基里德距離(squared-Euclidean distance)，式子為：
𝑲(𝑥𝑖, 𝑥𝑖′) = 𝑒𝑥𝑝 {−ℎ ×∑𝑝𝑗=1(𝑥𝑖𝑗 − 𝑥𝑖′𝑗)2
𝑝 } (14)
其中𝑥𝑖代表品系 i 基因型，𝑥𝑖𝑗為第 i 品系第 j 個分子標誌的基因型值，共有 p 個分 子標誌。h 為帶寬參數 (bandwidth parameter) 或稱為平滑參數 (smoothing
parameter)，影響輸入空間中距離與特徵空間中距離的關係，在此設為 0.5，表示
doi:10.6342/NTU201701806
15
兩者間相似；若為 2.5 則相關性小 (Lorenz et al. 2011; Perez and de los Campos 2014)。另外特徵空間中的個體效應變方事前分布為反卡方分布𝜎𝑢2~χ−2(𝑑𝑓, 𝑆)，
同樣以預設的自由度為 5 及分子標誌效應變方解釋外表型變方 50%之條件求得尺 度參數 S，並以吉布斯抽樣法估計效應，設定同 BRR 所述 (Perez and de los Campos 2014)。
2.5.6 Random Forests
RF 是一種整體學習方法，並非單一的迴歸模型，而是利用自助重抽總合法 (bootstrap aggregating, bagging) 建立新的訓練資料，再隨機抽取一部份分子標誌 資料來建造決策樹，本研究使用 Liaw 和 Wiener (2002) 發表的 randomForest 套 件。參數設定依照預設，節點大小最小值為 5，即至少由 5 個變數組成；每個節 點以包含𝑚𝑡𝑟𝑦個變數的隨機預測變數子集合來運算確立分支的模型，預設𝑚𝑡𝑟𝑦 =
𝑝
3，p 為分子標誌數目；設定建造的決策樹數目𝑛𝑡𝑟𝑒𝑒為 1000 (Liaw and Wiener 2002)。