该代码主要用到剪切板和不同软件的灵活运用。
SNPs <- read.table("clipboard",header=F)
refere <- read.table("clipboard",header=T)
res <- refere[which(rownames(refere) %in% SNPs[,1]),]
write.table(res,"clipboard",sep="\t",quote=F,col.names=F,row.names=F)
S2_44435558
S3_151630215
S3_144885149
S4_181062027
S7_141142136
S9_76535431
S9_115952576
S10_137204131
2. 运行代码的第一行。
3. 用EmEditor打开基因型文件（HMP格式），选择制表符分隔，复制前两列。
4. 运行第三行到最后的代码。
5. 在需要的地方粘贴皆可。
模型：
y = Xβ + Zu + ε代码：
# sommer的GS预测
library(sommer) # 加载sommer包，第一次安装请先运行install.packages("sommer")
data(DT_wheat)
DT <- DT_wheat
GT <- GT_wheat
colnames(DT) <- paste0("X",1:ncol(DT)) # 列名为X1~X4
DT <- as.data.frame(DT) # 转换为数据框
DT$id <- as.factor(rownames(DT)) # 把行号作为列并作为因子
# select environment 1
rownames(GT) <- rownames(DT) # 给基因型赋值相同的行号（名称）
K <- A.mat(GT) # additive relationship matrix
colnames(K) <- rownames(K) <- rownames(DT) # K的行列相同，都赋值名称
# GBLUP pedigree-based approach
set.seed(12345)
y.trn <- DT
vv <- sample(rownames(DT),round(nrow(DT)/5)) # 抽取20%作为验证群体（被预测）
y.trn[vv,"X1"] <- NA
head(y.trn)
## GBLUP
ans <- mmer(X1~1,
random=~vs(id,Gu=K),
rcov=~units,
data=y.trn, verbose = FALSE) # kinship based
ans$U$`u:id`$X1 <- as.data.frame(ans$U$`u:id`$X1)
rownames(ans$U$`u:id`$X1) <- gsub("id","",rownames(ans$U$`u:id`$X1))
cor(ans$U$`u:id`$X1[vv,],DT[vv,"X1"], use="complete")
## rrBLUP
ans2 <- mmer(X1~1,
random=~vs(list(GT)),
rcov=~units,
data=y.trn, verbose = FALSE) # kinship based
u <- GT %*% as.matrix(ans2$U$`u:GT`$X1) # BLUPs for individuals
rownames(u) <- rownames(GT)
cor(u[vv,],DT[vv,"X1"]) # same correlation
数据使用的是CIMMYT的全球小麦项目组提供的数据，共有599个个体。y.trn表示籽粒产量，每一列表示一个环境，共有4个环境，例子中只计算了80%预测20%的结果。
杂交种的预测使用GCA效应和SCA效应。数据来自1254年测定的玉米单交试验。123个dent和86个flint，用测定的1254材料预测未测定的9324个材料。例子中只计算了第一个性状GY。
模型：
y = Xβ + Zu1 + Zu2 + ZuS + ε代码：
library(sommer)
data(DT_technow)
DT <- DT_technow # 表型数据
Md <- Md_technow # dent马齿基因型矩阵
Mf <- Mf_technow # flint硬粒基因型矩阵
Ad <- Ad_technow # dent的加性效应矩阵
Af <- Af_technow # flint的加性效应矩阵
y.trn <- DT
vv1 <- which(!is.na(DT$GY)) # 获得有数据的材料行号
vv2 <- sample(vv1, 100) # 在有数据的材料行号中随机选择100个材料
y.trn[vv2,"GY"] <- NA # 将抽取的100个材料设为NA
anss2 <- mmer(GY~1, # 固定效应
random=~vs(dent,Gu=Ad) + vs(flint,Gu=Af), # 随机效应，并指定方差协方差矩阵
rcov=~units,
data=y.trn, verbose = FALSE)
summary(anss2)$varcomp
zu1 <- model.matrix(~dent-1,y.trn) %*% anss2$U$`u:dent`$GY
zu2 <- model.matrix(~flint-1,y.trn) %*% anss2$U$`u:flint`$GY
u <- zu1+zu2+anss2$Beta[1,"Estimate"]
cor(u[vv2,], DT$GY[vv2])
VarComp VarCompSE Zratio Constraint
u:dent.GY-GY 16.06423 2.5737578 6.241548 Positive
u:flint.GY-GY 11.42070 2.1591718 5.289390 Positive
units.GY-GY 16.81801 0.7689509 21.871368 Positive
> cor(u[vv2,], DT$GY[vv2])
[1] 0.922328预测时，所有有数据的个体作为建模群体，这里随机减去100个个体的测定值用于验证，即建模群体是有数据的个体-100。zu1是模型中的Zu1，Z是model.matrix(~dent-1,y.trn)，即dent的设计矩阵；u1是anss2$U$u:dent$GY是GY性状dent的BLUP值。相应的，zu2是flint的设计矩阵×GY性状flint的BLUP。
anns2$Beta[1,”Estimate”]是固定效应的BLUP。
最后，100个个体的测定值和这100个个体的BLUP值的相关系数作为预测结果。
上面的例子只使用了GCA进行预测，因为SCA在这里面没有对预测精度起到多大的影响，相反，由于使用了很大的矩阵，导致计算时间明显增加。
直接遗传效应（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表现型只由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表现型由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遗传力又叫遗传率，分为广义遗传力（board heritability）和狭义遗传力（narrow sense heritability）。广义遗传力一般用H2表示，分子是所有的遗传方差，分母是总方差（包括遗传方差、环境方差和误差方差）。相应的，狭义遗传力用h2表示，分子是加性效应的遗传方差，分母是总方差。
模型为：
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由于该数据主要模拟显性效应，因此，加性方差非常小，导致了h2很小。
本例中有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" （1）预测的背景
在混合线性模型中，y是n×1的观测值向量，混合线性模型可以写成：
这里τ是固定效应的t×1的向量；X是n×t的设计矩阵，它将观测值与固定效应联系起来；u是q×1的随机效应向量；Z是n×q的矩阵，它联系观测值与随机效应；e是n×1的残差误差向量。W和β分别表示组合的设计矩阵和效应向量，它假设：
G和R是随机效应和残差的协方差矩阵，函数的参数是γ和φ。协方差矩阵的数据应该写作：
var(y)=σ2(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-1DT。因为方差参数是未知的，因此，我们用它们的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。
方差协方差矩阵，多组变量组成的以方差和协方差数值得到的矩阵。要获得方差-协方差矩阵，需要先知道什么是方差。
方差：是衡量随机变量或一组数据时离散程度的度量，即（样本中的每个值减去均值）平方和/自由度。
例如，有一组数据，用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)
对角线是每个变量自身的方差，其余部分是协方差。方差协方差矩阵就是这么简单。
sommer是Solving Mixed Model Equations in R的首字母缩写，译为在R中解混合线性方程。
sommer的诞生是为了给R语言用户提供一个高效、稳定的多变量混合模型求解工具。该包主要关注p>n（估计的效应多于观测值），其核心算法由C++编写。该包可以拟合包含随机效应的方差-协方差结构的混合模型，指定非齐次方差（一个变量X对另一个变量Y的影响可能因个体而异），还能获得固定/随机效应的BLUP、BLUE、残差、拟合值、方差等参数。
包的核心是解混合模型，利用一个接口调用the NR Direct-Inversion Newton-Raphson or Average Information algorithms（直接反演牛顿-拉夫逊或平均信息算法）。支持多变量模型。多变量或拓展的单变量形式为：
这里，yi是性状的表型向量；βi是固定效应向量；ui是个体的随机效应向量；ei是第i个性状的残差（i=1…t）。假设随机效应ui和ei服从正态分布且均值为0。X和Z分别是固定效应和随机效应的发生率矩阵。多元变量的响应分布和表型的方差协方差（V）是：
这里，K是第K个随机效应（u=1…k）的关系矩阵或协方差矩阵；H=I是残差项的单位矩阵或部分单位矩阵；σgi2和σei2分别表示第i个性状的遗传（或第k个任意随机项）和残差的方差；σgij和σeij分别是遗传（或第k个任意随机项）和残差的协方差（i=1…t, j=1…t）。算法用对数似然数优化。
这里||是矩阵的行列式，用牛顿优化算法更新REML的估计。
这里，θ是性状间随机效应方差分量和协方差分量的向量。H-1是第k个循环二阶导数Hessian矩阵的逆。dL/dσi2是关于方差协方差组分的似然一阶倒数向量。Newton-Raphson算法包含了Lee和Van Der Werf (2016)提出的关系矩阵的特征分析，提高了时间效率。方差分量线性组合的流行的标准差估计函数（例如遗传力和遗传相关）也被纳入该包。
线性模型的优势在于灵活的指定方差、协方差结构。一般来说，混合模型的方差结构可以看作是多重方差协方差结构的张量（kronecker）积。例如，一个多响应模型（例如2个性状），包含g个个体（例如100个）在e个处理（例如3个环境），个体作为随机效应的方差协方差可以表示为下列乘法模型：
T⨂G⨂A
T是性状间个体的协方差结构。
G是环境中个体的协方差结构。
A是一个方形矩阵，表示个体水平之间的协方差（任何已知的关系矩阵）。
上述T和G矩阵需要估计，而A矩阵已知。T和G叫做非结构化（unstructured, US）矩阵。
参考，未必是这个意思（结构化矩阵（二维表）是大小固定，取值的类型固定的矩阵，实际指的是（可能和原定义不完全符合）每一列（变量）的取值是固定类型，如年龄，只能取整数数字。非结构化矩阵，指每一列（性状）可能取不确定的类型，即值的类型不同，如长度可能取小数，病害可能取级数，开花期可能取时间等。）
协方差结构除了上述例子外，还有其他形式：
对角线（Diagonal，DIAG）协方差结构：实际除了方差外，其他都是0。
复合对称性（Compound symmetry，CS）协方差结构：貌似加入了互作。
一阶自回归（First order autoregressive，AR1）协方差结构：
前面提到过的非结构化的方差结构：
sommer可以拟合上述协方差结构。
（1）单变量齐次方差模型
表型数据的形式为：
这类模型指的是单个响应变量模型，其中兴趣变量（即基因型）需要分析与第二个随机效应（即环境）的相互作用，但假设基因型跨环境具有相同的方差分量。这就是所谓的复合对称（CS）模型。
library(sommer)
data("DT_example")
DT <- DT_example
ans1 <- mmer(Yield~Env,
random=~Name+Env:Name,
rcov=~units,
data=DT,
verbose=FALSE)
summary(ans1)
============================================================
Multivariate Linear Mixed Model fit by REML
********************** sommer 4.1 **********************
============================================================
logLik AIC BIC Method Converge
Value -20.14538 46.29075 55.95182 NR TRUE
============================================================
Variance-Covariance components:
VarComp VarCompSE Zratio Constraint
Name.Yield-Yield 3.682 1.691 2.177 Positive
Env:Name.Yield-Yield 5.173 1.495 3.460 Positive
units.Yield-Yield 4.366 0.647 6.748 Positive
============================================================
Fixed effects:
Trait Effect Estimate Std.Error t.value
1 Yield (Intercept) 16.496 0.6855 24.065
2 Yield EnvCA.2012 -5.777 0.7558 -7.643
3 Yield EnvCA.2013 -6.380 0.7960 -8.015
============================================================
Groups and observations:
Yield
Name 41
Env:Name 123
============================================================
Use the '$' sign to access results and parametersmmer函数是sommer中拟合单变量或多变量混合线性模型的函数。
logLik（Log-likelihood），说明书中没有写，应该是似然函数取对数的最大值，即最大似然对数函数。
AIC（Akaike information criterion）， 赤池信息量标准。是评估统计模型的复杂度和衡量统计模型“拟合”资料之优良性（英语：Goodness of Fit）的一种标准，是由日本统计学家赤池弘次创立和发展的。赤池信息量准则建立在信息熵的概念基础上。
一般而言，模型越复杂，模型能解释的点的误差就越小，当然，这是在拟合建模群体时的情况。但是，用另一个测试集来测试，此时，越复杂的点可能误差越大，因为在拟合的时候，模型会把真实情况当作真是情况调整模型，而另外一组真实值，其误差还是误差，此时，复杂的模型很难拟合的很好。从下图看，似然值越小，说明拟合越好。在train的角度，似然值是一直下降的。在test的角度上，有下降、上升的趋势。在AIC角度上也有类似的趋势，因此，AIC作为选择模型的指标，最小值可以得到最优的模型。AIC=2k-2ln(L)，这里k是参数的个数，L是似然数。解读过来就是参数不能太少，也不能太多，找到一个平衡点。
BIC（Bayesian information criterion），贝叶斯信息量准则。是在有限集合中进行模型选择的准则：BIC最低的模型是最好的。该准则部分基于似然函数并与赤池信息量准则（AIC）紧密相关。
BIC=2k×ln(n)-2ln(L)，n是样本数。除了考虑模型参数外，还考虑的样本数。
拟合模型时，增加参数可提高似然，但如此下去可能导致过拟合。BIC与AIC都致力于向模型中引入关于参数数量的惩罚项；其中，BIC中的惩罚项会大于AIC中的惩罚项。
Method，方差估计的方法或算法，包括NR（Direct-inversion Newton-Raphson）和AI（Average Information）。
units是随机误差项。
zratio貌似是z检验比例？未确定。
（2）单变量非齐次方差模型
在多环境实验中，假设跨地点的遗传方差和残差齐次不太合理。因此，指定一个一般的遗传组分和地点特异性遗传方差是一种解决途径。这需要CS+DIAG模型（非齐次CS模型）。
library(sommer)
data("DT_example")
DT <- DT_example
ans2 <- mmer(Yield~Env,
random=~Name+vs(ds(Env),Name),
rcov=~ vs(ds(Env),units),
data=DT,verbose=FALSE)
============================================================
Multivariate Linear Mixed Model fit by REML
********************** sommer 4.1 **********************
============================================================
logLik AIC BIC Method Converge
Value -15.42983 36.85965 46.52072 NR TRUE
============================================================
Variance-Covariance components:
VarComp VarCompSE Zratio Constraint
Name.Yield-Yield 2.963 1.496 1.980 Positive
CA.2011:Name.Yield-Yield 10.146 4.507 2.251 Positive
CA.2012:Name.Yield-Yield 1.878 1.870 1.004 Positive
CA.2013:Name.Yield-Yield 6.629 2.503 2.649 Positive
CA.2011:units.Yield-Yield 4.942 1.525 3.242 Positive
CA.2012:units.Yield-Yield 5.725 1.312 4.363 Positive
CA.2013:units.Yield-Yield 2.560 0.640 4.000 Positive
============================================================
Fixed effects:
Trait Effect Estimate Std.Error t.value
1 Yield (Intercept) 16.508 0.8268 19.965
2 Yield EnvCA.2012 -5.817 0.8575 -6.783
3 Yield EnvCA.2013 -6.412 0.9356 -6.854
============================================================
Groups and observations:
Yield
Name 41
CA.2011:Name 41
CA.2012:Name 41
CA.2013:Name 41
============================================================
Use the '$' sign to access results and parameters如你所见，at和diag函数都能用于表示每个环境的基因型有不同的方差。残差也一样。不同的是，at可以指定用来指定特别的水平或不同方差的环境。
（3）非结构化模型
除了上面的 CS+DIAG模型假设外，还假设一个特定的随机效应（如环境）与另一个随机效应（如基因型）有协方差结构。可以用us()函数获得。us是非结构化（unstructured model）的缩写。
library(sommer)
data("DT_example")
DT <- DT_example
ans3 <- mmer(Yield~Env,
random=~vs(us(Env),Name),
rcov=~ vs(us(Env),units),
data=DT,verbose=FALSE)
summary(ans3)
===================================================================
Multivariate Linear Mixed Model fit by REML
************************** sommer 4.1 **************************
===================================================================
logLik AIC BIC Method Converge
Value -11.49971 28.99943 38.66049 NR TRUE
===================================================================
Variance-Covariance components:
VarComp VarCompSE Zratio Constraint
CA.2011:Name.Yield-Yield 15.665 5.4207 2.890 Positive
CA.2012:CA.2011:Name.Yield-Yield 6.110 2.4851 2.459 Unconstr
CA.2012:Name.Yield-Yield 4.530 1.8208 2.488 Positive
CA.2013:CA.2011:Name.Yield-Yield 6.384 3.0659 2.082 Unconstr
CA.2013:CA.2012:Name.Yield-Yield 0.393 1.5234 0.258 Unconstr
CA.2013:Name.Yield-Yield 8.597 2.4838 3.461 Positive
CA.2011:units.Yield-Yield 4.970 1.5323 3.243 Positive
CA.2012:CA.2011:units.Yield-Yield 4.087 0.0000 Inf Unconstr
CA.2012:units.Yield-Yield 5.673 1.3008 4.361 Positive
CA.2013:CA.2011:units.Yield-Yield 4.087 0.0000 Inf Unconstr
CA.2013:CA.2012:units.Yield-Yield 4.087 0.0000 Inf Unconstr
CA.2013:units.Yield-Yield 2.557 0.6393 4.000 Positive
===================================================================
Fixed effects:
Trait Effect Estimate Std.Error t.value
1 Yield (Intercept) 16.331 0.8137 20.070
2 Yield EnvCA.2012 -5.696 0.7404 -7.693
3 Yield EnvCA.2013 -6.271 0.8191 -7.656
===================================================================
Groups and observations:
Yield
CA.2011:Name 41
CA.2012:CA.2011:Name 82
CA.2012:Name 41
CA.2013:CA.2011:Name 82
CA.2013:CA.2012:Name 82
CA.2013:Name 41
===================================================================
Use the '$' sign to access results and parameters正如我们所见，us(Env)表明，基因型(name)在环境（Env）中有协方差结构。
（4）多变量齐次方差模型
多响应变量模型由某些变量隐藏的相关性决定，从预测角度看有作用。在sommer中，指定多变量模型，响应变量需要在响应中使用cbind()函数；而模型的随机部分使用us(trait)，diag(trait)或at(trait)。
library(sommer)
data("DT_example")
DT <- DT_example
DT$EnvName <- paste(DT$Env,DT$Name)
ans4 <- mmer(cbind(Yield,Weight)~Env,
random=~vs(Name,Gtc=unsm(2))+vs(EnvName,Gtc=unsm(2)),
rcov=~ vs(units,Gtc=unsm(2)),
data=DT,verbose=FALSE)
============================================================
Multivariate Linear Mixed Model fit by REML
********************** sommer 4.1 **********************
============================================================
logLik AIC BIC Method Converge
Value 167.0252 -322.0505 -298.5695 NR TRUE
============================================================
Variance-Covariance components:
VarComp VarCompSE Zratio Constraint
u:Name.Yield-Yield 3.7089 1.68117 2.206 Positive
u:Name.Yield-Weight 0.9071 0.37944 2.391 Unconstr
u:Name.Weight-Weight 0.2243 0.08775 2.557 Positive
u:EnvName.Yield-Yield 5.0921 1.47879 3.443 Positive
u:EnvName.Yield-Weight 1.0269 0.30767 3.338 Unconstr
u:EnvName.Weight-Weight 0.2101 0.06661 3.154 Positive
u:units.Yield-Yield 4.3837 0.64941 6.750 Positive
u:units.Yield-Weight 0.9077 0.14145 6.417 Unconstr
u:units.Weight-Weight 0.2280 0.03377 6.751 Positive
============================================================
Fixed effects:
Trait Effect Estimate Std.Error t.value
1 Yield (Intercept) 16.4093 0.6783 24.191
2 Weight (Intercept) 0.9806 0.1497 6.550
3 Yield EnvCA.2012 -5.6844 0.7474 -7.606
4 Weight EnvCA.2012 -1.1846 0.1593 -7.439
5 Yield EnvCA.2013 -6.2952 0.7850 -8.019
6 Weight EnvCA.2013 -1.3559 0.1681 -8.065
============================================================
Groups and observations:
Yield Weight
u:Name 41 41
u:EnvName 94 94
============================================================
Use the '$' sign to access results and parameters可能注意到了，我们在随机效应后面增加了us(trait)，这表明应该在多变量模型中假设结构。随机效应（如name）后面的diag(trait)表明性状模型（Yield and Weight）没有协方差组分，不应该被估计，而us(trait)假设随机效应有协方差组分（如随机效应Name的产量和重量的协方差）。残差部分相同。
（5）多变量非齐次方差模型
该模型拓展了单变量非齐次方差模型。这是一个CS+DIAG多变量模型。
library(sommer)
data("DT_example")
DT <- DT_example
DT$EnvName <- paste(DT$Env,DT$Name)
ans5 <- mmer(cbind(Yield,Weight)~Env,
random=~vs(Name,Gtc=unsm(2))+vs(ds(Env),Name,Gtc=unsm(2)),
rcov=~ vs(ds(Env),units,Gtc=unsm(2)),
data=DT,verbose=FALSE)
summary(ans5)
=============================================================
Multivariate Linear Mixed Model fit by REML
********************** sommer 4.1 **********************
=============================================================
logLik AIC BIC Method Converge
Value 177.8154 -343.6308 -320.1497 NR TRUE
=============================================================
Variance-Covariance components:
VarComp VarCompSE Zratio Constraint
u:Name.Yield-Yield 3.31936 1.45269 2.2850 Positive
u:Name.Yield-Weight 0.79393 0.32621 2.4338 Unconstr
u:Name.Weight-Weight 0.19085 0.07503 2.5438 Positive
CA.2011:Name.Yield-Yield 8.70657 4.01470 2.1687 Positive
CA.2011:Name.Yield-Weight 1.77892 0.83926 2.1196 Unconstr
CA.2011:Name.Weight-Weight 0.35966 0.17903 2.0089 Positive
CA.2012:Name.Yield-Yield 2.57109 1.94951 1.3188 Positive
CA.2012:Name.Yield-Weight 0.33245 0.39840 0.8345 Unconstr
CA.2012:Name.Weight-Weight 0.03842 0.08595 0.4470 Positive
CA.2013:Name.Yield-Yield 5.46908 2.16307 2.5284 Positive
CA.2013:Name.Yield-Weight 1.34713 0.50479 2.6687 Unconstr
CA.2013:Name.Weight-Weight 0.32902 0.12208 2.6952 Positive
CA.2011:units.Yield-Yield 4.93852 1.52318 3.2422 Positive
CA.2011:units.Yield-Weight 0.99447 0.32150 3.0932 Unconstr
CA.2011:units.Weight-Weight 0.23982 0.07394 3.2433 Positive
CA.2012:units.Yield-Yield 5.73887 1.31533 4.3631 Positive
CA.2012:units.Yield-Weight 1.28009 0.30157 4.2448 Unconstr
CA.2012:units.Weight-Weight 0.31806 0.07286 4.3652 Positive
CA.2013:units.Yield-Yield 2.56127 0.63993 4.0024 Positive
CA.2013:units.Yield-Weight 0.44569 0.12645 3.5246 Unconstr
CA.2013:units.Weight-Weight 0.12232 0.03057 4.0009 Positive
=============================================================
Fixed effects:
Trait Effect Estimate Std.Error t.value
1 Yield (Intercept) 16.4243 0.7891 20.815
2 Weight (Intercept) 0.9866 0.1683 5.863
3 Yield EnvCA.2012 -5.7339 0.8266 -6.937
4 Weight EnvCA.2012 -1.1998 0.1698 -7.066
5 Yield EnvCA.2013 -6.3128 0.8757 -7.209
6 Weight EnvCA.2013 -1.3621 0.1915 -7.114
=============================================================
Groups and observations:
Yield Weight
u:Name 41 41
CA.2011:Name 41 41
CA.2012:Name 41 41
CA.2013:Name 41 41
=============================================================
Use the '$' sign to access results and parameters（6）多变量非结构化方差模型
模型是单变量非结构化方差模型的拓展，是一个us多变量模型。
library(sommer)
data("DT_example")
DT <- DT_example
DT$EnvName <- paste(DT$Env,DT$Name)
ans6 <- mmer(cbind(Yield,Weight)~Env,
random=~vs(us(Env),Name,Gtc=unsm(2)),
rcov=~ vs(ds(Env),units,Gtc=unsm(2)),
data=DT,verbose=FALSE)
====================================================================
Multivariate Linear Mixed Model fit by REML
************************** sommer 4.1 **************************
====================================================================
logLik AIC BIC Method Converge
Value 181.7945 -351.5889 -328.1079 NR TRUE
====================================================================
Variance-Covariance components:
VarComp VarCompSE Zratio Constraint
CA.2011:Name.Yield-Yield 15.6451 5.35692 2.921 Positive
CA.2011:Name.Yield-Weight 3.3586 1.14633 2.930 Unconstr
CA.2011:Name.Weight-Weight 0.7182 0.24871 2.888 Positive
CA.2012:CA.2011:Name.Yield-Yield 6.5289 2.48615 2.626 Positive
CA.2012:CA.2011:Name.Yield-Weight 1.3505 0.52388 2.578 Unconstr
CA.2012:CA.2011:Name.Weight-Weight 0.2842 0.11259 2.524 Positive
CA.2012:Name.Yield-Yield 4.7893 1.86183 2.572 Positive
CA.2012:Name.Yield-Weight 0.8640 0.38377 2.251 Unconstr
CA.2012:Name.Weight-Weight 0.1693 0.08354 2.027 Positive
CA.2013:CA.2011:Name.Yield-Yield 5.9934 2.93830 2.040 Positive
CA.2013:CA.2011:Name.Yield-Weight 1.4232 0.64973 2.190 Unconstr
CA.2013:CA.2011:Name.Weight-Weight 0.3379 0.14680 2.302 Positive
CA.2013:CA.2012:Name.Yield-Yield 2.0987 1.44034 1.457 Positive
CA.2013:CA.2012:Name.Yield-Weight 0.5240 0.32356 1.619 Unconstr
CA.2013:CA.2012:Name.Weight-Weight 0.1342 0.07572 1.772 Positive
CA.2013:Name.Yield-Yield 8.6257 2.47811 3.481 Positive
CA.2013:Name.Yield-Weight 2.1048 0.58748 3.583 Unconstr
CA.2013:Name.Weight-Weight 0.5125 0.14285 3.588 Positive
CA.2011:units.Yield-Yield 4.9516 1.52694 3.243 Positive
CA.2011:units.Yield-Weight 0.9993 0.32286 3.095 Unconstr
CA.2011:units.Weight-Weight 0.2411 0.07432 3.244 Positive
CA.2012:units.Yield-Yield 5.7790 1.32423 4.364 Positive
CA.2012:units.Yield-Weight 1.2914 0.30408 4.247 Unconstr
CA.2012:units.Weight-Weight 0.3212 0.07356 4.366 Positive
CA.2013:units.Yield-Yield 2.5567 0.63883 4.002 Positive
CA.2013:units.Yield-Weight 0.4452 0.12631 3.524 Unconstr
CA.2013:units.Weight-Weight 0.1223 0.03056 4.001 Positive
====================================================================
Fixed effects:
Trait Effect Estimate Std.Error t.value
1 Yield (Intercept) 16.3342 0.8254 19.790
2 Weight (Intercept) 0.9677 0.1770 5.466
3 Yield EnvCA.2012 -5.6637 0.7449 -7.604
4 Weight EnvCA.2012 -1.1855 0.1604 -7.390
5 Yield EnvCA.2013 -6.2153 0.8340 -7.453
6 Weight EnvCA.2013 -1.3406 0.1806 -7.425
====================================================================
Groups and observations:
Yield Weight
CA.2011:Name 41 41
CA.2012:CA.2011:Name 82 82
CA.2012:Name 41 41
CA.2013:CA.2011:Name 82 82
CA.2013:CA.2012:Name 82 82
CA.2013:Name 41 41
====================================================================
Use the '$' sign to access results and parameters任何随机效应都可以指定不同的结构。
（7）方差模型特殊函数的细节
方差模型主函数vs()和辅助函数
sommer的函数vs()允许构建复杂的方差模型并传递给mmer()函数。vs()函数的格式如下：
random=~vs(..., Gu, Gti, Gtc)
这表示，vs()函数反映了矩阵中每个随机效应可能有特别的方差结构。
var(u)=T⊗E⊗…⊗A
这里，…表示vs()函数的参数，用来指定随机效应方差矩阵的keronecker products乘积。辅助函数ds()、us()、cs()、at()可用于定义这些方差结构。因为，随机变量x（如个体）的方差模型可能需要更加灵活，且参数不应只有一个：
random=~x~被读作预测，~x表示用x预测。
例如，如果个体在不同的时间点或环境下测试，我们可以假设在不同的环境-时间点组合中，个体间有不同的方差和协方差组分。例如：
var(u)=T⊗E⊗S⊗A
可以在vs()函数中指定：
random=~vs(us(e),us(s),x,Gu=A,Gtc=T)
这里，e是数据框中环境的列向量，s是数据框中时间点的列向量，x是数据框中个体的向量，A是个体间已知的方差协方差矩阵（通常是一个单位矩阵；默认不指定），T是协方差矩阵且行和列的数量等于性状的数量。
为随机效应构建协方差模型的辅助函数是+ds()对角线协方差结构+us()非结构化协方差+at()是指定水平的非齐次方差结构+cs()自定义协方差结构。
ds()指定对角线（DIAG）协方差结构
对角线协方差结构类似下面的形式：
考虑一个随机效应的例子，g为个体，在e个环境中测定，模型则为：
random=~vs(ds(e),g)
其意义是用向量g与对角线方差矩阵e的张量积做预测。
us()指定非结构化协方差
非结构化协方差表示如下：
考虑相同的例子，个体在环境中的模型为：
random=~vs(us(e),g)
at()用来指定特殊水平的非齐次方差
第二随机效应（方差和协方差）的个对角线协方差结构是这样个样子的：
还是之前的例子，g是个体，有ABC三个处理（环境），则模型可能是：
random=~vs(at(e,c("A","B")),g)
这里方差组分只为g拟合水平A和B。
cs()指定水平的方差协方差结构
第二个随机变量的自定义协方差结构可能形式如下：
依然用之前的例子，g是个体，e是三个处理（环境）A、B、C，模型可能是：
random=~vs(cs(e,mm),g)
这里mm表明为g估计的方差和协方差组分。
（8）方差组分的约束规则（Gtc参数）
sommer的一大优势是在多性状框架下能够非常灵活的指定方差协方差结构。sommer 3.7版后可以非常容易的通过vs()和Gtc参数实现。Gtc通过以下填充数字，用于随机效应方差协方差组分的约束矩阵的预测。
0：不顾及参数
1：有约束的估计
2：无约束的估计
3：不估计，但是在Gti中提供固定值
unsm()用来快速指定非结构化的约束矩阵，fixm()用于固定值约束，fcm()用于固定效应约束。
考虑到4个性状的（y1，y2，y3，y4）多性状模型，一个随机效应（u）和1个固定效应x
fixed=cbind(y1,y2,y3,y4)~x
random=~vs(u,Gtc=?)
4×4的方差协方差组分估计的约束可以被估计为：
（a）非结构化（方差组分必须是正值，协方差可正可负）
random=~vs(u,Gtc=unsm(4))
unsm(4)
[,1] [,2] [,3] [,4]
[1,] 1 2 2 2
[2,] 2 1 2 2
[3,] 2 2 1 2
[4,] 2 2 2 1（b）非协方差（任意组分方差或协方差可以正或负）
random=~vs(u,Gtc=uncm(4))
uncm(4)
[,1] [,2] [,3] [,4]
[1,] 2 2 2 2
[2,] 2 2 2 2
[3,] 2 2 2 2
[4,] 2 2 2 2（c）固定值（方差或协方差组分表明3被作为固定值）
random=~vs(u,Gtc=fixm(4),Gti=mm)
fixm(4)
[,1] [,2] [,3] [,4]
[1,] 3 3 3 3
[2,] 3 3 3 3
[3,] 3 3 3 3
[4,] 3 3 3 3这里mm是4×4的矩阵，该矩阵具有方差分量的初始值。
（d）考虑固定效应
fixed=cbind(y1,y2,y3,y4)~vs(x,Gtc=fcm(c(1,0,1,0)))
fcm(c(1,0,1,0))
[,1] [,2]
[1,] 1 0
[2,] 0 0
[3,] 0 1
[4,] 0 0这里1和0表示性状，固定效应1被估计0不估计。
（9）特殊模型的特殊函数
随机回归模型
拟合随机回归模型，用户可以使用leg()函数拟合Legendre polynomials多项式。这可以用其他协方差结构，如ds()和us()等。
library(orthopolynom)
data(DT_legendre)
DT <- DT_legendre
mRR2<-mmer(Y~ 1 + Xf,
random=~ vs(us(leg(X,1)),SUBJECT),
rcov=~vs(units),
data=DT, verbose = FALSE)
summary(mRR2)$varcomp
VarComp VarCompSE Zratio Constraint
leg0:SUBJECT.Y-Y 2.5782969 0.6717242 3.838326 Positive
leg1:leg0:SUBJECT.Y-Y 0.4765431 0.2394975 1.989763 Unconstr
leg1:SUBJECT.Y-Y 0.3497299 0.2183229 1.601893 Positive
u:units.Y-Y 2.6912226 0.3825197 7.035513 Positive这里，协方差X用于解释主体的轨迹，结合了一个非结构化的协方差矩阵。细节可以查看理论。
GWAS模型
尽管基因组相关分析可以考虑通过各种方法，混合模型仍然是最流行的，我们采用这个方法。最流行和经典的两个方法，第一个通过标记混合建模，获得标记效应，提供p值-log10。第二个假设所有标记的遗传方差组分相似，因此，方差组分只被估计一次，在一般线性模型中使用表型方差矩阵的倒数（V-inverse）测试所有标记，b=(X’V-X)-XV-y。此时的GWAS更快更高效。sommer提供直接的拓展，GWAS函数可以拟合以上两种方法。库中已经有很多类似的方法。sommer中获得标记效应的一般线性模型的形式是：
b=(X’V-X)X’V-y
X=ZMi
这里，b是标记效应（维度1×mt），y是响应变量（单或多变量）（维度1×nt），V是表型方差矩阵的逆（维度nt×nt），Z是随机效应选择的关系矩阵，用来执行GWAS（维度nt×ut），Mi是标记矩阵（M参数）第i列（维度u×m）。
t是性状数量，n是观测值数量，m是标记数量，u是随机效应水平数量。如果P3D=TRUE，计算一次V矩阵并用于所有的标记测试；FALSE，每个标记使用REML估计。
这里我们展示一个简单的GWAS模型，单个性状：
data(DT_cpdata)
DT <- DT_cpdata
GT <- GT_cpdata
MP <- MP_cpdata
#### create the variance-covariance matrix
A <- A.mat(GT) # additive relationship matrix
#### look at the data and fit the model
head(DT,3)
id Row Col Year color Yield FruitAver Firmness Rowf Colf
P003 P003 3 1 2014 0.10075269 154.67 41.93 588.917 3 1
P004 P004 4 1 2014 0.13891940 186.77 58.79 640.031 4 1
P005 P005 5 1 2014 0.08681502 80.21 48.16 671.523 5 1
head(MP,3)
Locus Position Chrom
1 scaffold_77830_839 0 1
2 scaffold_39187_895 0 1
3 scaffold_50439_2379 0 1
GT[1:3,1:4]
scaffold_50439_2381 scaffold_39344_153 uneak_3436043 uneak_2632033
P003 0 0 0 1
P004 0 0 0 1
P005 0 -1 0 1
mix1 <- GWAS(color~1,
random=~vs(id,Gu=A)
+ Rowf + Colf,
rcov=~units,
data=DT,
M=GT, gTerm = "u:id",
verbose = FALSE)
## Performing GWAS evaluation
ms <- as.data.frame(mix1$scores)
ms$Locus <- rownames(ms)
MP2 <- merge(MP,ms,by="Locus",all.x = TRUE);
manhattan(MP2, pch=20,cex=.5, PVCN = "color")
需要注意，标记矩阵M必须补缺（不允许缺失值）。确保M矩阵中的行数等于g项，如id有300个个体，M矩阵为有300×m，m是标记数量。
叠加模型[overlay()函数]
overlay()叠加不同随机效应的矩阵，估计叠加项的单个方差组分。
data("DT_halfdiallel")
DT <- DT_halfdiallel
head(DT)
rep geno male female sugar
1 1 12 1 2 13.950509
2 2 12 1 2 9.756918
3 1 13 1 3 13.906355
4 2 13 1 3 9.119455
5 1 14 1 4 5.174483
6 2 14 1 4 8.452221
DT$femalef <- as.factor(DT$female)
DT$malef <- as.factor(DT$male)
DT$genof <- as.factor(DT$geno)
#### model using overlay
modh <- mmer(sugar~1,
random=~vs(overlay(femalef,malef))
+ genof,
data=DT,verbose = FALSE)
这里，famalef和malef随机效应被叠加，变成了单个随机效应，具有相同的方差组分。
空间模型（使用二维条线）
在田间实验设计中，使用CPdata展示3维条线的使用，用于适应田间设计的空间效应。早期的各种试验，可用的种子很少，因此，需要使用非重复的设计。实验设计例如增强设计和部分重复设计（p-rep）现在非常流行。
适应空间趋势，田间空间协方差矩阵被提出（即自回归残差；arl）。不幸的是，这些协方差矩阵使得模型非常不稳定。最近，其他研究组提出了使用2维空间去克服上述问题，空间项建模更强大。
下面的例子，我们假设没有重复的群体，row和range信息可用，允许拟合2维条线模型。
data(DT_cpdata)
DT <- DT_cpdata
GT <- GT_cpdata
MP <- MP_cpdata
### mimic two fields
A <- A.mat(GT)
mix <- mmer(Yield~1,
random=~vs(id, Gu=A) +
vs(Rowf) +
vs(Colf) +
vs(spl2D(Row,Col)),
rcov=~vs(units),
data=DT,verbose = FALSE)
summary(mix)
============================================================
Multivariate Linear Mixed Model fit by REML
********************** sommer 4.1 **********************
============================================================
logLik AIC BIC Method Converge
Value -151.2011 304.4021 308.2938 NR TRUE
============================================================
Variance-Covariance components:
VarComp VarCompSE Zratio Constraint
u:id.Yield-Yield 783.4 319.3 2.4536 Positive
u:Rowf.Yield-Yield 814.7 390.5 2.0863 Positive
u:Colf.Yield-Yield 182.2 129.7 1.4053 Positive
u:Row.Yield-Yield 513.6 694.7 0.7393 Positive
u:units.Yield-Yield 2922.6 294.1 9.9368 Positive
============================================================
Fixed effects:
Trait Effect Estimate Std.Error t.value
1 Yield (Intercept) 132.1 8.791 15.03
============================================================
Groups and observations:
Yield
u:id 363
u:Rowf 13
u:Colf 36
u:Row 168
============================================================
Use the '$' sign to access results and parametersGT矩阵作为随机效应被封装在一个列表中，在vs()中使用。
（10）全基因组选择
你可以用rrBLUP和GBLUP模型拟合。本例子在群体中指定个体。基本形式为：
利用已知的基因型效应方差协方差矩阵作为加性关系矩阵（A），使用A.met函数建立所有个体的联系，并预测未测定个体的BLUP。本例中，标记矩阵是随机效应，标记的关系也可以被指定，但是在这里假设是一个对角线矩阵。
data(DT_wheat)
DT <- DT_wheat
GT <- GT_wheat
colnames(DT) <- paste0("X",1:ncol(DT))
DT <- as.data.frame(DT);DT$id <- as.factor(rownames(DT))
# select environment 1
rownames(GT) <- rownames(DT)
K <- A.mat(GT) # additive relationship matrix
colnames(K) <- rownames(K) <- rownames(DT)
# GBLUP pedigree-based approach
set.seed(12345)
y.trn <- DT
vv <- sample(rownames(DT),round(nrow(DT)/5))
y.trn[vv,"X1"] <- NA
## GBLUP
ans <- mmer(X1~1,
random=~vs(id,Gu=K),
rcov=~units,
data=y.trn,verbose = FALSE) # kinship based
ans$U$`u:id`$X1 <- as.data.frame(ans$U$`u:id`$X1)
rownames(ans$U$`u:id`$X1) <- gsub("id","",rownames(ans$U$`u:id`$X1))
cor(ans$U$`u:id`$X1[vv,],DT[vv,"X1"], use="complete")
[1] 0.5737594
## rrBLUP
ans2 <- mmer(X1~1,
random=~vs(list(GT)),
rcov=~units,
data=y.trn,verbose = FALSE) # kinship based
u <- GT %*% as.matrix(ans2$U$`u:GT`$X1) # BLUPs for individuals
rownames(u) <- rownames(GT)
cor(u[vv,],DT[vv,"X1"]) # same correlation
# the same can be applied in multi-response models in GBLUP or rrBLUP
[1] 0.5737594~1的1表示截距，基本公式为y=C+ε，即常数。
（11）似然比测试
似然比测试（LRT）是调查随机效应或特定方差协方差组分显著性的方式。
（11.1）方差组分显著性测试
例如，想象研究的模型在加入了空间内核效应后改善了，那么模型可能是：
data(DT_cpdata)
DT <- DT_cpdata
GT <- GT_cpdata
MP <- MP_cpdata
### mimic two fields
A <- A.mat(GT)
mix1 <- mmer(Yield~1,
random=~vs(id, Gu=A) +
vs(Rowf) +
vs(Colf),
rcov=~vs(units),
data=DT, verbose = FALSE)
带有空间内核的模型为：
mix2 <- mmer(Yield~1,
random=~vs(id, Gu=A) +
vs(Rowf) +
vs(Colf) +
vs(spl2D(Row,Col)),
rcov=~vs(units),
data=DT,verbose = FALSE)
似然比测试，检测第二个模型：
lrt <- anova(mix1, mix2)
Likelihood ratio test for mixed models
==============================================================
Df AIC BIC loLik Chisq ChiDf PrChisq
mod2 8 304.4021 308.2938 -151.2011
mod1 7 305.0477 308.9393 -151.5238 0.64554 1 0.42171
==============================================================
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1似然值增加了，不显著。
（11.2）测试协方差组分的显著性
有些时候，研究人员对协方差结构是否感相关感兴趣。假设我们有两个多性状模型，1）拟合非协方差（独立），2）拟合群体中产量和颜色的遗传协方差：
data(DT_example)
DT <- DT_example
DT$EnvName <- paste(DT$Env,DT$Name)
modelBase <- mmer(cbind(Yield, Weight) ~ Env,
random= ~ vs(Name, Gtc=diag(2)), # here is diag()
rcov= ~ vs(units, Gtc=unsm(2)),
data=DT,verbose = FALSE)
modelCov <- mmer(cbind(Yield, Weight) ~ Env,
random= ~ vs(us(Env),Name, Gtc=unsm(2)), # here is unsm()
rcov= ~ vs(ds(Env),units, Gtc=unsm(2)),
data=DT,verbose = FALSE)
lrt <- anova(modelBase, modelCov)
Likelihood ratio test for mixed models
==============================================================
Df AIC BIC loLik Chisq ChiDf PrChisq
mod2 45 -351.5889 -328.1079 181.7945
mod1 23 -253.4383 -229.9573 132.7192 98.15058 22 1.3552747104066e-11 ***
==============================================================
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1基因型的协方差显著改善了拟合，卡方分布小于0.05。带有协方差的模型更适合。
直接安装R语言一直是老大难的问题，Ubuntu等系统一般会内置R，但是版本较低，不能满足使用需要，升级通常也是失败的。因此，尝试升级或重新安装最新版本的R语言。经过一番折腾，发现此方法行不通。
替代方案：
以下是安装失败过程：
首先，我使用的是Ubuntu 20.04 TLS版本。
cd /etc/apt/
sudo vi sources.list2. 进入vim编辑模式。接下来，按键盘上的[i]键，进入插入编辑模式。将光标移动到最后一行末尾，回车，插入一行，复制下面的代码，然后按鼠标右键。
deb https://cloud.r-project.org/bin/linux/ubuntu focal-cran40/
# deb-src https://cloud.r-project.org/bin/linux/ubuntu focal-cran40/3. 然后按[ESC]键退出编辑模式，依次输入[:wq!]，回车，即可保存退出。
4. 接下来，输入下面代码升级源信息。
sudo apt-get update5. 若出现如下错误，则需要获得key。
Err:2 https://cran.itam.mx/bin/linux/ubuntu groovy-cran40/ InRelease
The following signatures couldn’t be verified because the public key is not available: NO_PUBKEY 51716619E084DAB9
6. 使用下面代码获得key，注意替换后面的一串数字。
sudo apt-key adv --keyserver hkp://keyserver.ubuntu.com:80 --recv-keys 51716619E084DAB9完成后，重复步骤4。
7. 接下来安装最新版本的R。
sudo apt-get install r-base8. 若出现下列错误，则按照提示操作。
E: Unmet dependencies. Try ‘apt –fix-broken install’ with no packages (or specify a solution).
9. 根据提示，输入下列代码。
sudo apt --fix-broken install10. 重复第7步，若仍然安不上，则先卸载当前R。
sudo apt-get autoremove r-base-core数据框是R语言最重要的数据表现形式，数以万计的数据通过形成一个巨大的“表”来标准化数据结构，再从中抽取部分或全部数据用于分析，这个“表”叫做数据框。数据框和矩阵非常类似，不同点在于，数据框可以指定每一列的数据类型，而整个矩阵只能有且只有一种数据类型。
下面，我们尝试录入数据类型。
数据表4-1
| 姓名 | 性别 | 英语成绩 | 评定 |
| 张三 | 男 | 95 | 优秀 |
| 李四 | 女 | 84 | 良好 |
| 王五 | 男 | 61 | 合格 |
从表4-1中不难发现，数据并非都是数值型，还有字符型的列，比如性别和评定。此时，用数据框就更为合适，因为数据框可以为每一列分配不同的数据类型。数据框在使用上更加符合人类的习惯，因为数据并不只有数值型，还是有很多其他类型，如定性类：姓名、性别、邮政编码、民族、态度（同意、不同意）、分级等。这些可以用字符串（string）来表示，也可以用因子（factor）表示。
我们可以通过下面的方法构建数据框：
names <- c("张三", "李四", "王五")
gender <- c("男", "女", "男")
English_scores <- c(95, 84, 61)
evaluation <- c("优", "良", "合格")
dataframe <- data.frame(names=names,gender=gender,English_scores=English_scores,evaluation=evaluation)
dataframe
class(dataframe)
上面列举了一种通过向量获得数据框的方法，当然还可以先将向量转换成矩阵，再转换成数据框。
dataframe2 <- cbind(names,gender,English_scores,evaluation)
dataframe2 <- as.data.frame(dataframe2)
dataframe2
当只有纸制记录的原始数据，需要录入的时候，比较麻烦，可以通过下面的方式录入的R语言的数据框。
dataframe <- data.frame(names=c("张三", "李四", "王五"),
gender=c("男", "女", "男"),
English_scores=c(95, 84, 61),
evaluation=c("优", "良", "合格"))
dataframe
在数据量大的时候，通过上述方式逐个录入是很容易出错的，最好是先录入到Excel中，在读取到R程序，后面会有章节专门介绍。
R语言的逻辑值（也称布尔值，Boolean）共有两种，TRUE和FALSE，可以简写为T和F，注意都是大写的。它们也可以用数值1和0表示。
一般而言，判断运算得到的是逻辑值。例如，1>2，会返回FALSE；5>=1，会返回TRUE；10==10，会返回TRUE，这里要注意，判断运算的相等需要写2个等号，因为一个等号已经做赋值符号用了等同于'<-‘。
在R语言中，不仅数值可以计算逻辑值，逻辑值和字符串也可以计算逻辑值。你可以尝试使用TRUE==TRUE来计算逻辑结果；也可以利用”a”==”A”来计算逻辑结果；还可以计算TRUE==1。
不等于也是一个判断符号，但是R语言中的不等于符号并不是≠，而是!=。例如，4!=5，会返回TRUE。!的作用是取反，!=的意思是对相等取反，就是不等于的意思。同理，!可以与逻辑值连用，起到对逻辑值取反的作用。如：!TRUE的结果是FALSE；!FALSE的结果是TRUE。!可以和不等式连用，像!5>7得到的结果是TRUE，这里的计算顺序是先算5>7，结果为FALSE，再取反，得到TRUE。可以看出，计算的优先顺序，不等式计算优先级要高于取反运算。
在实际应用中，我们常常会遇到多个条件判断的情形。看这个例子，一群学生考研，英语需要高于40，政治需要高于60，专业课需要高于70分才能录取。此时，就需要多个判断运算同时起作用。
假如：小红的英语、政治和专业课成绩分别为xiaohong <- c(41,75,90)，用判断语句判断小红是否考上的判断是：
xiaohong[1] > 40 & xiaohong[2] > 60 & xiaohong[3] > 70
不难看出，多个条件并列用&符号链接，&表示与，必须所有条件都为TRUE才为TRUE；相应的，|符合表示或，只要有一个条件为TRUE即为TRUE。
当一个逻辑值和一个逻辑值向量并列时，&符号会将单个逻辑值复制为向量的个数，然后一一进行运算。
TRUE & c(TRUE,FALSE,FALSE)
这可能并不是我们想要的结果，如果我们只想对向量中的第一个逻辑值进行运算，可以使用严格并列符号&&和严格或符号||。这样，向量中只有第一个逻辑值会和前面的逻辑值发生运算。
TRUE || c(TRUE, FALSE ,FALSE)
判断函数，isTRUE()和isFALSE()，判断逻辑值的判断函数，说起来很拗口。顾名思义，isTRUE的参数为TRUE时，结果为TRUE，否则为FALSE；isFALSE的参数为FALSE是，结果为TRUE，否则为FALSE。自己可以练习一下。
identical()函数可以判断两个对象是否完全相等。完全相等则为TRUE，不完全相等则为FALSE。传入的两个参数顺序可以颠倒，可以使用变量。
identical(c(1,2,3),c(1,2,3))
xor()函数叫做亦或判断函数，如果两个参数一个是TRUE一个是FALSE则返回TRUE，否则返回FALSE。在该函数中，FALSE可以用0表示，0以外的其他数字都相当于TRUE。
xor(0,-1)
xor(1,-1)
xor(TRUE,FALSE)
逻辑值在筛选数据上具有巨大的作用。当学校的英语成绩为14，34，56，34，78，68，87，69，88，100时，要快速挑出不及格的分数，可以用which函数和逻辑运算。
scores <- c(14, 34, 56, 34, 78, 68, 87, 69, 88, 100)
which(scores < 60)
scores[which(scores < 60)]
当统计全校的分数时，数据量已经很大，很难一眼看出是否有不及格的同学，这个时候可以使用any()函数：只要有一个满足条件即为TRUE。any(scores < 60)返回TRUE，说明数据中有不及格的分数。
scores <- c(14, 34, 56, 34, 78, 68, 87, 69, 88, 100)
any(scores < 60)
同理，我们可以用all()函数判断是否所有人都及格了。很明显，并不是所有人都及格了。
all(scores > 60)
当需要加载R包的时候，用requrie()函数可以返回逻辑值，加载成功为TRUE，加载失败为FALSE。
1. 下列返回FALSE的是哪个选项？
A：9 >= 10
B：7 == 7
C：0 > -36
D：6 < 8
2. 下列哪个选项返回TRUE？
A：9 >= 10
B：57 < 8
C：7 == 9
D：-6 > -7
3. 下列哪个选项返回FALSE？
A：9 < 10
B：!FALSE
C：!(0 >= -1)
D：7 != 8
4. 请说明逻辑与&和严格逻辑与&&的区别？
5. 5>8||6!=8&&4>3.9的结果？
6. 下列结果为TRUE的是？
A：TRUE && FALSE || 10 >= 5 && 1 < 8
B：FALSE || TRUE && FALSE
C：99.99 > 100 || 28 < 2.3 || 8 != 8.0
D：TRUE && 15 < 15 && 6 >= 6
7. 下列结果为FALSE的是？
A：FALSE && 12 >= 12 || 1 >= 4 || 80 <= 49.5
B：6 >= -3 && !(9 > 10) && !(!TRUE)
C：FALSE || TRUE && 8 != 1 || 26 > 3
D：!(8 > 7) || 19 == 19.0 && 7.8 >= 7.79
8. 下列哪个结果是TRUE？
A：isTRUE(!TRUE)
B：isTRUE(3)
C：!isTRUE(4 < 3)
D：!isTRUE(8 != 5)
E：isTRUE(NA)
9. 下列哪个选项的结果为TRUE？
A：!identical(7, 7)
B：identical(4, 3.1)
C：identical(5 > 4, 3 < 3.1)
D：identical(‘hello’, ‘Hello’)
10. xor(12,-1)的结果是什么？
11. 下列选项为FALSE的是哪个？
A：xor(!isTRUE(TRUE), 6 > -1)
B：xor(identical(xor, ‘xor’), 7 == 7.0)
C：xor(!!TRUE, !!FALSE)
D：xor(4 >= 9, 8 != 8.0)
12. 请说明any()函数的作用？
13. 请说明all()函数的作用？