DMU-多性狀動物模型-學習筆記4
多性狀動物模型
本次主要是演示如何使用DMU分析多性狀動物模型.
資料使用learnasreml包中的資料
learnasreml是我編寫的輔助學習asreml的R包, 裡面有相關的資料和程式碼, 這裡我們用其中的animalmodel.dat和animalmodel.ped的資料.
如果沒有軟體包, 首先安裝:
setwd("d:/dmu-test/") library(devtools) # install_github("dengfei2013/learnasreml") library(learnasreml) data("animalmodel.dat") data("animalmodel.ped") dat = animalmodel.dat ped = animalmodel.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,"animal-model.txt",sep = " ",col.names = F) fwrite(dmuped,"animal-ped.txt",sep = " ",col.names = F)
看一下資料:
> summary(dat) ANIMAL MOTHER BYEAR SEX BWT TARSUS 1 : 1 96 : 8 998 : 53 1:470 Min. : 0.000 Min. : 0.00 2 : 1 541 : 8 994 : 47 2:614 1st Qu.: 2.730 1st Qu.: 0.00 3 : 1 581 : 8 983 : 45 Median : 6.385 Median :16.27 5 : 1 584 : 8 987 : 45 Mean : 5.802 Mean :12.93 6 : 1 1302 : 8 991 : 45 3rd Qu.: 8.660 3rd Qu.:21.94 7 : 1 12 : 7 997 : 44 Max. :15.350 Max. :39.66 (Other):1078 (Other):1037 (Other):805 > 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, MOTHER, BYEAR, SEX
有變數2個: 分別是BWT和TARSUS
缺失值為0
系譜中,
有三列資料, 無出生時間一列, 缺失值為0
需要做的處理
- 系譜增加第四列出生時間, 因為資料都是數字, 沒有字串, 不需要轉化
- 在儲存資料時, 去掉行頭
- 編輯DIR檔案
編寫DIR檔案
想要分析的模型:
觀測值: BWT(第五列), TARSUS (第六列)
固定因子: BYEAR和SEX(第三列, 第四列)
隨機因子: ID
所以這裡編寫DIR
第一部分, 是註釋, 這裡所寫的東西會輸出到結果檔案, 基本上就是模型的解釋, 這部分沒有強制要求, 可以省略
$COMMENT
Model
y: BWT TARSUS
fixed: BYEAR + SEX
random: ANIMAL
第二部分是分析方法, 預設是AI
$ANALYSE 1 1 0 0
第三部分是定義因子數和變數數, 以及檔案位置:
$DATA ASCII (4,2,0) d:/dmu-test/animal-model.txt
上面的意思是, 資料是ASCII格式, 有4個固子, 2個變數, 缺失值用0表示, 資料的絕對路徑是: d:/dmu-test/animal-model.txt
第四部分, 定義變數名, 也是為了方便結果輸出, 相當於資料的行頭名
$VARIABLE
ANIMAL MOTHER BYEAR SEX
BWT TARSUS
第五部分, 有6行, 定義模型
整體來說是:
第一行: 兩性狀 # 2
第二行: 1性狀無吸收 # 0
第三行: 2性狀無吸收 # 0
第四行: 1性狀, 是由3個因子, 兩個固定因子:3,4, 一個隨機因子:1 # 1 0 3 3 4 1
第五行: 2性狀, 是由3個因子, 兩個固定因子:3,4, 一個隨機因子:1 # 2 0 3 3 4 1
第六行: 性狀1, 1個隨機因子 # 1 1
第七行: 性狀2, 1個隨機因子 # 1 1
第八行: 性狀1,無迴歸 # 0
第九行: 性狀2,無迴歸 # 0
第十行: 殘差0
$MODEL
2
0
0
1 0 3 3 4 1
2 0 3 3 4 1
1 1
1 1
0
0
0
第六部分: 指定系譜
$VAR_STR 1 PED 2 ASCII d:/dmu-test/animal-ped.txt
注意, 如果想要輸出BLUP值, 定義:$DMUAI
$DMUAI
10
1D-7
1D-6
1
完整DIR檔案, 放入model.txt中, 然後重新命名為: mul-animalmodel.DIR
$COMMENT
Model
y: BWT TARSUS
fixed: BYEAR + SEX
random: ANIMAL
$ANALYSE 1 1 0 0
$DATA ASCII (4,2,0) d:/dmu-test/animal-model.txt
$VARIABLE
ANIMAL MOTHER BYEAR SEX
BWT TARSUS
$MODEL
2
0
0
1 0 3 3 4 1
2 0 3 3 4 1
1 1
1 1
0
0
0
$VAR_STR 1 PED 2 ASCII d:/dmu-test/animal-ped.txt
$DMUAI
10
1D-7
1D-6
1
執行DIR檔案
這裡執行的run_dmuai.bat, 將DMU安裝路徑下的檔案run_dmuai.bat拷貝到d:/dmu-test資料夾, 在終端cmd介面鍵入:
run_dmuai.bat mul_animalmodel
執行結果:
D:\dmu-test>run_dmuai.bat mul_animalmodel
D:\dmu-test>Echo OFF
Starting DMU using mul_animalmodel.DIR as directive file
D:\dmu-test>
檢視結果
在檔案*lst中有估算的方差組分, 結果如下:
Eval Criterion !!Delta!! !!Gradient!! Parameters
---- --------- --------- ------------ |--------------------------------------------------------
1 12028.8 0.6039 4.096 | 1.5877 0.73966E-01 1.8936 1.4327
| 0.12929 1.9136
2 7774.73 0.9673 6.170 | 2.1162 0.31777 3.4204 1.7356
| 0.49187 3.5631
3 5909.74 1.510 8.930 | 2.3621 0.82080 5.6988 1.9228
| 1.1827 6.3310
4 5161.67 1.984 10.91 | 2.4806 1.4486 8.3095 2.1217
| 2.1515 10.257
5 4917.50 1.785 9.047 | 2.5638 1.8830 10.081 2.3066
| 3.0591 14.120
6 4867.84 0.7835 3.603 | 2.5821 1.9975 10.541 2.3927
| 3.4651 15.932
7 4864.20 0.8472E-01 0.3898 | 2.5817 2.0041 10.586 2.4033
| 3.5105 16.129
8 4864.17 0.1682E-02 0.4107E-02 | 2.5819 2.0049 10.590 2.4033
| 3.5102 16.128
9 4864.17 0.4621E-04 0.6606E-04 | 2.5819 2.0049 10.590 2.4032
| 3.5102 16.128
10 4864.17 0.7679E-05 0.1041E-04 | 2.5819 2.0049 10.590 2.4032
| 3.5102 16.128
11 4864.17 0.1192E-05 0.1748E-05 | 2.5819 2.0049 10.590 2.4032
| 3.5102 16.128
12 4864.17 0.1937E-06 0.3123E-06 | 2.5819 2.0049 10.590 2.4032
| 3.5102 16.128
可以看到模型收斂
方差組分為:
Estimated (co)-variance components
----------------------------------
Parameter vector for L at convergence
Asymptotic SE based on AI-information matrix
No Parameter Asymp. S.E.
1 2.58189 0.437110
2 2.00491 0.857216
3 10.5895 2.68116
4 2.40324 0.348455
5 3.51022 0.727723
6 16.1280 2.36436
遺傳力需要手動計算, 這裡還沒有找到解決方案.
對比asreml的結果:
程式碼:
library(asreml)
dat[dat$BWT==0,]$BWT=NA
dat[dat$TARSUS==0,]$TARSUS=NA
ainv = asreml.Ainverse(ped)$ginv
mod2 = asreml(cbind(BWT,TARSUS) ~ trait + trait:(BYEAR + SEX),
random = ~ us(trait):ped(ANIMAL), rcov = ~ units:us(trait),ginverse = list(ANIMAL=ainv),data=dat)
summary(mod2)$varcomp
方差組分:
> summary(mod2)$varcomp
gamma component std.error z.ratio constraint
trait:ped(ANIMAL)!trait.BWT:BWT 2.581883 2.581883 0.4371085 5.906732 Positive
trait:ped(ANIMAL)!trait.TARSUS:BWT 2.004949 2.004949 0.8572152 2.338910 Positive
trait:ped(ANIMAL)!trait.TARSUS:TARSUS 10.589430 10.589430 2.6811944 3.949520 Positive
R!variance 1.000000 1.000000 NA NA Fixed
R!trait.BWT:BWT 2.403246 2.403246 0.3484542 6.896879 Positive
R!trait.TARSUS:BWT 3.510189 3.510189 0.7277219 4.823531 Positive
R!trait.TARSUS:TARSUS 16.128117 16.128117 2.3643446 6.821390 Positive
DMU和asreml比較
兩者方差組分一致.