如何使用 R 在单个新数据框中输出多个模型的残差？

Question

我想针对多个不同的因变量为一组静态自变量运行多个回归模型，并将残差输出到一个看起来像...的新文件中。

SampleID     site_residual1    site_residual2    site_residual3
F001         0.003             0.988             0.776
F001         0.002             0.876             0.665
F002         0.134             0.234             0.786
...

我一直在使用以下代码来获取单个剩余输出，但未能成功实现将贯穿我所有站点的循环。

infile = sprintf("/path/siteinput.txt.gz")

infile 看起来像...

SampleID     site1  site2   site3   etc...
F001         0.003  0.988   0.776   etc...
F001         0.002  0.876   0.665   etc...
F002         0.134  0.234   0.786   etc...
...

...

pheno = read.table("/path/pheno_covar.txt", header=T, sep="\t")

现象看起来...

SampleID     indep1 indep2  indep3  chip1   etc...
F001         0.003  0.988   0.776   2       etc...
F001         0.002  0.876   0.665   2       etc...
F002         0.134  0.234   0.786   1       etc...
...

...

residfile = sprintf("/path/test_resid_out.txt")

library(lme4)

beta = read.table(infile, header=T, sep="\t")

merged = merge(beta, pheno, by="SampleID")

site<-merged$site1
chip <- as.factor(merged$chip1)

model1 <- lmer (formula= site ~ indep1 +indep2 + indep3 + (1|chip), data=merged)

print(summary(model1))
print(resid(model1))

site1_resid = resid(model1, na.action=na.exclude)

residout<-(data.frame(SampleID, site1_resid))
write.table(residout, file=residfile, sep="\t", col.names=TRUE, row.names=FALSE, quote=FALSE)

我的输出看起来像......

SampleID    site1_resid
F001        0.0110177454696274
F002        0.0923483180517723
F003        0.103686493563883
F004        -0.106193404096636
F005        -0.124621172636435
....

...

所以，我真的在寻找一种方法来为我的“infile”中的每个站点运行 model1 并将所有残差输出到一个新文件中。另外，我希望列标题包含“站点”的原始名称。我确实有一些缺失的信息（所有协变量都是完整的，但是某些 ID 缺少一些站点）。

任何建议将不胜感激。

Answer 1

借助 magrittr 管道 (%>%) 使其更易于阅读（虽然不是必需的）：

library(magrittr)
names(beta) %>% 
  setdiff("SampleID") %>% 
  setNames(., .) %>% 
  lapply(function(x) {
    model <- lmer(data = merged, formula = paste(x, "~ indep1 +indep2 + indep3 + (1|chip)"))
    # print(summary(model))
    # print(resid(model))
    resid(model, na.action=na.exclude)
  }) %>% 
  c(list(SampleID = merged$SampleID), .) %>% 
  do.call(what = "data.frame")

（顺便说一句，我担心你有重复的SampleIDs。这是有意的吗？如果是，你确定要SampleIDmerge()吗？你宁愿不要cbind(beta, pheno[, - 1, drop = FALSE]) ?)