r - 带有 2 向方差分析的 p 值的汇总表(平均值 + 标准误差)

Question

我正在尝试制作一个表格，用于输出我们通常通过 2-way anova 分析的大型研究的汇总统计数据，查看两个变量的主要影响以及相互作用。
我想要一种快速运行统计数据的方法，并以易于阅读的格式输出它们，如果格式很好，那就更好了。
我已经能够获得 2 路 anova 输出，而且我还使用了 gtsummary 包和 tbl_summary做一张 table 。但是，我无法弄清楚如何按 1 个以上的变量进行分组。我的解决方案是创建一个结合两个自变量的新变量，只是为了将数据分成正确的组。
下面的可重现示例。
我想知道是否有一种方法可以用我拥有的均值(sem)输出制作表格，但要获得我的 2 路方差分析的结果(也在下面粘贴)。在这个巨大的例子中，我想要一列表示“性别”的主要影响的 P 值，下一列表示“登船”的主要影响的 p 值，然后是交互的 p 值。
有什么想法吗？

library(titanic)
library(tidyverse)
library(gtsummary)
library(plotrix) #has a std.error function


##I really want to look at a 2-way anova, looking for the p-value for Sex, Embarked, and their interaction.
#This code just allows me to make a table with the 4 columns I want, but of course it now won't do the correct stats.
df <- titanic_train %>%
  filter(Embarked != "C" &  Embarked != "") %>%
  mutate(grp = paste(Sex, Embarked)) #add a new column that combines Sex & Pclass

#code to make my table 
  
table1 <- df %>%  
  select(grp, Age, Fare, Survived) %>%
  tbl_summary(
    by = grp,  ##can't figure out a way to put 2 variables here (Sex & Embarked)
    missing = "ifany", 
    statistic = all_continuous() ~ "{mean} ({std.error})",
    digits = all_continuous() ~ 1) %>% #this puts 1 decimal place for all values
   modify_header(stat_by = md("**{level}**<br>N =  {n}")) %>%
  bold_labels() %>%
  modify_spanning_header(all_stat_cols() ~ "**These are the Columns I Want**") %>%
  add_p(test = everything() ~ "aov",  ##This is a 1-way ANOVA, but I need 2 variables
  )

table1

#these are the p-values I want in my table:
two_way_anova_age <- aov(Age ~ Sex * Embarked, data = df)
summary(two_way_anova_age)

two_way_anova_fare <- aov(Fare ~ Sex * Embarked, data = df)
summary(two_way_anova_fare)

two_way_anova_surv <- aov(Survived ~ Sex * Embarked, data = df)
summary(two_way_anova_surv)

Answer 1

以下是将结果合并到 gtsummary 表中的方法。

library(gtsummary)
library(titanic)
library(tidyverse)
library(plotrix) #has a std.error function

packageVersion("gtsummary")
#> [1] '1.4.0'

# create smaller version of the dataset
df <- 
  titanic_train %>%
  select(Sex, Embarked, Age, Fare) %>%
  filter(Embarked != "") # deleting empty Embarked status

# first, write a little function to get the 2-way ANOVA p-values in a table
# function to get 2-way ANOVA p-values in tibble
twoway_p <- function(variable) {
  paste(variable, "~ Sex * Embarked") %>%
    as.formula() %>%
    aov(data = df) %>% 
    broom::tidy() %>%
    select(term, p.value) %>%
    filter(complete.cases(.)) %>%
    pivot_wider(names_from = term, values_from = p.value) %>%
    mutate(
      variable = .env$variable,
      row_type = "label"
    )
}

# add all results to a single table (will be merged with gtsummary table in next step)
twoway_results <-
  bind_rows(
    twoway_p("Age"),
    twoway_p("Fare")
  )
twoway_results
#> # A tibble: 2 x 5
#>            Sex Embarked `Sex:Embarked` variable row_type
#>          <dbl>    <dbl>          <dbl> <chr>    <chr>   
#> 1 0.00823      3.97e- 1         0.611  Age      label   
#> 2 0.0000000191 4.27e-16         0.0958 Fare     label


tbl <-
  # first build a stratified `tbl_summary()` table to get summary stats by two variables
  df %>%
  tbl_strata(
    strata =  Sex,
    .tbl_fun =
      ~.x %>%
      tbl_summary(
        by = Embarked,
        missing = "no",
        statistic = all_continuous() ~ "{mean} ({std.error})",
        digits = everything() ~ 1
      ) %>%
      modify_header(all_stat_cols() ~ "**{level}**")
  ) %>%
  # merge the 2way ANOVA results into tbl_summary table
  modify_table_body(
    ~.x %>%
      left_join(
        twoway_results,
        by = c("variable", "row_type")
      )
  ) %>%
  # by default the new columns are hidden, add a header to unhide them
  modify_header(list(
    Sex ~ "**Sex**", 
    Embarked ~ "**Embarked**", 
    `Sex:Embarked` ~ "**Sex * Embarked**"
  )) %>%
  # adding spanning header to analysis results
  modify_spanning_header(c(Sex, Embarked, `Sex:Embarked`) ~ "**Two-way ANOVA p-values**") %>%
  # format the p-values with a pvalue formatting function
  modify_fmt_fun(c(Sex, Embarked, `Sex:Embarked`) ~ style_pvalue) %>%
  # update the footnote to be nicer looking
  modify_footnote(all_stat_cols() ~ "Mean (SE)")

创建于 2021-03-27 由 reprex package (v1.0.0)