如何在 R 中使用 glm 循环多次曝光和结果以及不同的模型？答案

【问题标题】：How to loop multiple exposures and outcomes as well as different models with glm in R?如何在 R 中使用 glm 循环多次曝光和结果以及不同的模型？
【发布时间】：2020-11-12 08:19:53
【问题描述】：

下面的代码目前针对每个结果的每次曝光（每个结果 3 次曝光）运行未经调整的 glm，并将结果导出到列表中。对于每次曝光，我需要 3 个模型： model 1：未调整（我们目前有），model 2：针对 cov1 进行了调整，model 3：针对 cov1、cov2 和 cov3 进行了调整

我将如何在这段代码中实现不同的模型？

amino_df <- data.frame(y = rbinom(100, 1, 0.5), y2 = rbinom(100, 1, 0.3), y3 = rbinom(100, 1, 0.2), y4 = rbinom(100, 1, 0.22),
                       exp1 = rnorm(100), exp2 = rnorm(100), exp3 = rnorm(100),
                       cov1 = rnorm(100), cov2 = rnorm(100), cov3 = rnorm(100))

exp <- c("exp1", "exp2", "exp3")
y <- c("y", "y2","y3","y4")
cov <- c("cov1", "cov2", "cov3")

obs_results <- replicate(length(y), data.frame())  

for(j in seq_along(y)){
  for (i in seq_along(exp)){
    mod <- as.formula(paste(y[j], "~", exp[i]))
    glmmodel <- glm(formula = mod, family = binomial, data = amino_df)
    
    obs_results[[j]][i,1] <- names(coef(glmmodel))[2]
    obs_results[[j]][i,2] <- exp(glmmodel$coefficients[2])
    obs_results[[j]][i,3] <- summary(glmmodel)$coefficients[2,2]
    obs_results[[j]][i,4] <- summary(glmmodel)$coefficients[2,4]
    obs_results[[j]][i,5] <- exp(confint.default(glmmodel)[2,1])
    obs_results[[j]][i,6] <- exp(confint.default(glmmodel)[2,2])
  }
  colnames(obs_results[[j]]) <- c("exposure","OR", "SE", "P_value", "95_CI_LOW","95_CI_HIGH")
}
names(obs_results) <- y

obs_df <- do.call("rbind", lapply(obs_results, as.data.frame))

编辑 - 我现在有一个解决方案：

还有一个问题，下面的代码是否可以调整为包含针对不同曝光的不同模型？所以对于 exp1，调整所有 3 个缺点：cov1、cov2、cov3，但对于 exp2，只调整 cov1、cov2？只有 exp3 cov2 和 cov1？

amino_df <- data.frame(y = rbinom(100, 1, 0.5), y2 = rbinom(100, 1, 0.3), 
                       y3 = rbinom(100, 1, 0.2), y4 = rbinom(100, 1, 0.22),
                       exp1 = rnorm(100), exp2 = rnorm(100), exp3 = rnorm(100),
                       cov1 = rnorm(100), cov2 = rnorm(100), cov3 = rnorm(100))

exp <- c("exp1", "exp2", "exp3")
y <- c("y", "y2","y3","y4")
model <- c("", "+ cov1", "+ cov1 + cov2 + cov3")

obs_df <- lapply(y, function(j){
    lapply(exp, function(i){
        lapply(model, function(h){
            mod = as.formula(paste(j, "~", i, h))
            glmmodel = glm(formula = mod, family = binomial, data = amino_df)
            
            obs_results = data.frame(
                outcome = j,
                exposure = names(coef(glmmodel))[2], 
                covariate = h,
                OR = exp(glmmodel$coefficients[2]), 
                SE = summary(glmmodel)$coefficients[2,2], 
                P_value = summary(glmmodel)$coefficients[2,4], 
                `95_CI_LOW` = exp(confint.default(glmmodel)[2,1]), 
                `95_CI_HIGH` = exp(confint.default(glmmodel)[2,2])
            )
            return(obs_results)
        }) %>% bind_rows
    }) %>% bind_rows
}) %>% bind_rows %>% `colnames<-`(gsub("X95","95",colnames(.))) %>% `rownames<-`(NULL)

head(obs_df)

【问题讨论】：

标签： r for-loop glm

【解决方案1】：

就像你一开始指定exp和y一样，你可以指定不同的模型类型。

这是一种使用 lapply() 代替 for 循环的方法：

amino_df <- data.frame(y = rbinom(100, 1, 0.5), y2 = rbinom(100, 1, 0.3), 
                       y3 = rbinom(100, 1, 0.2), y4 = rbinom(100, 1, 0.22),
                       exp1 = rnorm(100), exp2 = rnorm(100), exp3 = rnorm(100),
                       cov1 = rnorm(100), cov2 = rnorm(100), cov3 = rnorm(100))

exp <- c("exp1", "exp2", "exp3")
y <- c("y", "y2","y3","y4")
model <- c("", "+ cov1", "+ cov1 + cov2 + cov3")

obs_df <- lapply(y, function(j){
    lapply(exp, function(i){
        lapply(model, function(h){
            mod = as.formula(paste(j, "~", i, h))
            glmmodel = glm(formula = mod, family = binomial, data = amino_df)
            
            obs_results = data.frame(
                outcome = j,
                exposure = names(coef(glmmodel))[2], 
                covariate = h,
                OR = exp(glmmodel$coefficients[2]), 
                SE = summary(glmmodel)$coefficients[2,2], 
                P_value = summary(glmmodel)$coefficients[2,4], 
                `95_CI_LOW` = exp(confint.default(glmmodel)[2,1]), 
                `95_CI_HIGH` = exp(confint.default(glmmodel)[2,2])
            )
            return(obs_results)
        }) %>% bind_rows
    }) %>% bind_rows
}) %>% bind_rows %>% `colnames<-`(gsub("X95","95",colnames(.))) %>% `rownames<-`(NULL)

head(obs_df)
#  outcome exposure            covariate        OR        SE   P_value 95_CI_LOW 95_CI_HIGH
#1       y     exp1                      0.9425290 0.2125285 0.7806305 0.6214270   1.429550
#2       y     exp1               + cov1 0.9356460 0.2138513 0.7557639 0.6152917   1.422794
#3       y     exp1 + cov1 + cov2 + cov3 0.9638427 0.2174432 0.8655098 0.6293876   1.476027
#4       y     exp2                      1.3297429 0.1865916 0.1266809 0.9224452   1.916879
#5       y     exp2               + cov1 1.3300740 0.1866225 0.1264124 0.9226190   1.917473
#6       y     exp2 + cov1 + cov2 + cov3 1.3558196 0.1903111 0.1097054 0.9337031   1.968770

我在末尾添加了gsub("X95","95",colnames(.))，因为在创建新数据框时，默认情况下，以数字开头的列名（即“95_CI_LOW”、“95_CI_HIGH”）会在开头插入一个“X”；此代码将其删除。

补充

如果在您的模型中使用不同的协变量对不同的曝光进行了唯一调整，则可以改为执行以下操作。最简单的解决方案是通过上面的代码运行所有可能的曝光+协变量组合，然后过滤obs_df（带有filter()）以仅选择您想要的分析。但是，如果您使用大型数据集，这意味着运行时间会不必要地延长。

更直接的方法是在model 中具体输入要包含哪些曝光+协变量组合并删除lapply(exp) 函数（并相应地编辑核心函数）：

model <- c("exp1 + cov1 + cov2 + cov3", "exp2 + cov1 + cov2", "exp3 + cov1")

obs_df <- lapply(y, function(j){
    lapply(model, function(h){
        mod = as.formula(paste(j, "~", h))
        glmmodel = glm(formula = mod, family = binomial, data = amino_df)
            
        obs_results = data.frame(
            outcome = j,
            exposure = names(coef(glmmodel))[2], 
            covariate = gsub(names(coef(glmmodel))[2],"",h), # gsub to remove exposure from covariate(s)
            OR = exp(glmmodel$coefficients[2]), 
            SE = summary(glmmodel)$coefficients[2,2], 
            P_value = summary(glmmodel)$coefficients[2,4], 
            `95_CI_LOW` = exp(confint.default(glmmodel)[2,1]), 
            `95_CI_HIGH` = exp(confint.default(glmmodel)[2,2])
        )
        return(obs_results)
    }) %>% bind_rows
}) %>% bind_rows %>% `colnames<-`(gsub("X95","95",colnames(.))) %>% `rownames<-`(NULL)

【讨论】：

谢谢！有没有办法针对不同的曝光调整不同的模型？所以对于 exp1，调整 cov1、cov2、cov3，但对于 exp2，只调整 cov1、cov2？
有办法的。选项 1) 运行所有可能的曝光+协变量组合（如上面的代码），然后使用 filter() 选择您想要的——这将是对代码的最简单修改，但运行不需要的分析会花费不必要的时间。选项 2) 更改 model 以仅包含您想要的曝光+协变量组合，然后从 mod 公式中删除 lapply(exp, ...) 循环和“i”。
@KurtTaylor，我用补充更新了答案以回答您的后续问题
非常感谢！！在模拟数据集上效果很好。知道为什么在使用我的真实数据集时会得到这个，而 CHD_any 是我唯一的结果吗？ Error in str2lang(x) : <text>:2:0: unexpected end of input 1: CHD_any ~
对不起，我以前没有看到过这个错误。可能要检查的事情： CHD_any 是数字类还是因子类？您的其他变量是数字/整数还是因子/字符类？有错别字吗？如果您要使用变量运行单个 glm() 模型，您还会收到错误吗？

【解决方案2】：

我建议收集您希望改变的不同组件模型之间的数据框，并相应地构建：

library(tidyverse)

y <- c("y", "y2","y3","y4")
exp <- c("exp1", "exp2", "exp3")
cov <- list(character(), "cov1", c("cov1", "cov2", "cov3"))

# each covariate for each exposure
models1 <- crossing(outcome = y, exposure = exp, covariates = cov)
models1
#> # A tibble: 36 x 3
#>    outcome exposure covariates
#>    <chr>   <chr>    <list>    
#>  1 y       exp1     <chr [0]> 
#>  2 y       exp1     <chr [1]> 
#>  3 y       exp1     <chr [3]> 
#>  4 y       exp2     <chr [0]> 
#>  5 y       exp2     <chr [1]> 
#>  6 y       exp2     <chr [3]> 
#>  7 y       exp3     <chr [0]> 
#>  8 y       exp3     <chr [1]> 
#>  9 y       exp3     <chr [3]> 
#> 10 y2      exp1     <chr [0]> 
#> # ... with 26 more rows

# covariates specific per exposure
models2 <- crossing(outcome = y, nesting(exposure = exp, covariates = cov))
models2
#> # A tibble: 12 x 3
#>    outcome exposure covariates
#>    <chr>   <chr>    <list>    
#>  1 y       exp1     <chr [0]> 
#>  2 y       exp2     <chr [1]> 
#>  3 y       exp3     <chr [3]> 
#>  4 y2      exp1     <chr [0]> 
#>  5 y2      exp2     <chr [1]> 
#>  6 y2      exp3     <chr [3]> 
#>  7 y3      exp1     <chr [0]> 
#>  8 y3      exp2     <chr [1]> 
#>  9 y3      exp3     <chr [3]> 
#> 10 y4      exp1     <chr [0]> 
#> 11 y4      exp2     <chr [1]> 
#> 12 y4      exp3     <chr [3]>

然后把你的模型拟合和总结放到一个函数中这些组件：

fit_model <- function(outcome, exposure, covariates) {
  formula = reformulate(c(exposure, covariates), outcome)
  glmmodel = glm(formula = formula, family = binomial, data = amino_df)
  
  # using data.frame would not handle the covariate list column properly
  obs_results = tibble(
    outcome = outcome,
    exposure = names(coef(glmmodel))[2], 
    covariate = list(covariates),
    OR = exp(glmmodel$coefficients[2]), 
    SE = summary(glmmodel)$coefficients[2,2], 
    P_value = summary(glmmodel)$coefficients[2,4], 
    `95_CI_LOW` = exp(confint.default(glmmodel)[2,1]), 
    `95_CI_HIGH` = exp(confint.default(glmmodel)[2,2])
  )
  
  return(obs_results)
}

有了这些，您可以使用pmap() 为每一行拟合模型在您的规范数据框中：

amino_df <- data.frame(y = rbinom(100, 1, 0.5), y2 = rbinom(100, 1, 0.3), 
                       y3 = rbinom(100, 1, 0.2), y4 = rbinom(100, 1, 0.22),
                       exp1 = rnorm(100), exp2 = rnorm(100), exp3 = rnorm(100),
                       cov1 = rnorm(100), cov2 = rnorm(100), cov3 = rnorm(100))

# each covariate for each exposure
pmap_df(models1, fit_model)
#> # A tibble: 36 x 8
#>    outcome exposure covariate    OR    SE P_value `95_CI_LOW` `95_CI_HIGH`
#>    <chr>   <chr>    <list>    <dbl> <dbl>   <dbl>       <dbl>        <dbl>
#>  1 y       exp1     <chr [0]> 1.01  0.191  0.944        0.697         1.47
#>  2 y       exp1     <chr [1]> 1.01  0.191  0.947        0.697         1.47
#>  3 y       exp1     <chr [3]> 0.990 0.194  0.960        0.677         1.45
#>  4 y       exp2     <chr [0]> 1.26  0.215  0.281        0.827         1.92
#>  5 y       exp2     <chr [1]> 1.29  0.220  0.244        0.840         1.99
#>  6 y       exp2     <chr [3]> 1.31  0.222  0.229        0.845         2.02
#>  7 y       exp3     <chr [0]> 1.43  0.216  0.0969       0.937         2.19
#>  8 y       exp3     <chr [1]> 1.43  0.217  0.101        0.933         2.18
#>  9 y       exp3     <chr [3]> 1.36  0.221  0.166        0.881         2.09
#> 10 y2      exp1     <chr [0]> 1.55  0.230  0.0580       0.985         2.43
#> # ... with 26 more rows

# covariates specific per exposure
pmap_df(models2, fit_model)
#> # A tibble: 12 x 8
#>    outcome exposure covariate    OR    SE P_value `95_CI_LOW` `95_CI_HIGH`
#>    <chr>   <chr>    <list>    <dbl> <dbl>   <dbl>       <dbl>        <dbl>
#>  1 y       exp1     <chr [0]> 1.01  0.191  0.944        0.697         1.47
#>  2 y       exp2     <chr [1]> 1.29  0.220  0.244        0.840         1.99
#>  3 y       exp3     <chr [3]> 1.36  0.221  0.166        0.881         2.09
#>  4 y2      exp1     <chr [0]> 1.55  0.230  0.0580       0.985         2.43
#>  5 y2      exp2     <chr [1]> 0.717 0.249  0.182        0.441         1.17
#>  6 y2      exp3     <chr [3]> 0.999 0.241  0.996        0.622         1.60
#>  7 y3      exp1     <chr [0]> 1.21  0.243  0.442        0.749         1.94
#>  8 y3      exp2     <chr [1]> 0.822 0.267  0.463        0.487         1.39
#>  9 y3      exp3     <chr [3]> 1.56  0.269  0.0980       0.921         2.64
#> 10 y4      exp1     <chr [0]> 1.12  0.224  0.601        0.725         1.74
#> 11 y4      exp2     <chr [1]> 0.721 0.255  0.200        0.437         1.19
#> 12 y4      exp3     <chr [3]> 0.767 0.252  0.291        0.468         1.26

【讨论】：