博文
random regression model之ASReml V4篇
|||
为何使用random regression(RR) model？举个简单的例子，当测量的次数为10时，直接使用US矩阵需要初始值为10*11/2=55个，而用pol(Age,2)来拟合的话，所需初始值为=(2+1)* (2+2)/2=6个，大大降低了参数数量，也使得ASReml运行更容易收敛。
现举个简单的例子，某物种，目标性状为壳长，分4个时间测定。对于不同时间测定时，ASReml分析时，一般不直接使用时间age，而是新增一个参数assess，其值为1，2，。。。，简单理解为测定次数，与时间对应。
对于测定时间或次数不超过4时，一般直接使用US model，反之，则使用RR model，后者可大大降低了模型参数，减少运算量，并使得模型容易收敛。RR model无法直接获得G矩阵，需要通过运行的结果，通过矩阵运算来获得，在ASReml win版本中，可以通过!pin功能来实现。
下面附上US model和RR model的分析代码。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 | !NOGRAPH !workspace 500 !rename !ARGS 6 #1 2 3 4 5 random regression analysis for longitudinal data #animal,dam,sire,family_id,gender,assess,ag3,shell_length #800001,947,408,201102081,F,1,6,26.77 #800001,947,408,201102081,F,2,18,NA #800001,947,408,201102081,F,3,24,NA #800001,947,408,201102081,F,4,30,80.17 animal !P # 800001 dam !I #!p # 947 sire !I #!P # 408 family_id !I # 201102081 gender !A # F assess * # 4 age # 30 shell_length # 80.17 # Check/Correct these field definitions. fam2011_ped.csv !skip 1 !repeat !MAKE !ALPHA !Groups 9 fam2011_trait_rr2.csv !SKIP 1 !maxit 30 !extra 5 !MVINCLUDE !continue !AISING !FCON !dopart $1 !part 1 # us model shell_length !SIGMAPAR ~ assess -assess.gender !r us(assess).animal !f mv residual id(8001).diag(assess) !part 6 # RR model, reduced fit2 shell_length !SIGMAPAR ~ assess -assess.gender !r us(pol(age,1)).animal !f mv residual id(8001).diag(assess) !pin !DEFINE F V_6_6 5 6*-2 7*1 # -1-1 -1*-1 # 8 V6 F V_618 5 6*-1 # -1+0 -1*0 # 9 V618 F V1818 5 # 0+0 0*0 # 10 V18 F V_624 5 6*-0.5 7*-.5 # -1+.5 -1*.5 # 11 V624 F V1824 5 6*0.5 # 0+.5 0*.5 # 12 V1824 F V2424 5 6 7*0.25 # .5+.5 .5*.5 # 13 V24 F V_630 5 7*-1 # -1+1 -1*1 # 14 V630 F V1830 5 6 # 0+1 0*1 # 15 V1830 F V2430 5 6*1.5 7*0.5 # .5+1 .5*1 # 16 V2430 F V3030 5 6*2 7 # 1+1 1*1 # 17 V30 F P6 1 8 # 18 F P18 2 10 # 19 F P24 3 13 # 20 F P30 4 17 # 21 H H2_06 8 18 H H2_18 10 19 H H2_24 13 20 H H2_30 17 21 R R_618 8:10 R R_624 8 11 13 R R_630 8 14 17 R R_1824 10 12 13 R R_1830 10 15 17 R R_2430 13 16 17 |
模型6--RR model的运行结果.asr如下：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 | - - - Results from analysis of shell_length - - - Akaike Information Criterion 93331.35 (assuming 7 parameters). Bayesian Information Criterion 93388.65 Model_Term Sigma Sigma Sigma/SE % C id(8001).diag(assess) 32004 effects assess DIAG_V 1 2.37043 2.37043 11.96 0 P assess DIAG_V 2 10.9314 10.9314 45.60 0 P assess DIAG_V 3 7.39985 7.39985 38.83 0 P assess DIAG_V 4 10.3899 10.3899 35.53 0 P us(pol(age,1)).animal 16430 effects pol(age,1) US_V 1 1 11.0905 11.0905 44.51 0 P pol(age,1) US_C 2 1 5.67169 5.67169 30.13 0 P pol(age,1) US_V 2 2 6.07165 6.07165 26.88 0 P animal NRM 8215 Covariance/Variance/Correlation Matrix US us(pol(age,1)).anima 11.09 0.6912 5.672 6.072 |
模型6--RR model的运行结果.res如下：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 | === === === === Residual statistics for RRmodel6.asr === === === === Notice: ASReml will have merged any design points closer than 0.002400 Use !POLPOINTS 100000 on the data file line to increase the resolution. pol(age,1) has 2 levels 6.00000 1.00000 -1.00000 7.20000 1.00000 -0.90000 8.40000 1.00000 -0.80000 9.60000 1.00000 -0.70000 10.80000 1.00000 -0.60000 12.00000 1.00000 -0.50000 13.20000 1.00000 -0.40000 14.40000 1.00000 -0.30000 15.60000 1.00000 -0.20000 16.80000 1.00000 -0.10000 18.00000 1.00000 -0.00000 19.20000 1.00000 0.10000 20.40000 1.00000 0.20000 21.60000 1.00000 0.30000 22.80000 1.00000 0.40000 24.00000 1.00000 0.50000 25.20000 1.00000 0.60000 26.40000 1.00000 0.70000 27.60000 1.00000 0.80000 28.80000 1.00000 0.90000 30.00000 1.00000 1.00000 |
结合.asr中的随机回归系数方差协方差矩阵与 .res中的随机回归系数设计矩阵来计算G矩阵。然后通过!pin来实现计算。
结果在.pvc中，具体如下：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 | ASReml 4.0 [01 Jan 2013] random regression analysis for lon RRModel6.pvc created 22 Jul 2014 17:11:55.609 - - - Results from analysis of shell_length - - - id(8001).diag(assess) 32004 effects 1 assess 2.37043 0.198196 2 assess 10.9314 0.239724 3 assess 7.39985 0.190570 4 assess 10.3899 0.292426 us(pol(age,1)).animal 16430 effects 5 pol(age 11.0905 0.249169 6 pol(age 5.67169 0.188241 7 pol(age 6.07165 0.225880 animal NRM 8215 8 V_6_6 5 5.8188 0.38711 9 V_618 5 5.4188 0.22224 10 V1818 5 11.091 0.24916 11 V_624 5 5.2188 0.21212 12 V1824 5 13.926 0.30823 13 V2424 5 18.280 0.40702 14 V_630 5 5.0189 0.26674 15 V1830 5 16.762 0.38163 16 V2430 5 22.634 0.51650 17 V3030 5 28.506 0.66620 18 P6 1 8.1892 0.23169 19 P18 2 22.022 0.34331 20 P24 3 25.680 0.43052 21 P30 4 38.895 0.65974 H2_6 = V_6_6 5 8/P6 1 18= 0.7105 0.0300 H2_18 = V1818 5 10/P18 2 19= 0.5036 0.0079 H2_24 = V2424 5 13/P24 3 20= 0.7118 0.0073 H2_30 = V3030 5 17/P30 4 21= 0.7329 0.0080 R_618 2 1 = V_618 9/SQR[V_6_6 8*V1818 10]= 0.6745 0.0128 R_624 = V_624 5/SQR[V_6_6 5*V2424 5]= 0.5060 0.0203 R_630 = V_630 5/SQR[V_6_6 5*V3030 5]= 0.3897 0.0255 R_1824 = V1824 5/SQR[V1818 5*V2424 5]= 0.9781 0.0018 R_1830 = V1830 5/SQR[V1818 5*V3030 5]= 0.9427 0.0048 R_2430 = V2430 5/SQR[V2424 5*V3030 5]= 0.9915 0.0007 Notice: The parameter estimates are followed by their approximate standard errors. |
RR model计算G矩阵，在ASReml win版本中，有点不方便。如果用R软件来计算，只需将构建随机回归系数设计矩阵Z和随机回归方差协方差矩阵Q，既可通过G=Z*Q*Z‘来获得G矩阵，不过要获得标准误SE，有点麻烦。
https://blog.sciencenet.cn/blog-1114360-830095.html
上一篇：ASReml-R与协变量模型分析
下一篇：AAfun2.0版本已完成
