Create ANOVA summary tables for colleagues

Jason.C · April 23, 2020, 9:45pm

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))
           diff    lwr     upr p adj
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")  %>% 
	tab_header(
		title = md("**BDI**"),
		subtitle = md("ANOVA PostHoc")) %>% 
	tab_stubhead(label = "Comparisons")

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.

Any suggestions, advice, or leads would be greatly appreciated.

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

FJCC · April 24, 2020, 2:15am

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")
head(Models)
#> # 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)}

Jason.C · April 24, 2020, 2:11pm

Thank you!
I really appreciate your time and will try your suggestions and reply later today.

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?

system · May 15, 2020, 2:11pm

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