比较模型系数的百分比变化答案

【问题标题】：Comparing percent change of model coefficients比较模型系数的百分比变化
【发布时间】：2018-06-01 14:42:29
【问题描述】：

我正在完成来自 Hosmer-Lemeshow 的有目的的模型构建的第 3 步，它建议比较完整模型 [Iris.mod1] 和简化模型 [Iris.mod2] 之间系数的百分比变化。如果可能的话，我想自动化这一步。

现在我有以下代码：

#Make species a binomial DV
iris = subset(iris, iris$Species != 'virginica')
iris$Species = as.numeric(ifelse(iris$Species == 'setosa', 1, 0))

#Build models
Iris.mod1 = glm(Species~Sepal.Length+Sepal.Width+Petal.Length+Petal.Width, 
data = iris, family = binomial())

Iris.mod2 = glm(Species~Sepal.Length+Petal.Length, data = iris, family = 
binomial())

我实际使用的数据集有大约 93 个变量和 170 万行。但我只是在这个例子中使用了虹膜数据。

#Try to see if any coefficients changed by > 20%
paste(names(which((summary(Iris.mod1)$coefficients[2: 
(nrow(summary(Iris.mod1)$coefficients)),1] - 
(summary(Iris.mod2)$coefficients[2: 
(nrow(summary(Iris.mod2)$coefficients)),1]/
(summary(Iris.mod1)$coefficients[2:nrow(summary(Iris.mod1)$coefficients)),1] 
> 0.2 == TRUE)))))

但是，这段代码充满了错误，我迷失在括号的海洋中。

有没有一种有效的方法来确定哪些变量的系数变化超过 20%？

提前谢谢你。

【问题讨论】：

很高兴看到改进 - 我会在接受答案的建议中添加一个警告。这样做很重要，但不要太快。人们通常建议等待 12-24 小时，看看是否有其他更好的答案出现。（除了例外，如果答案看起来很完美，请立即接受。）

标签： r logistic-regression

【解决方案1】：

broom 包非常适合制作模型系数和项的数据框。我们可以使用它来获得可行的格式：

library(broom)
m_list = list(m1 = Iris.mod1, m2 = Iris.mod2)
t_list = lapply(m_list, tidy)

library(dplyr)
library(tidyr)
bind_rows(t_list, .id = "mod") %>%
    select(term, estimate, mod) %>%
    spread(key = mod, value = estimate) %>%
    mutate(p_change = (m2 - m1) / m1 * 100,
           p_change_gt_20 = p_change > 20)
#           term         m1        m2 p_change p_change_gt_20
# 1  (Intercept)  -6.556265 -65.84266 904.2709           TRUE
# 2 Petal.Length -19.053588 -49.04616 157.4117           TRUE
# 3  Petal.Width -25.032928        NA       NA             NA
# 4 Sepal.Length   9.878866  37.56141 280.2199           TRUE
# 5  Sepal.Width   7.417640        NA       NA             NA

【讨论】：