下面的代码当前对每个结果的每个暴露运行未调整的glm(每个结果3个暴露(，并将结果导出到列表中。每次曝光，我需要3个模型：模型1：未调整(我们目前有(，模式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个CON：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)
就像一开始指定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"(得到一个"0"；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)
我建议收集您希望更改的不同组件在模型之间转换为数据帧，并相应地构建：
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