简介

R包sommer内置了C++，运算速度还是比较快的，功能也很丰富，可求解各种复杂模型。语法相比于lme4包也要好懂一些。
建议查看文档：vignette("v1.sommer.quick.start")

image.png

混合线性模型关键在于协方差结构的建立，有以下几类：

• 复合对称（Compound Symmetry，CS），所有方差相等，所有协方差也相等，对应于单变量方法。但是对于不同尺度的变量是无意义的.
• 方差组分（Variance Components），每个方差都不相同，并且全部协方差等于0。如果变量完全独立，并且彼此测量尺度不同，才是有意义的模式。
• 非结构化（Unstructured），没有模式，每个方差和协方差完全不同，彼此间没有关系，对应于多变量方法。协方差结构有很多种，只有在特定的统计条件下才有意义。

调用模型并不难，难的是在理解的基础上如何随心所欲地应用。

GS示例代码

• 预处理

library(sommer)

data(DT_wheat)
DT <- DT_wheat
GT <- GT_wheat
dim(GT)
GT[1:10,1:10]

colnames(DT) <- paste0("X",1:ncol(DT))
DT <- as.data.frame(DT)
DT$id <- as.factor(rownames(DT))
rownames(GT) <- rownames(DT)

# check NA
which(is.na(GT))
which(is.na(DT))

set.seed(12345)
y.trn <- DT

#制造1/5的缺失值，作为验证集
vv <- sample(rownames(DT),round(nrow(DT)/5))
y.trn[vv,"X1"] <- NA
y.trn[1:5,1:5]

• GBLUP

#######-----------GBLUP--------------------------------------
# GBLUP pedigree-based approach
K <- A.mat(GT) # additive relationship matrix
K[1:5,1:5]
colnames(K) <- rownames(K) <- rownames(DT)
#test first trait X1
system.time(
  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")

image.png

• RRBLUP

#######-----------rrBLUP--------------------------------------
system.time(
  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

image.png

两者结果相差不大（如果去掉随机种子，循环运行的结果相差还是很大的），运算时间相差比较大。

Ref: 协方差矩阵，协方差结构