簡介
R包sommer內置了C++，運算速度還是比較快的，功能也很豐富，可求解各種復雜模型。語法相比於lme4包也要好懂一些。
建議查看文檔：vignette("v1.sommer.quick.start")
混合線性模型關鍵在於協方差結構的建立，有以下幾類：
- 復合對稱（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")
- 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
兩者結果相差不大（如果去掉隨機種子，循環運行的結果相差還是很大的），運算時間相差比較大。
Ref: 協方差矩陣，協方差結構