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

Loops through variable lists - embedded loops

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

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

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’…
    }
}

希望能帮到你。

谢谢。

乔希

试试这个。使用 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 package (v0.3.0)

于 2020-03-28 创建