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可以做的高级模型

如样条函数的随机回归模型，如协方差函数模型等，快来免费申请试用吧！