时间戳

首页 > 解决方案 > 循环遍历变量列表 - 嵌入式循环

r - 循环遍历变量列表 - 嵌入式循环

问题描述

我是一个长期的 STATA 用户，终于尝试掌握 R，以便我可以改进我的图形。在以下代码中，我使用 5 个暴露变量 x1 - x5 将 GAM 拟合到结果变量 y1，然后绘制 x1 的预测。

我想要做的是有两个循环，一个嵌入到另一个循环中，以便第一个循环迭代 5 个结果，为每个结果拟合 GAM，然后在嵌入的第二个循环中，它迭代 5 个曝光，绘制预测对于每个。结果将是来自五个 GAM 中的每一个的 25 个 5 个变量的图。在真实数据库中，变量没有编号，因此它必须将变量名称作为字符串循环。

y1.gam <- mgcv::gam(y1~s(x1,bs="cr",fx=TRUE)+
                   s(x2,bs="cr",fx=TRUE)+
                   s(x3,bs="cr",fx=TRUE)+
                   s(x4,bs="cr",fx=TRUE)+
                   s(x5,bs="cr",fx=TRUE)+
                  family = poisson(link = "log"),
                  data = data)
y1.x1.plot <- plotGAM(gamFit = y1.gam , smooth.cov = "x1", groupCovs = NULL,
       plotCI=TRUE, orderedAsFactor = FALSE)

如果它有帮助，这就是它在 STATA 中的运行方式：

global outcome y1 y2 y3 y4 y5
global exposure x1 x2 x3 x4 x5

foreach v of varlist $outcome {
    gam `v’ $exposure, …
    foreach w of varlist $exposure{
         plot `w’…
    }
}

希望你能帮忙。

谢谢。

乔什

标签： rloopsggplot2regressiongam

解决方案

尝试这个。用作mtcars示例数据集，即使它肯定不是此类模型的合适示例数据集。但是，我希望它足以展示一般方法。循环的关键是substitute用于为估计步骤设置公式对象。结果是一个包含模型、绘图、公式、...的列表

vars_outcome <- c("mpg", "disp")
vars_exposure <- c("hp", "qsec")

# Grid of outcome and exposure variables
vars_grid <- expand.grid(out = vars_outcome, exp = vars_exposure, stringsAsFactors = FALSE)
# Init list for formulas, models, plots
mods <- list(out = vars_grid$out, exp = vars_grid$exp, fmla = list(), mod = list(), mod = list())

for (i in seq_len(nrow(vars_grid))) {
  # Set up the formula
  mods$fmla[[i]] <- substitute(out ~ s(exp, bs="cr",fx=TRUE), list(out = as.name(mods$out[[i]]), exp = as.name(mods$exp[[i]])))
  # Estimate Model
  mods$mod[[i]] <- mgcv::gam(mods$fmla[[i]], family = poisson(link = "log"), data = mtcars)
  # Plot Model
  mods$plt[[i]] <- voxel::plotGAM(gamFit = mods$mod[[i]] , smooth.cov = mods$exp[[i]], groupCovs = NULL, plotCI=TRUE, orderedAsFactor = FALSE)

  # Create a "variable" containing the plot 
  assign(paste(mods$out[[i]], mods$exp[[i]], sep = "_"), mods$plt[[i]])

}

## Name the list with the plots
names(mods$plt) <- paste(mods$out, mods$exp, sep = "_")

mods$fmla[[1]]
#> mpg ~ s(hp, bs = "cr", fx = TRUE)
mods$mod[[1]]
#> 
#> Family: poisson 
#> Link function: log 
#> 
#> Formula:
#> mpg ~ s(hp, bs = "cr", fx = TRUE)
#> attr(,".Environment")
#> <environment: R_GlobalEnv>
#> 
#> Estimated degrees of freedom:
#> 9  total = 10 
#> 
#> UBRE score: -0.1518201
mods$plt[[1]]

由reprex 包（v0.3.0）于 2020 年 3 月 28 日创建

推荐阅读