1. 程式人生 > >DMU-多性狀動物模型-學習筆記4

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比較

兩者方差組分一致.

在這裡插入圖片描述