分类：基础统计

sommer包拟合直接和间接遗传效应

直接遗传效应（direct genetic effects，DGE）：个体的表现型直接受基因型影响。

间接遗传效应（indirect genetic effects，IGE）：个体的表现除了受本身基因的影响，还间接受到相邻个体（社会伙伴）的基因影响。它属于一种环境影响。间接遗传效应的存在改变了基因型-表型的关系，在非直观的条件下改变进化的进程。

使用CIMMYT小麦项目组提供的数据。数据有2个block，98个个体。该数据用于计算邻接材料对焦点材料的效应。

表型数据：

     block   focal neighbour trait
3878    b2 id_4568   id_1995   506
2038    b2 id_9637   id_9637   390
3994    b2 id_1720   id_4571   251
1642    b2  id_755   id_1771   438
1505    b2 id_6115   id_4568   486
1813    b2 id_1536    id_784   421
1459    b2 id_6115   id_5764   632
3345    b2 id_7489   id_4669   212
1680    b2  id_755   id_4960   569
1293    b2 id_4669   id_3613   652
2063    b2  id_845    id_845   629
2051    b2 id_4960   id_4960   490
1412    b2 id_5229   id_5746   891
1760    b2 id_1995   id_9779   182
3526    b2 id_8509   id_6115   455
2014    b2 id_6180   id_6180   702
1138    b2 id_5963   id_2734   351
3282    b2 id_1436   id_5963   500
2026    b2 id_5042   id_5042   332
1492    b2 id_6115   id_1802   587

基因型数据：

             id_5440      id_4568     id_5890     id_1309     id_4640      id_4007
id_5440  1.956696361 -0.061769158 -0.01634397 -0.06447941 -0.04040527 -0.004130619
id_4568 -0.061769158  2.002057918  0.02756001 -0.06931443  0.01471098  0.008653623
id_5890 -0.016343974  0.027560012  1.94589038  0.10144781  0.02089325  0.006026260
id_1309 -0.064479410 -0.069314435  0.10144781  1.95258925 -0.01528626  0.016297550
id_4640 -0.040405271  0.014710976  0.02089325 -0.01528626  1.92236611 -0.016776517
id_4007 -0.004130619  0.008653623  0.00602626  0.01629755 -0.01677652  1.953613106

（1）直接遗传效应

表现型只由focal决定。

library(sommer)
data(DT_ige)
DT <- DT_ige
## Direct genetic effects model
modDGE <- mmer(trait ~ block,   # 固定效应为block
random = ~ focal,   # 随机效应为focal
rcov = ~ units,
data = DT, verbose=FALSE)
summary(modDGE)$varcomp

                   VarComp VarCompSE    Zratio Constraint
focal.trait-trait 19894.45 3118.3474  6.379806   Positive
units.trait-trait 10134.22  477.9483 21.203584   Positive

（2）间接遗传效应

表现型由focal和neighbour共同决定。这里使用DGE的方式计算focal和neighbour的方差组分。

data(DT_ige)
DT <- DT_ige
## Indirect genetic effects model
modDGE <- mmer(trait ~ block,
random = ~ focal + neighbour,
rcov = ~ units,
data = DT, verbose=FALSE)
summary(modDGE)$varcomp

                        VarComp VarCompSE    Zratio Constraint
focal.trait-trait     20550.511 3148.6833  6.526700   Positive
neighbour.trait-trait  2926.704  607.4191  4.818261   Positive
units.trait-trait      7301.084  363.8236 20.067649   Positive

此外，还可以用IGE的方式计算间接遗传效应。gvs()用来计算基因型的协方差矩阵。

data(DT_ige)
DT <- DT_ige
### Indirect genetic effects model
modIGE <- mmer(trait ~ block,
               random = ~ gvs(focal, neighbour),
               rcov = ~ units,
               data = DT, verbose=FALSE)
summary(modIGE)$varcomp

                                  VarComp VarCompSE    Zratio Constraint
focal:focal.trait-trait         21014.516 3212.3586  6.541772   Positive
focal:neighbour.trait-trait     -7469.401 1246.1105 -5.994173   Unconstr
neighbour:neighbour.trait-trait  2964.707  576.9991  5.138149   Positive
units.trait-trait                7297.715  357.8869 20.391120   Positive

另外，可以通过gvs()的Gu参数提供方差协方差矩阵列表，这里Af=An都是个体的加性效应矩阵。

data(DT_ige)
DT <- DT_ige
Af <- A_ige
An <- A_ige
### Indirect genetic effects model
modIGE <- mmer(trait ~ block,
               random = ~ gvs(focal, neighbour, Gu=list(Af,An)),
               rcov = ~ units,
               data = DT, verbose=FALSE)
summary(modIGE)$varcomp

                                  VarComp VarCompSE    Zratio Constraint
focal:focal.trait-trait         27806.797 4162.7014  6.679988   Positive
focal:neighbour.trait-trait     -9901.351 1532.8048 -6.459630   Unconstr
neighbour:neighbour.trait-trait  3638.534  611.4065  5.951089   Positive
units.trait-trait                7409.998  359.9827 20.584320   Positive

2021年4月16日

sommer包计算双列杂交的遗传力

遗传力又叫遗传率，分为广义遗传力（board heritability）和狭义遗传力（narrow sense heritability）。广义遗传力一般用H²表示，分子是所有的遗传方差，分母是总方差（包括遗传方差、环境方差和误差方差）。相应的，狭义遗传力用h²表示，分子是加性效应的遗传方差，分母是总方差。

（1）完全双列杂交设计

模型为：

y = Xβ + Zu1 + Zu2 + ZuS + ε

这里，X是固定效应的发生率矩阵，β是固定效应向量，Z是随机效应的发生率矩阵，u1和u2分别是杂种优势群GCA1和GCA2的随机效应向量，uS是SCA的随机效应向量，ε是残差向量。GCA和SCA的BLUP可以用于预测杂交。

计算代码为：

library(sommer)   # 加载sommer包，第一次安装请先运行install.packages("sommer")
data(DT_cornhybrids)   # 加载数据集
DT <- DT_cornhybrids   # 杂交数据矩阵赋值给DT变量
DTi <- DTi_cornhybrids   # 亲本信息矩阵赋值给DTi变量，计算遗传力的时候用不到
GT <- GT_cornhybrids   # 数值化的基因型矩阵赋值给GT变量，计算遗传力时用不到
modFD <- mmer(Yield~Location,   # 用地点做固定效应
              random=~GCA1+GCA2+SCA,   # GCA1、GCA2和SCA做随机效应
              rcov=~units,   # 
              data=DT, verbose = FALSE)   # 数据集使用DT，不显示迭代过程
(suma <- summary(modFD)$varcomp)   # 获得方差组分并显示
Vgca <- sum(suma[1:2,1])   # 两个GCA方差的和
Vsca <- suma[3,1]   # SCA方差的
Ve <- suma[4,1]   # 环境方差的
Va = 4*Vgca   # 4个环境×GCA方差 ？？？为什么不直接用GCA而要×4？
Vd = 4*Vsca   # 4个环境×SCA方差 ？？？为什么不直接用SCA？
Vg <- Va + Vd   # 遗传方差=加性方差+显性方差
(H2 <- Vg / (Vg + (Ve)) )   # 广义遗传力=遗传方差/总方差，Ve外面的括号是多余的
(h2 <- Va / (Vg + (Ve)) )   # 狭义遗传力=加性方差/总方差

本例中，使用来自2个杂种优势群的40个自交系，每个杂种优势群有20个自交系。因此，共有20×20=400个可能的组合。该模拟数据在4个地点鉴定了100个系的表型数据，每个地点只有1次重复，共计4×100=400条数据。

下面是杂交数据DT的前100行（1600×6）：

head(GT,100)

    Location GCA1   GCA2         SCA    Yield PlantHeight
1          1 A258 AS5707 A258:AS5707       NA          NA
2          1 A258     B2     A258:B2       NA          NA
3          1 A258    B99    A258:B99       NA          NA
4          1 A258   LH51   A258:LH51       NA          NA
5          1 A258   Mo44   A258:Mo44       NA          NA
6          1 A258  NC320  A258:NC320       NA          NA
7          1 A258  Oh40B  A258:Oh40B 129.0082    1.632108
8          1 A258 W117Ht A258:W117Ht 139.5254    1.613334
9          1 A258 W182BN A258:W182BN 122.8658    1.688060
10         1 A258  W611S  A258:W611S       NA          NA
11         1 A258   A641   A258:A641 166.6094    2.010321
12         1 A258    B10    A258:B10 159.3310    1.947324
13         1 A258   B119   A258:B119       NA          NA
14         1 A258   B14A   A258:B14A       NA          NA
15         1 A258  CM174  A258:CM174       NA          NA
16         1 A258    H91    A258:H91       NA          NA
17         1 A258    N7A    A258:N7A       NA          NA
18         1 A258  PHG86  A258:PHG86       NA          NA
19         1 A258   R226   A258:R226       NA          NA
20         1 A258  W610S  A258:W610S 147.6366    1.564184
21         1 B118 AS5707 B118:AS5707       NA          NA
22         1 B118     B2     B118:B2 126.9928    2.027177
23         1 B118    B99    B118:B99       NA          NA
24         1 B118   LH51   B118:LH51       NA          NA
25         1 B118   Mo44   B118:Mo44 135.6515    1.862859
26         1 B118  NC320  B118:NC320       NA          NA
27         1 B118  Oh40B  B118:Oh40B       NA          NA
28         1 B118 W117Ht B118:W117Ht       NA          NA
29         1 B118 W182BN B118:W182BN       NA          NA
30         1 B118  W611S  B118:W611S 130.1698    1.629009
31         1 B118   A641   B118:A641 156.6514    2.130744
32         1 B118    B10    B118:B10       NA          NA
33         1 B118   B119   B118:B119       NA          NA
34         1 B118   B14A   B118:B14A       NA          NA
35         1 B118  CM174  B118:CM174       NA          NA
36         1 B118    H91    B118:H91       NA          NA
37         1 B118    N7A    B118:N7A 145.3416    1.525327
38         1 B118  PHG86  B118:PHG86       NA          NA
39         1 B118   R226   B118:R226       NA          NA
40         1 B118  W610S  B118:W610S 165.0516    1.603079
41         1  B97 AS5707  B97:AS5707 119.5537    2.031154
42         1  B97     B2      B97:B2       NA          NA
43         1  B97    B99     B97:B99       NA          NA
44         1  B97   LH51    B97:LH51 121.0262    1.728354
45         1  B97   Mo44    B97:Mo44 141.8828    2.007985
46         1  B97  NC320   B97:NC320       NA          NA
47         1  B97  Oh40B   B97:Oh40B 121.4198    1.555613
48         1  B97 W117Ht  B97:W117Ht 141.3195    1.525166
49         1  B97 W182BN  B97:W182BN 129.2736    1.558707
50         1  B97  W611S   B97:W611S       NA          NA
51         1  B97   A641    B97:A641       NA          NA
52         1  B97    B10     B97:B10       NA          NA
53         1  B97   B119    B97:B119       NA          NA
54         1  B97   B14A    B97:B14A 146.8906    1.884599
55         1  B97  CM174   B97:CM174       NA          NA
56         1  B97    H91     B97:H91       NA          NA
57         1  B97    N7A     B97:N7A       NA          NA
58         1  B97  PHG86   B97:PHG86       NA          NA
59         1  B97   R226    B97:R226       NA          NA
60         1  B97  W610S   B97:W610S       NA          NA
61         1 C102 AS5707 C102:AS5707       NA          NA
62         1 C102     B2     C102:B2 134.6043    1.918804
63         1 C102    B99    C102:B99       NA          NA
64         1 C102   LH51   C102:LH51 148.6455    2.105199
65         1 C102   Mo44   C102:Mo44       NA          NA
66         1 C102  NC320  C102:NC320       NA          NA
67         1 C102  Oh40B  C102:Oh40B       NA          NA
68         1 C102 W117Ht C102:W117Ht       NA          NA
69         1 C102 W182BN C102:W182BN       NA          NA
70         1 C102  W611S  C102:W611S       NA          NA
71         1 C102   A641   C102:A641       NA          NA
72         1 C102    B10    C102:B10       NA          NA
73         1 C102   B119   C102:B119       NA          NA
74         1 C102   B14A   C102:B14A       NA          NA
75         1 C102  CM174  C102:CM174       NA          NA
76         1 C102    H91    C102:H91       NA          NA
77         1 C102    N7A    C102:N7A       NA          NA
78         1 C102  PHG86  C102:PHG86       NA          NA
79         1 C102   R226   C102:R226 136.8265    1.684973
80         1 C102  W610S  C102:W610S       NA          NA
81         1 LH61 AS5707 LH61:AS5707       NA          NA
82         1 LH61     B2     LH61:B2       NA          NA
83         1 LH61    B99    LH61:B99       NA          NA
84         1 LH61   LH51   LH61:LH51       NA          NA
85         1 LH61   Mo44   LH61:Mo44       NA          NA
86         1 LH61  NC320  LH61:NC320       NA          NA
87         1 LH61  Oh40B  LH61:Oh40B       NA          NA
88         1 LH61 W117Ht LH61:W117Ht 126.6400    1.607918
89         1 LH61 W182BN LH61:W182BN       NA          NA
90         1 LH61  W611S  LH61:W611S       NA          NA
91         1 LH61   A641   LH61:A641       NA          NA
92         1 LH61    B10    LH61:B10       NA          NA
93         1 LH61   B119   LH61:B119       NA          NA
94         1 LH61   B14A   LH61:B14A       NA          NA
95         1 LH61  CM174  LH61:CM174       NA          NA
96         1 LH61    H91    LH61:H91       NA          NA
97         1 LH61    N7A    LH61:N7A       NA          NA
98         1 LH61  PHG86  LH61:PHG86 138.8772    1.652411
99         1 LH61   R226   LH61:R226       NA          NA
100        1 LH61  W610S  LH61:W610S       NA          NA

下面是亲本信息数据DTi（40×4）：

DTi

   Genotype Group     Yield PlantHeight
1      A258   NSS  88.20679   0.9272210
2      A634   SSS  93.31265   1.0226329
3      A641   SSS 109.00473   1.0095320
4      A680   SSS  98.81964   1.1096407
5    AS5707   NSS  80.65517   1.3279382
6       B10   SSS  95.84411   1.3262020
7      B105   SSS 100.65840   1.1560672
8      B118   NSS 111.36915   1.4949915
9      B119   SSS  86.68451   0.9925757
10     B121   SSS 104.95121   0.6704880
11     B14A   SSS 114.69768   1.0709695
12       B2   NSS 108.69982   1.1506838
13      B84   SSS  92.18648   1.2417499
14      B97   NSS  84.54403   1.2251734
15      B99   NSS  88.14189   1.3141365
16     C102   NSS  99.86213   1.3150211
17    CM174   SSS 107.16554   1.0458131
18    H105W   SSS  94.94959   0.9217209
19      H91   SSS 105.75978   0.8894477
20     LH51   NSS  81.85276   1.2395905
21     LH61   NSS 109.39862   0.9673939
22      LP5   SSS 103.10739   0.9972820
23     Mo44   NSS 108.44876   0.8871609
24      N7A   SSS 103.19176   1.4401722
25    NC262   NSS  99.52945   0.7187014
26    NC320   NSS 104.42666   1.2892791
27    NC344   NSS  96.68368   0.8644697
28    Oh40B   NSS 102.63953   0.9376076
29     Pa91   NSS 124.03104   1.1982550
30    PHG71   SSS  83.70858   1.1174897
31    PHG86   SSS  90.85945   1.2136497
32    PHW17   SSS 104.01616   1.1334628
33     R226   SSS  95.31983   1.3587364
34    Tzi25   SSS  91.65843   1.0268473
35   W117Ht   NSS 116.60487   0.9949743
36    W153R   NSS 113.74265   1.5680403
37   W182BN   NSS  89.71262   1.4102404
38      W23   NSS  93.20116   1.2686134
39    W610S   SSS 106.17348   0.9484871
40    W611S   NSS  92.02427   1.1937574

基因型矩阵前6行和前6列（40×40）：

GT[1:6,1:6]

              A258       A634        A641        A680      AS5707         B10
A258    2.23285528 -0.3504778 -0.04756856 -0.32239362 -0.07776163 -0.01257374
A634   -0.35047780  1.4529169  0.45203869 -0.02293680 -0.43538636  0.19984929
A641   -0.04756856  0.4520387  1.96940221 -0.09896791 -0.28059417  0.02019641
A680   -0.32239362 -0.0229368 -0.09896791  1.65221984 -0.33095920  0.12259252
AS5707 -0.07776163 -0.4353864 -0.28059417 -0.33095920  2.36536453 -0.18705982
B10    -0.01257374  0.1998493  0.02019641  0.12259252 -0.18705982  1.78689265

modFD <- mmer(Yield~Location,
              random=~GCA1+GCA2+SCA,
              rcov=~units,
              data=DT, verbose = FALSE)
(suma <- summary(modFD)$varcomp)

                     VarComp VarCompSE     Zratio Constraint
GCA1.Yield-Yield    0.000000  16.50337  0.0000000   Positive
GCA2.Yield-Yield    7.412226  18.94200  0.3913116   Positive
SCA.Yield-Yield   187.560303  41.59428  4.5092817   Positive
units.Yield-Yield 221.142463  18.14716 12.1860656   Positive

Vgca <- sum(suma[1,1])
Vsca <- suma[2,1]
Ve <- suma[3,1]
Va = 2*Vgca
Vd = 2*Vsca
Vg <- Va + Vd
(H2 <- Vg / (Vg + Ve) )
(h2 <- Va / (Vg + (Ve)) )

> (H2 <- Vg / (Vg + Ve) )
[1] 0.7790856
> (h2 <- Va / (Vg + (Ve)) )
[1] 0.02961832

由于该数据主要模拟显性效应，因此，加性方差非常小，导致了h²很小。

（2）不考虑自交的半双列杂交

本例中有7个亲本，共有7×(7-1)/2=21个杂交组合。实验设计是2个重复的完全随机设计（completely randomized disign，CRD）。双亲的GCA和SCA的方差组分可以被估计，

模型为：

y = Xβ + Zug + Zus + ε

代码为：

data("DT_halfdiallel")   # 加载数据集
DT <- DT_halfdiallel   # 将数据集赋值给DT变量
DT$femalef <- as.factor(DT$female)   # 增加femalef列，作为因子
DT$malef <- as.factor(DT$male)   # 增加malef列，作为因子
DT$genof <- as.factor(DT$geno)   # 增加genof列，作为因子
#### model using overlay
modh <- mmer(sugar~1,   # 斜率为固定效应，预测糖
             random=~vs(overlay(femalef,malef))   # 亲本的设计矩阵+基因型向量为随机效应
             + genof,
             data=DT, verbose = FALSE)   # 数据集是DT，不显示迭代
(suma <- summary(modh)$varcomp)
summary(modh)$varcomp
Vgca <- sum(suma[1,1])
Vsca <- suma[2,1]
Ve <- suma[3,1]
Va = 2*Vgca
Vd = 2*Vsca
Vg <- Va + Vd
(H2 <- Vg / (Vg + Ve) )
(h2 <- Va / (Vg + (Ve)) )

                       VarComp VarCompSE   Zratio Constraint
u:femalef.sugar-sugar 5.507899 3.5741151 1.541052   Positive
genof.sugar-sugar     1.815784 1.3629575 1.332238   Positive
units.sugar-sugar     3.117538 0.9626094 3.238632   Positive

这里，vs()是构造方差协方差矩阵的主函数，里面包含了overlay()辅助函数。overlay()函数的作用是构造综合考虑两个变量的发生率矩阵，对于每个因子（本例中共7个），无论是父本还是母本，只要出现即为1，否则为0。

本例中发生率矩阵如下：

overlay(DT$femalef,DT$malef)

   1 2 3 4 5 6 7
1  1 1 0 0 0 0 0
2  1 1 0 0 0 0 0
3  1 0 1 0 0 0 0
4  1 0 1 0 0 0 0
5  1 0 0 1 0 0 0
6  1 0 0 1 0 0 0
7  1 0 0 0 1 0 0
8  1 0 0 0 1 0 0
9  1 0 0 0 0 1 0
10 1 0 0 0 0 1 0
11 1 0 0 0 0 0 1
12 1 0 0 0 0 0 1
13 0 1 1 0 0 0 0
14 0 1 1 0 0 0 0
15 0 1 0 1 0 0 0
16 0 1 0 1 0 0 0
17 0 1 0 0 1 0 0
18 0 1 0 0 1 0 0
19 0 1 0 0 0 1 0
20 0 1 0 0 0 1 0
21 0 1 0 0 0 0 1
22 0 1 0 0 0 0 1
23 0 0 1 1 0 0 0
24 0 0 1 1 0 0 0
25 0 0 1 0 1 0 0
26 0 0 1 0 1 0 0
27 0 0 1 0 0 1 0
28 0 0 1 0 0 1 0
29 0 0 1 0 0 0 1
30 0 0 1 0 0 0 1
31 0 0 0 1 1 0 0
32 0 0 0 1 1 0 0
33 0 0 0 1 0 1 0
34 0 0 0 1 0 1 0
35 0 0 0 1 0 0 1
36 0 0 0 1 0 0 1
37 0 0 0 0 1 1 0
38 0 0 0 0 1 1 0
39 0 0 0 0 1 0 1
40 0 0 0 0 1 0 1
41 0 0 0 0 0 1 1
42 0 0 0 0 0 1 1
attr(,"variables")
[1] "DT$femalef" "DT$malef"

2021年4月15日

图解双列杂交（diallel design）
双列杂交是玉米生产上经常应用的组配实验设计，用来获取最好的单交组合。
双列杂交设计有4种：
（1）完全双列杂交（full diallel）
```
	A1	A2	A3
A1	A1A1	A1A2	A1A3
A2	A2A1	A2A2	A2A3
A3	A3A1	A3A2	A3A3
```
一组数据两两组合，没有例外的情况就是完全双列杂交。其组合数为n^2，如本例n=3，则共有9个组合。
在玉米育种中，本身和本身杂交相当于自交，所以这种方式在玉米上不常用。
（2）不考虑自交的完全双列杂交（full diallel without parents）
```
	A1	A2	A3
A1		A1A2	A1A3
A2	A2A1		A2A3
A3	A3A1	A3A2	
```
这种情况去掉了自交（对角线）的情况。组合数为n(n-1)，本例子中为3*(3-1)=6。由于正反交的效果几乎一致，因此这种方法也不常见。
（3）半双列杂交（half diallel）
```
	A1	A2	A3
A1	A1A1	A1A2	A1A3
A2		A2A2	A2A3
A3			A3A3
```
由于这种组合的矩阵，左下三角和右上三角是对称的，因此只保留一个。组合数公式为：n(n+1)/2。本例3*(3+1)/2=6。该方式解决了正反交的问题。但是，却没有考虑自交问题，在玉米育种上也不常见。
（4）不考虑自交的半双列杂交（half diallel without parents）
```
	A1	A2	A3
A1		A1A2	A1A3
A2			A2A3
A3			
```
这种方式既考虑到了，重复的组合，又去掉了自交。组合数公式为：n(n-1)/2。本例中为3*(3-1)/2=3。这种组配方案是玉米育种中比较常见的方案。
（5）不同杂种优势群的双列杂交
实际上，还有第5种情况，即考虑到不同杂种优势群的双列杂交。组合公式为：n×m，n和m分别为两个杂种优势群中个体的数量。本例子中，3×3=9种组合。这种情况没有自交，因此也不必考虑。这种组配设计也可以叫做Line by Tester设计，材料少时，两种组配设计一致，但是通常情况下，Tester的数量不会很多。
```
	B1	B2	B3
A1	A1B1	A1B2	A1B3
A2	A2B1	A2B2	A2B3
A3	A3B1	A3B2	A3B3
```
2021年4月15日
sommer包的预测函数
（1）预测的背景
在混合线性模型中，y是n×1的观测值向量，混合线性模型可以写成：
这里τ是固定效应的t×1的向量；X是n×t的设计矩阵，它将观测值与固定效应联系起来；u是q×1的随机效应向量；Z是n×q的矩阵，它联系观测值与随机效应；e是n×1的残差误差向量。W和β分别表示组合的设计矩阵和效应向量，它假设：
G和R是随机效应和残差的协方差矩阵，函数的参数是γ和φ。协方差矩阵的数据应该写作：
var(y)=σ²(ZGZ’+R)
协方差的参数γ和φ 用有限的最大似然或最大似然估计。
已知D、 γ和φ 的线性组合的BLUP，然后有：
这里
是混合模型公式的解：
这些也可以写作：
这里的τbar是τ最佳线性无偏估计（BLUE），ubar是u的最佳线性无偏预测（BLUP），这些用协方差var(βbar-β)=C^-1。考虑到，为了更加清晰，使用能形成置信区间预测误差方差（Prediction error variance，PEV或var(βbar-β)，而不是通常感兴趣的方差估计者var(βbar)。对于线性组合Dβbar，PEV是DC^-1D^T。因为方差参数是未知的，因此，我们用它们的REML估计替换掉未知方差参数，并使用经验值。
sommer包中的预测函数用混合模型公式和线性组合D矩阵构建C^-1矩阵用于预测。对于均值，它使用线性组合的D矩阵乘以需要固定和随机效应的向量Xτbar+Zubar或Dβbar，这里D矩阵是来自X和/或Z矩阵的固定和/或随机效应的线性组合。
注意，由于这里打不出来βbar，因此用βbar代替。τbar和ubar同理。
（2）预测均值
sommer用predict()函数计算mmer函数中固定和随机效应的均值和标准差。使用yatesoats数据集你和固定和随机效应。
```
data(DT_yatesoats)
DT <- DT_yatesoats
m3 <- mmer(fixed=Y ~ V + N + V:N,
           random = ~ B + B:MP,
           rcov=~units,
           data = DT, verbose=FALSE)
summary(m3)$varcomp
```
```
           VarComp VarCompSE   Zratio Constraint
B.Y-Y     214.4477 168.62790 1.271722   Positive
B:MP.Y-Y  106.0508  67.83280 1.563415   Positive
units.Y-Y 177.0883  37.34293 4.742217   Positive
```
现在，该模型可以和分类参数一起使用，获得分类的均值。例如，模型包括了固定公式项”V“用于品种，”N“用于指定氮肥处理，”V:N“表示品种和氮肥互作。分类参数可以用于指定需要的均值。下面的例子，氮肥处理的均值：
```
p0 <- predict.mmer(object=m3, classify = "N")
```
```
iteration    LogLik     wall    cpu(sec)   restrained
    1      -0.745921   16:3:8      0           0
    2      -0.745921   16:3:8      0           0
    3      -0.745921   16:3:8      0           0
    4      -0.745921   16:3:8      0           0
```
```
p0$pvals
```
```
  trait   N predicted.value standard.error
1     Y   0        79.38889       9.006796
2     Y 0.2        98.88889       9.006796
3     Y 0.4       114.22222       9.006796
4     Y 0.6       123.38889       9.006796
```
结束语
记住，sommer使用直接反演算法（direct inversion，DI），对于大数据集来说可能非常慢。该包主要针p>n的情况（随机效应比观测值多）且模型使用密集的协方差结构。例如，适用于协方差结构密集且重复较少的实验。也适用于随机效应较多的基因组问题。对于重复很多的情况，如200个个体，10次重复，有2000条记录时，asreml的速度要高于sommer。
2021年4月1日
R语言计算方差协方差矩阵
方差协方差矩阵，多组变量组成的以方差和协方差数值得到的矩阵。要获得方差-协方差矩阵，需要先知道什么是方差。
方差：是衡量随机变量或一组数据时离散程度的度量，即（样本中的每个值减去均值）平方和/自由度。
例如，有一组数据，用R手动计算其方差。
a <- c(2, 2, 1, 5, 4, 1, 2, 5, 2, 2)
sum((a-mean(a))^2)/(length(a)-1) #手动计算方差
var(a) #函数计算方差
可以对比一下，函数计算的方差与手动计算的方差是否一致。
协方差：是用于衡量两个随机变量的联合变化程度。而方差是协方差的一种特殊情况，即变量与自身的协方差。
如果我们又有一组数据，用R计算一下协方差。
b <- c(2, 2, 1, 4, 4, 3, 5, 1, 3, 5)
sum((a-mean(a))*(b-mean(b)))/(length(a)-1)
var(a,b)
可自行查看结果是否一致。
方差-协方差矩阵，是用不同数据得到的协方差矩阵。
如果要得到a，b的方差-协方差矩阵，直接使用var()就可以。
matrixNoHeader <- cbind(a,b)
var(matrixNoHeader)
协方差矩阵行列数量相等，等于变量的个数。
如果是多个数据
a <- c(2, 2, 1, 5, 4, 1, 2, 5, 2, 2) b <- c(2, 2, 1, 4, 4, 3, 5, 1, 3, 5) c <- c(1, 5, 5, 5, 3, 3, 2, 1, 4, 2) d <- c(4, 2, 5, 1, 5, 5, 4, 3, 5, 1)
matrixNoHeader <- cbind(a,b,c,d) var(matrixNoHeader)
对角线是每个变量自身的方差，其余部分是协方差。方差协方差矩阵就是这么简单。
2021年3月23日
两个杂种优势群双列杂交设计的杂交种数量
两个杂种优势群，利用双列杂交组配设计，理论上有j×k种组合（j是A群材料数，k是B群材料数）。
若无杂种优势群的划分，则组合数为n(n+1)/2，n是群体中的材料总数。
2020年12月31日
中心化和标准化的R语言函数scale()
中心化(Zero-centered或者Mean-subtraction)是一个过程，是将一组数据（通常是R语言数据框中的列）的每个值减去该组数据的均值，是的，平均数为零。
标准化(Standardization或Normalization)是在中心化的基础上，一组数据中的每一个中心化过的值除以改组数据的标准差。标准化的元素，均值为0且方差为1。原始数据不同维度的特征（单位）不一致时，通过标准化消除特征之间的差异。
```
scale(x, center = TRUE, scale = FALSE)   # 中心化
scale(x, center = TRUE, scale = TRUE)   # 标准化
```
2020年12月19日
矩阵的阿达玛乘积（Hadamard product）
两个行、列数目均相等的矩阵，每个对应位置上的元素相乘，得到另一个矩阵，最终的矩阵叫做两个矩阵的阿达玛乘积。该乘法有别与矩阵乘法，因此乘号使用#。
2020年12月19日
逆矩阵（Invers Matrix）
在文献中，模型通常会以矩阵相乘的形式给出。在矩阵乘法中，常常见到矩阵的逆矩阵R^-1。
若有两个矩阵，矩阵A×矩阵B=单位矩阵（左上到右下对角线的元素均为1，其余部分均为0），则A和B叫做互逆矩阵，某矩阵的逆矩阵用该矩阵的-1次方表示。
当一个矩阵A×一个向量x=另一个向量b时，可以用逆矩阵求出前面的向量。因此，文献中矩阵的逆是用来就系数的过程。证明如下：
2020年12月18日
关系矩阵（Incidence Matrix）
关系矩阵表示的是两类对象的关系。假设第一类对象是X，第二类对象是Y。每个X占一行，每个Y占一列。如果X中的某一个项（条目，entry）x与Y中某一个项（条目）有关系，则赋值为1，否则为0。以此类推，所有的X与Y的关系可以形成一个矩阵，则这个矩阵为Incidence Matrix。
2020年12月18日