f_aov: Perform multiple 'aov()' functions with optional data...

In rfriend: Provides Batch Functions and Visualisation for Basic Statistical Procedures

Perform multiple aov() functions with optional data transformation, inspection and Post Hoc test.

Description

Performs an Analysis of Variance (ANOVA) on a given dataset with options for (Box-Cox) transformations, normality tests, and post hoc analysis. The omnibus table is computed with Type II Sums of Squares via Anova, which is order-invariant for the main effects in unbalanced designs (default summary(aov()) uses Type I SS, where the main-effect p-values depend on the order in which terms appear in the formula). Type II also aligns with the model-based emmeans post hoc tests, so the omnibus table and the pairwise comparisons cannot tell mismatched stories on unbalanced data. Several response parameters can be analysed in sequence and the generated output can be in various formats ('Word', 'pdf', 'Excel').

Usage

f_aov(
  formula,
  data = NULL,
  norm_plots = TRUE,
  interaction_plots = TRUE,
  contrast_plots = FALSE,
  ANCOVA = FALSE,
  transformation = TRUE,
  force_transformation = NULL,
  force_aov = FALSE,
  alpha = 0.05,
  adjust = "sidak",
  anova_type = 2,
  intro_text = TRUE,
  close_generated_files = FALSE,
  open_generated_files = interactive(),
  output_type = "default",
  save_as = NULL,
  save_in_wdir = FALSE,
  ...
)

Arguments

formula	A formula specifying the model to be fitted. More response variables can be added using - or + (e.g., response1 + response2 ~ predictor) to do a sequential aov() for each response parameter.
data	A data frame containing the variables in the model.
norm_plots	Logical. If TRUE, diagnostic residual plots are included in the output files. Default is TRUE.
interaction_plots	Logical. If TRUE, estimated means / interaction plots are included in the output files after the post hoc table. Default is TRUE.
contrast_plots	Logical. If TRUE, a contrast forest plot is added for each categorical post hoc term: one row per pairwise comparison, showing the estimated difference between two levels with its confidence interval and a reference line at zero. A CI that excludes zero indicates a significant difference; because the interval is on the difference itself, this "excludes zero" reading is exact (it is the same information the compact-letter display encodes, but it also shows the direction and magnitude of each difference). Default FALSE because the number of pairwise contrasts grows quickly with the number of factor levels (k levels give k(k-1)/2 contrasts); turn it on when you want the detailed pairwise picture. Main-effect and interaction contrast plots are kept separate: main-effect plots are stored as out$y1$contrast_plot_<term> (e.g. contrast_plot_treatment) while interaction cell-contrast plots are stored as out$y1$interaction_contrast_plot_<term> (e.g. interaction_contrast_plot_a_b). Contrast CIs use the same adjust method as the post hoc p-values, so figure and table agree.
ANCOVA	Logical. If TRUE, prevents automatic conversion of predictors to factors, allowing for Analysis of Covariance (ANCOVA). Default is FALSE.
transformation	Logical or character string. If TRUE, or if "boxcox" applies a f_boxcox() transformation if residuals are not normal. If "bestnormalize", applies f_bestNormalize() transformation. If FALSE no transformation will be applied. Default is TRUE.
force_transformation	Character string. A vector containing the names of response variables that should be transformed regardless of the normality test. Default is NULL
force_aov	Logical. If TRUE, runs the ANOVA even when at least one cell has n = 1 (saturated model). By default (FALSE), such responses are skipped with a warning because F-statistics and p-values are undefined for saturated models. Set to TRUE only for diagnostic purposes – results should not be reported or interpreted as valid. Default is FALSE.
alpha	Numeric. Significance level for ANOVA, post hoc tests, and Shapiro-Wilk test. Default is 0.05.
adjust	Character string specifying the method used to adjust p-values for multiple comparisons. Available methods include: "tukey" Tukey's Honest Significant Difference method, appropriate for all pairwise comparisons. Controls family-wise error rate. "sidak" Sidak correction that controls the family-wise error rate. Less conservative than Bonferroni. "bonferroni" Conservative adjustment that multiplies p-values by the number of comparisons. "none" No adjustment. Equivalent to Fisher's LSD method. "fdr" False Discovery Rate adjustment, controls the expected proportion of false positives among significant results. Default is "sidak".
anova_type	Integer, either 2 or 3. Sums of Squares type for the omnibus ANOVA table computed via Anova. 2 (Default) Type II. Order-invariant in unbalanced designs (drug * dose and dose * drug give the same main-effect p-values), respects the marginality principle (each main effect is tested after all other main effects, ignoring interactions containing it), and is safe with R's default treatment contrasts. Recommended for most unbalanced designs and consistent with the emmeans-based post hoc tests. 3 Type III. Also order-invariant, but tests each term after all other terms including higher-order interactions. Type III is the SPSS / SAS default. For its main-effect rows to be interpretable as effects averaged across the other factors, the model must be fitted with orthogonal (sum / effect / Helmert / polynomial) contrasts. When anova_type = 3 and the user has not supplied their own contrasts via ..., f_aov automatically installs contr.sum / contr.poly for the duration of the call (the previous options("contrasts") is restored on exit). Note that under Type III, when an interaction is significant the main effect rows are conditional on the interaction and should be interpreted with care – the cell means table that f_aov reports automatically when an interaction is significant remains the appropriate summary.
intro_text	Logical. If TRUE, includes a short explanation about ANOVA assumptions in the output file. Default is TRUE.
close_generated_files	Logical. Closes open Excel or Word (NOT pdf) files before writing, depending on the output format. Works on Windows (taskkill), macOS (pkill) and Linux (pkill/soffice). Default FALSE. WARNING: Always save your work before using this option!!
open_generated_files	Logical. Whether to open the generated output files after creation. Defaults to TRUE in an interactive R session and FALSE otherwise (e.g. in scripts or automated pipelines). Set to TRUE or FALSE to override this behaviour explicitly.
output_type	Character string specifying the output format. Default is "default". "default": Returns the object and lets R decide whether to print; auto-prints if unassigned, silent if assigned to a variable. Use print(result) or plot(result) to display the returned object. "console": Forces immediate printing to the console regardless of object assignment. "pdf", "word", "excel": Saves results to a file of the corresponding format. See save_as, save_in_wdir, and open_generated_files for file path and opening behavior. "rmd": Stores the raw markdown string inside the returned object for use in R Markdown documents.
save_as	Character string specifying the output file path (without extension). If a full path is provided, output is saved to that location. If only a filename is given, the file is saved in tempdir(). If only a directory is specified (providing an existing directory with trailing slash), the file is named "dataname_aov_output" in that directory. If an extension is provided the output format specified with option "output_type" will be overruled. Defaults to file.path(tempdir(), "dataname_summary.pdf").
save_in_wdir	Logical. If TRUE, saves the file in the working directory. Default is FALSE, this avoid unintended changes to the global environment. If save_as location is specified save_in_wdir is overwritten by save_as.
...	Additional arguments forwarded to aov. The arguments subset, na.action, and weights are handled specially: when supplied, they are applied via model.frame so that the n=1 cell check, Shapiro-Wilk test, Levene test, optional transformations, residual diagnostics, and emmeans post hoc tests all see the exact same row set as aov() itself. Any other aov() arguments (e.g. contrasts, projections, qr, contrasts.arg) are passed through unchanged.

Details

The function performs the following steps:

• Check if all specified variables are present in the data.

• Ensure that the response variable is numeric.

• Fit the model with aov and compute the omnibus ANOVA table with Type II Sums of Squares via Anova. Type II is used (instead of the default Type I from summary(aov())) because Type I main-effect SS depend on the order of terms in the formula in unbalanced designs, whereas the emmeans-based post hoc tests are model-based and therefore order-invariant. Pairing Type I with emmeans can produce mismatched stories between the omnibus and post hoc tables. Type II keeps both order-invariant and is safe with R's default treatment contrasts (unlike Type III, which would require sum / effect contrasts to be interpretable for main effects).

• Check normality of residuals using the Shapiro-Wilk test.

• If residuals are not normal and transformation = TRUE apply a data transformation.

• If significant differences are found in ANOVA, proceed with post hoc tests using estimated marginal means from emmeans() and Sidak adjustment (or another option of adjust =.

Effect and interaction plots. When interaction_plots = TRUE, an estimated marginal means plot (estimate \pm 95% CI, with jittered raw data and compact-letter-display labels) is added after the post hoc table for each categorical predictor. For a significant categorical interaction, interaction plots are drawn instead: a two-way interaction uses the x-axis plus colour (both orientations), while three- and four-way interactions add facet panels for the remaining factor(s), with one plot per choice of x-axis factor. Interactions involving five or more categorical factors are not plotted (a warning is issued); consult the post hoc cell-means table instead. When the response was transformed (Box-Cox or bestNormalize), the plotted estimates are back-transformed to the original scale (medians). The plots themselves are kept clean for publication (data, axes, and legend only); the descriptive label and explanatory caption are emitted as text above and below each figure in the report. All effect and interaction plots are ggplot2 objects and are stored in the returned object (e.g. out$y1$effect_plot_treatment, out$y1$interaction_plot_a_b_1) so they can be retrieved and customised afterwards. Matches f_glm.

More response variables can be added using - or + (e.g., response1 + response2 ~ predictor) to do a sequential aov() for each response parameter captured in one output file.

Outputs can be generated in multiple formats ("pdf", "word", "excel" and "rmd") as specified by output_type. The function also closes any open 'Word' files to avoid conflicts when generating 'Word' documents. If output_type = "rmd" is used it is adviced to use it in a chunk with {r, echo=FALSE, results='asis'}

*Non-significant ANOVA results*: When the overall F-test is not significant, f_aov still reports the estimated marginal means table, but with all pairwise comparison letters replaced by *"ns"*. The numeric estimates (and their confidence intervals) are provided because they are often needed for manuscript tables, especially when the response was back-transformed from a Box-Cox or bestNormalize scale - the raw descriptive means and the emmeans values can differ, and it is the emmeans values that correspond to the actual model. The *"ns"* labels signal that pairwise differences should not be interpreted.

This function requires [Pandoc](https://github.com/jgm/pandoc/releases/tag) (version 1.12.3 or higher), a universal document converter.

• Windows: Install Pandoc and ensure the installation folder.
  (e.g., "C:/Users/your_username/AppData/Local/Pandoc") is added to your system PATH.

• macOS: If using Homebrew, Pandoc is typically installed in "/usr/local/bin". Alternatively, download the .pkg installer and verify that the binary's location is in your PATH.

• Linux: Install Pandoc through your distribution's package manager (commonly installed in "/usr/bin" or "/usr/local/bin") or manually, and ensure the directory containing Pandoc is in your PATH.

• If Pandoc is not found, this function may not work as intended.

Value

An object of class 'f_aov' containing the fitted model (aov_test), the Type II omnibus ANOVA table from Anova (aov_summary), normality and homogeneity diagnostics, optional transformation results, and the emmeans post hoc tests. Using the option output_type, it can also generate output as R Markdown, 'Word', 'pdf', or 'Excel' files. Includes print and plot methods for 'f_aov' objects.

Multiple Testing Across Response Variables

When several response variables are analysed in a single call (e.g. y1 + y2 + y3 ~ treatment), each ANOVA is an independent null-hypothesis test at level alpha. The post hoc adjustments (adjust = "sidak", "tukey", etc.) only control the family-wise error rate within one ANOVA (across pairwise group comparisons for that response). They do not protect against the inflation of Type I error across the set of responses.

Practical implication: With k independent response variables all tested at \alpha = 0.05, the probability of obtaining at least one false positive is 1 - (1 - 0.05)^k, which reaches ~40% for k = 10.

When this matters: The risk is highest in exploratory studies where many responses are screened simultaneously without a clear a priori hypothesis for each one. It is less of a concern when each response is a pre-specified primary outcome with its own biological rationale.

Possible remedies:

• Bonferroni correction across responses: use alpha = 0.05 / k where k is the number of response variables. Conservative but simple.

• False Discovery Rate (FDR): apply p.adjust(p_values, method = "fdr") to the vector of per-response ANOVA p-values after the fact.

• MANOVA: if the responses are correlated and you want a single omnibus test across all of them, use manova() before interpreting individual ANOVAs.

• Pre-registration: declare primary vs. exploratory responses before data collection to justify differential correction thresholds.

Author(s)

Sander H. van Delden [e-mail]

Examples

# The left hand side contains two response variables,
# so two aov's will be conducted, i.e. "Sepal.Width"
# and "Sepal.Length" in response to the explanatory variable: "Species".
f_aov_out <- f_aov(Sepal.Width + Sepal.Length ~ Species,
                   data = iris,
                   # Save output in MS Word file (Default is console)
                   output_type = "word",
                   # Do bestNormalize transformation for non-normal residual (Default is boxcox)
                   transformation = "bestnormalize"
                   )

# Print output to the console.
print(f_aov_out)

# Plot residual plots.
plot(f_aov_out)


#To print rmd output set chunck option to results = 'asis' and use cat().
f_aov_rmd_out <- f_aov(Sepal.Width ~ Species, data = iris, output_type = "rmd")
cat(f_aov_rmd_out$rmd)

rfriend documentation built on July 7, 2026, 1:06 a.m.