Goals for today

Learn how to compare means of a continuous outcome from several independent groups to address the effect group on the dependent variable (modeling a continuous DV with a categorical IV). Extend this model to groups that differ on two factors (2x2 design).

Step 1 - Get organized

Decision process for comparing several independent groups (from Field textbook Chapter 12)
Decision process for comparing several independent groups (from Field textbook Chapter 12)

Step 2 - Import data and check for outliers, normality

Let’s get going!

  1. Make a new code chunk and use readr::read_csv() to read in the data. Make sure that NA values are handled the way you want (click on the tibble in the Environment window pane to take a quick look).
  2. Then make sure the columns that contain nominal vals are stored as class factor, using forcats::as_factor(), and use forcats::fct_relevel() to reorder the levels of emoreg_cond so that “control” is level 1, “suppress” is level 2, and “reappraise” is level 3 ( you can verify the level order of a factor by running levels(emo_tib$emoreg_cond))
    • emoreg_cond should be a factor (with 3 levels)
    • direct_cond should be a factor (with 2 levels)
  3. Then get descriptives for the DV card_out_react_prep and the covariate prepIOS_cent, organized by just the first IV emoreg_cond - use psych::describeBy() or dplyr::summarise()
    • there are no missing cases in this data (this is a subset of the original dataset) but you should be mindful of missing cases when working with your own data
    • note: if you choose to use psych::describeBy() you will get an error if you only ask for descriptives for a single column (so be sure to send both card_out_react_prep and prepIOS_cent to the function - see code below)
  4. We are going to start by considering only 1 of the categorical IVs, emoreg_cond. Group the data by just this first IV and check the distribution shape (for a normal-ish shape) and potential outliers for the DV (histogram and boxplot)
    • tip: to split a histogram into a separate plot for each group, use facet_wrap(~emoreg_cond) as an added layer to the ggplot code (see example solution)
# 1. first import the data
emo_tib <- readr::read_csv("data/oveis-2020-recoded-data.csv", na = "NA")
# 2. now make sure the columns we want as factors are treated that way, using forcats::as_factor()
#while you're at it, reorder the levels of 
emo_tib <- emo_tib %>% dplyr::mutate(
  emoreg_cond = forcats::as_factor(emoreg_cond) %>% 
    forcats::fct_relevel("control","suppress","reappraise"),
  direct_cond = forcats::as_factor(direct_cond) 
)

# 3. descriptives by emoreg_cond group
emo_tib %>% select(card_out_react_prep, prepIOS_cent, emoreg_cond) %>%
  psych::describeBy(card_out_react_prep + prepIOS_cent ~ emoreg_cond, data = ., mat = TRUE) %>% 
  knitr::kable(caption = "Cardiac Output and partner closeness by Emotion Regulation Condition", digits = 3) %>% 
  kableExtra::kable_classic(lightable_options = "hover")  

# 4. look at the data
p1 <- emo_tib %>% ggplot( aes(x = card_out_react_prep)) + 
  geom_histogram(position = "identity", bins = 15) +  
  facet_wrap( ~ emoreg_cond) +
  theme_classic()
p1
p2 <- emo_tib %>% ggplot( aes(x = emoreg_cond, y = card_out_react_prep)) + 
  geom_boxplot(position = "identity") +  
  theme_classic()
p2
Cardiac Output and partner closeness by Emotion Regulation Condition
item group1 vars n mean sd median trimmed mad min max range skew kurtosis se
card_out_react_prep1 1 control 1 81 -0.258 0.854 -0.294 -0.264 0.741 -2.857 1.892 4.749 -0.050 0.363 0.095
card_out_react_prep2 2 suppress 1 75 -0.131 0.907 -0.232 -0.144 1.001 -2.012 2.337 4.350 0.183 -0.523 0.105
card_out_react_prep3 3 reappraise 1 78 0.284 0.700 0.243 0.250 0.575 -1.195 2.563 3.757 0.562 0.701 0.079
prepIOS_cent1 4 control 2 81 -0.321 1.694 -0.457 -0.319 2.965 -3.457 2.543 6.000 -0.088 -1.187 0.188
prepIOS_cent2 5 suppress 2 75 0.156 1.423 -0.457 0.117 1.483 -2.457 2.543 5.000 0.190 -0.994 0.164
prepIOS_cent3 6 reappraise 2 78 0.184 1.512 0.543 0.230 1.483 -3.457 2.543 6.000 -0.211 -0.577 0.171

 

Step 3 - ANOVA using lm()

Step 3.1: specify the model

The distributions look normal-ish, and there is a decent sample size, so we can go ahead with the ANOVA model, remembering to examine the model residuals for sources of bias (e.g., extreme values, heteroskedasticity).
ANOVA and linear regression are both based on the General Linear Model and use the same underlying math, so we will use the same lm() function that we used in the regression activity a few weeks ago:
1. use the lm() function to build a linear model with 1 categorical predictor (lm(formula= card_out_react_prep ~ emoreg_cond)), store the model in a variable called card_lm
2. then use the anova() function on the resulting fitted model to partition the variance in the same way we did with SPSS to get the omnibus F-stat 3. and then we can use effect_size::eta_squared(model, partial = False, ci = .95) to get effect size η2 (eta-squared; since there is only one predictor this is the same as partial η2, and the same as the overall R2 for the model). You could also use omega_squared() in place of eta_squared(), which is a similar effect size measure that is adjusted to generalize to a population more accurately (it will always be lower than eta-squared).
Look at the solution below and make it work for you, then let’s think about the output.
Note: To be sure that you always get two-sided confidence intervals, set the argment alternative = "two.sided" in the call to the eta_squared function. The same issue applies to many functions in the effectsize package (e.g., cohens_d())

card_lm <- emo_tib %>% 
  lm(formula = card_out_react_prep ~ emoreg_cond)
anova(card_lm)
effectsize::eta_squared(card_lm, ci = .95, partial = FALSE, 
                        alternative = "two.sided") %>% 
  kableExtra::kbl(caption = "Effect size", digits = 4) %>% 
  kableExtra::kable_styling(full_width = FALSE)
## Analysis of Variance Table
## 
## Response: card_out_react_prep
##              Df  Sum Sq Mean Sq F value    Pr(>F)    
## emoreg_cond   2  12.674  6.3370  9.3238 0.0001276 ***
## Residuals   231 157.000  0.6797                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Effect size
Parameter Eta2 CI CI_low CI_high
emoreg_cond 0.0747 0.95 0.0196 0.1433
Understanding the output of anova():
  1. The Analysis of Variance Table gives the F-statistic, which is equal to the mean square between groups divided by the mean square within groups. In our regression examples we referred to this value more generally as the ratio of variance explained by the model to leftover error (residual variance).
  2. The p-value for the F-stat tells you the probability of an F-statistic at least that large (mean differences that large) under the null hypothesis that there is no difference between emo_reg group means. When you report the F-stat you include the model (between groups) and residual (within groups) degrees of freedom. E.g., F(2,231) = 9.324, p = .0001
  3. The eta-squared (η2) effect size can be interpreted as proportion of variance in the DV explained by the model. Since there is only one IV, partial eta-squared for the one predictor is the same as eta-squared. When we talk about an effect of a factor with many levels this is often called the “omnibus effect” to refer to the effect across all levels of the factor.
  4. The F-stat and effect size do not tell you anything about the direction of the group effect, so we need to think more about the model. Let’s look at contrasts of our categorical predictor (referred to as planned contrasts).

Step 3.2 - Planned contrasts

When we used lm() to run a linear regression (a few weeks ago), we examined the output with the summary() function rather than the anova() function, and it gave us some information about the model (e.g., coefficients for predictors) that we didn’t see above in the model description we got from anova(). In a code chunk or in your console, run summary(card_lm). You should get something like this:

summary(card_lm)
## 
## Call:
## lm(formula = card_out_react_prep ~ emoreg_cond, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.59904 -0.53027 -0.03875  0.47554  2.46870 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            -0.2578     0.0916  -2.814  0.00531 ** 
## emoreg_condsuppress     0.1264     0.1321   0.957  0.33975    
## emoreg_condreappraise   0.5422     0.1308   4.146 4.76e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8244 on 231 degrees of freedom
## Multiple R-squared:  0.0747, Adjusted R-squared:  0.06668 
## F-statistic: 9.324 on 2 and 231 DF,  p-value: 0.0001276
Notice the following in the output:
  1. The “Multiple R-squared” is the same as the eta-squared effect size you just calculated (proportion of variance explained by the model).
  2. The overall F-stat and p-value are the same as what you got from the anova() function.
  3. But now you see a Coefficients table with rows for the intercept and 2 predictors (emoreg_condsuppress and emoreg_condreappraise)
    • The intercept has the same meaning that it had in our linear regression examples. It is the predicted value of the DV when the values of all predictors are zero. But what does it mean for a categorical predictor to have a value of zero?
    • The answer is that it depends on how the categorical variable is “coded”. To estimate the lm() model, R recodes the 3-level categorical variable into 2 numerical columns (always the number of levels minus 1), and uses those numerical “contrasts” as predictors in a linear regression model ( DV = intercept + B1(contrast1) + B2(contrast2)). By default, R will use a “treatment” coding scheme (also called “dummy” coding) where the second level (“suppress”) gets a value of 1 in contrast1 (0 in contrast2) and the third level gets a value of 1 in contrast2 (0 in contrast1).
    • Run the command contrasts(emo_tib$emoreg_cond) to see how R has coded the categorical variable. Notice that first level gets a value of 0 in both contrasts- this is called the “reference level”. This means that the intercept is the mean of the reference level (you can go back to your descriptives above and verify that the mean cardiac output for the “control” group is -.26, same as the intercept in our model output) - this is true because the intercept is the predicted value when the two predictors (contrasts) are both zero.
    • The two contrasts in the model summary() output get named by the factor name and the contrast name put together. The coefficient for each contrast represents (just like in linear regression) the predicted change in DV for a one unit change in the IV. So for contrast1 (emoreg_condsuppress) a one unit change is the difference between the reference level (“control”) and the second level (“suppress”). This means that adding the coefficient (.126) to the intercept (-.258) should give you the mean of the “suppress” group (-.132) - the math checks out. Likewise, the coefficient for contrast2 (emoreg_condreappraise) represents the difference between the reference group (“control”) and the third level (“reappraise”).
    • The t-value and Pr(>|t|) (p-value) for each contrast have the same meaning as in linear regression, and allow you to make inferences based on the contrasts. The p-value represents the probability of getting a t-stat at least that far from 0 under the null hypothesis that the true contrast coefficient is 0.
    • So, for contrast1 (emoreg_condsuppress), the p-value of .3398 means it is not surprising to observe a contrast coefficient that far from 0 under the null hypothesis. The contrast represents the difference between “control” and “suppress” groups, so we would fail to reject the null hypothesis that those two groups have the same mean cardiac output in this paradigm.
    • For contrast2 (emoreg_condreappraise), the p-value of .000048 (4.76e-05) means it is unlikely to observe a contrast coefficient that far from 0 under the null, so we have reason to reject the null hypothesis that the “control” and “reappraise” groups do not differ.
    • Contrasts from this “treatment”/“dummy” contrast coding scheme are not orthogonal. A consequence of non-orthogonal contrasts is that they are not independent (e.g., [suppress vs control] and [reappraise vs control] effects may both be driven by a low control group mean.

Step 3.3 - other contrast coding schemes

Notice that the contrasts above don’t give us a comparison between “suppress” and “reappraise” groups. In planned contrasts, we are limited to k-1 contrasts, where k is the number of levels of the categorical predictor. The default (“treatment”/“dummy”) coding scheme allows us to test a reference level vs. each other level, but other coding schemes allow us to make other comparisons if they make more sense for our purposes. Check out this page for an overview of contrast coding (from the UCLA Institute for Digital Research and Education) to learn how to use the contrasts() function to set contrasts for a categorical variable.

Polynomial contrast coding

The default treatment/dummy coding scheme makes sense if we hypothesized that both “suppress” and “reappraise” instructions would influence cardiac output. But if the groups were ordered (e.g., “control”, “suppress”, and “reappraise” represent increasing levels of regulation), then we might instead want to test for a linear trend (increasing cardiac output between level 1 and level 3) and quadratic trend (level 2 is lower or higher than the extremes, like a U or inverted U shape) by including a contrast for each. If we had 4 levels we could include a cubic contrast. This contrast coding scheme is referred to as polynomial contrast coding. Use the code below to see how to set the contrasts to test for a linear and quadratic trend (using the contr.poly() function).

# set contrasts using contr.poly()
contrasts(emo_tib$emoreg_cond) = contr.poly(3) #3 corresponds to number of levels
# view the contrast coding scheme
contrasts(emo_tib$emoreg_cond)
# re-estimate the model
card_lm <- emo_tib %>% 
  lm(formula = card_out_react_prep ~ emoreg_cond)
# output the model summary
summary(card_lm)
##                       .L         .Q
## control    -7.071068e-01  0.4082483
## suppress   -9.073800e-17 -0.8164966
## reappraise  7.071068e-01  0.4082483
## 
## Call:
## lm(formula = card_out_react_prep ~ emoreg_cond, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.59904 -0.53027 -0.03875  0.47554  2.46870 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.03493    0.05392  -0.648    0.518    
## emoreg_cond.L  0.38338    0.09248   4.146 4.76e-05 ***
## emoreg_cond.Q  0.11816    0.09430   1.253    0.211    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8244 on 231 degrees of freedom
## Multiple R-squared:  0.0747, Adjusted R-squared:  0.06668 
## F-statistic: 9.324 on 2 and 231 DF,  p-value: 0.0001276
Look at the output
  1. The overall model stats have not changed (R-squared, F, p-value are all the same as before)
  2. emoreg_cond.L is the linear contrast. The t-stat and low p-value mean that a linear contrast coefficient this far from 0 is unlikely under the null hypothesis that there is no linear (increasing or decreasing) trend in cardiac output across levels of emotion regulation.
  3. emoreg_cond.Q is the quadratic contrast. If we had hypothesized that the influence of emotion regulation on cardiac output was U-shaped then we would have expected a significant quadratic contrast coefficient, which is not the case here.

Step 3.4- check assumptions

We are using the same modeling framework that we used in linear regression, so the assumptions concerning residuals (normality and homogeneity) are the same, and we can check them the same way. Use the plot() function to get the default plots of residuals, like we did in the linear regression activity.
We have a categorical predictor with 3 levels, so you should notice that there are only 3 possible values for model fitted values, which correspond to the 3 group means.
If you are concerned about assumption violations, particularly homogeneity of variance, you can run the Kruskal-Wallis non-parametric test using kruskal.test(formula = card_out_react_prep ~ emoreg_cond, data = emo_tib). The Kruskal-Wallis test, like other nonparametric tests we have used before, is based on ranking the data.

plot(card_lm)
kruskal.test(formula = card_out_react_prep ~ emoreg_cond, data = emo_tib)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  card_out_react_prep by emoreg_cond
## Kruskal-Wallis chi-squared = 18.521, df = 2, p-value = 9.51e-05

 

Step 3.5 - Now let’s describe the key results so far

How would we write up what we found? Let’s focus on the original model with the default treatment coding scheme. We found that there was a significant effect of emotion regulation instruction on cardiac output, and a significant difference between the “reappraise” and “control” conditions, but no significant difference between “control” and “suppress” conditions. We could write that up like this:

Emotion regulation condition impacted cardiac output during partner interaction (F(2,231) = 9.324, p = .0001, η2 = .075). Planned contrasts showed that participants in the reappraisal condition exhibited greater cardiac output (M = 0.284, SD = 0.760) compared to the control condition (M = -0.258, SD = 0.854, t(231) = 4.146, p < .0001), but cardiac output in the suppression condition (M = -.131, SD = .907) did not significantly differ from the control condition (t(231) = .957, p = .340).

Of course, if we had a specific hypothesis that the reappraisal condition differed from the suppression condition we could be more complete and test the difference between all pairs of groups. We will get into that when we do post hoc tests in factorial ANOVA. The idea behind planned contrasts is to conduct limited tests for specific a priori hypotheses.

Step 4 - ANCOVA - independent groups with a continuous covariate

Now that we have seen how to model a continous DV with a categorical IV, we can extend the model. The variable prepIOS_cent is closeness ratings from each participant, rating how connected they felt to the partner they interacted with. By including prepIOS_cent as a covariate, we can see what effect a the emotion regulation variable has, adjusting for the effect of the covariate (described in full in Field textbook section 13.2).
Based on the code you wrote in Step 3.1, include prepIOS_cent as a predictor in the model. Simply add the term in the lm() model formula. But there’s one more thing. When there was just one predictor in the model we used the anova() function to get the omnibus F-stat for the overall emotion regulation effect. But the anova() function uses sequential (type 1) sums of squares, where variance is assigned to predictors in the order they appear in the model. So the F-stat will depend on the order of the variables in the model. We don’t want that and it is more common to use type 3 sums of squares (equivalent to type 2 if there is no interaction term - see a good explanation here). We can get type 3 SS using the car package Anova() function, but first we need to make sure that our contrast coding for the categorical variable is orthogonal (otherwise we may get uninterpretable results). A simple way to set orthogonal contrasts for a categorical variable is with the contr.helmert() function we can discuss this at the end of the activity. For now, use the code below to set the contrasts, fit the model, and get the ANOVA table using car::Anova(type=3)

#set an orthogonal coding scheme for emoreg_cond
contrasts(emo_tib$emoreg_cond) <- contr.helmert(3)
card_cov_lm <- emo_tib %>% lm(formula = card_out_react_prep ~  emoreg_cond + prepIOS_cent)
car::Anova(card_cov_lm, type=3)
## Anova Table (Type III tests)
## 
## Response: card_out_react_prep
##               Sum Sq  Df F value   Pr(>F)    
## (Intercept)    0.289   1  0.4250 0.515108    
## emoreg_cond   11.896   2  8.7527 0.000217 ***
## prepIOS_cent   0.703   1  1.0341 0.310256    
## Residuals    156.297 230                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 

Understanding the output

The anova table now gives us an F-statistic for each predictor. The F-stat for the emotion regulation condition now tells us about the difference in group means, adjusting for closeness. The F-stat for closeness (prepIOS_cent) tells us about the covariance with the DV. You can then examine planned contrasts (with the summary() function), which are adjusted for the covariate, or compare adjusted means using post-hoc tests that we will discuss in the final step.

Step 5 - Factorial Anova

In the data description there is actually a second factor we haven’t considered yet: direct_cond with two levels. This variable indicates whether a participant received instructions to regulate their emotions directly (“self”), or whether they were exposed indirectly (“other”) by interacting with a partner who had received emotion regulation instructions. The design is crossed, so there are a total of 6 (3*2) groups. Let’s examine the full design now.
We can hypothesize that there is an interaction of emoreg_cond and direct_cond such that the effect of emoreg_cond on cardiac output is greater when the instructions are received directly (“self” condition) compared to indirectly (“other” condition).

Step 5.1 - plot the means in a clustered bar chart

Use ggplot and map emoreg_cond to x, and direct_cond to fill. Use geom_col() to make bars for the group means, then add error bars (representing 95% CI) with geom_errorbar(). We did a similar chart in Step 3.1 of the data visualization activity, but it’s been a while, so take a peek at the solution if it helps.

card_out_summary <- emo_tib %>% 
  group_by(emoreg_cond,direct_cond) %>% 
  dplyr::summarise(
    mean_card_out = mean(card_out_react_prep),
    ci.low = mean_cl_normal(card_out_react_prep)$ymin,
    ci.upp = mean_cl_normal(card_out_react_prep)$ymax,
  ) %>% ungroup()

# add ymin and ymax to the aesthetics list, which are needed for geom_errorbar
p5 <- card_out_summary %>% 
  ggplot(aes(x = emoreg_cond, y = mean_card_out, fill = direct_cond,
             ymin=ci.low, ymax=ci.upp)) + 
    geom_col(position = position_dodge(.9)) + 
    geom_errorbar(width=.2,position = position_dodge(.9)) +
    theme_classic() + labs(title="Cardiac output means by condition")
p5

Step 5.2 - Factorial ANOVA with afex::aov_4()

Now that you have a picture of the pattern of means, let’s run the full model (DV predicted by IV1 + IV2 + IV1*IV2). We will use the afex::aov_4() function, which takes arguments formatted similarly to the lm() and other modeling functions. One important difference is that in the formula we need to specify an identifier variable in the formula.
Try it now (look at the solution).

card_afx <- afex::aov_4(card_out_react_prep ~ emoreg_cond*direct_cond + (1|Sbj_ID), data = emo_tib)
## Contrasts set to contr.sum for the following variables: emoreg_cond, direct_cond
card_afx
## Anova Table (Type 3 tests)
## 
## Response: card_out_react_prep
##                    Effect     df  MSE        F  ges p.value
## 1             emoreg_cond 2, 228 0.68 9.07 *** .074   <.001
## 2             direct_cond 1, 228 0.68     0.77 .003    .381
## 3 emoreg_cond:direct_cond 2, 228 0.68     0.71 .006    .493
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

 

Check out the output
  1. Notice the message “Contrasts set to contr.sum …” – this is telling us that the function set the contrasts for the categorical variables. This is done to ensure that the contrasts are orthogonal, it’s nice that the function takes care of it for us but it means that whatever contrast coding we might have set before is overridden when the function runs.
  2. We get an F-stat for each of the three predictors (emoreg_cond, direct_cond and their interaction emoreg_cond:direct_cond). If we had found a significant interaction, we would infer that the effect of emotion regulation depended on whether the instructions were received directly or just by interacting with a partner who received emotion regulation instructions.
  3. We do find a significant main effect of emotion regulation condition (consistent with what we saw earlier). The column “ges” gives a “generalized eta-squared” effect size, which is similar to the eta-squared effect size that we saw for the effect of emotion regulation in the model with just one predictor.

Step 5.3 - Simple effects

  • Although there is no interaction, we might want to confirm that there is a significant emotion regulation effect within each level of the direct_cond factor. It is interesting if there is an effect just by interacting with a person who is regulating emotions, so we should confirm that.
  • We can run simple effects tests, which just means testing the effect of one factor at each level of the second factor
  • let’s test the effect of emoreg_cond at each level of direct_cond now. We use the emmeans::joint_tests() function and specify our model and the name of the variable (in quotes) that we will test across:
emmeans::joint_tests(card_afx, "direct_cond")
## direct_cond = self:
##  model term  df1 df2 F.ratio p.value
##  emoreg_cond   2 228   4.582  0.0112
## 
## direct_cond = partner:
##  model term  df1 df2 F.ratio p.value
##  emoreg_cond   2 228   5.360  0.0053

 

This function basically runs a one way anova at each level of direct_cond. We can see that there is a significant effect of emotion regulation in both the “self” and “other” conditions. We might write up the findings like this:
A 3 (emotion regulation condition) X 2 (emotion regulation instruction received by oneself or by partner) factorial ANOVA revealed a main effect of emotion regulation on cardiac output during the task (F(2,228) = 9.07, p < .001). The main effect of receiving instructions to oneself versus one’s partner was not significant (F(2,228) = .77, p = .381) and the interaction was not significant (F(2,228) = .71, p = .493). The effect of emotion regulation was significant in both groups, whether emotion regulation instructions were received by oneself (F(2,228) = 4.58, p = .011) or by one’s partner (F(2,228) = 5.36, p = .005).

We would then follow up with more specific information about the means as we did above. If you are interested, try running post-hoc tests with emmeans() like in the code below. If you want to test whether the “reappraise” condition differs from the “suppress” condition then this is one way you could address that.

Extra stuff

Check out the code below for easy plotting with afex and post-hoc tests to compare all pairs of groups while correcting for multiple comparisons.

# Easy plot of the means with points overlayed (must estimate model with afex first):
afex::afex_plot(card_afx, "emoreg_cond", "direct_cond")

# marginal means and post hoc comparisons between all pairs of groups
emmeans::emmeans(card_afx, pairwise ~ emoreg_cond*direct_cond, adjust = "tukey")

# or, compare levels of emoreg_cond collapsing across groups 
emmeans::emmeans(card_afx, pairwise ~ emoreg_cond, adjust = "tukey")
## NOTE: Results may be misleading due to involvement in interactions
## $emmeans
##  emoreg_cond direct_cond  emmean    SE  df lower.CL upper.CL
##  control     self        -0.1295 0.134 228  -0.3935   0.1344
##  suppress    self        -0.1663 0.132 228  -0.4268   0.0942
##  reappraise  self         0.3421 0.132 228   0.0816   0.6027
##  control     partner     -0.3711 0.126 228  -0.6193  -0.1230
##  suppress    partner     -0.0936 0.138 228  -0.3648   0.1776
##  reappraise  partner      0.2267 0.132 228  -0.0339   0.4872
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                              estimate    SE  df t.ratio p.value
##  control self - suppress self            0.0368 0.188 228   0.195  1.0000
##  control self - reappraise self         -0.4717 0.188 228  -2.506  0.1267
##  control self - control partner          0.2416 0.184 228   1.314  0.7769
##  control self - suppress partner        -0.0359 0.192 228  -0.187  1.0000
##  control self - reappraise partner      -0.3562 0.188 228  -1.892  0.4096
##  suppress self - reappraise self        -0.5085 0.187 228  -2.719  0.0754
##  suppress self - control partner         0.2048 0.183 228   1.122  0.8720
##  suppress self - suppress partner       -0.0727 0.191 228  -0.381  0.9989
##  suppress self - reappraise partner     -0.3930 0.187 228  -2.101  0.2902
##  reappraise self - control partner       0.7133 0.183 228   3.906  0.0017
##  reappraise self - suppress partner      0.4357 0.191 228   2.283  0.2052
##  reappraise self - reappraise partner    0.1155 0.187 228   0.618  0.9897
##  control partner - suppress partner     -0.2775 0.187 228  -1.488  0.6724
##  control partner - reappraise partner   -0.5978 0.183 228  -3.274  0.0153
##  suppress partner - reappraise partner  -0.3202 0.191 228  -1.678  0.5478
## 
## P value adjustment: tukey method for comparing a family of 6 estimates 
## 
## $emmeans
##  emoreg_cond emmean     SE  df lower.CL upper.CL
##  control     -0.250 0.0919 228   -0.431  -0.0692
##  suppress    -0.130 0.0954 228   -0.318   0.0581
##  reappraise   0.284 0.0935 228    0.100   0.4686
## 
## Results are averaged over the levels of: direct_cond 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast              estimate    SE  df t.ratio p.value
##  control - suppress      -0.120 0.132 228  -0.909  0.6355
##  control - reappraise    -0.535 0.131 228  -4.078  0.0002
##  suppress - reappraise   -0.414 0.134 228  -3.102  0.0061
## 
## Results are averaged over the levels of: direct_cond 
## P value adjustment: tukey method for comparing a family of 3 estimates

That’s all for this activity!

References