如何解决如何在R中使用glm循环多次曝光和结果以及不同模型？
下面的代码当前针对每个结果的每次曝光运行未调整的glm（每个结果3次曝光），并将结果导出到列表中。对于每次曝光，我需要3个模型： 模型1 ：未经调整（我们目前拥有），模型2 ：针对cov1进行了调整，模型3 ：针对cov1，cov2和cov3进行了调整
如何在此代码中实现不同的模型？
amino_df <- data.frame(y = rbinom(100,1,0.5),y2 = rbinom(100,0.3),y3 = rbinom(100,0.2),y4 = rbinom(100,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,"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,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,0.3),y3 = rbinom(100,0.2),y4 = rbinom(100,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",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,data = amino_df)
obs_results = data.frame(
outcome = j,covariate = gsub(names(coef(glmmodel))[2],"",h),# gsub to remove exposure from covariate(s)
OR = exp(glmmodel$coefficients[2]),2])
)
return(obs_results)
}) %>% bind_rows
}) %>% bind_rows %>% `colnames<-`(gsub("X95",colnames(.))) %>% `rownames<-`(NULL)
我建议收集您希望更改的不同组件 在模型之间建立数据框，并相应地构建模型：
library(tidyverse)
y <- c("y","y4")
exp <- c("exp1","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,data = amino_df)
# using data.frame would not handle the covariate list column properly
obs_results = tibble(
outcome = outcome,covariate = list(covariates),2])
)
return(obs_results)
}
有了这些参数后，您可以使用pmap()为每一行拟合模型
在您的规范数据框中：
amino_df <- data.frame(y = rbinom(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
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。