library(prmisc)
library(effectsize)
library(afex)
#> Warning in check_dep_version(): ABI version mismatch:
#> lme4 was built with Matrix ABI version 1
#> Current Matrix ABI version is 2
#> Please re-install lme4 from source or restore original 'Matrix' packageprmisc is a small package that formats numbers and statistical test results according to APA style guidelines. The output can be included in Quarto or R Markdown documents. It covers some basic statistical tests (t-test, ANOVA, correlation, \(\chi^2\) test, Mann-Whitney-U test) and some basic number printing manipulations: formatting p-values, removing leading zeros for numbers that cannot be greater than one, and others.
Packages with a more rich set of functionality are papaja and apa. However, we find ourselves returning to prmisc because it offers some convenience functions not available in these packages
prmisc uses the output of the base R t.test() function, and optionally effectsize::cohens_d() to display the results of a t-test:
ttest <- t.test(1:10, y = c(7:20), var.equal = TRUE)
library("effectsize") # for Cohen's d
cohend <- cohens_d(1:10, c(7:20))
print_ttest(ttest, cohend)
#> [1] "$t(22) = -5.15$, $p < .001$, $d = -2.13$"We must include this as “inline code” to properly display it in the output of an R Markdown / Quarto file:
The results were significant, \(t(22) = -5.15\), \(p < .001\), \(d = -2.13\).
An example for paired data:
data(sleep) # ?sleep
tt <- t.test(sleep$extra[sleep$group == 1],
sleep$extra[sleep$group == 2], paired = TRUE)
cd <- cohens_d(sleep$extra[sleep$group == 1],
sleep$extra[sleep$group == 2], paired = TRUE)
#> For paired samples, 'repeated_measures_d()' provides more options.
print_ttest(tt, cd)
#> [1] "$t(9) = -4.06$, $p = .003$, $d_z = -1.28$"\(t(9) = -4.06\), \(p = .003\), \(d_z = -1.28\)
We can also print the confidence interval of the effect size:
print_ttest(tt, cd, confidence = TRUE)
#> [1] "$t(9) = -4.06$, $p = .003$, $d_z = -1.28$, $95\\%$ $CI$ $[-2.12$, $-0.41]$"\(t(9) = -4.06\), \(p = .003\), \(d_z = -1.28\), \(95\%\) \(CI\) \([-2.12\), \(-0.41]\)
The information about the CI is taken from the effectsize object:
cd <- cohens_d(sleep$extra[sleep$group == 1],
sleep$extra[sleep$group == 2], paired = TRUE, ci = .8)
#> For paired samples, 'repeated_measures_d()' provides more options.
print_ttest(tt, cd, confidence = TRUE)
#> [1] "$t(9) = -4.06$, $p = .003$, $d_z = -1.28$, $80\\%$ $CI$ $[-1.81$, $-0.70]$"\(t(9) = -4.06\), \(p = .003\), \(d_z = -1.28\), \(80\%\) \(CI\) \([-1.81\), \(-0.70]\)
The effect size object can also be left out:
\(t(9) = -4.06\), \(p = .003\)
Printing the results of a correlation including the significance test is straight forward. prmisc::print_cortest() uses the output of the base R function cor.test() as input.
x <- c(44.4, 45.9, 41.9, 53.3, 44.7, 44.1, 50.7, 45.2, 60.1)
y <- c( 2.6, 3.1, 2.5, 5.0, 3.6, 4.0, 5.2, 2.8, 3.8)
cor_results <- cor.test(x, y)
print_cortest(cor_results)
#> [1] "$r(7) = .57$, $p = .108$"\(r(7) = .57\), \(p = .108\)
prmisc can format ANOVA results as returned by the afex package.
library("afex")
# see ?aov_ez
data(md_12.1)
aov_results <- aov_ez("id", "rt", md_12.1, within = c("angle", "noise"))
print_anova(aov_results) # returns a list with all effects in this ANOVA
#> $angle
#> [1] "$F(1.92$, $17.31) = 40.72$, $p < .001$, $\\eta_G^2 = .39$"
#>
#> $noise
#> [1] "$F(1$, $9) = 33.77$, $p < .001$, $\\eta_G^2 = .39$"
#>
#> $`angle:noise`
#> [1] "$F(1.81$, $16.27) = 45.31$, $p < .001$, $\\eta_G^2 = .19$"\(F(1.92\), \(17.31) = 40.72\), \(p < .001\), \(\eta_G^2 = .39\)
Via the package papaja (papaja::apa_print(aov_results)), we would obtain a more rich output, including an ANOVA table. However, unlike papaja, prmisc can print the nonitalic \(\eta\), which would be required according to APA guidelines. because \(\eta\) is a greek symbol:
print_anova(aov_results, italic_eta = FALSE)
#> $angle
#> [1] "$F(1.92$, $17.31) = 40.72$, $p < .001$, $\\upeta_\\mathrm{G}^2 = .39$"
#>
#> $noise
#> [1] "$F(1$, $9) = 33.77$, $p < .001$, $\\upeta_\\mathrm{G}^2 = .39$"
#>
#> $`angle:noise`
#> [1] "$F(1.81$, $16.27) = 45.31$, $p < .001$, $\\upeta_\\mathrm{G}^2 = .19$"This only works with Latex/PDF output. To use it in R Markdown documents, we need to load the upgreek package in the YAML header, using
header-includes:
-\usepackageupgreekWe can also use the partial \(\eta^2\), or no effect size index:
print_anova(
aov_ez("id", "rt", md_12.1, within = c("angle", "noise"),
anova_table = list(es = "pes"))
)
#> $angle
#> [1] "$F(1.92$, $17.31) = 40.72$, $p < .001$, $\\eta_p^2 = .82$"
#>
#> $noise
#> [1] "$F(1$, $9) = 33.77$, $p < .001$, $\\eta_p^2 = .79$"
#>
#> $`angle:noise`
#> [1] "$F(1.81$, $16.27) = 45.31$, $p < .001$, $\\eta_p^2 = .83$"
print_anova(
aov_ez("id", "rt", md_12.1, within = c("angle", "noise"),
anova_table = list(es = "none"))
)
#> $angle
#> [1] "$F(1.92$, $17.31) = 40.72$, $p < .001$"
#>
#> $noise
#> [1] "$F(1$, $9) = 33.77$, $p < .001$"
#>
#> $`angle:noise`
#> [1] "$F(1.81$, $16.27) = 45.31$, $p < .001$"Unlike the other functions that display statistical tests, prmisc::print_chi2() can perform the statistical test by itself. By default, prmisc::print_chi2() does not use a continuity correction.
x <- matrix(c(12, 5, 7, 7), ncol = 2)
print_chi2(x) # does not use continuity correction by default
#> [1] "$\\chi^2(1, N = 31) = 1.37$, $p = .242$, $\\phi = .21$"\(\chi^2(1, N = 31) = 1.37\), \(p = .242\), \(\phi = .21\)
print_chi2(x, correct = TRUE) # use continuity correction
#> [1] "$\\chi^2(1, N = 31) = 0.64$, $p = .423$, $\\phi = .14$"\(\chi^2(1, N = 31) = 0.64\), \(p = .423\), \(\phi = .14\)
The \(\phi\) index is displayed as measure of effect size if the input is a 2x2 contingency table. Using the es argument can prevent this:
Like for \(\eta\) in the output of print_anova(), we can also use the nonitalic symbol for \(\chi\) and \(\phi\), which would be required by APA guidelines (because they are greek symbols). Again, this only works for Latex/PDF output and we would need the upgreek package.
prmisc includes several functions to format numbers. We find these particularly useful for writing papers according to APA style, for example when we create custom tables with numeric results.
An important function is force_decimals(), which ensures that a specified number of decimals is shown in the output.
Note that the function round() will not produce the same results as force_decimals() because unfortunately trailing zeros are removed in the R Markdown / Quarto output.
force_decimals(): 1.23, 0.87, 2.35, 1.20round(): 1.23, 0.87, 2.35, 1.2This is because round() returns a numeric vector and force_decimals() (as all other functions in prmisc) returns a string.
Small numbers that are close to zero are rounded to zero by default. This behaviour can be controlled using the argument round_zero:
force_decimals(c(0.004, 0.001, 0.0005, 0.02))
#> [1] "0.00" "0.00" "0.00" "0.02"
force_decimals(c(0.004, 0.001, 0.0005, 0.02), round_zero = FALSE)
#> [1] "< 0.01" "< 0.01" "< 0.01" "0.02"Sometimes we want to leave integers intact (without decimals), which can be accomplished via force_or_cut():
Compare to:
force_decimals(c(1:3, 1.23456, 0.873, 2.3456), decimals = 2)
#> [1] "1.00" "2.00" "3.00" "1.23" "0.87" "2.35"If we do not want the zero, e.g., for p-values or correlation coefficients, use decimals_only():
We can format a p-value, the default is 3 decimals:
Or we format several p-values with one function call:
format_p(c(0.3123, 0.001, 0.00001, 0.19))
#> [1] "$p = .312$" "$p = .001$" "$p < .001$" "$p = .190$"While all other number printing functions have two decimals by default, format_p() has three decimals by default.
We often want to just print the mean and standard deviation of a set of observations, we can do this with print_mean_sd():
Inline: \((M = 5.84\), \(SD = 0.83)\)