DMU-單性狀重複力模型-學習筆記3
阿新 • • 發佈：2018-12-22
單性狀重複力模型
本次主要是演示如何使用DMU分析單性狀重複力模型.
重複力模型和動物模型的區別:
不是所有的性狀都可以分析重複力模型, 首先重複力模型是動物模型的拓展, 它適合一個個體多個觀測值的情況.
- 比如豬的產仔數, 一個母豬有多個胎次
- 比如雞的產蛋, 不同時間段, 雞都有產蛋量
- 牛的產奶量, 不同的測定日, 產奶量不同
- 豬的飼料消耗, 也是重複測量的資料
只有這樣的資料才可以將永久環境效應剖分出來.
重複力是遺傳力的上限:
教科書上這樣說, 這句話怎麼理解呢?
首先, 我們認為
遺傳力為
這裡的Vg如果有重複測量的資料, 可以剖分為可以遺傳的部分, 和不可以遺傳的部分(永久環境效應), 那麼遺傳力的計算公式為:
重複力的計算公式為:
當Vpe為0時, 重複力=遺傳力, 當Vpe>0時, 重複力>遺傳力, 所以說重複力是遺傳力的上限!
資料使用learnasreml包中的資料
learnasreml是我編寫的輔助學習asreml的R包, 裡面有相關的資料和程式碼, 這裡我們用其中的repeatmodel.dat和repeatmodel.ped的資料.
如果沒有軟體包, 首先安裝:
setwd("d:/dmu-test/")
library(devtools)
# install_github("dengfei2013/learnasreml")
library(learnasreml)
data("repeatmodel.dat")
data("repeatmodel.ped")
dat = repeatmodel.dat
ped = repeatmodel.ped
summary(dat)
summary(ped)
看一下資料:
> summary(dat) ANIMAL BYEAR AGE YEAR LAYDATE 1 : 5 1000 : 109 2:308 1004 : 79 Min. : 0.00 3 : 5 1001 : 98 3:322 1005 : 78 1st Qu.:20.00 9 : 5 999 : 86 4:339 1003 : 69 Median :24.00 17 : 5 1002 : 85 5:315 1006 : 64 Mean :23.54 42 : 5 987 : 70 6:323 1002 : 60 3rd Qu.:27.00 50 : 5 989 : 66 988 : 54 Max. :41.00 (Other):1577 (Other):1093 (Other):1203 > summary(ped) ID FATHER MOTHER Min. : 1 Min. : 0.0 Min. : 0.0 1st Qu.: 328 1st Qu.: 0.0 1st Qu.: 135.0 Median : 655 Median : 0.0 Median : 538.0 Mean : 655 Mean : 261.5 Mean : 547.4 3rd Qu.: 982 3rd Qu.: 458.0 3rd Qu.: 932.0 Max. :1309 Max. :1304.0 Max. :1306.0
資料介紹:
- 有因子4個: 分別是 ANIMAL BYEAR AGE YEAR
- 有變數1個: LAYDATE
- 缺失值用0表示
- 系譜中有三列資料, 無出生時間一列, 缺失值為0
需要做的處理
- 系譜增加第四列出生時間, 因為資料都是數字, 沒有字串, 不需要轉化
- 在儲存資料時, 去掉行頭
- 編輯DIR檔案
使用R語言清洗資料, 並儲存資料到D盤dmu-test
dat = repeatmodel.dat
ped = repeatmodel.ped
summary(dat)
summary(ped)
dmuped = ped
dmuped$Birth = 2018
head(dat)
library(data.table)
# write.table(dat,"animal-model.txt",row.names = F,col.names = F)
fwrite(dat,"repeat-model.txt",sep = " ",col.names = F)
fwrite(dmuped,"repeat-ped.txt",sep = " ",col.names = F)
編寫DIR檔案
想要分析的模型:
觀測值: LAYDATE (第四列)
固定因子: 第二列, 第三列, 第四列
隨機因子: ID, ide(ID)
所以這裡編寫DIR
第一部分, 是註釋, 這裡所寫的東西會輸出到結果檔案, 基本上就是模型的解釋, 這部分沒有強制要求, 可以省略
$COMMENT
Estimate variance components for BWT
Model
y: LAYDATE
fixed: BYEAR + AGE + YEAR
random: ANIMAL +ide(ANIMAL)
第二部分是分析方法, 預設是AI
$ANALYSE 1 1 0 0
第三部分是定義因子數和變數書, 以及檔案位置:
$DATA ASCII (4,1,0) d:/dmu-test/repeat-model.txt
第四部分, 定義變數名, 也是為了方便結果輸出, 相當於資料的行頭名
$VARIABLE
ANIMAL BYEAR AGE YEAR
BWT
第五部分, 有6行, 定義模型
整體來說是:
第一行: 單性狀 # 1
第二行: 無吸收 # 0
第三行: 主要定義y變數, 固定因子, 隨機因子
- 分析的是第一個變數 # 1
- 無權重考慮 # 0
- 共五個因子(固定+隨機, 固定寫前面, 隨機寫後面) # 5
- 第一個固定因子是第二列因子 #2 #BYEAR
- 第二個固定一致是第三列因子 #3 #AGE
- 第三個固定因子是第四列 #4 #YEAR
- 第四個隨機因子是第一列 #1 #ANIMAL
- 第五個隨機因子是第一列 #1 #ANIMAL
- 所以, 5個因子, 三個固定因子:2,3,4, 兩個隨機因子:1,1 #1 0 5 2 3 4 1 1
第四行: 有兩個隨機因子, 他們的關係是獨立的, 所以是2 1 2
1
0
1 0 5 2 3 4 1 1
2 1 2
0
0
第六部分: 指定系譜
$VAR_STR 1 PED 2 ASCII d:/dmu-test/repeat-ped.txt
注意, 如果想要輸出BLUP值, 定義:$DMUAI
$DMUAI
10
1D-7
1D-6
1
完整DIR檔案, 放入model.txt中, 然後重新命名為: Uni_repeatmodel.DIR
$COMMENT
Estimate variance components for BWT
Model
y: LAYDATE
fixed: BYEAR + AGE + YEAR
random: ANIMAL +ide(ANIMAL)
$ANALYSE 1 1 0 0
$DATA ASCII (4,1,0) d:/dmu-test/repeat-model.txt
$VARIABLE
ANIMAL BYEAR AGE YEAR
BWT
$MODEL
1
0
1 0 5 2 3 4 1 1
2 1 2
0
0
$VAR_STR 1 PED 2 ASCII d:/dmu-test/repeat-ped.txt
$DMUAI
10
1D-7
1D-6
1
執行DIR檔案
這裡執行的run_dmuai.bat, 將DMU安裝路徑下的檔案run_dmuai.bat拷貝到d:/dmu-test資料夾, 在終端cmd介面鍵入:
run_dmuai.bat Uni_repeatmodel
執行結果:
D:\dmu-test>run_dmuai.bat Uni_repeatmodel
D:\dmu-test>Echo OFF
Starting DMU using Uni_repeatmodel.DIR as directive file
檢視結果
在檔案*lst中有估算的方差組分, 結果如下:
SUMMARY OF MINIMIZATION PROCESS
Eval Criterion !!Delta!! !!Gradient!! Parameters
---- --------- --------- ------------ |------------------------------------------
1 12629.2 0.8574 4.330 | 1.8100 1.8898 1.8705
2 8234.59 1.370 6.822 | 2.9917 3.3812 3.2879
3 6444.28 1.776 8.529 | 4.2397 5.4642 5.1761
4 5857.47 1.566 6.869 | 4.9013 7.4662 6.8832
5 5736.36 0.6798 2.497 | 4.9407 8.3737 7.6324
6 5727.01 0.7387E-01 0.2311 | 4.9325 8.4634 7.7233
7 5726.91 0.1399E-02 0.1596E-02 | 4.9341 8.4621 7.7245
8 5726.91 0.7706E-04 0.5119E-04 | 4.9340 8.4622 7.7245
9 5726.91 0.4564E-05 0.2610E-05 | 4.9340 8.4622 7.7245
10 5726.91 0.2695E-06 0.1558E-06 | 4.9340 8.4622 7.72
可以看到模型收斂
方差組分為:
Estimated (co)-variance components
----------------------------------
Parameter vector for L at convergence
Asymptotic SE based on AI-information matrix
No Parameter Asymp. S.E.
1 4.93404 1.76364
2 8.46217 1.63818
3 7.72445 0.329943
遺傳力需要手動計算, 這裡還沒有找到解決方案.
對比asreml的結果:
程式碼:
library(asreml)
head(dat)
dat[dat$LAYDATE==0,]$LAYDATE=NA
ainv = asreml.Ainverse(ped)$ginv
mod = asreml(LAYDATE ~ BYEAR + AGE + YEAR, random = ~ ped(ANIMAL)+ide(ANIMAL), ginverse = list(ANIMAL=ainv),data=dat)
summary(mod)$varcomp
pin(mod,h2 ~ V1/(V1+V2+V3))
方差組分:
> summary(mod)$varcomp
gamma component std.error z.ratio constraint
ped(ANIMAL)!ped 0.6387559 4.934041 1.7636385 2.797649 Positive
ide(ANIMAL)!id 1.0955038 8.462169 1.6381812 5.165588 Positive
R!variance 1.0000000 7.724454 0.3299432 23.411466 Positive
遺傳力:
> pin(mod,h2 ~ V1/(V1+V2+V3))
Estimate SE
h2 0.233612 0.07907261
DMU和asreml比較
兩者方差組分一致.
