统计模型的漂亮摘要答案

【问题标题】：Pretty summaries for statistical models统计模型的漂亮摘要
【发布时间】：2019-07-25 05:01:47
【问题描述】：

我正在寻找一种漂亮的方式来查看 R 中的统计模型摘要。在以下示例中，我想查看 cyl_6 或 cyl.6 而不是 cyl6。我该怎么做？

library(dplyr)
library(broom)

mean_mpg <- mean(mtcars$mpg)

# creating a new variable that shows that Miles/(US) gallon is greater than the mean or not

mtcars <-
  mtcars %>%
  mutate(mpg_cat = ifelse(mpg > mean_mpg, 1,0))

mtcars$cyl <- as.factor(mtcars$cyl)

model <-
mtcars %>%
  select (cyl,vs, am, mpg_cat) %>%

  glm(formula = mpg_cat ~ .,
      data = ., family = "binomial")


tidy(model)

【问题讨论】：

这与“整洁”或“漂亮”的统计摘要没有任何关系。您只是试图操纵字符串。最简单的方法是使用正则表达式。您应该能够根据您的需要调整this，例如sub("(\\d+)", "_\\1", term)。

标签： r glm broom

【解决方案1】：

我可以想到一种方法来做到这一点，但它非常笨拙：在运行模型之前更改 cyl 的 contrasts 属性（以及您想要包含的任何其他因素）：

mtcars$cyl <- as.factor(mtcars$cyl)
cont = contrasts(mtcars$cyl)
colnames(cont) = paste0("_", colnames(cont))
contrasts(mtcars$cyl) = cont

model <-
    mtcars %>%
    select (cyl,vs, am, mpg_cat) %>%

    glm(formula = mpg_cat ~ .,
        data = ., family = "binomial")

tidy(model)

输出：

# A tibble: 5 x 5
  term        estimate std.error  statistic p.value
  <chr>          <dbl>     <dbl>      <dbl>   <dbl>
1 (Intercept)   22.9      24034.  0.000953    0.999
2 cyl_6        -22.4      12326. -0.00182     0.999
3 cyl_8        -44.5      23246. -0.00191     0.998
4 vs            -1.59     13641. -0.000117    1.000
5 am             0.201    13641.  0.0000147   1.000

如果您希望默认使用此行为，我想您可以编写 contr.treatment 的修改版本，根据需要设置列名，然后使用 options(contrasts = ...) 将其设置为默认值？我还没有测试过它是否有效。

【讨论】：

如何在dplyr中写出这三行代码？
我不确定有没有好的方法 - 我不知道有一种优雅的方法可以将属性设置为dplyr 链的一部分，并且您需要在cyl 专栏。有关想法，请参阅 stackoverflow.com/questions/25662859/…。
谢谢。我的第二个问题是：我怎样才能编写一个函数来为我的因子变量做呢？

【解决方案2】：

例如，在管道中使用sub。
我首先简化模型代码。

model <-
  mtcars %>%
  mutate(mpg_cat = as.integer(mpg > mean(mpg)),
         cyl = factor(cyl)) %>%
  select (cyl,vs, am, mpg_cat) %>%
  glm(formula = mpg_cat ~ .,
      data = ., family = "binomial")

现在是应用正则表达式的问题：

"^cyl" 匹配字符串开头的 "cyl"。

管道是

model %>%
  tidy() %>%
  mutate(term = sub("^cyl", "cyl_", term))
## A tibble: 5 x 5
#  term        estimate std.error  statistic p.value
#  <chr>          <dbl>     <dbl>      <dbl>   <dbl>
#1 (Intercept)   22.9      24034.  0.000953    0.999
#2 cyl_6        -22.4      12326. -0.00182     0.999
#3 cyl_8        -44.5      23246. -0.00191     0.998
#4 vs            -1.59     13641. -0.000117    1.000
#5 am             0.201    13641.  0.0000147   1.000

【讨论】：

有道理，但在我的真实模型中，我有大约 20 个变量级字符串，我不想为 glm 函数的每次调用操作术语列。