目录

• 简介
• GS示例代码

简介

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: 协方差矩阵，协方差结构