``` r
library(tidyverse)

gss_clean <- gss_cat %>%
  filter(race != "Other", marital != "No answer") %>%
  mutate(relig = recode(relig, Protestant = "Protestant", .default = "Not Protestant"),
         marital = recode(marital, Married = "Married", .default = "Not married"),
         race = fct_drop(race))

gss_all_3 <- gss_clean %>%
  group_by(marital, race, relig) %>%
  summarize(avg_age = mean(age, na.rm = TRUE),
            sd_age = sd(age, na.rm = TRUE),
            N = n())

gss_race_marriage <- gss_clean %>%
  group_by(marital, race) %>%
  summarize(avg_age = mean(age, na.rm = TRUE),
            sd_age = sd(age, na.rm = TRUE),
            N = n())

gss_marriage <- gss_clean %>%
  group_by(marital) %>%
  summarize(avg_age = mean(age, na.rm = TRUE),
            sd_age = sd(age, na.rm = TRUE),
            N = n())

bind_rows(gss_all_3, gss_race_marriage, gss_marriage) %>%
  ungroup() %>%
  arrange(marital, race, relig)

#> # A tibble: 14 x 6
#>    marital     race  relig          avg_age sd_age     N
#>    <fct>       <fct> <fct>            <dbl>  <dbl> <int>
#>  1 Not married Black Not Protestant    38.6   15.5   672
#>  2 Not married Black Protestant        44.8   17.1  1586
#>  3 Not married Black <NA>              42.9   16.9  2258
#>  4 Not married White Not Protestant    44.4   18.7  4389
#>  5 Not married White Protestant        51.7   19.8  3677
#>  6 Not married White <NA>              47.7   19.5  8066
#>  7 Not married <NA>  <NA>              46.7   19.1 10324
#>  8 Married     Black Not Protestant    42.8   11.6   185
#>  9 Married     Black Protestant        47.4   13.7   684
#> 10 Married     Black <NA>              46.4   13.4   869
#> 11 Married     White Not Protestant    47.8   14.6  3810
#> 12 Married     White Protestant        51.3   15.6  4506
#> 13 Married     White <NA>              49.7   15.2  8316
#> 14 Married     <NA>  <NA>              49.4   15.1  9185
```