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

我想为一组静态针对多个不同的因变量和输出残差到一个新文件，看起来像自变量的运行多元回归模型...如何使用R在单个新数据框中输出来自多个模型的残差？

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）。

任何意见，将不胜感激。

来源

2017-03-15 JustGenetics

也许看看'broom'包是为这种情况设计的：https：// cran。 r-project.org/web/packages/broom/vignettes/broom.html –

谢谢您的建议。扫帚看起来对于向原始数据框添加残差看起来很有用，但我真的很有兴趣为每个结果变量创建一个具有残差列的新数据框。我在Broom包中没有看到这个功能，但也许我错过了这个地方？ – JustGenetics

随着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")

（在一个侧面说明，我很担心，你有重复SampleID S是它的目的是什么？如果是，你确定要merge()SampleID？你宁愿做cbind(beta, pheno[, - 1, drop = FALSE])？）

来源

2017-03-15 23:30:34

SampleIDs不应包含重复项，但会产生1：1合并。但是，它们在两个文件中的排序方式不同。 – JustGenetics

谢谢，@ a-p-o-m。当我应用你的建议时，我得到以下错误...'Data.frame中的错误（SampleID = c（1L，2L，3L，4L，5L，6L，7L，8L，9L，10L，：参数意味着不同的数字行：2684，2676 调用：％>％... withVisible - > - > do.call - > data.frame Execution halted'。我假设这是来自某些网站的缺失值，但我是不清楚为什么'na.action = na.exclude'没有处理这个问题，欢迎和赞赏任何额外的建议 – JustGenetics

谢谢@Apom我可以在'na.action = na.exclude '放入模型声明中而不是'resid'。 – JustGenetics

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

回答

相关问题