来自 GAM 的卡方值的子集导致 R答案

【问题标题】：Subset by Chi Sq Values from GAM results in R来自 GAM 的卡方值的子集导致 R
【发布时间】：2021-01-27 20:25:44
【问题描述】：

我在 R 中运行了一个 for 循环，它为 200 个不同的随机数据组合（200 个不同的 set.seed 值）生成二项式 GAM 模型。

for-loop 和 GAM 运行良好，它们将模型存储在适当的列表 model[[i]] 中，每个列表元素代表特定 set.seed 迭代的模型。

我可以在单个列表元素（例如model[[5]]）上运行summary()，然后得到如下结果：

Approximate significance of smooth terms:
                      edf Ref.df Chi.sq  p-value    
s(Petal.Width)  1.133e+00      9  5.414 0.008701 ** 
s(Sepal.Length) 8.643e-01      9  2.338 0.067509 .  

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.361   Deviance explained = 32.7%
-REML = 83.084  Scale est. = 1         n = 160

由于我的 model 列表中有 200 个这样的元素，我想知道是否有一种快速的方法可以计算出 Chi.sq 值等于多少次（总共 200 个） Petal.Width 变量为 0。基本上，由于我使用bs = "cs" 设置了 GAM，因此 Chi.sq 等于 0 的次数表示变量选择过程从模型中删除该变量的频率。

这是我用于 for 循环构建模型的代码的清理版本：


a <- c(1:200)
model <- list()


for (i in a) {
  
#In my version there's some extra code here that subsets iris based on i 
  
  model[[i]] <- gam(Species~ s(Petal.Width, bs="cs") + 
                      s(Sepal.Length, bs="cs"),
                    data = iris,
                    family = binomial,
                    method = "REML")

}

【问题讨论】：

您能提供您的模型的代码吗？我使用了broom 包中的tidy() 命令，它将模型结果放在一个小标题中以收集输出，然后您可以将列放在一起。
@brianavery 进行了一些编辑 - 这就是您要找的吗？
有帮助，谢谢，但是您使用哪个软件包用于 gam -- gam 或 mgcv？
@brianavery 我正在使用mgcv

标签： r subset chi-squared

【解决方案1】：

我发现 broom 包中的 tidy() 函数在这种情况下很有帮助。
这是您的模型（我将迭代次数减少到 10 次以使其运行更快，
它还使用mgcv::gam() 发出警告）

a <- c(1:10)
model <- list()

for (i in a) {
  model[[i]] <- gam(Species~ s(Petal.Width, bs="cs") + 
                      s(Sepal.Length, bs="cs"),
                    data = iris,
                    family = binomial,
                    method = "REML")
}

如果您需要或想要所有模型的列表，那么第二个循环将从每个模型中提取 statistic 并将它们组合成一个数据框：

results_df <- data.frame(term=c("Petal.Width","Sepal.Length"))
for(m in model){
  results_df <- bind_cols(results_df, model=tidy(m)$statistic)
}
results_df

带输出：

          term    model...2    model...3    model...4    model...5    model...6
1  Petal.Width 1.304372e-09 1.304372e-09 1.304372e-09 1.304372e-09 1.304372e-09
2 Sepal.Length 3.844738e-16 3.844738e-16 3.844738e-16 3.844738e-16 3.844738e-16
     model...7    model...8    model...9   model...10   model...11
1 1.304372e-09 1.304372e-09 1.304372e-09 1.304372e-09 1.304372e-09
2 3.844738e-16 3.844738e-16 3.844738e-16 3.844738e-16 3.844738e-16

您可以重命名列以使其更好，等等。
如果您只想保留这些，您也可以将建模和统计数据结合起来。

【讨论】：

哎呀 - 当我简化这个问题的代码时，我减少了变量的数量，所以忽略之前的评论（现已删除）。
但是，我现在有 results_df 数据框，但其中的统计数据与我手动提取 summary(model[[1]]) 时获得的 Chi.sq 值不匹配

【解决方案2】：

使用来自@brian avery 和this question 的信息，我想出了以下内容，可以准确回答我最初的问题。我确信有一种更有效的方法可以使用管道来完成其中的一些工作，但这对我有用。

a <- c(1:200)
model <- list()
sum <- list ()
Chisq <- list()


for (i in a) {
  model[[i]] <- gam(Species~ s(Petal.Width, bs="cs") + 
                      s(Sepal.Length, bs="cs"),
                    data = iris,
                    family = binomial,
                    method = "REML")

  sum[[i]] <- summary(model[[i]])
  Chisq[[i]] <- sum[[i]]$chi.sq

}


results_df <- data.frame(term=c("Petal.Width","Sepal.Length"))

for(i in a){
  results_df <- bind_cols(results_df, Chisq[[i]])
}

rowSums(results_df<0.001)

最后一行计算了变量最终 Chi.sq 值低于 0.001 的次数。它回答了mgcv 包自动将该变量从模型中缩小了多少次（共 200 次）的问题。

【讨论】：