Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Degree of Freedom is different between rstatix and afex packages #21

Open
jooyoungseo opened this issue Dec 27, 2019 · 2 comments
Open

Degree of Freedom is different between rstatix and afex packages #21

jooyoungseo opened this issue Dec 27, 2019 · 2 comments

Comments

@jooyoungseo
Copy link

I have found that df (Degree of Freedom) reported by rstatix package is different from what is returned by afex package.

For example, you could see the difference in the following. F(2, 18) in rstatix package; F(1.38, 12.42) in afex package.

Would you mind addressing this issue using the reproducible code below?:

# Load packages:
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(rstatix))
suppressPackageStartupMessages(library(afex))
library(papaja)

# Prepare data:
data(selfesteem, package = "datarium")

# Transform data into a long format:
df <- selfesteem %>%
  tidyr::gather(key = "time", value = "score", -id) %>%
  rstatix::convert_as_factor(id, time)

# Anova results from `rstatix` package:
res.aov1 <- rstatix::anova_test(data = df, dv = score, wid = id, within = time)
res.aov1
#> ANOVA Table (type III tests)
#> 
#> $ANOVA
#>   Effect DFn DFd      F        p p<.05   ges
#> 1   time   2  18 55.469 2.01e-08     * 0.829
#> 
#> $`Mauchly's Test for Sphericity`
#>   Effect     W     p p<.05
#> 1   time 0.551 0.092      
#> 
#> $`Sphericity Corrections`
#>   Effect  GGe      DF[GG]    p[GG] p[GG]<.05   HFe      DF[HF]    p[HF]
#> 1   time 0.69 1.38, 12.42 2.16e-06         * 0.774 1.55, 13.94 6.03e-07
#>   p[HF]<.05
#> 1         *
rstatix::get_test_label(res.aov1, description = NULL, detailed = TRUE, type = "text")
#> [1] "F(2,18) = 55.47, p = <0.0001, eta2[g] = 0.829"

# Anova results from `afex` package:
res.aov2 <- afex::aov_ez(data = df, dv = "score", id = "id", within = "time")
res.aov2
#> Anova Table (Type 3 tests)
#> 
#> Response: score
#>   Effect          df  MSE         F ges p.value
#> 1   time 1.38, 12.42 1.34 55.47 *** .83  <.0001
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
#> 
#> Sphericity correction method: GG
papaja::apa_print(res.aov2)$full_result$time
#> [1] "$F(1.38, 12.42) = 55.47$, $\\mathit{MSE} = 1.34$, $p < .001$, $\\hat{\\eta}^2_G = .829$"

Created on 2019-12-27 by the reprex package (v0.3.0.9001)

Session info 
sessioninfo::session_info()
#> - Session info ---------------------------------------------------------------
#>  setting  value                       
#>  version  R version 3.6.1 (2019-07-05)
#>  os       Windows 10 x64              
#>  system   x86_64, mingw32             
#>  ui       RTerm                       
#>  language (EN)                        
#>  collate  English_United States.1252  
#>  ctype    English_United States.1252  
#>  tz       America/New_York            
#>  date     2019-12-27                  
#> 
#> - Packages -------------------------------------------------------------------
#>  package      * version     date       lib
#>  abind          1.4-5       2016-07-21 [1]
#>  afex         * 0.25-1      2019-12-26 [1]
#>  assertthat     0.2.1       2019-03-21 [1]
#>  backports      1.1.5       2019-10-02 [1]
#>  boot           1.3-23      2019-07-05 [1]
#>  broom          0.5.3.9000  2019-12-19 [1]
#>  car            3.0-6       2019-12-23 [1]
#>  carData        3.0-3       2019-11-16 [1]
#>  cellranger     1.1.0       2016-07-27 [1]
#>  cli            2.0.0       2019-12-09 [1]
#>  coda           0.19-3      2019-07-05 [1]
#>  codetools      0.2-16      2018-12-24 [1]
#>  colorspace     1.4-1       2019-03-18 [1]
#>  crayon         1.3.4       2017-09-16 [1]
#>  curl           4.3         2019-12-02 [1]
#>  data.table     1.12.8      2019-12-09 [1]
#>  DBI            1.0.0       2018-05-02 [1]
#>  dbplyr         1.4.2       2019-11-08 [1]
#>  digest         0.6.23      2019-11-23 [1]
#>  dplyr        * 0.8.3       2019-07-04 [1]
#>  ellipsis       0.3.0       2019-09-20 [1]
#>  emmeans        1.4.3.01    2019-11-28 [1]
#>  estimability   1.3         2018-02-11 [1]
#>  evaluate       0.14        2019-05-28 [1]
#>  fansi          0.4.0       2018-10-05 [1]
#>  forcats      * 0.4.0       2019-02-17 [1]
#>  foreign        0.8-72      2019-08-02 [1]
#>  fs             1.3.1       2019-05-06 [1]
#>  generics       0.0.2       2018-11-29 [1]
#>  ggplot2      * 3.2.1.9000  2019-12-23 [1]
#>  glue           1.3.1.9000  2019-12-19 [1]
#>  gtable         0.3.0       2019-03-25 [1]
#>  haven          2.2.0.9000  2019-11-08 [1]
#>  highr          0.8         2019-03-20 [1]
#>  hms            0.5.2       2019-10-30 [1]
#>  htmltools      0.4.0       2019-10-04 [1]
#>  httr           1.4.1       2019-08-05 [1]
#>  jsonlite       1.6         2018-12-07 [1]
#>  knitr          1.26.1      2019-12-05 [1]
#>  lattice        0.20-38     2018-11-04 [1]
#>  lifecycle      0.1.0       2019-08-01 [1]
#>  lme4         * 1.1-21      2019-03-05 [1]
#>  lmerTest       3.1-1       2019-12-13 [1]
#>  lubridate      1.7.4       2018-04-11 [1]
#>  magrittr       1.5         2014-11-22 [1]
#>  MASS           7.3-51.4    2019-03-31 [1]
#>  Matrix       * 1.2-18      2019-11-27 [1]
#>  minqa          1.2.4       2014-10-09 [1]
#>  modelr         0.1.5       2019-08-08 [1]
#>  multcomp       1.4-10      2019-03-05 [1]
#>  munsell        0.5.0       2018-06-12 [1]
#>  mvtnorm        1.0-11      2019-06-19 [1]
#>  nlme           3.1-142     2019-11-07 [1]
#>  nloptr         1.2.1       2018-10-03 [1]
#>  numDeriv       2016.8-1.1  2019-06-06 [1]
#>  openxlsx       4.1.4       2019-12-06 [1]
#>  papaja       * 0.1.0.9842  2019-12-19 [1]
#>  pillar         1.4.3       2019-12-20 [1]
#>  pkgconfig      2.0.3       2019-09-22 [1]
#>  plyr           1.8.5       2019-12-10 [1]
#>  purrr        * 0.3.3       2019-10-18 [1]
#>  R6             2.4.1       2019-11-12 [1]
#>  Rcpp           1.0.3       2019-11-08 [1]
#>  readr        * 1.3.1.9000  2019-12-23 [1]
#>  readxl         1.3.1.9000  2019-12-17 [1]
#>  reprex         0.3.0.9001  2019-11-13 [1]
#>  reshape2       1.4.3       2017-12-11 [1]
#>  rio            0.5.16      2018-11-26 [1]
#>  rlang          0.4.2.9000  2019-12-26 [1]
#>  rmarkdown      2.0.5       2019-12-20 [1]
#>  rstatix      * 0.3.1.999   2019-12-23 [1]
#>  rvest          0.3.5.9000  2019-11-09 [1]
#>  sandwich       2.5-1       2019-04-06 [1]
#>  scales         1.1.0       2019-11-18 [1]
#>  sessioninfo    1.1.1       2018-11-05 [1]
#>  stringi        1.4.3       2019-03-12 [1]
#>  stringr      * 1.4.0.9000  2019-11-11 [1]
#>  styler         1.2.0.9000  2019-11-13 [1]
#>  survival       3.1-8       2019-12-03 [1]
#>  TH.data        1.0-10      2019-01-21 [1]
#>  tibble       * 2.1.3       2019-06-06 [1]
#>  tidyr        * 1.0.0.9000  2019-12-09 [1]
#>  tidyselect     0.2.99.9000 2019-12-19 [1]
#>  tidyverse    * 1.3.0.9000  2019-12-17 [1]
#>  vctrs          0.2.99.9000 2019-12-26 [1]
#>  withr          2.1.2       2018-03-15 [1]
#>  xfun           0.11.3      2019-11-22 [1]
#>  xml2           1.2.2.9000  2019-12-19 [1]
#>  xtable         1.8-4       2019-04-21 [1]
#>  yaml           2.2.0       2018-07-25 [1]
#>  zip            2.0.4       2019-09-01 [1]
#>  zoo            1.8-6       2019-05-28 [1]
#>  source                               
#>  CRAN (R 3.6.0)                       
#>  Github (singmann/afex@23af68c)       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  Github (tidymodels/broom@3c922d5)    
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.0)                       
#>  Github (hadley/dbplyr@99d4db4)       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.2)                       
#>  CRAN (R 3.6.0)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  Github (hadley/ggplot2@d05e437)      
#>  Github (tidyverse/glue@b9ffe6c)      
#>  CRAN (R 3.6.1)                       
#>  Github (tidyverse/haven@690f9d7)     
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  Github (yihui/knitr@33d69c3)         
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.0)                       
#>  CRAN (R 3.6.2)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.0)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.0)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.0)                       
#>  CRAN (R 3.6.0)                       
#>  CRAN (R 3.6.1)                       
#>  Github (crsh/papaja@081b13f)         
#>  CRAN (R 3.6.2)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  Github (hadley/readr@4bb6ab2)        
#>  Github (hadley/readxl@641ce4e)       
#>  Github (tidyverse/reprex@27aa69a)    
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.0)                       
#>  Github (r-lib/rlang@ce4f717)         
#>  Github (rstudio/rmarkdown@7f51e23)   
#>  Github (kassambara/rstatix@7c03d07)  
#>  Github (hadley/rvest@a5200f2)        
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.0)                       
#>  Github (hadley/stringr@80aaaac)      
#>  Github (r-lib/styler@a8acde5)        
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.0)                       
#>  CRAN (R 3.6.1)                       
#>  Github (hadley/tidyr@885d1af)        
#>  Github (tidyverse/tidyselect@6b8f176)
#>  Github (tidyverse/tidyverse@1d7f9b7) 
#>  Github (r-lib/vctrs@fa48942)         
#>  CRAN (R 3.6.1)                       
#>  Github (yihui/xfun@4fcd6ff)          
#>  Github (r-lib/xml2@9fa8a7b)          
#>  CRAN (R 3.6.0)                       
#>  CRAN (R 3.6.0)                       
#>  CRAN (R 3.6.1)                       
#>  CRAN (R 3.6.1)                       
#> 
#> [1] C:/Program Files/R/R-3.6.1/library
@kassambara
Copy link
Owner

There is no difference between rstatix and afex in computing the repeated measures ANOVA. However, I agree with you that there is a small difference between the two packages when printing the results.

Let's clarify that...

As you know, one of the assumption of the repeated measure ANOVA is the assumption of sphericity: the variance of the differences between groups should be equal. This can be checked using the Mauchly’s test of sphericity, which is automatically reported when using the R function rstatix::anova_test() . Read more at Repeated Measures ANOVA in R using the rstatix package.

As you can see in the rstatix output, the Mauchly's Test for Sphericity is not significant (p = 0.09), so we can assume the sphericity and there is no need to apply corrections to the degrees of freedom.

$`Mauchly's Test for Sphericity`
  Effect     W     p p<.05
1   time 0.551 0.092

Note also that, the Sphericity corrections results are available in the rstatix output and should be considered in case we could not maintain the sphericity assumption. Two common corrections used in the literature are provided: Greenhouse-Geisser epsilon (GGe), and Huynh-Feldt epsilon (HFe) and their corresponding p-values. You cand see that the df for GG correction is DF[GG] = 1.38, 12.42 like in the afex output.

$`Sphericity Corrections`
  Effect  GGe      DF[GG]    p[GG] p[GG]<.05   HFe      DF[HF]    p[HF] p[HF]<.05
1   time 0.69 1.38, 12.42 2.16e-06         * 0.774 1.55, 13.94 6.03e-07         *

The default behaviour of the afex package is to systematically apply the GG correction when printing the repeated measures ANOVA output no mather whether the Mauchly's test is significant or not.

In rstatix, the correction is automatically applied only when the sphericity can't be assumed. We follow the same philosophy as the popular commercial software SPSS by printing detailed results. If users want to force the sphericity correction in the rstatix package, they can go as follow:

library(rstatix)

# Data preparation
data(selfesteem, package = "datarium")
df <- selfesteem %>%
  tidyr::gather(key = "time", value = "score", -id) %>%
  rstatix::convert_as_factor(id, time)

# Anova results from `rstatix` package:
res.aov1 <- rstatix::anova_test(data = df, dv = score, wid = id, within = time)

# Display ANOVA table with GG correction, default is "auto
get_anova_table(res.aov1, correction = "GG")
ANOVA Table (type III tests)

  Effect  DFn  DFd    F        p p<.05   ges
1   time 1.38 12.4 55.5 2.16e-06     * 0.829

# Print test abel
rstatix::get_test_label(
  res.aov1, correction = "GG", 
  description = NULL, detailed = TRUE, type = "text"
  )
[1] "F(1.38,12.42) = 55.47, p = <0.0001, eta2[g] = 0.829"

@jooyoungseo
Copy link
Author

Your explanation is tremendously instrumental, @kassambara!

I really appreciate your help!

Now everything is very clear.

This could be another question; however, would there be any way to print the results with text formattings (e.g., italics and bold)?

Currently, rstatix::get_test_label(..., type = "expression") cannot be used in body text area. type = "text", as you know, prints without any formatting adjustments.

I hope that I can use rstatix::get_test_label() in a similar fashion to papaja::apa_print().

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@kassambara @jooyoungseo