R中glm的输出不包括变量名

Question

我正在编写一个 R 脚本，使用 glm.nb

测试 150,000 个遗传标记与连续变量的关联

我写了以下内容来做到这一点：

fhandle<-file("ichip_nb_model.csv","a")
for (i in seq(7, ncol(ICHPdt), 1)) {
     glmmod<-glm.nb(OverllTot0 ~ EurAdmix + Sex + DisDurMonths + BMI + Group + SmokingStatus + eval(parse(text = paste("ICHPdt$", colnames(ICHPdt)[i], sep=""))), data=covfiledt)
     writeLines(capture.output(coef(summary(glmmod))), con=fhandle)
     writeLines(colnames(ICHPdt)[i], con=fhandle)
 }

然而，这会导致一个问题。第一个 writeLines 语句中写入的输出不包含列名，而是包含整个 eval 表达式（请参阅下面的当前输出部分）。但我不想要 eval 表达式，我想要 评估为 的内容，因为列 header 是正在测试的遗传标记的名称。

因此，作为 stop-gap，我添加了第二个 writeLines 语句，但我更愿意找到一种解决方案，它只会导致显示变体的名称。脚本生成的输出格式如下：

当前输出

    Estimate      Std. Error    z value
        (Intercept)                                                         -0.960341597 0.898711395 -1.0685762
        EurAdmix                                                             2.055048065 1.132148532  1.8151753
        Sex                                                                  0.783616302 0.369298081  2.1219073
        DisDurMonths                                                        -0.013458018 0.002786449 -4.8298098
        BMI                                                                 -0.008077163 0.012892372 -0.6265071
        Group                                                               -0.059876340 0.288834615 -0.2073032
        SmokingStatus                                                       -0.089029296 0.185598042 -0.4796888
        eval(parse(text = paste("ICHPdt$", colnames(ICHPdt)[i], sep = ""))) -0.108334664 0.169464866 -0.6392751
                                                                                Pr(>|z|)
        (Intercept)                                                         2.852607e-01
        EurAdmix                                                            6.949697e-02
        Sex                                                                 3.384552e-02
        DisDurMonths                                                        1.366635e-06
        BMI                                                                 5.309824e-01
        Group                                                               8.357731e-01
        SmokingStatus                                                   

6.314487e-01
    eval(parse(text = paste("ICHPdt$", colnames(ICHPdt)[i], sep = ""))) 5.226440e-01

期望的输出：

        Estimate      Std. Error    z value
            (Intercept)                                                         -0.960341597 0.898711395 -1.0685762
            EurAdmix                                                             2.055048065 1.132148532  1.8151753
            Sex                                                                  0.783616302 0.369298081  2.1219073
            DisDurMonths                                                        -0.013458018 0.002786449 -4.8298098
            BMI                                                                 -0.008077163 0.012892372 -0.6265071
            Group                                                               -0.059876340 0.288834615 -0.2073032
            SmokingStatus                                                       -0.089029296 0.185598042 -0.4796888
            ICHPdt$rs728931
-0.108334664 0.169464866 -0.6392751


        Pr(>|z|)
                (Intercept)                                                         2.852607e-01
                EurAdmix                                                            6.949697e-02
                Sex                                                                 3.384552e-02
                DisDurMonths                                                        1.366635e-06
                BMI                                                                 5.309824e-01
                Group                                                               8.357731e-01
                SmokingStatus                                                       6.314487e-01
                ICHPdt$rs728931
5.226440e-01

Answer 1

几乎不需要使用 eval(parse())。在这种情况下，将公式构建为字符串然后将其发送到调用会更容易。我假设您正在使用 MASS 包中的 glm.nb() 函数。这是动态构建公式的示例

library(MASS)

other <- data.frame(matrix(runif(nrow(quine)*3), nrow=nrow(quine)))    

lapply(names(other), function(x) {
    ff<-as.formula(paste0("Days ~ Sex/(Age + Eth*Lrn) + other$", x))
    glm.nb(ff, data = quine)
})

虽然实际上，从不同的数据集中获取一些值通常不是一个好主意。公式中有 $ 通常是一个不好的迹象。您可能会考虑将额外的协变量数据合并到 data.frame 中。在这里，我还展示了使用 bquote()

构建公式的另一种方法

lapply(names(other), function(x) {
    glm.nb(bquote(Days ~ Sex/(Age + Eth*Lrn) + .(as.name(x))), 
        data = cbind(quine, other))
})

Answer 2

summary.glm(fit)$coefficients 的结果就是您想要的结果，将第 8 个行名的名称更改为所需的字符值很简单。也许使用 write.table 会节省一些 capture.output 强加的笨拙（当然，除非你真的想要使用 'Pr(>|t|)' 名称值的环绕：

fhandle<-file("ichip_nb_model.csv","a")
for (i in seq(7, ncol(ICHPdt), 1)) {
     glmmod<-glm.nb(OverllTot0 ~ EurAdmix + Sex + DisDurMonths + BMI + Group + SmokingStatus + eval(parse(text = paste("ICHPdt$", colnames(ICHPdt)[i], sep=""))), data=covfiledt)
     summ <- coef(summary(glmmod)))
     rownames(summ)[8] <- paste0("ICHPdt$", colnames(ICHPdt)[i])
     write.table( round(summ, 4) file=fhandle)
 }

R中glm的输出不包括变量名

output from glm in R does not include variable name

r

glm