{"id":10875,"date":"2019-11-29T08:37:37","date_gmt":"2019-11-29T06:37:37","guid":{"rendered":"https:\/\/www.datanovia.com\/en\/?post_type=dt_lessons&#038;p=10875"},"modified":"2019-11-29T08:41:02","modified_gmt":"2019-11-29T06:41:02","slug":"mixed-anova-in-r","status":"publish","type":"dt_lessons","link":"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/","title":{"rendered":"Mixed ANOVA in R"},"content":{"rendered":"<div id=\"rdoc\">\n<p><strong>Mixed ANOVA<\/strong> is used to compare the means of groups cross-classified by two different types of factor variables, including:<\/p>\n<ul>\n<li><strong>between-subjects factors<\/strong>, which have independent categories (e.g., gender: male\/female)<\/li>\n<li><strong>within-subjects factors<\/strong>, which have related categories also known as repeated measures (e.g., time: before\/after treatment).<\/li>\n<\/ul>\n<p>The mixed ANOVA test is also referred as <em>mixed design ANOVA<\/em> and <em>mixed measures ANOVA<\/em>.<\/p>\n<p>This chapter describes different types of mixed ANOVA, including:<\/p>\n<ul>\n<li><strong>two-way mixed ANOVA<\/strong>, used to compare the means of groups cross-classified by two independent categorical variables, including one between-subjects and one within-subjects factors.<\/li>\n<li><strong>three-way mixed ANOVA<\/strong>, used to evaluate if there is a three-way interaction between three independent variables, including between-subjects and within-subjects factors. You can have two different designs for three-way mixed ANOVA:\n<ol style=\"list-style-type: decimal;\">\n<li>one between-subjects factor and two within-subjects factors<\/li>\n<li>two between-subjects factor and one within-subjects factor<\/li>\n<\/ol>\n<\/li>\n<\/ul>\n<p>You will learn how to:<\/p>\n<ul>\n<li><strong>Compute and interpret the different mixed ANOVA tests in R<\/strong>.<\/li>\n<li><strong>Check mixed ANOVA test assumptions<\/strong><\/li>\n<li><strong>Perform post-hoc tests<\/strong>, multiple pairwise comparisons between groups to identify which groups are different<\/li>\n<li><strong>Visualize the data<\/strong> using box plots, add ANOVA and pairwise comparisons p-values to the plot<\/li>\n<\/ul>\n<p>Contents:<\/p>\n<div id=\"TOC\">\n<ul>\n<li><a href=\"#assumptions\">Assumptions<\/a><\/li>\n<li><a href=\"#prerequisites\">Prerequisites<\/a><\/li>\n<li><a href=\"#two-way-mixed\">Two-way mixed ANOVA<\/a>\n<ul>\n<li><a href=\"#data-preparation\">Data preparation<\/a><\/li>\n<li><a href=\"#summary-statistics\">Summary statistics<\/a><\/li>\n<li><a href=\"#visualization\">Visualization<\/a><\/li>\n<li><a href=\"#check-assumptions\">Check assumptions<\/a><\/li>\n<li><a href=\"#computation\">Computation<\/a><\/li>\n<li><a href=\"#post-hoc-tests\">Post-hoc tests<\/a><\/li>\n<li><a href=\"#report\">Report<\/a><\/li>\n<\/ul>\n<\/li>\n<li><a href=\"#three-way-bbw-b\">Three-way mixed ANOVA: 2 between- and 1 within-subjects factors<\/a>\n<ul>\n<li><a href=\"#data-preparation-1\">Data preparation<\/a><\/li>\n<li><a href=\"#summary-statistics-1\">Summary statistics<\/a><\/li>\n<li><a href=\"#visualization-1\">Visualization<\/a><\/li>\n<li><a href=\"#check-assumptions-1\">Check assumptions<\/a><\/li>\n<li><a href=\"#computation-1\">Computation<\/a><\/li>\n<li><a href=\"#post-hoc-tests-1\">Post-hoc tests<\/a><\/li>\n<li><a href=\"#report-1\">Report<\/a><\/li>\n<\/ul>\n<\/li>\n<li><a href=\"#three-way-bww-b\">Three-way Mixed ANOVA: 1 between- and 2 within-subjects factors<\/a>\n<ul>\n<li><a href=\"#data-preparation-2\">Data preparation<\/a><\/li>\n<li><a href=\"#summary-statistics-2\">Summary statistics<\/a><\/li>\n<li><a href=\"#visualization-2\">Visualization<\/a><\/li>\n<li><a href=\"#check-assumptions-2\">Check assumptions<\/a><\/li>\n<li><a href=\"#computation-2\">Computation<\/a><\/li>\n<li><a href=\"#post-hoc-tests-2\">Post-hoc tests<\/a><\/li>\n<li><a href=\"#report-2\">Report<\/a><\/li>\n<\/ul>\n<\/li>\n<li><a href=\"#summary\">Summary<\/a><\/li>\n<\/ul>\n<\/div>\n<div class='dt-sc-hr-invisible-medium  '><\/div>\n<div class='dt-sc-ico-content type1'><div class='custom-icon' ><a href='https:\/\/www.datanovia.com\/en\/product\/practical-statistics-in-r-for-comparing-groups-numerical-variables\/' target='_blank'><span class='fa fa-book'><\/span><\/a><\/div><h4><a href='https:\/\/www.datanovia.com\/en\/product\/practical-statistics-in-r-for-comparing-groups-numerical-variables\/' target='_blank'> Related Book <\/a><\/h4>Practical Statistics in R II - Comparing Groups: Numerical Variables<\/div>\n<div class='dt-sc-hr-invisible-medium  '><\/div>\n<div id=\"assumptions\" class=\"section level2\">\n<h2>Assumptions<\/h2>\n<p>The mixed ANOVA makes the following assumptions about the data:<\/p>\n<ul>\n<li><strong>No significant outliers<\/strong> in any cell of the design. This can be checked by visualizing the data using box plot methods and by using the function <code>identify_outliers()<\/code> [rstatix package].<\/li>\n<li><strong>Normality<\/strong>: the outcome (or dependent) variable should be approximately normally distributed in each cell of the design. This can be checked using the <strong>Shapiro-Wilk normality test<\/strong> (<code>shapiro_test()<\/code> [rstatix]) or by visual inspection using <strong>QQ plot<\/strong> (<code>ggqqplot()<\/code> [ggpubr package]).<\/li>\n<li><strong>Homogeneity of variances<\/strong>: the variance of the outcome variable should be equal between the groups of the between-subjects factors. This can be assessed using the <strong>Levene\u2019s test for equality of variances<\/strong> (<code>levene_test()<\/code> [rstatix]).<\/li>\n<li><strong>Assumption of sphericity<\/strong>: the variance of the differences between within-subjects groups should be equal. This can be checked using the <strong>Mauchly\u2019s test of sphericity<\/strong>, which is automatically reported when using the <code>anova_test()<\/code> R function.<\/li>\n<li><strong>Homogeneity of covariances<\/strong> tested by Box\u2019s M. The covariance matrices should be equal across the cells formed by the between-subjects factors.<\/li>\n<\/ul>\n<p>Before computing mixed ANOVA test, you need to perform some preliminary tests to check if the assumptions are met.<\/p>\n<\/div>\n<div id=\"prerequisites\" class=\"section level2\">\n<h2>Prerequisites<\/h2>\n<p>Make sure that you have installed the following R packages:<\/p>\n<ul>\n<li><code>tidyverse<\/code> for data manipulation and visualization<\/li>\n<li><code>ggpubr<\/code> for creating easily publication ready plots<\/li>\n<li><code>rstatix<\/code> provides pipe-friendly R functions for easy statistical analyses<\/li>\n<li><code>datarium<\/code>: contains required data sets for this chapter<\/li>\n<\/ul>\n<p>Start by loading the following R packages:<\/p>\n<pre class=\"r\"><code>library(tidyverse)\r\nlibrary(ggpubr)\r\nlibrary(rstatix)<\/code><\/pre>\n<p>Key R functions:<\/p>\n<ul>\n<li><code>anova_test()<\/code> [rstatix package], a wrapper around <code>car::Anova()<\/code> for making easy the computation of repeated measures ANOVA. Key arguments for performing repeated measures ANOVA:\n<ul>\n<li><code>data<\/code>: data frame<\/li>\n<li><code>dv<\/code>: (numeric) the dependent (or outcome) variable name.<\/li>\n<li><code>wid<\/code>: variable name specifying the case\/sample identifier.<\/li>\n<li><code>between<\/code>: between-subjects factor or grouping variable<\/li>\n<li><code>within<\/code>: within-subjects factor or grouping variable<\/li>\n<\/ul>\n<\/li>\n<li><code>get_anova_table()<\/code> [rstatix package]. Extracts the ANOVA table from the output of <code>anova_test()<\/code>. It returns ANOVA table that is automatically corrected for eventual deviation from the sphericity assumption. The default is to apply automatically the Greenhouse-Geisser sphericity correction to only within-subject factors violating the sphericity assumption (i.e., Mauchly\u2019s test p-value is significant, p &lt;= 0.05). Read more in Chapter @ref(mauchly-s-test-of-sphericity-in-r).<\/li>\n<\/ul>\n<\/div>\n<div id=\"two-way-mixed\" class=\"section level2\">\n<h2>Two-way mixed ANOVA<\/h2>\n<div id=\"data-preparation\" class=\"section level3\">\n<h3>Data preparation<\/h3>\n<p>We\u2019ll use the <code>anxiety<\/code> dataset [in the datarium package], which contains the anxiety score, measured at three time points (t1, t2 and t3), of three groups of individuals practicing physical exercises at different levels (grp1: basal, grp2: moderate and grp3: high)<\/p>\n<div class=\"block\">\n<p>Two-way mixed ANOVA can be used to evaluate if there is interaction between group and time in explaining the anxiety score.<\/p>\n<\/div>\n<p>Load and show one random row by group:<\/p>\n<pre class=\"r\"><code># Wide format\r\nset.seed(123)\r\ndata(\"anxiety\", package = \"datarium\")\r\nanxiety %&gt;% sample_n_by(group, size = 1)<\/code><\/pre>\n<pre><code>## # A tibble: 3 x 5\r\n##   id    group    t1    t2    t3\r\n##   &lt;fct&gt; &lt;fct&gt; &lt;dbl&gt; &lt;dbl&gt; &lt;dbl&gt;\r\n## 1 5     grp1   16.5  15.8  15.7\r\n## 2 27    grp2   17.8  17.7  16.9\r\n## 3 37    grp3   17.1  15.6  14.3<\/code><\/pre>\n<pre class=\"r\"><code># Gather the columns t1, t2 and t3 into long format.\r\n# Convert id and time into factor variables\r\nanxiety &lt;- anxiety %&gt;%\r\n  gather(key = \"time\", value = \"score\", t1, t2, t3) %&gt;%\r\n  convert_as_factor(id, time)\r\n# Inspect some random rows of the data by groups\r\nset.seed(123)\r\nanxiety %&gt;% sample_n_by(group, time, size = 1)<\/code><\/pre>\n<pre><code>## # A tibble: 9 x 4\r\n##   id    group time  score\r\n##   &lt;fct&gt; &lt;fct&gt; &lt;fct&gt; &lt;dbl&gt;\r\n## 1 5     grp1  t1     16.5\r\n## 2 12    grp1  t2     17.7\r\n## 3 7     grp1  t3     16.5\r\n## 4 29    grp2  t1     18.4\r\n## 5 30    grp2  t2     18.9\r\n## 6 16    grp2  t3     12.7\r\n## # \u2026 with 3 more rows<\/code><\/pre>\n<\/div>\n<div id=\"summary-statistics\" class=\"section level3\">\n<h3>Summary statistics<\/h3>\n<p>Group the data by <code>time<\/code> and <code>group<\/code>, and then compute some summary statistics of the <code>score<\/code> variable: mean and sd (standard deviation)<\/p>\n<pre class=\"r\"><code>anxiety %&gt;%\r\n  group_by(time, group) %&gt;%\r\n  get_summary_stats(score, type = \"mean_sd\")<\/code><\/pre>\n<pre><code>## # A tibble: 9 x 6\r\n##   group time  variable     n  mean    sd\r\n##   &lt;fct&gt; &lt;fct&gt; &lt;chr&gt;    &lt;dbl&gt; &lt;dbl&gt; &lt;dbl&gt;\r\n## 1 grp1  t1    score       15  17.1  1.63\r\n## 2 grp2  t1    score       15  16.6  1.57\r\n## 3 grp3  t1    score       15  17.0  1.32\r\n## 4 grp1  t2    score       15  16.9  1.70\r\n## 5 grp2  t2    score       15  16.5  1.70\r\n## 6 grp3  t2    score       15  15.0  1.39\r\n## # \u2026 with 3 more rows<\/code><\/pre>\n<\/div>\n<div id=\"visualization\" class=\"section level3\">\n<h3>Visualization<\/h3>\n<p>Create a box plots:<\/p>\n<pre class=\"r\"><code>bxp &lt;- ggboxplot(\r\n  anxiety, x = \"time\", y = \"score\",\r\n  color = \"group\", palette = \"jco\"\r\n  )\r\nbxp<\/code><\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/dn-tutorials\/r-statistics-2-comparing-groups-means\/figures\/047-mixed-anova-two-way-boxplots-1.png\" width=\"384\" \/><\/p>\n<\/div>\n<div id=\"check-assumptions\" class=\"section level3\">\n<h3>Check assumptions<\/h3>\n<div id=\"outliers\" class=\"section level4\">\n<h4>Outliers<\/h4>\n<p>Outliers can be easily identified using box plot methods, implemented in the R function <code>identify_outliers()<\/code> [rstatix package].<\/p>\n<pre class=\"r\"><code>anxiety %&gt;%\r\n  group_by(time, group) %&gt;%\r\n  identify_outliers(score)<\/code><\/pre>\n<pre><code>## [1] group      time       id         score      is.outlier is.extreme\r\n## &lt;0 rows&gt; (or 0-length row.names)<\/code><\/pre>\n<div class=\"success\">\n<p>There were no extreme outliers.<\/p>\n<\/div>\n<div class=\"warning\">\n<p>Note that, in the situation where you have extreme outliers, this can be due to: 1) data entry errors, measurement errors or unusual values.<\/p>\n<p>Yo can include the outlier in the analysis anyway if you do not believe the result will be substantially affected. This can be evaluated by comparing the result of the ANOVA with and without the outlier.<\/p>\n<p>It\u2019s also possible to keep the outliers in the data and perform robust ANOVA test using the WRS2 package.<\/p>\n<\/div>\n<\/div>\n<div id=\"normality-assumption\" class=\"section level4\">\n<h4>Normality assumption<\/h4>\n<p>The normality assumption can be checked by computing Shapiro-Wilk test for each combinations of factor levels. If the data is normally distributed, the p-value should be greater than 0.05.<\/p>\n<pre class=\"r\"><code>anxiety %&gt;%\r\n  group_by(time, group) %&gt;%\r\n  shapiro_test(score)<\/code><\/pre>\n<pre><code>## # A tibble: 9 x 5\r\n##   group time  variable statistic     p\r\n##   &lt;fct&gt; &lt;fct&gt; &lt;chr&gt;        &lt;dbl&gt; &lt;dbl&gt;\r\n## 1 grp1  t1    score        0.964 0.769\r\n## 2 grp2  t1    score        0.977 0.949\r\n## 3 grp3  t1    score        0.954 0.588\r\n## 4 grp1  t2    score        0.956 0.624\r\n## 5 grp2  t2    score        0.935 0.328\r\n## 6 grp3  t2    score        0.952 0.558\r\n## # \u2026 with 3 more rows<\/code><\/pre>\n<div class=\"success\">\n<p>The score were normally distributed (p &gt; 0.05) for each cell, as assessed by Shapiro-Wilk\u2019s test of normality.<\/p>\n<\/div>\n<p>Note that, if your sample size is greater than 50, the normal QQ plot is preferred because at larger sample sizes the Shapiro-Wilk test becomes very sensitive even to a minor deviation from normality.<\/p>\n<p>QQ plot draws the correlation between a given data and the normal distribution.<\/p>\n<pre class=\"r\"><code>ggqqplot(anxiety, \"score\", ggtheme = theme_bw()) +\r\n  facet_grid(time ~ group)<\/code><\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/dn-tutorials\/r-statistics-2-comparing-groups-means\/figures\/047-mixed-anova-two-way-qq-plot-1.png\" width=\"576\" \/><\/p>\n<div class=\"success\">\n<p>All the points fall approximately along the reference line, for each cell. So we can assume normality of the data.<\/p>\n<\/div>\n<div class=\"warning\">\n<p>In the situation where the assumptions are not met, you could consider running the two-way repeated measures ANOVA on the transformed or performing a robust ANOVA test using the WRS2 R package.<\/p>\n<\/div>\n<\/div>\n<div id=\"homogneity-of-variance-assumption\" class=\"section level4\">\n<h4>Homogneity of variance assumption<\/h4>\n<p>The homogeneity of variance assumption of the between-subject factor (<code>group<\/code>) can be checked using the Levene\u2019s test. The test is performed at each level of <code>time<\/code> variable:<\/p>\n<pre class=\"r\"><code>anxiety %&gt;%\r\n  group_by(time) %&gt;%\r\n  levene_test(score ~ group)<\/code><\/pre>\n<pre><code>## # A tibble: 3 x 5\r\n##   time    df1   df2 statistic     p\r\n##   &lt;fct&gt; &lt;int&gt; &lt;int&gt;     &lt;dbl&gt; &lt;dbl&gt;\r\n## 1 t1        2    42     0.176 0.839\r\n## 2 t2        2    42     0.249 0.781\r\n## 3 t3        2    42     0.335 0.717<\/code><\/pre>\n<div class=\"success\">\n<p>There was homogeneity of variances, as assessed by Levene\u2019s test (p &gt; 0.05).<\/p>\n<\/div>\n<div class=\"warning\">\n<p>Note that, if you do not have homogeneity of variances, you can try to transform the outcome (dependent) variable to correct for the unequal variances.<\/p>\n<p>It\u2019s also possible to perform robust ANOVA test using the WRS2 R package.<\/p>\n<\/div>\n<\/div>\n<div id=\"homogeneity-of-covariances-assumption\" class=\"section level4\">\n<h4>Homogeneity of covariances assumption<\/h4>\n<p>The homogeneity of covariances of the between-subject factor (<code>group<\/code>) can be evaluated using the <strong>Box\u2019s M-test<\/strong> implemented in the <code>rstatix<\/code> package. If this test is statistically significant (i.e., p &lt; 0.001), you do not have equal covariances, but if the test is not statistically significant (i.e., p &gt; 0.001), you have equal covariances and you have not violated the assumption of homogeneity of covariances.<\/p>\n<div class=\"warning\">\n<p>Note that, the Box\u2019s M is highly sensitive, so unless p &lt; 0.001 and your sample sizes are unequal, ignore it. However, if significant and you have unequal sample sizes, the test is not robust (<a href=\"https:\/\/en.wikiversity.org\/wiki\/Box%27s_M\">https:\/\/en.wikiversity.org\/wiki\/Box%27s_M<\/a>, Tabachnick &amp; Fidell, 2001).<\/p>\n<\/div>\n<p>Compute Box\u2019s M-test:<\/p>\n<pre class=\"r\"><code>box_m(anxiety[, \"score\", drop = FALSE], anxiety$group)<\/code><\/pre>\n<pre><code>## # A tibble: 1 x 4\r\n##   statistic p.value parameter method                                             \r\n##       &lt;dbl&gt;   &lt;dbl&gt;     &lt;dbl&gt; &lt;chr&gt;                                              \r\n## 1      1.93   0.381         2 Box's M-test for Homogeneity of Covariance Matrices<\/code><\/pre>\n<div class=\"success\">\n<p>There was homogeneity of covariances, as assessed by Box\u2019s test of equality of covariance matrices (p &gt; 0.001).<\/p>\n<\/div>\n<div class=\"warning\">\n<p>Note that, if you do not have homogeneity of covariances, you could consider separating your analyses into distinct repeated measures ANOVAs for each group. Alternatively, you could omit the interpretation of the interaction term.<\/p>\n<p>Unfortunately, it is difficult to remedy a failure of this assumption. Often, a mixed ANOVA is run anyway and the violation noted.<\/p>\n<\/div>\n<\/div>\n<div id=\"assumption-of-sphericity\" class=\"section level4\">\n<h4>Assumption of sphericity<\/h4>\n<p>As mentioned in previous sections, the assumption of sphericity will be automatically checked during the computation of the ANOVA test using the R function <code>anova_test()<\/code> [rstatix package]. The Mauchly\u2019s test is internally used to assess the sphericity assumption.<\/p>\n<p>By using the function <code>get_anova_table()<\/code> [rstatix] to extract the ANOVA table, the Greenhouse-Geisser sphericity correction is automatically applied to factors violating the sphericity assumption.<\/p>\n<\/div>\n<\/div>\n<div id=\"computation\" class=\"section level3\">\n<h3>Computation<\/h3>\n<pre class=\"r\"><code># Two-way mixed ANOVA test\r\nres.aov &lt;- anova_test(\r\n  data = anxiety, dv = score, wid = id,\r\n  between = group, within = time\r\n  )\r\nget_anova_table(res.aov)<\/code><\/pre>\n<pre><code>## ANOVA Table (type II tests)\r\n## \r\n##       Effect DFn DFd      F        p p&lt;.05   ges\r\n## 1      group   2  42   4.35 1.90e-02     * 0.168\r\n## 2       time   2  84 394.91 1.91e-43     * 0.179\r\n## 3 group:time   4  84 110.19 1.38e-32     * 0.108<\/code><\/pre>\n<div class=\"success\">\n<p>From the output above, it can be seen that, there is a statistically significant two-way interactions between group and time on anxiety score, F(4, 84) = 110.18, p &lt; 0.0001.<\/p>\n<\/div>\n<\/div>\n<div id=\"post-hoc-tests\" class=\"section level3\">\n<h3>Post-hoc tests<\/h3>\n<p>A <strong>significant two-way interaction<\/strong> indicates that the impact that one factor has on the outcome variable depends on the level of the other factor (and vice versa). So, you can decompose a significant two-way interaction into:<\/p>\n<ul>\n<li><strong>Simple main effect<\/strong>: run one-way model of the first variable (factor A) at each level of the second variable (factor B),<\/li>\n<li><strong>Simple pairwise comparisons<\/strong>: if the simple main effect is significant, run multiple pairwise comparisons to determine which groups are different.<\/li>\n<\/ul>\n<p>For a <strong>non-significant two-way interaction<\/strong>, you need to determine whether you have any statistically significant <strong>main effects<\/strong> from the ANOVA output.<\/p>\n<div id=\"procedure-for-a-significant-two-way-interaction\" class=\"section level4\">\n<h4>Procedure for a significant two-way interaction<\/h4>\n<p><strong>Simple main effect of group variable<\/strong>. In our example, we\u2019ll investigate the effect of the between-subject factor <code>group<\/code> on anxiety score at every <code>time<\/code> point.<\/p>\n<pre class=\"r\"><code># Effect of group at each time point\r\none.way &lt;- anxiety %&gt;%\r\n  group_by(time) %&gt;%\r\n  anova_test(dv = score, wid = id, between = group) %&gt;%\r\n  get_anova_table() %&gt;%\r\n  adjust_pvalue(method = \"bonferroni\")\r\none.way<\/code><\/pre>\n<pre><code>## # A tibble: 3 x 9\r\n##   time  Effect   DFn   DFd      F         p `p&lt;.05`   ges     p.adj\r\n##   &lt;fct&gt; &lt;chr&gt;  &lt;dbl&gt; &lt;dbl&gt;  &lt;dbl&gt;     &lt;dbl&gt; &lt;chr&gt;   &lt;dbl&gt;     &lt;dbl&gt;\r\n## 1 t1    group      2    42  0.365 0.696     \"\"      0.017 1        \r\n## 2 t2    group      2    42  5.84  0.006     *       0.218 0.018    \r\n## 3 t3    group      2    42 13.8   0.0000248 *       0.396 0.0000744<\/code><\/pre>\n<pre class=\"r\"><code># Pairwise comparisons between group levels\r\npwc &lt;- anxiety %&gt;%\r\n  group_by(time) %&gt;%\r\n  pairwise_t_test(score ~ group, p.adjust.method = \"bonferroni\")\r\npwc<\/code><\/pre>\n<pre><code>## # A tibble: 9 x 10\r\n##   time  .y.   group1 group2    n1    n2       p p.signif   p.adj p.adj.signif\r\n## * &lt;fct&gt; &lt;chr&gt; &lt;chr&gt;  &lt;chr&gt;  &lt;int&gt; &lt;int&gt;   &lt;dbl&gt; &lt;chr&gt;      &lt;dbl&gt; &lt;chr&gt;       \r\n## 1 t1    score grp1   grp2      15    15 0.43    ns       1       ns          \r\n## 2 t1    score grp1   grp3      15    15 0.895   ns       1       ns          \r\n## 3 t1    score grp2   grp3      15    15 0.51    ns       1       ns          \r\n## 4 t2    score grp1   grp2      15    15 0.435   ns       1       ns          \r\n## 5 t2    score grp1   grp3      15    15 0.00212 **       0.00636 **          \r\n## 6 t2    score grp2   grp3      15    15 0.0169  *        0.0507  ns          \r\n## # \u2026 with 3 more rows<\/code><\/pre>\n<div class=\"success\">\n<p>Considering the Bonferroni adjusted p-value (p.adj), it can be seen that the simple main effect of group was significant at t2 (p = 0.018) and t3 (p &lt; 0.0001) but not at t1 (p = 1).<\/p>\n<p>Pairwise comparisons show that the mean anxiety score was significantly different in grp1 vs grp3 comparison at t2 (p = 0.0063); in grp1 vs grp3 (p &lt; 0.0001) and in grp2 vs grp3 (p = 0.0013) at t3.<\/p>\n<\/div>\n<p><strong>Simple main effects of time variable<\/strong>. It\u2019s also possible to perform the same analyze for the within-subject <code>time<\/code> variable at each level of <code>group<\/code> as shown in the following R code. You don\u2019t necessarily need to do this analysis.<\/p>\n<pre class=\"r\"><code># Effect of time at each level of exercises group\r\none.way2 &lt;- anxiety %&gt;%\r\n  group_by(group) %&gt;%\r\n  anova_test(dv = score, wid = id, within = time) %&gt;%\r\n  get_anova_table() %&gt;%\r\n  adjust_pvalue(method = \"bonferroni\")\r\none.way2<\/code><\/pre>\n<pre><code>## # A tibble: 3 x 9\r\n##   group Effect   DFn   DFd     F        p `p&lt;.05`   ges    p.adj\r\n##   &lt;fct&gt; &lt;chr&gt;  &lt;dbl&gt; &lt;dbl&gt; &lt;dbl&gt;    &lt;dbl&gt; &lt;chr&gt;   &lt;dbl&gt;    &lt;dbl&gt;\r\n## 1 grp1  time       2    28  14.8 4.05e- 5 *       0.024 1.21e- 4\r\n## 2 grp2  time       2    28  77.5 3.88e-12 *       0.086 1.16e-11\r\n## 3 grp3  time       2    28 490.  1.64e-22 *       0.531 4.92e-22<\/code><\/pre>\n<pre class=\"r\"><code># Pairwise comparisons between time points at each group levels\r\n# Paired t-test is used because we have repeated measures by time\r\npwc2 &lt;- anxiety %&gt;%\r\n  group_by(group) %&gt;%\r\n  pairwise_t_test(\r\n    score ~ time, paired = TRUE, \r\n    p.adjust.method = \"bonferroni\"\r\n    ) %&gt;%\r\n  select(-df, -statistic, -p) # Remove details\r\npwc2<\/code><\/pre>\n<pre><code>## # A tibble: 9 x 8\r\n##   group .y.   group1 group2    n1    n2        p.adj p.adj.signif\r\n## * &lt;fct&gt; &lt;chr&gt; &lt;chr&gt;  &lt;chr&gt;  &lt;int&gt; &lt;int&gt;        &lt;dbl&gt; &lt;chr&gt;       \r\n## 1 grp1  score t1     t2        15    15 0.194        ns          \r\n## 2 grp1  score t1     t3        15    15 0.002        **          \r\n## 3 grp1  score t2     t3        15    15 0.006        **          \r\n## 4 grp2  score t1     t2        15    15 0.268        ns          \r\n## 5 grp2  score t1     t3        15    15 0.000000151  ****        \r\n## 6 grp2  score t2     t3        15    15 0.0000000612 ****        \r\n## # \u2026 with 3 more rows<\/code><\/pre>\n<div class=\"success\">\n<p>There was a statistically significant effect of time on anxiety score for each of the three groups. Using pairwise paired t-test comparisons, it can be seen that for grp1 and grp2, the mean anxiety score was not statistically significantly different between t1 and t2 time points.<\/p>\n<p>The pairwise comparisons t1 vs t3 and t2 vs t3 were statistically significantly different for all groups.<\/p>\n<\/div>\n<\/div>\n<div id=\"procedure-for-non-significant-two-way-interaction\" class=\"section level4\">\n<h4>Procedure for non-significant two-way interaction<\/h4>\n<p>If the interaction is not significant, you need to interpret the main effects for each of the two variables: <code>group<\/code> and `time. A significant main effect can be followed up with pairwise comparisons.<\/p>\n<div class=\"success\">\n<p>In our example, there was a statistically significant main effects of group (F(2, 42) = 4.35, p = 0.02) and time (F(2, 84) = 394.91, p &lt; 0.0001) on the anxiety score.<\/p>\n<\/div>\n<p>Perform multiple pairwise paired t-tests for the <code>time<\/code> variable, ignoring group. P-values are adjusted using the Bonferroni multiple testing correction method.<\/p>\n<pre class=\"r\"><code>anxiety %&gt;%\r\n  pairwise_t_test(\r\n    score ~ time, paired = TRUE, \r\n    p.adjust.method = \"bonferroni\"\r\n  )<\/code><\/pre>\n<div class=\"success\">\n<p>All pairwise comparisons are significant.<\/p>\n<\/div>\n<p>You can perform a similar analysis for the <code>group<\/code> variable.<\/p>\n<pre class=\"r\"><code>anxiety %&gt;%\r\n  pairwise_t_test(\r\n    score ~ group, \r\n    p.adjust.method = \"bonferroni\"\r\n  )<\/code><\/pre>\n<div class=\"success\">\n<p>All pairwise comparisons are significant except grp1 vs grp2.<\/p>\n<\/div>\n<\/div>\n<\/div>\n<div id=\"report\" class=\"section level3\">\n<h3>Report<\/h3>\n<p>There was a statistically significant interaction between exercises group and time in explaining the anxiety score, F(4, 84) = 110.19, p &lt; 0.0001.<\/p>\n<p>Considering the Bonferroni adjusted p-value, the simple main effect of exercises group was significant at t2 (p = 0.018) and t3 (p &lt; 0.0001) but not at t1 (p = 1).<\/p>\n<p>Pairwise comparisons show that the mean anxiety score was significantly different in grp1 vs grp3 comparison at t2 (p = 0.0063); in grp1 vs grp3 (p &lt; 0.0001) and in grp2 vs grp3 (p = 0.0013) at t3.<\/p>\n<div class=\"warning\">\n<p>Note that, for the plot below, we only need the pairwise comparison results for t2 and t3 but not for t1 (because the simple main effect of exercises group was not significant at this time point). We\u2019ll filter the comparison results accordingly.<\/p>\n<\/div>\n<pre class=\"r\"><code># Visualization: boxplots with p-values\r\npwc &lt;- pwc %&gt;% add_xy_position(x = \"time\")\r\npwc.filtered &lt;- pwc %&gt;% filter(time != \"t1\")\r\nbxp + \r\n  stat_pvalue_manual(pwc.filtered, tip.length = 0, hide.ns = TRUE) +\r\n  labs(\r\n    subtitle = get_test_label(res.aov, detailed = TRUE),\r\n    caption = get_pwc_label(pwc)\r\n  )<\/code><\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/dn-tutorials\/r-statistics-2-comparing-groups-means\/figures\/047-mixed-anova-two-way-mixed-boxplots-with-p-values-1.png\" width=\"576\" \/><\/p>\n<\/div>\n<\/div>\n<div id=\"three-way-bbw-b\" class=\"section level2\">\n<h2>Three-way mixed ANOVA: 2 between- and 1 within-subjects factors<\/h2>\n<p>This section describes how to compute the three-way mixed ANOVA, in R, for a situation where you have <strong>two between-subjects factors and one within-subjects factor<\/strong>.<\/p>\n<p>This setting is for investigating group differences over time (i.e., the within-subjects factor) where groups are formed by the combination of two between-subjects factors. For example, you might want to understand how performance score changes over time (e.g., 0, 4 and 8 months) depending on gender (i.e., male\/female) and stress (low, moderate and high stress).<\/p>\n<div id=\"data-preparation-1\" class=\"section level3\">\n<h3>Data preparation<\/h3>\n<p>We\u2019ll use <code>performance<\/code> dataset [datarium package] containing the performance score measures of participants at two time points. The aim of this study is to evaluate the effect of gender and stress on performance score.<\/p>\n<p>The data contains the following variables:<\/p>\n<ol style=\"list-style-type: decimal;\">\n<li>Performance score (outcome or dependent variable) measured at two time points, <code>t1<\/code> and <code>t2<\/code>.<\/li>\n<li>Two between-subjects factors: <code>gender<\/code> (levels: male and female) and <code>stress<\/code> (low, moderate, high)<\/li>\n<\/ol>\n<ol style=\"list-style-type: decimal;\" start=\"3\">\n<li>One within-subjects factor, <code>time<\/code>, which has two time points: <code>t1<\/code> and <code>t2<\/code>.<\/li>\n<\/ol>\n<p>Load and inspect the data by showing one random row by group:<\/p>\n<pre class=\"r\"><code># Load and inspect the data\r\n# Wide format\r\nset.seed(123)\r\ndata(\"performance\", package = \"datarium\")\r\nperformance %&gt;% sample_n_by(gender, stress, size = 1)<\/code><\/pre>\n<pre><code>## # A tibble: 6 x 5\r\n##      id gender stress      t1    t2\r\n##   &lt;int&gt; &lt;fct&gt;  &lt;fct&gt;    &lt;dbl&gt; &lt;dbl&gt;\r\n## 1     3 male   low       5.63  5.47\r\n## 2    18 male   moderate  5.57  5.78\r\n## 3    25 male   high      5.48  5.74\r\n## 4    39 female low       5.50  5.66\r\n## 5    50 female moderate  5.96  5.32\r\n## 6    51 female high      5.59  5.06<\/code><\/pre>\n<pre class=\"r\"><code># Gather the columns t1, t2 and t3 into long format.\r\n# Convert id and time into factor variables\r\nperformance &lt;- performance %&gt;%\r\n  gather(key = \"time\", value = \"score\", t1, t2) %&gt;%\r\n  convert_as_factor(id, time)\r\n# Inspect some random rows of the data by groups\r\nset.seed(123)\r\nperformance %&gt;% sample_n_by(gender, stress, time, size = 1)<\/code><\/pre>\n<pre><code>## # A tibble: 12 x 5\r\n##   id    gender stress   time  score\r\n##   &lt;fct&gt; &lt;fct&gt;  &lt;fct&gt;    &lt;fct&gt; &lt;dbl&gt;\r\n## 1 3     male   low      t1     5.63\r\n## 2 8     male   low      t2     5.92\r\n## 3 15    male   moderate t1     5.96\r\n## 4 19    male   moderate t2     5.76\r\n## 5 30    male   high     t1     5.38\r\n## 6 21    male   high     t2     5.64\r\n## # \u2026 with 6 more rows<\/code><\/pre>\n<\/div>\n<div id=\"summary-statistics-1\" class=\"section level3\">\n<h3>Summary statistics<\/h3>\n<p>Group the data by <code>gender<\/code>, <code>stress<\/code> and <code>time<\/code>, and then compute some summary statistics of the <code>score<\/code> variable: mean and sd (standard deviation)<\/p>\n<pre class=\"r\"><code>performance %&gt;%\r\n  group_by(gender, stress, time ) %&gt;%\r\n  get_summary_stats(score, type = \"mean_sd\")<\/code><\/pre>\n<pre><code>## # A tibble: 12 x 7\r\n##   gender stress   time  variable     n  mean    sd\r\n##   &lt;fct&gt;  &lt;fct&gt;    &lt;fct&gt; &lt;chr&gt;    &lt;dbl&gt; &lt;dbl&gt; &lt;dbl&gt;\r\n## 1 male   low      t1    score       10  5.72 0.19 \r\n## 2 male   low      t2    score       10  5.70 0.143\r\n## 3 male   moderate t1    score       10  5.72 0.193\r\n## 4 male   moderate t2    score       10  5.77 0.155\r\n## 5 male   high     t1    score       10  5.48 0.121\r\n## 6 male   high     t2    score       10  5.64 0.195\r\n## # \u2026 with 6 more rows<\/code><\/pre>\n<\/div>\n<div id=\"visualization-1\" class=\"section level3\">\n<h3>Visualization<\/h3>\n<p>Create box plots of performance <code>score<\/code> by <code>gender<\/code> colored by <code>stress<\/code> levels and faceted by <code>time<\/code>:<\/p>\n<pre class=\"r\"><code>bxp &lt;- ggboxplot(\r\n  performance, x = \"gender\", y = \"score\",\r\n  color = \"stress\", palette = \"jco\",\r\n  facet.by =  \"time\"\r\n  )\r\nbxp<\/code><\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/dn-tutorials\/r-statistics-2-comparing-groups-means\/figures\/047-mixed-anova-three-way-bbw-boxplots-1.png\" width=\"576\" \/><\/p>\n<\/div>\n<div id=\"check-assumptions-1\" class=\"section level3\">\n<h3>Check assumptions<\/h3>\n<div id=\"outliers-1\" class=\"section level4\">\n<h4>Outliers<\/h4>\n<pre class=\"r\"><code>performance %&gt;%\r\n  group_by(gender, stress, time) %&gt;%\r\n  identify_outliers(score)<\/code><\/pre>\n<pre><code>## # A tibble: 1 x 7\r\n##   gender stress time  id    score is.outlier is.extreme\r\n##   &lt;fct&gt;  &lt;fct&gt;  &lt;fct&gt; &lt;fct&gt; &lt;dbl&gt; &lt;lgl&gt;      &lt;lgl&gt;     \r\n## 1 female low    t2    36     6.15 TRUE       FALSE<\/code><\/pre>\n<div class=\"success\">\n<p>There were no extreme outliers.<\/p>\n<\/div>\n<\/div>\n<div id=\"normality-assumption-1\" class=\"section level4\">\n<h4>Normality assumption<\/h4>\n<p>Compute Shapiro-Wilk test for each combinations of factor levels:<\/p>\n<pre class=\"r\"><code>performance %&gt;%\r\n  group_by(gender, stress, time ) %&gt;%\r\n  shapiro_test(score)<\/code><\/pre>\n<pre><code>## # A tibble: 12 x 6\r\n##   gender stress   time  variable statistic      p\r\n##   &lt;fct&gt;  &lt;fct&gt;    &lt;fct&gt; &lt;chr&gt;        &lt;dbl&gt;  &lt;dbl&gt;\r\n## 1 male   low      t1    score        0.942 0.574 \r\n## 2 male   low      t2    score        0.966 0.849 \r\n## 3 male   moderate t1    score        0.848 0.0547\r\n## 4 male   moderate t2    score        0.958 0.761 \r\n## 5 male   high     t1    score        0.915 0.319 \r\n## 6 male   high     t2    score        0.925 0.403 \r\n## # \u2026 with 6 more rows<\/code><\/pre>\n<div class=\"success\">\n<p>The score were normally distributed (p &gt; 0.05) for each cell, as assessed by Shapiro-Wilk\u2019s test of normality.<\/p>\n<\/div>\n<p>Create QQ plot for each cell of design:<\/p>\n<pre class=\"r\"><code>ggqqplot(performance, \"score\", ggtheme = theme_bw()) +\r\n  facet_grid(time ~ stress, labeller = \"label_both\")<\/code><\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/dn-tutorials\/r-statistics-2-comparing-groups-means\/figures\/047-mixed-anova-three-way-bbw-qq-plot-1.png\" width=\"480\" \/><\/p>\n<div class=\"success\">\n<p>All the points fall approximately along the reference line, for each cell. So we can assume normality of the data.<\/p>\n<\/div>\n<\/div>\n<div id=\"homogneity-of-variance-assumption-1\" class=\"section level4\">\n<h4>Homogneity of variance assumption<\/h4>\n<p>Compute the Levene\u2019s test at each level of the within-subjects factor, here <code>time<\/code> variable:<\/p>\n<pre class=\"r\"><code>performance %&gt;%\r\n  group_by(time) %&gt;%\r\n  levene_test(score ~ gender*stress)<\/code><\/pre>\n<pre><code>## # A tibble: 2 x 5\r\n##   time    df1   df2 statistic     p\r\n##   &lt;fct&gt; &lt;int&gt; &lt;int&gt;     &lt;dbl&gt; &lt;dbl&gt;\r\n## 1 t1        5    54     0.974 0.442\r\n## 2 t2        5    54     0.722 0.610<\/code><\/pre>\n<div class=\"success\">\n<p>There was homogeneity of variances, as assessed by Levene\u2019s test of homogeneity of variance (p &gt; .05).<\/p>\n<\/div>\n<\/div>\n<div id=\"assumption-of-sphericity-1\" class=\"section level4\">\n<h4>Assumption of sphericity<\/h4>\n<p>As mentioned in the two-way mixed ANOVA section, the Mauchly\u2019s test of sphericity and the sphericity corrections are internally done using the R function <code>anova_test()<\/code> and <code>get_anova_table()<\/code> [rstatix package].<\/p>\n<\/div>\n<\/div>\n<div id=\"computation-1\" class=\"section level3\">\n<h3>Computation<\/h3>\n<pre class=\"r\"><code>res.aov &lt;- anova_test(\r\n  data = performance, dv = score, wid = id,\r\n  within = time, between = c(gender, stress)\r\n  )\r\nget_anova_table(res.aov)<\/code><\/pre>\n<pre><code>## ANOVA Table (type II tests)\r\n## \r\n##               Effect DFn DFd      F        p p&lt;.05      ges\r\n## 1             gender   1  54  2.406 1.27e-01       0.023000\r\n## 2             stress   2  54 21.166 1.63e-07     * 0.288000\r\n## 3               time   1  54  0.063 8.03e-01       0.000564\r\n## 4      gender:stress   2  54  1.554 2.21e-01       0.029000\r\n## 5        gender:time   1  54  4.730 3.40e-02     * 0.041000\r\n## 6        stress:time   2  54  1.821 1.72e-01       0.032000\r\n## 7 gender:stress:time   2  54  6.101 4.00e-03     * 0.098000<\/code><\/pre>\n<div class=\"success\">\n<p>There was a statistically significant three-way interaction between time, gender, and stress F(2, 54) = 6.10, p = 0.004.<\/p>\n<\/div>\n<\/div>\n<div id=\"post-hoc-tests-1\" class=\"section level3\">\n<h3>Post-hoc tests<\/h3>\n<p><strong>If there is a significant three-way interaction effect<\/strong>, you can decompose it into:<\/p>\n<ul>\n<li><strong>Simple two-way interaction<\/strong>: run two-way interaction at each level of third variable,<\/li>\n<li><strong>Simple simple main effect<\/strong>: run one-way model at each level of second variable, and<\/li>\n<li><strong>simple simple pairwise comparisons<\/strong>: run pairwise or other post-hoc comparisons if necessary.<\/li>\n<\/ul>\n<p><strong>If you do not have a statistically significant three-way interaction<\/strong>, you need to determine whether you have any statistically significant two-way interaction from the ANOVA output. A significant two-way interaction can be followed up by a simple main effect analysis, which can be followed up by simple pairwise comparisons if significant.<\/p>\n<p>In this section we\u2019ll describe the procedure for a significant three-way interaction.<\/p>\n<div id=\"compute-simple-two-way-interaction\" class=\"section level4\">\n<h4>Compute simple two-way interaction<\/h4>\n<p>You are free to decide which two variables will form the simple two-way interactions and which variable will act as the third variable.<\/p>\n<p>In the following R code, we have considered the simple two-way interaction of <code>gender*stress<\/code> at each level of <code>time<\/code>.<\/p>\n<p>Group the data by <code>time<\/code> (the within-subject factors) and analyze the simple two-way interaction between <code>gender<\/code> and <code>stress<\/code>, which are the between-subjects factors.<\/p>\n<pre class=\"r\"><code># two-way interaction at each time levels\r\ntwo.way &lt;- performance %&gt;%\r\n  group_by(time) %&gt;%\r\n  anova_test(dv = score, wid = id, between = c(gender, stress))\r\ntwo.way<\/code><\/pre>\n<pre><code>## # A tibble: 6 x 8\r\n##   time  Effect          DFn   DFd      F          p `p&lt;.05`   ges\r\n##   &lt;fct&gt; &lt;chr&gt;         &lt;dbl&gt; &lt;dbl&gt;  &lt;dbl&gt;      &lt;dbl&gt; &lt;chr&gt;   &lt;dbl&gt;\r\n## 1 t1    gender            1    54  0.186 0.668      \"\"      0.003\r\n## 2 t1    stress            2    54 14.9   0.00000723 *       0.355\r\n## 3 t1    gender:stress     2    54  2.12  0.131      \"\"      0.073\r\n## 4 t2    gender            1    54  5.97  0.018      *       0.1  \r\n## 5 t2    stress            2    54  9.60  0.000271   *       0.262\r\n## 6 t2    gender:stress     2    54  4.95  0.011      *       0.155<\/code><\/pre>\n<div class=\"success\">\n<p>There was a statistically significant simple two-way interaction of gender and stress at t2, F(2, 54) = 4.95, p = 0.011, but not at t1, F(2, 54) = 2.12, p = 0.13.<\/p>\n<\/div>\n<div class=\"warning\">\n<p>Note that, statistical significance of a simple two-way interaction was accepted at a Bonferroni-adjusted alpha level of 0.025. This corresponds to the current level you declare statistical significance at (i.e., p &lt; 0.05) divided by the number of simple two-way interaction you are computing (i.e., 2).<\/p>\n<\/div>\n<\/div>\n<div id=\"compute-simple-simple-main-effects\" class=\"section level4\">\n<h4>Compute simple simple main effects<\/h4>\n<p>A statistically significant simple two-way interaction can be followed up with <strong>simple simple main effects<\/strong>.<\/p>\n<p>In our example, you could therefore investigate the effect of <code>stress<\/code> on the performance score at every level of <code>gender<\/code> or investigate the effect of <code>gender<\/code> at every level of <code>stress<\/code>.<\/p>\n<div class=\"warning\">\n<p>Note that, you will only need to do this for the simple two-way interaction for \u201ct2\u201d as this was the only simple two-way interaction that was statistically significant.<\/p>\n<\/div>\n<p>Group the data by <code>time<\/code> and <code>gender<\/code>, and analyze the simple main effect of <code>stress<\/code> on performance score:<\/p>\n<pre class=\"r\"><code>stress.effect &lt;- performance %&gt;%\r\n  group_by(time, gender) %&gt;%\r\n  anova_test(dv = score, wid = id, between = stress)\r\nstress.effect %&gt;% filter(time == \"t2\")<\/code><\/pre>\n<pre><code>## # A tibble: 2 x 9\r\n##   gender time  Effect   DFn   DFd     F        p `p&lt;.05`   ges\r\n##   &lt;fct&gt;  &lt;fct&gt; &lt;chr&gt;  &lt;dbl&gt; &lt;dbl&gt; &lt;dbl&gt;    &lt;dbl&gt; &lt;chr&gt;   &lt;dbl&gt;\r\n## 1 male   t2    stress     2    27  1.57 0.227    \"\"      0.104\r\n## 2 female t2    stress     2    27 10.5  0.000416 *       0.438<\/code><\/pre>\n<div class=\"warning\">\n<p>In the table above, we only need the results for <code>time = t2<\/code>. Statistical significance of a simple simple main effect was accepted at a Bonferroni-adjusted alpha level of 0.025, that is 0.05 divided y the number of simple simple main effects you are computing (i.e., 2).<\/p>\n<\/div>\n<div class=\"success\">\n<p>There was a statistically significant simple simple main effect of stress on the performance score for female at t2 time point, F(2, 27) = 10.5, p = 0.0004, but not for males, F(2, 27) = 1.57, p = 0.23.<\/p>\n<\/div>\n<\/div>\n<div id=\"compute-simple-simple-comparisons\" class=\"section level4\">\n<h4>Compute simple simple comparisons<\/h4>\n<p>A statistically significant simple simple main effect can be followed up by <strong>multiple pairwise comparisons<\/strong> to determine which group means are different.<\/p>\n<div class=\"warning\">\n<p>Note that, you will only need to concentrate on the pairwise comparison results for female, because the effect of stress was significant for female only in the previous section.<\/p>\n<\/div>\n<p>Group the data by <code>time<\/code> and <code>gender<\/code>, and perform pairwise comparisons between <code>stress<\/code> levels with Bonferroni adjustment:<\/p>\n<pre class=\"r\"><code># Fit pairwise comparisons\r\npwc &lt;- performance %&gt;%\r\n  group_by(time, gender) %&gt;%\r\n  pairwise_t_test(score ~ stress, p.adjust.method = \"bonferroni\") %&gt;%\r\n  select(-p, -p.signif) # Remove details\r\n# Focus on the results of \"female\" at t2\r\npwc %&gt;% filter(time == \"t2\", gender == \"female\")<\/code><\/pre>\n<pre><code>## # A tibble: 3 x 9\r\n##   gender time  .y.   group1   group2      n1    n2    p.adj p.adj.signif\r\n##   &lt;fct&gt;  &lt;fct&gt; &lt;chr&gt; &lt;chr&gt;    &lt;chr&gt;    &lt;int&gt; &lt;int&gt;    &lt;dbl&gt; &lt;chr&gt;       \r\n## 1 female t2    score low      moderate    10    10 0.323    ns          \r\n## 2 female t2    score low      high        10    10 0.000318 ***         \r\n## 3 female t2    score moderate high        10    10 0.0235   *<\/code><\/pre>\n<div class=\"success\">\n<p>For female, the mean performance score was statistically significantly different between low and high stress levels (p &lt; 0.001) and between moderate and high stress levels (p = 0.023).<\/p>\n<p>There was no significant difference between low and moderate stress groups (p = 0.32)<\/p>\n<\/div>\n<\/div>\n<\/div>\n<div id=\"report-1\" class=\"section level3\">\n<h3>Report<\/h3>\n<p>A three-way mixed ANOVA was performed to evaluate the effects of gender, stress and time on performance score.<\/p>\n<p>There were no extreme outliers, as assessed by box plot method. The data was normally distributed, as assessed by Shapiro-Wilk\u2019s test of normality (p &gt; 0.05). There was homogeneity of variances (p &gt; 0.05) as assessed by Levene\u2019s test of homogeneity of variances.<\/p>\n<p>There was a statistically significant three-way interaction between gender, stress and time, F(2, 54) = 6.10, p = 0.004.<\/p>\n<p>For the simple two-way interactions and simple simple main effects, a Bonferroni adjustment was applied leading to statistical significance being accepted at the p &lt; 0.025 level.<\/p>\n<p>There was a statistically significant simple two-way interaction between gender and stress at time point t2, F(2, 54) = 4.95, p = 0.011, but not at t1, F(2, 54) = 2.12, p = 0.13.<\/p>\n<p>There was a statistically significant simple simple main effect of stress on the performance score for female at t2 time point, F(2, 27) = 10.5, p = 0.0004, but not for males, F(2, 27) = 1.57, p = 0.23.<\/p>\n<p>All simple simple pairwise comparisons were run between the different stress groups for female at time point t2. A Bonferroni adjustment was applied.<\/p>\n<p>The mean performance score was statistically significantly different between low and high stress levels (p &lt; 0.001) and between moderate and high stress levels (p = 0.024). There was no significant difference between low and moderate stress groups (p = 0.32).<\/p>\n<pre class=\"r\"><code># Visualization: box plots with p-values\r\npwc &lt;- pwc %&gt;% add_xy_position(x = \"gender\") \r\npwc.filtered &lt;- pwc %&gt;% filter(time == \"t2\", gender == \"female\")\r\nbxp + \r\n  stat_pvalue_manual(pwc.filtered, tip.length = 0, hide.ns = TRUE) +\r\n  labs(\r\n    subtitle = get_test_label(res.aov, detailed = TRUE),\r\n    caption = get_pwc_label(pwc)\r\n  )<\/code><\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/dn-tutorials\/r-statistics-2-comparing-groups-means\/figures\/047-mixed-anova-three-way-bbw-boxplots-with-p-values-1.png\" width=\"576\" \/><\/p>\n<\/div>\n<\/div>\n<div id=\"three-way-bww-b\" class=\"section level2\">\n<h2>Three-way Mixed ANOVA: 1 between- and 2 within-subjects factors<\/h2>\n<p>This section describes how to compute the three-way mixed ANOVA, in R, for a situation where you have <strong>one between-subjects factor and two within-subjects factors<\/strong>. For example, you might want to understand how weight loss score differs in individuals doing exercises vs not doing exercises over three <code>time<\/code> points (t1, t2, t3) depending on participant diets (<code>diet:no<\/code> and <code>diet:yes<\/code>).<\/p>\n<div id=\"data-preparation-2\" class=\"section level3\">\n<h3>Data preparation<\/h3>\n<p>We\u2019ll use the <code>weightloss<\/code> dataset available in the datarium package. This dataset was originally created for three-way repeated measures ANOVA. However, for our example in this article, we\u2019ll modify slightly the data so that it corresponds to a three-way mixed design.<\/p>\n<p>A researcher wanted to assess the effect of <code>time<\/code> on weight loss score depending on <code>exercises<\/code> programs and <code>diet<\/code>.<\/p>\n<p>The weight loss score was measured in two different groups: a group of participants doing exercises (<code>exercises:yes<\/code>) and in another group not doing exercises (<code>excises:no<\/code>).<\/p>\n<p>Each participant was also enrolled in two trials: (1) no diet and (2) diet. The order of the trials was counterbalanced and sufficient time was allowed between trials to allow any effects of previous trials to have dissipated.<\/p>\n<p>Each trial lasted 9 weeks and the weight loss score was measured at the beginning of each trial (t1), at the midpoint of each trial (t2) and at the end of each trial (t3).<\/p>\n<p>In this study design, 24 individuals were recruited. Of these 24 participants, 12 belongs to the <code>exercises:no<\/code> group and 12 were in <code>exercises:yes<\/code> group. The 24 participants were enrolled in two successive trials (<code>diet:no<\/code> and <code>diet:yes<\/code>) and the weight loss <code>score<\/code> was repeatedly measured at three time points.<\/p>\n<p>In this setting, we have:<\/p>\n<ul>\n<li>one dependent (or outcome) variable: <code>score<\/code><\/li>\n<li>One between-subjects factor: <code>exercises<\/code><\/li>\n<li>two within-subjects factors: <code>diet<\/code> end <code>time<\/code><\/li>\n<\/ul>\n<div class=\"block\">\n<p>Three-way mixed ANOVA can be performed in order to determine whether there is a significant interaction between diet, exercises and time on the weight loss score.<\/p>\n<\/div>\n<p>Load the data and inspect some random rows by group:<\/p>\n<pre class=\"r\"><code># Load the original data\r\n# Wide format\r\ndata(\"weightloss\", package = \"datarium\")\r\n# Modify it to have three-way mixed design\r\nweightloss &lt;- weightloss %&gt;%\r\n  mutate(id = rep(1:24, 2)) # two trials\r\n# Show one random row by group\r\nset.seed(123)\r\nweightloss %&gt;% sample_n_by(diet, exercises, size = 1)<\/code><\/pre>\n<pre><code>## # A tibble: 4 x 6\r\n##      id diet  exercises    t1    t2    t3\r\n##   &lt;int&gt; &lt;fct&gt; &lt;fct&gt;     &lt;dbl&gt; &lt;dbl&gt; &lt;dbl&gt;\r\n## 1     4 no    no         11.1   9.5  11.1\r\n## 2    22 no    yes        10.2  11.8  17.4\r\n## 3     5 yes   no         11.6  13.4  13.9\r\n## 4    23 yes   yes        12.7  12.7  15.1<\/code><\/pre>\n<pre class=\"r\"><code># Gather the columns t1, t2 and t3 into long format.\r\n# Convert id and time into factor variables\r\nweightloss &lt;- weightloss %&gt;%\r\n  gather(key = \"time\", value = \"score\", t1, t2, t3) %&gt;%\r\n  convert_as_factor(id, time)\r\n# Inspect some random rows of the data by groups\r\nset.seed(123)\r\nweightloss %&gt;% sample_n_by(diet, exercises, time, size = 1)<\/code><\/pre>\n<pre><code>## # A tibble: 12 x 5\r\n##   id    diet  exercises time  score\r\n##   &lt;fct&gt; &lt;fct&gt; &lt;fct&gt;     &lt;fct&gt; &lt;dbl&gt;\r\n## 1 4     no    no        t1     11.1\r\n## 2 10    no    no        t2     10.7\r\n## 3 5     no    no        t3     12.3\r\n## 4 23    no    yes       t1     10.2\r\n## 5 24    no    yes       t2     13.2\r\n## 6 13    no    yes       t3     15.8\r\n## # \u2026 with 6 more rows<\/code><\/pre>\n<\/div>\n<div id=\"summary-statistics-2\" class=\"section level3\">\n<h3>Summary statistics<\/h3>\n<p>Group the data by <code>exercises<\/code>, <code>diet<\/code> and <code>time<\/code>, and then compute some summary statistics of the <code>score<\/code> variable: mean and sd (standard deviation)<\/p>\n<pre class=\"r\"><code>weightloss %&gt;%\r\n  group_by(exercises, diet, time) %&gt;%\r\n  get_summary_stats(score, type = \"mean_sd\")<\/code><\/pre>\n<pre><code>## # A tibble: 12 x 7\r\n##   diet  exercises time  variable     n  mean    sd\r\n##   &lt;fct&gt; &lt;fct&gt;     &lt;fct&gt; &lt;chr&gt;    &lt;dbl&gt; &lt;dbl&gt; &lt;dbl&gt;\r\n## 1 no    no        t1    score       12  10.9 0.868\r\n## 2 no    no        t2    score       12  11.6 1.30 \r\n## 3 no    no        t3    score       12  11.4 0.935\r\n## 4 yes   no        t1    score       12  11.7 0.938\r\n## 5 yes   no        t2    score       12  12.4 1.42 \r\n## 6 yes   no        t3    score       12  13.8 1.43 \r\n## # \u2026 with 6 more rows<\/code><\/pre>\n<\/div>\n<div id=\"visualization-2\" class=\"section level3\">\n<h3>Visualization<\/h3>\n<p>Create box plots of weightloss <code>score<\/code> by <code>exercises<\/code> groups, colored by <code>time<\/code> points and faceted by <code>diet<\/code> trials:<\/p>\n<pre class=\"r\"><code>bxp &lt;- ggboxplot(\r\n  weightloss, x = \"exercises\", y = \"score\",\r\n  color = \"time\", palette = \"jco\",\r\n  facet.by = \"diet\", short.panel.labs = FALSE\r\n  )\r\nbxp<\/code><\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/dn-tutorials\/r-statistics-2-comparing-groups-means\/figures\/047-mixed-anova-three-way-bww-boxplots-1.png\" width=\"576\" \/><\/p>\n<\/div>\n<div id=\"check-assumptions-2\" class=\"section level3\">\n<h3>Check assumptions<\/h3>\n<div id=\"outliers-2\" class=\"section level4\">\n<h4>Outliers<\/h4>\n<pre class=\"r\"><code>weightloss %&gt;%\r\n  group_by(diet, exercises, time) %&gt;%\r\n  identify_outliers(score)<\/code><\/pre>\n<pre><code>## # A tibble: 5 x 7\r\n##   diet  exercises time  id    score is.outlier is.extreme\r\n##   &lt;fct&gt; &lt;fct&gt;     &lt;fct&gt; &lt;fct&gt; &lt;dbl&gt; &lt;lgl&gt;      &lt;lgl&gt;     \r\n## 1 no    no        t3    2      13.2 TRUE       FALSE     \r\n## 2 yes   no        t1    1      10.2 TRUE       FALSE     \r\n## 3 yes   no        t1    3      13.2 TRUE       FALSE     \r\n## 4 yes   no        t1    4      10.2 TRUE       FALSE     \r\n## 5 yes   no        t2    10     15.3 TRUE       FALSE<\/code><\/pre>\n<div class=\"success\">\n<p>There were no extreme outliers.<\/p>\n<\/div>\n<\/div>\n<div id=\"normality-assumption-2\" class=\"section level4\">\n<h4>Normality assumption<\/h4>\n<p>Compute Shapiro-Wilk test for each combinations of factor levels:<\/p>\n<pre class=\"r\"><code>weightloss %&gt;%\r\n  group_by(diet, exercises, time) %&gt;%\r\n  shapiro_test(score)<\/code><\/pre>\n<pre><code>## # A tibble: 12 x 6\r\n##   diet  exercises time  variable statistic     p\r\n##   &lt;fct&gt; &lt;fct&gt;     &lt;fct&gt; &lt;chr&gt;        &lt;dbl&gt; &lt;dbl&gt;\r\n## 1 no    no        t1    score        0.917 0.264\r\n## 2 no    no        t2    score        0.957 0.743\r\n## 3 no    no        t3    score        0.965 0.851\r\n## 4 no    yes       t1    score        0.922 0.306\r\n## 5 no    yes       t2    score        0.912 0.229\r\n## 6 no    yes       t3    score        0.953 0.674\r\n## # \u2026 with 6 more rows<\/code><\/pre>\n<div class=\"success\">\n<p>The weight loss score was normally distributed (p &gt; 0.05), as assessed by Shapiro-Wilk\u2019s test of normality.<\/p>\n<\/div>\n<p>Create QQ plot for each cell of design:<\/p>\n<pre class=\"r\"><code>ggqqplot(weightloss, \"score\", ggtheme = theme_bw()) +\r\n  facet_grid(diet + exercises ~ time, labeller = \"label_both\")<\/code><\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/dn-tutorials\/r-statistics-2-comparing-groups-means\/figures\/047-mixed-anova-three-way-bww-qq-plot-1.png\" width=\"576\" \/><\/p>\n<div class=\"success\">\n<p>From the plot above, as all the points fall approximately along this reference line, we can assume normality.<\/p>\n<\/div>\n<\/div>\n<div id=\"homogneity-of-variance-assumption-2\" class=\"section level4\">\n<h4>Homogneity of variance assumption<\/h4>\n<p>Compute the Levene\u2019s test after grouping the data by <code>diet<\/code> and <code>time<\/code> categories:<\/p>\n<pre class=\"r\"><code>weightloss %&gt;%\r\n  group_by(diet, time) %&gt;%\r\n  levene_test(score ~ exercises)<\/code><\/pre>\n<pre><code>## # A tibble: 6 x 6\r\n##   diet  time    df1   df2 statistic      p\r\n##   &lt;fct&gt; &lt;fct&gt; &lt;int&gt; &lt;int&gt;     &lt;dbl&gt;  &lt;dbl&gt;\r\n## 1 no    t1        1    22    2.44   0.132 \r\n## 2 no    t2        1    22    0.691  0.415 \r\n## 3 no    t3        1    22    2.87   0.105 \r\n## 4 yes   t1        1    22    0.376  0.546 \r\n## 5 yes   t2        1    22    0.0574 0.813 \r\n## 6 yes   t3        1    22    5.14   0.0336<\/code><\/pre>\n<div class=\"success\">\n<p>There was homogeneity of variances for all cells (p &gt; 0.05), except for the condition diet:yes at time:t3 (p = 0.034), as assessed by Levene\u2019s test of homogeneity of variance.<\/p>\n<\/div>\n<div class=\"warning\">\n<p>Note that, if you do not have homogeneity of variances, you can try to transform the outcome (dependent) variable to correct for the unequal variances.<\/p>\n<p>If group sample sizes are (approximately) equal, run the three-way mixed ANOVA anyway because it is somewhat robust to heterogeneity of variance in these circumstances.<\/p>\n<p>It\u2019s also possible to perform robust ANOVA test using the WRS2 R package.<\/p>\n<\/div>\n<\/div>\n<div id=\"assumption-of-sphericity-2\" class=\"section level4\">\n<h4>Assumption of sphericity<\/h4>\n<p>As mentioned in the two-way mixed ANOVA section, the Mauchly\u2019s test of sphericity and the sphericity corrections are internally done using the R function <code>anova_test()<\/code> and <code>get_anova_table()<\/code> [rstatix package].<\/p>\n<\/div>\n<\/div>\n<div id=\"computation-2\" class=\"section level3\">\n<h3>Computation<\/h3>\n<pre class=\"r\"><code>res.aov &lt;- anova_test(\r\n  data = weightloss, dv = score, wid = id,\r\n  between = exercises, within = c(diet, time)\r\n  )\r\nget_anova_table(res.aov)<\/code><\/pre>\n<pre><code>## ANOVA Table (type II tests)\r\n## \r\n##                Effect DFn DFd      F        p p&lt;.05   ges\r\n## 1           exercises   1  22 38.771 2.88e-06     * 0.284\r\n## 2                diet   1  22  7.912 1.00e-02     * 0.028\r\n## 3                time   2  44 82.199 1.38e-15     * 0.541\r\n## 4      exercises:diet   1  22 51.698 3.31e-07     * 0.157\r\n## 5      exercises:time   2  44 26.222 3.18e-08     * 0.274\r\n## 6           diet:time   2  44  0.784 4.63e-01       0.013\r\n## 7 exercises:diet:time   2  44  9.966 2.69e-04     * 0.147<\/code><\/pre>\n<div class=\"success\">\n<p>From the output above, it can be seen that there is a statistically significant three-way interactions between exercises, diet and time F(2, 44) = 9.96, p = 0.00027.<\/p>\n<\/div>\n<div class=\"warning\">\n<p>Note that, if the three-way interaction is not statistically significant, you need to consult the two-way interactions in the output.<\/p>\n<p>In our example, there was a statistically significant two-way exercises*diet interaction (p &lt; 0.0001), and two-way exercises*time (p &lt; 0.0001). The two-way diet*time interaction was not statistically significant (p = 0.46).<\/p>\n<\/div>\n<\/div>\n<div id=\"post-hoc-tests-2\" class=\"section level3\">\n<h3>Post-hoc tests<\/h3>\n<p><strong>If there is a significant three-way interaction effect<\/strong>, you can decompose it into:<\/p>\n<ul>\n<li><strong>Simple two-way interaction<\/strong>: run two-way interaction at each level of third variable,<\/li>\n<li><strong>Simple simple main effect<\/strong>: run one-way model at each level of second variable, and<\/li>\n<li><strong>simple simple pairwise comparisons<\/strong>: run pairwise or other post-hoc comparisons if necessary.<\/li>\n<\/ul>\n<p><strong>If you do not have a statistically significant three-way interaction<\/strong>, you need to determine whether you have any statistically significant two-way interaction from the ANOVA output. You can follow up a significant two-way interaction by simple main effects analyses and pairwise comparisons between groups if necessary.<\/p>\n<p>In this section we\u2019ll describe the procedure for a significant three-way interaction.<\/p>\n<div id=\"compute-simple-two-way-interaction-1\" class=\"section level4\">\n<h4>Compute simple two-way interaction<\/h4>\n<p>In this example, we\u2019ll consider the <code>diet*time<\/code> interaction at each level of <code>exercises<\/code>. Group the data by <code>exercises<\/code> and analyze the simple two-way interaction between <code>diet<\/code> and <code>time<\/code>:<\/p>\n<pre class=\"r\"><code># Two-way ANOVA at each exercises group level\r\ntwo.way &lt;- weightloss %&gt;%\r\n  group_by(exercises) %&gt;%\r\n  anova_test(dv = score, wid = id, within = c(diet, time))\r\ntwo.way<\/code><\/pre>\n<pre><code>## # A tibble: 2 x 2\r\n##   exercises anova     \r\n##   &lt;fct&gt;     &lt;list&gt;    \r\n## 1 no        &lt;anov_tst&gt;\r\n## 2 yes       &lt;anov_tst&gt;<\/code><\/pre>\n<pre class=\"r\"><code># Extract anova table\r\nget_anova_table(two.way)<\/code><\/pre>\n<pre><code>## # A tibble: 6 x 8\r\n##   exercises Effect      DFn   DFd      F        p `p&lt;.05`   ges\r\n##   &lt;fct&gt;     &lt;chr&gt;     &lt;dbl&gt; &lt;dbl&gt;  &lt;dbl&gt;    &lt;dbl&gt; &lt;chr&gt;   &lt;dbl&gt;\r\n## 1 no        diet          1    11  56.4  1.18e- 5 *       0.262\r\n## 2 no        time          2    22   5.90 9.00e- 3 *       0.181\r\n## 3 no        diet:time     2    22   2.91 7.60e- 2 \"\"      0.09 \r\n## 4 yes       diet          1    11   8.60 1.40e- 2 *       0.066\r\n## 5 yes       time          2    22 148.   1.73e-13 *       0.746\r\n## 6 yes       diet:time     2    22   7.81 3.00e- 3 *       0.216<\/code><\/pre>\n<div class=\"success\">\n<p>There was a statistically significant simple two-way interaction between diet and time for exercises:yes group, F(2, 22) = 7.81, p = 0.0027, but not for exercises:no group, F(2, 22) = 2.91, p = 0.075.<\/p>\n<\/div>\n<div class=\"warning\">\n<p>Note that, statistical significance of a simple two-way interaction was accepted at a Bonferroni-adjusted alpha level of 0.025. This corresponds to the current level you declare statistical significance at (i.e., p &lt; 0.05) divided by the number of simple two-way interaction you are computing (i.e., 2).<\/p>\n<\/div>\n<\/div>\n<div id=\"compute-simple-simple-main-effect\" class=\"section level4\">\n<h4>Compute simple simple main effect<\/h4>\n<p>A statistically significant simple two-way interaction can be followed up with <strong>simple simple main effects<\/strong>.<\/p>\n<p>In our example, you could therefore investigate the effect of <code>time<\/code> on the weight loss score at every level of <code>diet<\/code> and\/or investigate the effect of <code>diet<\/code> at every level of <code>time<\/code>.<\/p>\n<div class=\"warning\">\n<p>Note that, you will only need to do this for the simple two-way interaction for \u201cexercises:yes\u201d group, as this was the only simple two-way interaction that was statistically significant.<\/p>\n<\/div>\n<p>Group the data by <code>exercises<\/code> and <code>diet<\/code>, and analyze the simple main effect of <code>time<\/code>:<\/p>\n<pre class=\"r\"><code>time.effect &lt;- weightloss %&gt;%\r\n  group_by(exercises, diet) %&gt;%\r\n  anova_test(dv = score, wid = id, within = time) %&gt;%\r\n  get_anova_table()\r\ntime.effect %&gt;% filter(exercises == \"yes\")<\/code><\/pre>\n<pre><code>## # A tibble: 2 x 9\r\n##   diet  exercises Effect   DFn   DFd     F        p `p&lt;.05`   ges\r\n##   &lt;fct&gt; &lt;fct&gt;     &lt;chr&gt;  &lt;dbl&gt; &lt;dbl&gt; &lt;dbl&gt;    &lt;dbl&gt; &lt;chr&gt;   &lt;dbl&gt;\r\n## 1 no    yes       time       2    22  78.8 9.30e-11 *       0.801\r\n## 2 yes   yes       time       2    22  30.9 4.06e- 7 *       0.655<\/code><\/pre>\n<div class=\"warning\">\n<p>In the table above, we only need the results for <code>exercises = yes<\/code>. Statistical significance of a simple simple main effect was accepted at a Bonferroni-adjusted alpha level of 0.025, that is 0.05 divided y the number of simple simple main effects you are computing (i.e., 2).<\/p>\n<\/div>\n<div class=\"success\">\n<p>The simple simple main effect of time on weight loss score was statistically significant under exercises condition for both diet:no (F(2,22) = 78.81, p &lt; 0.0001) and diet:yes (F(2, 22) = 30.92, p &lt; 0.0001) groups.<\/p>\n<\/div>\n<\/div>\n<div id=\"compute-simple-simple-comparisons-1\" class=\"section level4\">\n<h4>Compute simple simple comparisons<\/h4>\n<p>A statistically significant simple simple main effect can be followed up by <strong>multiple pairwise comparisons<\/strong> to determine which group means are different.<\/p>\n<div class=\"warning\">\n<p>Recall that, you will only need to concentrate on the pairwise comparison results for exercises:yes.<\/p>\n<\/div>\n<p>Group the data by <code>exercises<\/code> and <code>diet<\/code>, and perform pairwise comparisons between <code>time<\/code> points with Bonferroni adjustment. Paired t-test is used:<\/p>\n<pre class=\"r\"><code># compute pairwise comparisons\r\npwc &lt;- weightloss %&gt;%\r\n  group_by(exercises, diet) %&gt;%\r\n  pairwise_t_test(\r\n    score ~ time, paired = TRUE, \r\n    p.adjust.method = \"bonferroni\"\r\n    ) %&gt;%\r\n  select(-statistic, -df) # Remove details\r\n# Focus on the results of exercises:yes group\r\npwc %&gt;% filter(exercises == \"yes\") %&gt;%\r\n  select(-p)    # Remove p column<\/code><\/pre>\n<pre><code>## # A tibble: 6 x 9\r\n##   diet  exercises .y.   group1 group2    n1    n2        p.adj p.adj.signif\r\n##   &lt;fct&gt; &lt;fct&gt;     &lt;chr&gt; &lt;chr&gt;  &lt;chr&gt;  &lt;int&gt; &lt;int&gt;        &lt;dbl&gt; &lt;chr&gt;       \r\n## 1 no    yes       score t1     t2        12    12 0.000741     ***         \r\n## 2 no    yes       score t1     t3        12    12 0.0000000121 ****        \r\n## 3 no    yes       score t2     t3        12    12 0.000257     ***         \r\n## 4 yes   yes       score t1     t2        12    12 0.01         **          \r\n## 5 yes   yes       score t1     t3        12    12 0.00000124   ****        \r\n## 6 yes   yes       score t2     t3        12    12 0.02         *<\/code><\/pre>\n<div class=\"success\">\n<p>All simple simple pairwise comparisons were run between the different time points under exercises condition (i.e., exercises:yes) for diet:no and diet:yes trials. A Bonferroni adjustment was applied.<\/p>\n<p>The mean weight loss score was significantly different in all time point comparisons when exercises are performed (p &lt; 0.05).<\/p>\n<\/div>\n<\/div>\n<\/div>\n<div id=\"report-2\" class=\"section level3\">\n<h3>Report<\/h3>\n<p>A three-way mixed ANOVA was performed to evaluate the effects of diet, exercises and time on weight loss.<\/p>\n<p>There were no extreme outliers, as assessed by box plot method. The data was normally distributed, as assessed by Shapiro-Wilk\u2019s test of normality (p &gt; 0.05). There was homogeneity of variances (p &gt; 0.05) as assessed by Levene\u2019s test of homogeneity of variances. For the three-way interaction effect, Mauchly\u2019s test of sphericity indicated that the assumption of sphericity was met (p &gt; 0.05).<\/p>\n<p>There was a statistically significant three-way interaction between exercises, diet and time F(2, 44) = 9.96, p &lt; 0.001.<\/p>\n<p>For the simple two-way interactions and simple simple main effects, a Bonferroni adjustment was applied leading to statistical significance being accepted at the p &lt; 0.025 level.<\/p>\n<p>There was a statistically significant simple two-way interaction between diet and time for exercises:yes group, F(2, 22) = 7.81, p = 0.0027, but not for exercises:no group, F(2, 22) = 2.91, p = 0.075.<\/p>\n<p>The simple simple main effect of time on weight loss score was statistically significant under exercises condition for both diet:no (F(2,22) = 78.81, p &lt; 0.0001) and diet:yes (F(2, 22) = 30.92, p &lt; 0.0001) groups.<\/p>\n<p>All simple simple pairwise comparisons were run between the different time points under exercises condition (i.e., exercises:yes) for diet:no and diet:yes trials. A Bonferroni adjustment was applied. The mean weight loss score was significantly different in all time point comparisons when exercises are performed (p &lt; 0.05).<\/p>\n<pre class=\"r\"><code># Visualization: box plots with p-values\r\npwc &lt;- pwc %&gt;% add_xy_position(x = \"exercises\")\r\npwc.filtered &lt;- pwc %&gt;% filter(exercises == \"yes\")\r\nbxp + \r\n  stat_pvalue_manual(pwc.filtered, tip.length = 0, hide.ns = TRUE) +\r\n  labs(\r\n    subtitle = get_test_label(res.aov, detailed = TRUE),\r\n    caption = get_pwc_label(pwc)\r\n  )<\/code><\/pre>\n<p><img decoding=\"async\" src=\"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/dn-tutorials\/r-statistics-2-comparing-groups-means\/figures\/047-mixed-anova-three-way-bww-boxplot-with-p-values-1.png\" width=\"576\" \/><\/p>\n<\/div>\n<\/div>\n<div id=\"summary\" class=\"section level2\">\n<h2>Summary<\/h2>\n<p>This article describes how to compute and interpret mixed ANOVA in R. We also explain the assumptions made by mixed ANOVA tests and provide practical examples of R codes to check whether the test assumptions are met.<\/p>\n<\/div>\n<\/div>\n<p><!--end rdoc--><\/p>\n","protected":false},"excerpt":{"rendered":"<p>The Mixed ANOVA is used to compare the means of groups cross-classified by two different types of factor variables, including: i) between-subjects factors,  which have independent categories (e.g., gender: male\/female). ii) within-subjects factors, which have related categories also known as repeated measures (e.g., time: before\/after treatment).  This chapter describes how to compute and interpret the different mixed ANOVA tests in R.<\/p>\n","protected":false},"author":1,"featured_media":9015,"parent":0,"menu_order":0,"comment_status":"open","ping_status":"closed","template":"","class_list":["post-10875","dt_lessons","type-dt_lessons","status-publish","has-post-thumbnail","hentry"],"yoast_head":"<!-- This site is optimized with the Yoast SEO plugin v25.2 - https:\/\/yoast.com\/wordpress\/plugins\/seo\/ -->\n<title>Mixed ANOVA in R: The Ultimate Guide - Datanovia<\/title>\n<meta name=\"robots\" content=\"index, follow, max-snippet:-1, max-image-preview:large, max-video-preview:-1\" \/>\n<link rel=\"canonical\" href=\"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/\" \/>\n<meta property=\"og:locale\" content=\"en_US\" \/>\n<meta property=\"og:type\" content=\"article\" \/>\n<meta property=\"og:title\" content=\"Mixed ANOVA in R: The Ultimate Guide - Datanovia\" \/>\n<meta property=\"og:description\" content=\"The Mixed ANOVA is used to compare the means of groups cross-classified by two different types of factor variables, including: i) between-subjects factors, which have independent categories (e.g., gender: male\/female). ii) within-subjects factors, which have related categories also known as repeated measures (e.g., time: before\/after treatment). This chapter describes how to compute and interpret the different mixed ANOVA tests in R.\" \/>\n<meta property=\"og:url\" content=\"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/\" \/>\n<meta property=\"og:site_name\" content=\"Datanovia\" \/>\n<meta property=\"article:modified_time\" content=\"2019-11-29T06:41:02+00:00\" \/>\n<meta property=\"og:image\" content=\"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/2019\/05\/X39410210_704020863271194_6530535646589616128_n.jpg\" \/>\n\t<meta property=\"og:image:width\" content=\"1024\" \/>\n\t<meta property=\"og:image:height\" content=\"512\" \/>\n\t<meta property=\"og:image:type\" content=\"image\/jpeg\" \/>\n<meta name=\"twitter:card\" content=\"summary_large_image\" \/>\n<meta name=\"twitter:label1\" content=\"Est. reading time\" \/>\n\t<meta name=\"twitter:data1\" content=\"35 minutes\" \/>\n<script type=\"application\/ld+json\" class=\"yoast-schema-graph\">{\"@context\":\"https:\/\/schema.org\",\"@graph\":[{\"@type\":\"WebPage\",\"@id\":\"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/\",\"url\":\"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/\",\"name\":\"Mixed ANOVA in R: The Ultimate Guide - Datanovia\",\"isPartOf\":{\"@id\":\"https:\/\/www.datanovia.com\/en\/#website\"},\"primaryImageOfPage\":{\"@id\":\"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/#primaryimage\"},\"image\":{\"@id\":\"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/#primaryimage\"},\"thumbnailUrl\":\"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/2019\/05\/X39410210_704020863271194_6530535646589616128_n.jpg\",\"datePublished\":\"2019-11-29T06:37:37+00:00\",\"dateModified\":\"2019-11-29T06:41:02+00:00\",\"breadcrumb\":{\"@id\":\"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/#breadcrumb\"},\"inLanguage\":\"en-US\",\"potentialAction\":[{\"@type\":\"ReadAction\",\"target\":[\"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/\"]}]},{\"@type\":\"ImageObject\",\"inLanguage\":\"en-US\",\"@id\":\"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/#primaryimage\",\"url\":\"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/2019\/05\/X39410210_704020863271194_6530535646589616128_n.jpg\",\"contentUrl\":\"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/2019\/05\/X39410210_704020863271194_6530535646589616128_n.jpg\",\"width\":1024,\"height\":512},{\"@type\":\"BreadcrumbList\",\"@id\":\"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/#breadcrumb\",\"itemListElement\":[{\"@type\":\"ListItem\",\"position\":1,\"name\":\"Home\",\"item\":\"https:\/\/www.datanovia.com\/en\/\"},{\"@type\":\"ListItem\",\"position\":2,\"name\":\"Lessons\",\"item\":\"https:\/\/www.datanovia.com\/en\/lessons\/\"},{\"@type\":\"ListItem\",\"position\":3,\"name\":\"Mixed ANOVA in R\"}]},{\"@type\":\"WebSite\",\"@id\":\"https:\/\/www.datanovia.com\/en\/#website\",\"url\":\"https:\/\/www.datanovia.com\/en\/\",\"name\":\"Datanovia\",\"description\":\"Data Mining and Statistics for Decision Support\",\"publisher\":{\"@id\":\"https:\/\/www.datanovia.com\/en\/#organization\"},\"potentialAction\":[{\"@type\":\"SearchAction\",\"target\":{\"@type\":\"EntryPoint\",\"urlTemplate\":\"https:\/\/www.datanovia.com\/en\/?s={search_term_string}\"},\"query-input\":{\"@type\":\"PropertyValueSpecification\",\"valueRequired\":true,\"valueName\":\"search_term_string\"}}],\"inLanguage\":\"en-US\"},{\"@type\":\"Organization\",\"@id\":\"https:\/\/www.datanovia.com\/en\/#organization\",\"name\":\"Datanovia\",\"url\":\"https:\/\/www.datanovia.com\/en\/\",\"logo\":{\"@type\":\"ImageObject\",\"inLanguage\":\"en-US\",\"@id\":\"https:\/\/www.datanovia.com\/en\/#\/schema\/logo\/image\/\",\"url\":\"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/2018\/09\/datanovia-logo.png\",\"contentUrl\":\"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/2018\/09\/datanovia-logo.png\",\"width\":98,\"height\":99,\"caption\":\"Datanovia\"},\"image\":{\"@id\":\"https:\/\/www.datanovia.com\/en\/#\/schema\/logo\/image\/\"}}]}<\/script>\n<!-- \/ Yoast SEO plugin. -->","yoast_head_json":{"title":"Mixed ANOVA in R: The Ultimate Guide - Datanovia","robots":{"index":"index","follow":"follow","max-snippet":"max-snippet:-1","max-image-preview":"max-image-preview:large","max-video-preview":"max-video-preview:-1"},"canonical":"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/","og_locale":"en_US","og_type":"article","og_title":"Mixed ANOVA in R: The Ultimate Guide - Datanovia","og_description":"The Mixed ANOVA is used to compare the means of groups cross-classified by two different types of factor variables, including: i) between-subjects factors, which have independent categories (e.g., gender: male\/female). ii) within-subjects factors, which have related categories also known as repeated measures (e.g., time: before\/after treatment). This chapter describes how to compute and interpret the different mixed ANOVA tests in R.","og_url":"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/","og_site_name":"Datanovia","article_modified_time":"2019-11-29T06:41:02+00:00","og_image":[{"width":1024,"height":512,"url":"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/2019\/05\/X39410210_704020863271194_6530535646589616128_n.jpg","type":"image\/jpeg"}],"twitter_card":"summary_large_image","twitter_misc":{"Est. reading time":"35 minutes"},"schema":{"@context":"https:\/\/schema.org","@graph":[{"@type":"WebPage","@id":"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/","url":"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/","name":"Mixed ANOVA in R: The Ultimate Guide - Datanovia","isPartOf":{"@id":"https:\/\/www.datanovia.com\/en\/#website"},"primaryImageOfPage":{"@id":"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/#primaryimage"},"image":{"@id":"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/#primaryimage"},"thumbnailUrl":"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/2019\/05\/X39410210_704020863271194_6530535646589616128_n.jpg","datePublished":"2019-11-29T06:37:37+00:00","dateModified":"2019-11-29T06:41:02+00:00","breadcrumb":{"@id":"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/#breadcrumb"},"inLanguage":"en-US","potentialAction":[{"@type":"ReadAction","target":["https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/"]}]},{"@type":"ImageObject","inLanguage":"en-US","@id":"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/#primaryimage","url":"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/2019\/05\/X39410210_704020863271194_6530535646589616128_n.jpg","contentUrl":"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/2019\/05\/X39410210_704020863271194_6530535646589616128_n.jpg","width":1024,"height":512},{"@type":"BreadcrumbList","@id":"https:\/\/www.datanovia.com\/en\/lessons\/mixed-anova-in-r\/#breadcrumb","itemListElement":[{"@type":"ListItem","position":1,"name":"Home","item":"https:\/\/www.datanovia.com\/en\/"},{"@type":"ListItem","position":2,"name":"Lessons","item":"https:\/\/www.datanovia.com\/en\/lessons\/"},{"@type":"ListItem","position":3,"name":"Mixed ANOVA in R"}]},{"@type":"WebSite","@id":"https:\/\/www.datanovia.com\/en\/#website","url":"https:\/\/www.datanovia.com\/en\/","name":"Datanovia","description":"Data Mining and Statistics for Decision Support","publisher":{"@id":"https:\/\/www.datanovia.com\/en\/#organization"},"potentialAction":[{"@type":"SearchAction","target":{"@type":"EntryPoint","urlTemplate":"https:\/\/www.datanovia.com\/en\/?s={search_term_string}"},"query-input":{"@type":"PropertyValueSpecification","valueRequired":true,"valueName":"search_term_string"}}],"inLanguage":"en-US"},{"@type":"Organization","@id":"https:\/\/www.datanovia.com\/en\/#organization","name":"Datanovia","url":"https:\/\/www.datanovia.com\/en\/","logo":{"@type":"ImageObject","inLanguage":"en-US","@id":"https:\/\/www.datanovia.com\/en\/#\/schema\/logo\/image\/","url":"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/2018\/09\/datanovia-logo.png","contentUrl":"https:\/\/www.datanovia.com\/en\/wp-content\/uploads\/2018\/09\/datanovia-logo.png","width":98,"height":99,"caption":"Datanovia"},"image":{"@id":"https:\/\/www.datanovia.com\/en\/#\/schema\/logo\/image\/"}}]}},"multi-rating":{"mr_rating_results":[]},"_links":{"self":[{"href":"https:\/\/www.datanovia.com\/en\/wp-json\/wp\/v2\/dt_lessons\/10875","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.datanovia.com\/en\/wp-json\/wp\/v2\/dt_lessons"}],"about":[{"href":"https:\/\/www.datanovia.com\/en\/wp-json\/wp\/v2\/types\/dt_lessons"}],"author":[{"embeddable":true,"href":"https:\/\/www.datanovia.com\/en\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.datanovia.com\/en\/wp-json\/wp\/v2\/comments?post=10875"}],"version-history":[{"count":0,"href":"https:\/\/www.datanovia.com\/en\/wp-json\/wp\/v2\/dt_lessons\/10875\/revisions"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/www.datanovia.com\/en\/wp-json\/wp\/v2\/media\/9015"}],"wp:attachment":[{"href":"https:\/\/www.datanovia.com\/en\/wp-json\/wp\/v2\/media?parent=10875"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}