# Create ANOVA summary tables for colleagues

Greetings,
I’ve learned there are several ways to perform an ANOVA with PostHoc testing in R via RStudio. Just the other day, someone helped me construct the following (data sample at end of post):

``````# construct anova model for each column of dat1[, my_outcome_variables]
model_summaries <- dat1[,c("slpQual", "slpLat", "slpDur",
"slpEff", "slpDist","slpMeds",
"slpDayFcn","psqi_Global","ess_total",
"bdi_total", "mcgill_total",
"TIB")] %>%
purrr::map(~ aov(.x ~ group.factor, data = dat1))  %>%
# append the TukeyHSD CI estimates
purrr::map(function(x) {
list(
model = x,
tukey = TukeyHSD(x)
)
})
``````

To view the output I’ve done this:

``````#		TO ACCESS THE RESULTS
> summary(model_summaries\$bdi_total\$model)
Df Sum Sq Mean Sq F value   Pr(>F)
group.factor   2   3894  1947.1   32.84 1.29e-12 ***
Residuals    155   9191    59.3
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> print(round(model_summaries\$bdi_total\$tukey\$group.factor, 4))
CLBP-HC  6.8650 3.4255 10.3045 0e+00
FM-HC   12.8536 9.0869 16.6204 0e+00
FM-CLBP  5.9886 2.4199  9.5573 3e-04
``````

I’ve spent most of today working on a way to put this information into a table format for colleagues and have only partially succeeded with the second set of information. Specifically:

``````z_table1 <- print(signif(model_summaries\$bdi_total\$tukey\$group.factor, 4))
z_table1 <- as.data.frame(z_table1)
z_table1\$row <- c("CLBP-HC", "FM-HC", "FM-CLBP")

z_table1 %>%
gt(rowname_col = "row")  %>%
title = md("**BDI**"),
subtitle = md("ANOVA PostHoc")) %>%
``````

However, after trying to follow the examples of various packages, I cannot work out how to get the initial summary information into a table. Moreover, as I have 12 variables, there must be a more concise way of extracting this information other than creating 12 copies of the code and only changing the variable name in each one.

Jason the #rstatsnewbie

Data sample

``````dput(head(dat.SM, 30))
structure(list(group.factor = structure(c(1L, 1L, 1L, 1L, 1L,
1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
2L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L), .Label = c("HC", "CLBP",
"FM", "TRUE"), class = "factor"), slpQual = c(0, 1, 1, 1, 1,
0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 2, 0, 2, 1,
1, 2, 1, 2), slpLat = c(1, 1, 2, 1, 1, 1, 1, 1, 1, 2, 1, 2, 1,
1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 2, 0, 1, 1, 1, 2), slpDur = c(1,
0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0,
1, 0, 2, 1, 1, 0, 1, 3), slpEff = c(0, 0, 1, 0, 1, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 3
), slpDist = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), slpMeds = c(0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 1,
0, 0, 0, 0, 0, 3, 1), slpDayFcn = c(0, 0, 1, 0, 1, 0, 0, 1, 0,
1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
), psqi_Global = c(3, 3, 7, 3, 5, 4, 3, 4, 2, 6, 1, 4, 4, 3,
4, 1, 4, 5, 6, 4, 2, 2, 6, 1, 8, 3, 4, 5, 8, 12), ess_total = c(5,
5, 4, 7, 10, 2, 5, 5, 4, 9, 6, 10, 8, 9, 6, 8, 9, 3, 6, 10, 8,
10, 3, 0, 10, 6, 6, 6, 9, 6), bdi_total = c(0, 1, 1, 13, 5, 0,
1, 6, 0, 7, 0, 2, 0, 6, 0, 3, 1, 0, 0, 7, 12, 5, 7, 0, 2, 8,
1, 9, 12, 4), mcgill_total = c(0, 0, 0, 0, 0, 0, 9, 0, 0, 4,
0, 1, 0, 5, 0, 0, 0, 0, 0, 38, 34, 16, 0, 0, 0, 5, 2, 14, 0,
0), TIB = c(7, 7.5, 8.5, 8, 9.5, 8.5, 7.5, 8.5, 7.5, 7, 9, 8,
8, 8, 7, 8.5, 8, 8.5, 8.75, 8, 8, 8, 7.5, 8.25, 6, 6, 7.5, 9,
8.5, 9)), row.names = c(NA, 30L), class = "data.frame")
``````

For making a table of the summaries of models, how about starting with this roll up.

``````library(purrr)
dat1 <- structure(list(group.factor =
structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L,
2L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 1L),
.Label = c("HC", "CLBP", "FM", "TRUE"),
class = "factor"),
slpQual = c(0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1,
1, 0, 0, 2, 0, 2, 1, 1, 2, 1, 2),
slpLat = c(1, 1, 2, 1, 1, 1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 0, 1, 1, 1,
1, 0, 1, 1, 0, 2, 0, 1, 1, 1, 2),
slpDur = c(1,0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 0,
1, 1, 0, 1, 0, 2, 1, 1, 0, 1, 3),
slpEff = c(0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 1, 0, 0, 1, 1, 3),
slpDist = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
slpMeds = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,
0, 0, 1, 0, 0, 0, 0, 0, 3, 1),
slpDayFcn = c(0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
psqi_Global = c(3, 3, 7, 3, 5, 4, 3, 4, 2, 6, 1, 4, 4, 3, 4, 1, 4, 5, 6,
4, 2, 2, 6, 1, 8, 3, 4, 5, 8, 12),
ess_total = c(5,5, 4, 7, 10, 2, 5, 5, 4, 9, 6, 10, 8, 9, 6, 8, 9, 3, 6, 10,
8, 10, 3, 0, 10, 6, 6, 6, 9, 6),
bdi_total = c(0, 1, 1, 13, 5, 0, 1, 6, 0, 7, 0, 2, 0, 6, 0, 3, 1, 0, 0, 7,
12, 5, 7, 0, 2, 8, 1, 9, 12, 4),
mcgill_total = c(0, 0, 0, 0, 0, 0, 9, 0, 0, 4, 0, 1, 0, 5, 0, 0, 0, 0, 0,
38, 34, 16, 0, 0, 0, 5, 2, 14, 0, 0),
TIB = c(7, 7.5, 8.5, 8, 9.5, 8.5, 7.5, 8.5, 7.5, 7, 9, 8, 8, 8, 7, 8.5, 8,
8.5, 8.75, 8, 8, 8, 7.5, 8.25, 6, 6, 7.5, 9, 8.5, 9)),
row.names = c(NA, 30L), class = "data.frame")

model_summaries <- dat1[,c("slpQual", "slpLat", "slpDur",
"slpEff", "slpDist","slpMeds",
"slpDayFcn","psqi_Global","ess_total",
"bdi_total", "mcgill_total",
"TIB")] %>%
purrr::map(~ aov(.x ~ group.factor, data = dat1))  %>%
# append the TukeyHSD CI estimates
purrr::map(function(x) {
list(
model = x,
tukey = TukeyHSD(x)
)
})
GetModels <- function(l) l\$model
Models <- purrr::map(model_summaries, GetModels) %>%
purrr::map_dfr(broom::tidy, .id = "Name")
#> # A tibble: 6 x 7
#>   Name    term            df    sumsq  meansq statistic p.value
#>   <chr>   <chr>        <dbl>    <dbl>   <dbl>     <dbl>   <dbl>
#> 1 slpQual group.factor     1  0.327   0.327      0.773    0.387
#> 2 slpQual Residuals       28 11.8     0.423     NA       NA
#> 3 slpLat  group.factor     1  0.327   0.327      1.06     0.312
#> 4 slpLat  Residuals       28  8.64    0.309     NA       NA
#> 5 slpDur  group.factor     1  0.00667 0.00667    0.0125   0.912
#> 6 slpDur  Residuals       28 15.0     0.534     NA       NA
``````

Created on 2020-04-23 by the reprex package (v0.3.0)

Thank you!

Edit 1
This works pretty well in this table:

``````Models %>%
dplyr::select("term" , "df", "p.value")  %>%
gt() %>%
fmt_number(
columns = vars("p.value"),
decimals = 4
)
``````

Your suggestion raised several questions for me.

Q1: What do these actually do? I’ve played with it but can't figure it out.

Q2: How can I remove the row of "Residuals" or at least the "NA" ?

Q3: How can I incorporate the post-hoc tests in this function?

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.