Summary Table (mean + std.error) with p-values for 2-way anova

Question

I'm trying to make a table that outputs the summary statistics for a large study that we usually analyze by 2-way anova, looking at main effects of both variables as well as an interaction.

I'd like a way to run the stats quickly, and output them in a format that is easy to read, and if had nice formatting that would be even better.

I've been able to get both the 2-way anova output, and I've also used gtsummary package and tbl_summary to make a table. However, I can't figure out how to group by more than 1 variable. My solution has been to create a new variable that combines the two independent variables, just to split data into the right groups.

Reproducible Example Below.

I'm wondering if there's a way to make a table with the mean(sem) output as I have it, but to have the results of my 2-way anova (also pasted below). In this titanic example, I'd like a column for P-value for main effect of "Sex", next column for p-value main effect of "Embarked", then a p value for the interaction.

Any thoughts?

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

Here is how you can combine the results in a gtsummary table.

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)")

^{Created on 2021-03-27 by the reprex package (v1.0.0)}