R: 如何制作 MCMCglmm 循环

Question

我想通过从数据中提取系统发育信号来转换我的基因表达数据。我使用 R 包 MCMCglmm。我可以将 MCMCglmm 应用于表达式列之一:

require(ape)
library("MCMCglmm")
expr1 <- c(5,6,5, 11,12,13, 32,33,36)
expr2 <- c(1100,1212,1333, 32,33,36, 34, 38, 49)
expr3 <- c(32,33,36, 110,120,130, 320,330,360)
animal <- seq.int(9)
popGroup <- c(rep('A', 3),rep('B', 3), rep('C', 3))
data <- as.data.frame(cbind(expr1, expr2, expr3, animal, popGroup))
class(data$expr1)<-'integer'
class(data$expr3)<-'integer'
class(data$expr2)<-'integer'

# tree file content:
# (((1:2.0,(2:1.0,3:1.0):1.0):3.0,((4:1.0,5:1.0):1.0,6:2.0):3.0):3.0,(7:2.0,(8:1.0,9:1.0):1.0):6.0);
tree <- read.tree("tree.nwk")

prior<-list(R=list(V=1, nu=1), G=list(G1=list(V=1, nu=1)))
summary(MCMCglmm(expr1~popGroup-1, random=~animal, pedigree=tree, data=data, family="poisson", prior = prior))

但我有超过 20000 个这样的列。所以，我的想法是遍历所有这些:

for (i in 1:3) {
  M <- ( (colnames(data)[i]~popGroup-1, random=~animal, pedigree=tree, data=data, family="poisson", prior = prior))
  summM <- summary(M)
  statM <- summM$statistics[,1:2]
  print(statM)
}

问题在于在循环中定义响应变量。我尝试了很多方法，但都不起作用。

Answer 1

这是包 MCMCglmm 的作者 Jarrod Hadfield 的解决方案:

Unfortunately MCMCglmm works slightly differently than glm when defining response variables. What you could do is:

fixed<-as.formula(paste(colnames(data)[i], "~popGroup-1", sep=""))

MCMCglmm(fixed=fixed, ...)

脚本现在可以运行了。