ASReml|随机回归模型的示例代码
1、模拟数据来源
Mrode R A, Mrode R A, Thompson R. Linear models for the prediction of animal breeding values.[M]// Linear Models for the Prediction of Animal Breeding Values. 1996:121-134. Test day model: Chapter 9, Fixed Regression Model & Random Regression Model
2、什么是纵向数据
不同时间或不同年龄得到的重复度量值被称为纵向数据（longitudinal data），用于分析此类数据的模型应该能解释随时间或年龄而变化的平均数以及度量值间的协方差，并且能估计出相应参数。
育种中用于随机回归模型的方法，比如奶牛的测定日模型（Test day milk record）。
这里我们通过使用asreml软件，介绍一下书籍中的两个模型：
•固定回归模型
•随机回归模型
3、生成数据
# read in dataset
y <- c(17.0, 18.6, 24.0, 20.0, 20.0, 15.6, 16.0, 13.0, 8.2, 8.0,
23.0, 21.0, 18.0, 17.0, 16.2, 14.0, 14.2, 13.4, 11.8, 11.4,
10.4, 12.3, 13.2, 11.6, 8.4,
22.8, 22.4, 21.4, 18.8, 18.3, 16.2, 15.0,
22.2, 20.0, 21.0, 23.0, 16.8, 11.0, 13.0, 17.0, 13.0, 12.6)
# DIM: days in mild
DIM <- c(4, 38, 72, 106, 140, 174, 208, 242, 276, 310,
4, 38, 72, 106, 140, 174, 208, 242, 276, 310,
174, 208, 242, 276, 310,
106, 140, 174, 208, 242, 276, 310,
4, 38, 72, 106, 140, 174, 208, 242, 276, 310)
# HTD: herd-test-day
HTD <- c(1,2,3,4,5,6,7,8,9,10,
1,2,3,4,5,6,7,8,9,10,
6,7,8,9,10,
4,5,6,7,8,9,10,
1,2,3,4,5,6,7,8,9,10)
Animal <- c(rep(4, 10), rep(5, 10), rep(6, 5), rep(7, 7), rep(8, 10))
dat <- data.frame(Animal, DIM, HTD, y)
for (i in 1:3) dat[,i] <- as.factor(dat[,i])
head(dat)
str(dat)
# Pedigree
Animal <- 1:8
Sire <- c(0, 0, 0, 1, 3, 1, 3, 1)
Dam <- c(0, 0, 0, 2, 2, 5, 4, 7)
ped <- data.frame(Animal, Sire, Dam)
head(ped)
表型数据：
> dat
Animal DIM HTD y
1 4 4 1 17.0
2 4 38 2 18.6
3 4 72 3 24.0
4 4 106 4 20.0
5 4 140 5 20.0
6 4 174 6 15.6
7 4 208 7 16.0
8 4 242 8 13.0
9 4 276 9 8.2
10 4 310 10 8.0
11 5 4 1 23.0
12 5 38 2 21.0
13 5 72 3 18.0
14 5 106 4 17.0
15 5 140 5 16.2
16 5 174 6 14.0
17 5 208 7 14.2
18 5 242 8 13.4
19 5 276 9 11.8
20 5 310 10 11.4
21 6 174 6 10.4
22 6 208 7 12.3
23 6 242 8 13.2
24 6 276 9 11.6
25 6 310 10 8.4
26 7 106 4 22.8
27 7 140 5 22.4
28 7 174 6 21.4
29 7 208 7 18.8
30 7 242 8 18.3
31 7 276 9 16.2
32 7 310 10 15.0
33 8 4 1 22.2
34 8 38 2 20.0
35 8 72 3 21.0
36 8 106 4 23.0
37 8 140 5 16.8
38 8 174 6 11.0
39 8 208 7 13.0
40 8 242 8 17.0
41 8 276 9 13.0
42 8 310 10 12.6
系谱数据：
> ped
Animal Sire Dam
1 1 0 0
2 2 0 0
3 3 0 0
4 4 1 2
5 5 3 2
6 6 1 5
7 7 3 4
8 8 1 7
4、固定回归模型
由于书中给定了方差组分，所以这里，我们用asreml固定方差组分的方法，进行演示。
# 给定方差组分
Va = 5.521
Vpe = 8.47
Ve = 3.71
init = asreml(y ~ pol(DIM,4) + HTD,
random =~ vm(Animal,ainv) + ide(Animal),
residual =~ idv(units),start.values = T,
maxit=100,
data=dat)
vc = init$vparameters.table
vc[c(1,2,4),2] = c(Va,Vpe,Ve)
vc[c(1,2,4),3] = c(“F”,”F”,”F”)
mod1 = asreml(y ~ pol(DIM,4) + HTD,
random =~ vm(Animal,ainv) + ide(Animal),
residual =~ idv(units),G.param = vc,R.param = vc,
maxit=100,
data=dat)
summary(mod1)$varcomp
coef(mod1)
相关BLUE值和BLUP值：
> coef(mod1)
$fixed
effect
HTD_1 0.00000
HTD_2 0.00000
HTD_3 0.00000
HTD_4 0.00000
HTD_5 0.00000
HTD_6 -10.40514
HTD_7 -40.61782
HTD_8 -117.37303
HTD_9 -272.33896
HTD_10 -536.53441
pol(DIM, 4)_order0 0.00000
pol(DIM, 4)_order1 180.36904
pol(DIM, 4)_order2 159.06823
pol(DIM, 4)_order3 77.47144
pol(DIM, 4)_order4 17.91839
(Intercept) 102.07270
$random
effect
vm(Animal, ainv)_1 -0.30706603
vm(Animal, ainv)_2 -0.19331164
vm(Animal, ainv)_3 0.50037767
vm(Animal, ainv)_4 -0.03487732
vm(Animal, ainv)_5 -0.25509014
vm(Animal, ainv)_6 -0.73248534
vm(Animal, ainv)_7 1.14175101
vm(Animal, ainv)_8 0.34637219
ide(Animal)_1 0.00000000
ide(Animal)_2 0.00000000
ide(Animal)_3 0.00000000
ide(Animal)_4 -0.73389965
ide(Animal)_5 -0.56124921
ide(Animal)_6 -1.38504598
ide(Animal)_7 2.89795191
ide(Animal)_8 -0.2177570
书中的结果，两者基本一致：
5、随机回归模型
由于书中给定了方差组分，所以这里，我们用asreml固定方差组分的方法，进行演示。
Va = c(3.297, 0.594 ,0.921 ,-1.381, -0.289 ,1.005)
Vpe = c(6.872, -0.254 ,3.171, -1.101 ,0.167, 2.457 )
Ve = 3.710
init1 = asreml(y ~ pol(DIM,4) + HTD,
random =~ us(pol(DIM,2)):vm(Animal,ainv) +
us(pol(DIM,2)):ide(Animal),
residual =~ idv(units),start.values = T,
data=dat)
vc1 = init1$vparameters.table
vc1
vc1[1:14,2] = c(Va,Vpe,1,Ve)
vc1[1:14,3] = rep(“F”,14)
vc1
str(dat)
moded <- asreml(y ~ pol(DIM,4) + HTD,
random =~ us(pol(DIM,2)):vm(Animal,ainv) +
us(pol(DIM,2)):ide(Animal),
residual =~ idv(units),
G.param = vc1,
R.param = vc1,
data=dat)
summary(moded)$varcomp
coef(moded)
asreml计算的结果：
> coef(moded)
$fixed
effect
HTD_1 0.000000
HTD_2 0.000000
HTD_3 0.000000
HTD_4 0.000000
HTD_5 0.000000
HTD_6 -9.394095
HTD_7 -36.751883
HTD_8 -107.271590
HTD_9 -250.797553
HTD_10 -496.141083
pol(DIM, 4)_order0 0.000000
pol(DIM, 4)_order1 166.025789
pol(DIM, 4)_order2 146.982790
pol(DIM, 4)_order3 72.210875
pol(DIM, 4)_order4 16.894125
(Intercept) 95.191996
$random
effect
pol(DIM, 2)_order0:vm(Animal, ainv)_1 -0.12782765
pol(DIM, 2)_order0:vm(Animal, ainv)_2 -0.04213793
pol(DIM, 2)_order0:vm(Animal, ainv)_3 0.16996558
pol(DIM, 2)_order0:vm(Animal, ainv)_4 0.36683012
pol(DIM, 2)_order0:vm(Animal, ainv)_5 -0.43003701
pol(DIM, 2)_order0:vm(Animal, ainv)_6 -0.64729099
pol(DIM, 2)_order0:vm(Animal, ainv)_7 0.93231427
pol(DIM, 2)_order0:vm(Animal, ainv)_8 0.19096140
pol(DIM, 2)_order1:vm(Animal, ainv)_1 -0.04252153
pol(DIM, 2)_order1:vm(Animal, ainv)_2 -0.04581112
pol(DIM, 2)_order1:vm(Animal, ainv)_3 0.08833264
pol(DIM, 2)_order1:vm(Animal, ainv)_4 0.02529622
pol(DIM, 2)_order1:vm(Animal, ainv)_5 -0.09401289
pol(DIM, 2)_order1:vm(Animal, ainv)_6 -0.13115909
pol(DIM, 2)_order1:vm(Animal, ainv)_7 0.26042073
pol(DIM, 2)_order1:vm(Animal, ainv)_8 0.05985742
pol(DIM, 2)_order2:vm(Animal, ainv)_1 -0.01042919
pol(DIM, 2)_order2:vm(Animal, ainv)_2 -0.04994381
pol(DIM, 2)_order2:vm(Animal, ainv)_3 0.06037300
pol(DIM, 2)_order2:vm(Animal, ainv)_4 -0.39268200
pol(DIM, 2)_order2:vm(Animal, ainv)_5 0.31776628
pol(DIM, 2)_order2:vm(Animal, ainv)_6 0.31828446
pol(DIM, 2)_order2:vm(Animal, ainv)_7 -0.41833318
pol(DIM, 2)_order2:vm(Animal, ainv)_8 -0.02693079
pol(DIM, 2)_order0:ide(Animal)_1 0.00000000
pol(DIM, 2)_order0:ide(Animal)_2 0.00000000
pol(DIM, 2)_order0:ide(Animal)_3 0.00000000
pol(DIM, 2)_order0:ide(Animal)_4 -0.81848094
pol(DIM, 2)_order0:ide(Animal)_5 -0.39316974
pol(DIM, 2)_order0:ide(Animal)_6 -1.46405660
pol(DIM, 2)_order0:ide(Animal)_7 2.72622327
pol(DIM, 2)_order0:ide(Animal)_8 -0.05051598
pol(DIM, 2)_order1:ide(Animal)_1 0.00000000
pol(DIM, 2)_order1:ide(Animal)_2 0.00000000
pol(DIM, 2)_order1:ide(Animal)_3 0.00000000
pol(DIM, 2)_order1:ide(Animal)_4 -0.63999319
pol(DIM, 2)_order1:ide(Animal)_5 -0.07126467
pol(DIM, 2)_order1:ide(Animal)_6 0.09505531
pol(DIM, 2)_order1:ide(Animal)_7 0.56907161
pol(DIM, 2)_order1:ide(Animal)_8 0.04713094
pol(DIM, 2)_order2:ide(Animal)_1 0.00000000
pol(DIM, 2)_order2:ide(Animal)_2 0.00000000
pol(DIM, 2)_order2:ide(Animal)_3 0.00000000
pol(DIM, 2)_order2:ide(Animal)_4 -1.92020634
pol(DIM, 2)_order2:ide(Animal)_5 1.12195344
pol(DIM, 2)_order2:ide(Animal)_6 0.34869767
pol(DIM, 2)_order2:ide(Animal)_7 -0.61522760
pol(DIM, 2)_order2:ide(Animal)_8 1.06478283
这里，我们将个体4的回归系数进行提取：
•0阶：pol(DIM, 2)_order0:vm(Animal, ainv)_4 0.36683012
•1阶：pol(DIM, 2)_order1:vm(Animal, ainv)_4 0.02529622
•2阶：pol(DIM, 2)_order2:vm(Animal, ainv)_4 -0.39268200
书中的结果为
•0阶：0.3445
•1阶：0.0063
•2阶：-0.3164
结果差别不大。
由此数据，我们可以计算305天的育种值。
6、asreml可以做的高级模型
如样条函数的随机回归模型，如协方差函数模型等，快来免费申请试用吧！