apa_print for emmeans

[andreifoldes]

This is just a feature request, but I think there is great need for apa_print for the emmeans package.

[crsh]

Hi Andrei, absolutely! I have already started to implement support for emmeans but it currently does not work for all possible cases and is not covered by unit tests. I therefore have not released it on the master branch. If you install the development version (devtools::install_github("crsh/papaja@devel")) you can try it out; bug reports and suggestions are always welcome!

[crsh]

I'm currently focusing on implementing support for ANOVA models. As emmeans supports a lot of models, support for the output can only be added incrementally. I think the methods now support all follow-up tests for individual cells, pairwise comparisons, and almost all forms of custom contrasts.

I'll try to track what's left to do for ANOVA here:

• Add support for simple contrasts (see simple argument)
• Check out what broom can do for us here
• Add support for joint_tests() output
• Indicate families of tests in result table
• Indicate which p-values have been adjusted for multiple comparisons

[crsh]

I think it would be sensible to use broom here but the emmeans support is currently somewhat limited. I have submitted a PR to broom and have another one cued. So this issue is deferred until broom extends its emmeans-support. :/

[shirdekel]

What's the update on this? It seems that those PRs have been merged in broom.

My specific request (let me know if you'd prefer this in a new issue) is to be able to use apa_print() with emtrends(), as using it with emmeans() currently works fine (except for the CI NAs, I guess):

  library(papaja)
#> Loading required package: tinylabels
  library(emmeans)

  fiber.lm <- lm(strength ~ diameter * machine, data = fiber)

  emmeans_result <- emmeans(fiber.lm, ~machine)
#> NOTE: Results may be misleading due to involvement in interactions

  apa_print(emmeans_result)
#> $estimate
#> $estimate$A
#> [1] "$M = 40.22$, 95\\% CI NA"
#> 
#> $estimate$B
#> [1] "$M = 41.60$, 95\\% CI NA"
#> 
#> $estimate$C
#> [1] "$M = 38.53$, 95\\% CI NA"
#> 
#> 
#> $statistic
#> $statistic$A
#> [1] "$t(9) = 51.76$, $p < .001$"
#> 
#> $statistic$B
#> [1] "$t(9) = 48.50$, $p < .001$"
#> 
#> $statistic$C
#> [1] "$t(9) = 39.88$, $p < .001$"
#> 
#> 
#> $full_result
#> $full_result$A
#> [1] "$M = 40.22$, 95\\% CI NA, $t(9) = 51.76$, $p < .001$"
#> 
#> $full_result$B
#> [1] "$M = 41.60$, 95\\% CI NA, $t(9) = 48.50$, $p < .001$"
#> 
#> $full_result$C
#> [1] "$M = 38.53$, 95\\% CI NA, $t(9) = 39.88$, $p < .001$"
#> 
#> 
#> $table
#> A data.frame with 5 labelled columns:
#> 
#>   machine estimate       conf.int statistic p.value
#> A       A    40.22 [38.46, 41.98]     51.76  < .001
#> B       B    41.60 [39.66, 43.54]     48.50  < .001
#> C       C    38.53 [36.35, 40.72]     39.88  < .001
#> 
#> machine  : machine 
#> estimate : $M$ 
#> conf.int : 95\\% CI 
#> statistic: $t(9)$ 
#> p.value  : $p$ 
#> attr(,"class")
#> [1] "apa_results" "list"

  emtrends_result <- emtrends(fiber.lm, ~machine, var = "diameter")

  apa_print(emtrends_result)
#> Error in printnum.default(contrast_table$estimate, ...): The parameter 'x' is NULL. Please provide a value for 'x'

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

Session info

sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 4.0.2 (2020-06-22)
#>  os       macOS Mojave 10.14.6        
#>  system   x86_64, darwin17.0          
#>  ui       X11                         
#>  language (EN)                        
#>  collate  en_GB.UTF-8                 
#>  ctype    en_GB.UTF-8                 
#>  tz       Australia/Sydney            
#>  date     2021-03-11                  
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package      * version     date       lib source                            
#>  assertthat     0.2.1       2019-03-21 [1] CRAN (R 4.0.0)                    
#>  backports      1.2.1       2020-12-09 [1] CRAN (R 4.0.2)                    
#>  cli            2.3.1       2021-02-23 [1] CRAN (R 4.0.2)                    
#>  coda           0.19-4      2020-09-30 [1] CRAN (R 4.0.2)                    
#>  crayon         1.4.1       2021-02-08 [1] CRAN (R 4.0.2)                    
#>  debugme        1.1.0       2017-10-22 [1] CRAN (R 4.0.2)                    
#>  digest         0.6.27      2020-10-24 [1] CRAN (R 4.0.2)                    
#>  ellipsis       0.3.1       2020-05-15 [1] CRAN (R 4.0.0)                    
#>  emmeans      * 1.5.4-09004 2021-03-11 [1] Github (rvlenth/emmeans@7815617)  
#>  estimability   1.3         2018-02-11 [1] CRAN (R 4.0.0)                    
#>  evaluate       0.14        2019-05-28 [1] CRAN (R 4.0.0)                    
#>  fansi          0.4.2       2021-01-15 [1] CRAN (R 4.0.2)                    
#>  fs             1.5.0       2020-07-31 [1] CRAN (R 4.0.2)                    
#>  glue           1.4.2       2020-08-27 [1] CRAN (R 4.0.2)                    
#>  highr          0.8         2019-03-20 [1] CRAN (R 4.0.0)                    
#>  htmltools      0.5.1.1     2021-01-22 [1] CRAN (R 4.0.2)                    
#>  knitr          1.31        2021-01-27 [1] CRAN (R 4.0.2)                    
#>  lattice        0.20-41     2020-04-02 [1] CRAN (R 4.0.2)                    
#>  lifecycle      1.0.0       2021-02-15 [1] CRAN (R 4.0.2)                    
#>  magrittr       2.0.1       2020-11-17 [1] CRAN (R 4.0.2)                    
#>  mvtnorm        1.1-1       2020-06-09 [1] CRAN (R 4.0.2)                    
#>  papaja       * 0.1.0.9997  2021-03-11 [1] Github (crsh/papaja@d3b39cb)      
#>  pillar         1.5.1       2021-03-05 [1] CRAN (R 4.0.2)                    
#>  pkgconfig      2.0.3       2019-09-22 [1] CRAN (R 4.0.0)                    
#>  plyr           1.8.6       2020-03-03 [1] CRAN (R 4.0.0)                    
#>  purrr          0.3.4       2020-04-17 [1] CRAN (R 4.0.0)                    
#>  Rcpp           1.0.6       2021-01-15 [1] CRAN (R 4.0.2)                    
#>  reprex         1.0.0       2021-01-27 [1] CRAN (R 4.0.2)                    
#>  rlang          0.4.10      2020-12-30 [1] CRAN (R 4.0.2)                    
#>  rmarkdown      2.7.3       2021-03-11 [1] Github (rstudio/rmarkdown@61db7a9)
#>  sessioninfo    1.1.1       2018-11-05 [1] CRAN (R 4.0.0)                    
#>  stringi        1.5.3       2020-09-09 [1] CRAN (R 4.0.2)                    
#>  stringr        1.4.0       2019-02-10 [1] CRAN (R 4.0.0)                    
#>  styler         1.3.2       2020-02-23 [1] CRAN (R 4.0.0)                    
#>  tibble         3.1.0       2021-02-25 [1] CRAN (R 4.0.2)                    
#>  tinylabels   * 0.2.0       2021-02-06 [1] CRAN (R 4.0.2)                    
#>  utf8           1.1.4       2018-05-24 [1] CRAN (R 4.0.0)                    
#>  vctrs          0.3.6       2020-12-17 [1] CRAN (R 4.0.2)                    
#>  withr          2.4.1       2021-01-26 [1] CRAN (R 4.0.2)                    
#>  xfun           0.21        2021-02-10 [1] CRAN (R 4.0.2)                    
#>  xtable         1.8-4       2019-04-21 [1] CRAN (R 4.0.0)                    
#>  yaml           2.2.1       2020-02-01 [1] CRAN (R 4.0.2)                    
#> 
#> [1] /Library/Frameworks/R.framework/Versions/4.0/Resources/library

[Misch007]

Hey!

I got two bugs, which I think are independent from one another, but not sure!

1. The first one is already visible in the post above mine, but also is reproducible with my data (which you already know):
   When you ask for the $emmeans output, (not the $contrasts, this one works fine despite bug number 2), some random p-values (<.001) are returned, even though it should be purely descriptive and also is not returned by the emmeans() function itself (By the way, note the $t($\\infty$)$ issue I reported here Incorrect estimate label in apa_print() for emmeans interaction contrast #456)

participant<-c("vp1","vp2","vp3","vp4","vp5")
group<-c("intervention", "control","control","intervention","control")

df<-data.frame(participant,group) %>% 
  group_by(participant,group) %>% 
  summarise(session=c("t1","t2","t3","t4")) %>% 
  group_by(participant, group, session) %>% 
  summarise(task=c("a","b")) %>% 
  ungroup() %>% 
  mutate(errors=floor(runif(n=40,min=0,max=30)))

glmm <- mixed(errors~group*session*task+(1|participant), df, family='poisson', method='LRT')

emm.glmm<-emmeans(glmm, pairwise~session)

apa_print(emm.glmm$emmeans)$table 

$table
A data.frame with 5 labelled columns:

   session estimate     conf.int statistic p.value
t1      t1     2.42 [2.19, 2.66]     19.98  < .001
t2      t2     2.54 [2.30, 2.77]     21.19  < .001
t3      t3     2.29 [2.00, 2.59]     15.34  < .001
t4      t4     2.45 [2.21, 2.69]     20.12  < .001

session  : session 
estimate : $M$ 
conf.int : 95\\% CI 
statistic: $t($\\infty$)$ 
p.value  : $p$ 
attr(,"class")
[1] "apa_results" "list" 

2. The second one occurs when I subset the type argument in the emmeans() call, here I am getting different error messages, depending on whether I am using glmer or lmer, but for both subsets $emmeans and $contrasts:

a)

emm.glmm<-emmeans(glmm, pairwise~session, type="response")

apa_print(emm.glmm$emmeans)

Error in printnum.default(contrast_table$estimate, ...) : 
  The parameter 'x' is NULL. Please provide a value for 'x'

b)

lmm <- mixed(errors~group*session*task+(1|participant), df)

emm.lmm<-emmeans(lmm, pairwise~session, type="response")

apa_print(emm.lmm$emmeans)

Error in if (x@misc$inv.lbl == "ratio" && !is.null(x@misc$predict.type) &&  : 
  missing value where TRUE/FALSE needed

I hope this helps and thanks a lot for your work!
Mischa

[Misch007]

Hey @shirdekel,

I found a kinda work-around to use apa_print() with emtrends() for the time being, maybe this will be helpful for you as it was for me:

Before passing the result of emtrends() to apa_print() you have to use the function summary(). Its $contrasts output is then sort of recognized by apa_print() with the confidence intervals and the label for the estimate missing, which you can retrieve with the confint() call and rename using variable_label(). It does not work for the $emtrends output, though.

participant<-c("vp1","vp2","vp3","vp4","vp5")
group<-c("intervention", "control","control","intervention","control")

df<-data.frame(participant,group) %>% 
  group_by(participant,group) %>% 
  summarise(session=c("t1","t2","t3","t4")) %>% 
  ungroup() %>% 
  mutate(errors=floor(runif(n=20,min=0,max=30)),
         variable=floor(runif(n=20,min=-100,max=100)))

lmm <- lmer(errors~group*session*variable+(1|participant), df)

contrasts<-emtrends(lmm, pairwise~session|group, 'variable')$contrasts

apa.contrasts<-apa_print(summary(contrasts))
variable_label(apa.contrasts$table$estimate)<-"my label"
confint(contrasts)

Looking forwards for the professionals solving this one ;)
Mischa

[crsh]

Hi Mischa, thanks for sharing your workaround and for your patience. I think I will find time to tackle this issue soonish; it's on the very top of my priority list.