DMU-單性狀動物模型-母體效應--學習筆記5
單性狀動物模型-母體效應
本次主要是演示如何使用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(第五列)
固定因子: BYEAR和SEX(第三列, 第四列)
隨機因子: ID + MOTHER
所以這裡編寫DIR
第一部分, 是註釋, 這裡所寫的東西會輸出到結果檔案, 基本上就是模型的解釋, 這部分沒有強制要求, 可以省略
$COMMENT
Model
y: BWT TARSUS
fixed: BYEAR + SEX
random: ANIMAL+MOTHER
第二部分是分析方法, 預設是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
第五部分, 有6行, 定義模型
整體來說是:
第一行: 單性狀 # 1
第二行: 1性狀無吸收 # 0
第三行: 1個性狀, 是由3個因子, 兩個固定因子:3,4, 一個隨機因子:1 # 1 0 3 3 4 1
第四行: 2個隨機因子, 他們的關係是獨立, 沒有協方差, 1 2 # 2 1 2
注意, 如果兩個隨機因子之間, 是有協方差的, 那麼寫作 2 1 1, 即表示隨機因子有: 加性+ 母體 + 加性:母體協方差
第五行: 無迴歸項 # 0
第六行: 無約束 # 0
$MODEL
1
0
1 0 4 3 4 1 2
2 1 2
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中, 然後重新命名為: Uni_animalmodel-maternal.DIR
$COMMENT
Estimate variance components for BWT
Model
y: BWT
fixed: BYEAR + SEX
random: ANIMAL + MOTHER
$ANALYSE 1 1 0 0
$DATA ASCII (4,2,0) d:/dmu-test/animal-model.txt
$VARIABLE
ANIMAL MOTHER BYEAR SEX
BWT TARSUS
$MODEL
1
0
1 0 4 3 4 1 2
2 1 2
0
0
$VAR_STR 1 PED 2 ASCII d:/dmu-test/animal-ped.txt
$SOLUTION
$DMUAI
10
1D-7
1D-6
1
執行DIR檔案
這裡執行的run_dmuai.bat, 將DMU安裝路徑下的檔案run_dmuai.bat拷貝到d:/dmu-test資料夾, 在終端cmd介面鍵入:
run_dmuai.bat Uni_animalmodel-maternal
執行結果:
D:\dmu-test>run_dmuai.bat Uni_animalmodel-maternal
D:\dmu-test>Echo OFF
Starting DMU using Uni_animalmodel-maternal.DIR as directive file
D:\dmu-test>
檢視結果
在檔案*lst中有估算的方差組分, 結果如下:
SUMMARY OF MINIMIZATION PROCESS
Eval Criterion !!Delta!! !!Gradient!! Parameters
---- --------- --------- ------------ |------------------------------------------
1 2311.63 0.4281 3.796 | 1.6003 1.1819 1.3955
2 2170.63 0.3174 2.343 | 2.1014 1.1469 1.6188
3 2150.86 0.1004 0.5529 | 2.2655 1.1053 1.6588
4 2150.13 0.7153E-02 0.2589E-01 | 2.2777 1.1040 1.6570
5 2150.13 0.8765E-04 0.8250E-04 | 2.2778 1.1040 1.6569
6 2150.13 0.8343E-06 0.1020E-05 | 2.2778 1.1040 1.6569
7 2150.13 0.5385E-07 0.4235E-07 | 2.2778 1.1040 1.6569
可以看到模型收斂
方差組分為:
Estimated (co)-variance components
----------------------------------
Parameter vector for L at convergence
Asymptotic SE based on AI-information matrix
No Parameter Asymp. S.E.
1 2.27778 0.497101
2 1.10404 0.239802
3 1.65690 0.373448
Asymp. correlation matrix of parameter vector
可以看到:
加性方差組分為: 2.2778
母體效應方差組分為: 1.10404
殘差方差組分為: 1.6569
對比asreml的結果:
程式碼:
library(asreml)
head(dat)
dat[dat$BWT==0,]$BWT=NA
ainv = asreml.Ainverse(ped)$ginv
mod = asreml(BWT ~ BYEAR + SEX, random = ~ ped(ANIMAL) + MOTHER, ginverse = list(ANIMAL=ainv),data=dat)
summary(mod)$varcomp
pin(mod,h2 ~ V2/(V1+V2+V3))
方差組分:
> summary(mod)$varcomp
gamma component std.error z.ratio constraint
MOTHER!MOTHER.var 0.666325 1.104035 0.2398025 4.603936 Positive
ped(ANIMAL)!ped 1.374724 2.277783 0.4971012 4.582131 Positive
R!variance 1.000000 1.656902 0.3734477 4.436772 Positive
遺傳力:
> pin(mod,h2 ~ V2/(V1+V2+V3))
Estimate SE
h2 0.4520558 0.08842455
DMU和asreml比較
兩者方差組分一致.
