R 中 GAM 的卡方值子集结果

Question

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

for 循环和 GAMs 运行很好，它们将模型存储在适当的列表中，model[[i]]，每个列表元素代表某个 set.seed 的模型迭代。

我可以运行 summary() 在单个列表元素上（例如 model[[5]]）并得到如下内容：

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 个这样的元素，我想知道是否有快速的方法来计算 [=29] 有多少次（总共 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")

}

Answer 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

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

Answer 2

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

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 次）的问题。

R 中 GAM 的卡方值子集结果

Subset by Chi Sq Values from GAM results in R

r

subset

chi-squared