You will need the Rmd file for today's lab - can download it [here](https://uopsych.github.io/psy612/labs/lab-8/lab-8.Rmd).


Be sure to have the following packages installed and loaded:

```{r lab-8-1, message=FALSE}
library(tidyverse) # for plotting and data wrangling
library(rio) # for importing data
library(psych) # for descriptives
library(lsr) # for eta squared functions
library(emmeans) # for marginal means and simple effects
library(sjPlot) # for plotting model results 
library(apaTables) # for tables of means
library(car) # for testing model assumptions
library(broom) # for tidying model output

# suppress scientific notation
options(scipen = 999)
```


# Purpose

Factorial ANOVA (aka multiple regression with categorical predictors) refers to a special case of the general linear model in which there is an interaction of two or more categorical variables (i.e. *factors*). A factorial design is used when there is an interest in how two or more variables (or factors) affect some outcomes variable. Rather than conduct separate one-way ANOVAs for each factor, they are all included in one analysis. Today we will review how to run factorial ANOVA models in R and how to interpret and visualize the results. 


# Research scenario

Based on subjects' self-reports of rejection sensitivity (N = 80), a researcher divides subjects into two equal groups (low RS and high RS).  Whereas half of the subjects in each group interact with a partner who displays a happy emotional expression during the interaction, the other half of the subjects in each group interact with a partner who displays a neutral emotional expression during the interaction.  After the interaction, subjects are asked to rate the statement, "My interaction partner likes me", on a scale from 1 (strongly disagree) to 7 (strongly agree).

**Factor 1: Rejection Sensitivity**

* Low
* High

**Factor 2: Partner's Emotional Expression**

* Neutral
* Happy

**Dependent Variable: Perceived Liking**



# Hypothesis Testing

With a factorial ANOVA, three different null hypotheses can be tested regarding whether there is 1) a main effect of factor 1, 2) a main effect of factor 2, and 3) an interaction between factor 1 and factor 2.

Imagine a table with means like the following:

|         |High | Low | Row Means |
| ------- |-----| ----- | ----- |
| Neutral |x.xx   | x.xx  | x.xx  |
| Happy   |x.xx   | x.xx  | x.xx  |
| Column Means  |x.xx | x.xx  | x.xx |

**Null Hypothesis 1** (Main Effect of Rejection Sensitivity):
$$H_0: \mu_{Low} = \mu_{High}$$

**Null Hypothesis 2** (Main Effect of Partner's Emotional Expression):
$$H_0: \mu_{Neutral} = \mu_{Happy}$$
**Null Hypothesis 3** (Interaction Effect):

H0: There is no interaction between rejection sensitivity and partner's emotional expression.


# Import & inspect the Data

```{r lab-8-2}
reject <- import("https://raw.githubusercontent.com/uopsych/psy612/master/labs/lab-8/data/reject.csv")
```

```{r}
str(reject)
head(reject)
```

* It looks like `rs` and `partner` are both being read in as character variables. Let's go ahead and change those to factors.

```{r lab-8-7}
reject <- reject %>% 
  mutate(rs = as.factor(rs),
         partner = as.factor(partner))
```

* Check the structure again. 

```{r}
str(reject)
```


Notice that by default R orders factor levels alphabetically. In our case, this means that `High` will be the reference group of rejection sensitivity and `Happy` will be the reference group of interaction partner's emotional expression. However, it might be more intuitive to have `Low` and `Neutral` be the reference groups, respectively. 


* To accomplish this, we can simply re-order the levels of our factor variables with `fct_relevel()`.
```{r}
reject <- reject %>% 
  mutate(rs = fct_relevel(rs, c("Low", "High")), 
         partner = fct_relevel(partner, c("Neutral", "Happy")))

str(reject) # Notice the re-ordering of levels
```


* We can also see how the factors are being coded using contrasts().
```{r}
contrasts(reject$rs)
contrasts(reject$partner)
```

Notice that each are being dummy coded. It's very important to be aware of how the factors in your regression model are being coded to be able to accurately interpret the regression coefficients in your output.

Speaking of regression, note we are going to conduct the analysis using the lm() function. The following is how you could write the factorial ANOVA as the regression equation that we're estimating:

$$Liking_i = \beta_0 + \beta_1RejectionSensitivity_i + \beta_2PartnerEmotion_i + \beta_3RejectionXPartner_i + \epsilon_i$$


# Descriptive Statistics

Obtaining sample size per condition.

```{r}
reject %>%
  group_by(rs, partner) %>% 
  summarize(n = n())
```


## Table of Means

The results of a factorial ANOVA are often represented in a table of means, which contains the means of each combination of conditions (the cell means) and means for the levels of each variable ignoring the other variables (the marginal means).

### Table of Means using `apa.2way.table()`

* The apa.2way.table() function from the apaTables package is a very convenient way to get our cell means and marginal means. This function works for any type of 2-way ANOVA, regardless of the number of levels of your factors, e.g. it would work for a 3 x 3 ANOVA. All you need to do is indicate what the IV's (aka factors) and DV are and specify `show.marginal.means = TRUE`. 

```{r}
table_of_means <- apa.2way.table(iv1 = rs, 
               iv2 = partner, 
               dv = liking, 
               data = reject,
               show.marginal.means = TRUE)
table_of_means
```


We can also get the output of the apa.2way.table() function exported to a word document.
```{r}
apa.2way.table(iv1 = rs,
               iv2 = partner,
               dv = liking,
               data = reject,
               show.marginal.means = TRUE,
               filename = "table_of_means.doc")
```

### Table of Means using `dplyr`

The following can be used to obtain cell means.
```{r}
reject %>%
  group_by(rs, partner) %>% # Group by both rs and partner
  summarize(mean = mean(liking, na.rm = TRUE), # Obtain the mean for each condition
            sd = sd(liking, na.rm = TRUE)) %>%  # Obtain the SD for each condition
knitr::kable(digits = c(NA, NA, 2, 2, 2), # Specifying number of digits to report
             caption = "Cell Means & SD")
```


The following can be used to obtain marginal means for rejection sensitivity.
```{r}
# Marginal Means for Rejection Sensitivity
reject %>%
  group_by(rs) %>% # Just group by rejection sensitivity
  summarize(mean = mean(liking, na.rm = TRUE),
            sd = sd(liking, na.rm = TRUE)) %>% 
knitr::kable(digits = c(NA, 2, 2, 2),
             caption = "Marginal Means & SD for Rejection Sensitivity")
```


And the marginal means for partner's emotional expression.
```{r}
# Marginal Means for Partner's Emotional Expression
reject %>%
  group_by(partner) %>% # Just group by partner
  summarize(mean = mean(liking, na.rm = TRUE),
            sd = sd(liking, na.rm = TRUE)) %>% 
knitr::kable(digits = c(NA, 2, 2, 2),
             caption = "Marginal Means & SD for Partner's Emotional Expression")
```


The following can be used to obtain the grand mean (the overall mean on the DV).
```{r}
reject %>% # note that we don't need group_by for this one
  summarize(mean = mean(liking, na.rm = TRUE),
            sd = sd(liking, na.rm = TRUE)) %>% 
knitr::kable(digits = c(2, 2, 2),
             caption = "Grand Mean")
```


>**Question:** Which means are being compared in the main effect of rejection sensitivity?

>**Question:** Which means are being compared in the main effect of interaction partner?

>**Question:** Which means are involved in the interaction?


# Conduct the Factorial ANOVA 

## Using Regression: `lm`()

Factorial ANOVA is the method by which we can examine whether two (or more) categorical IVs have joint effects on a continuous outcome of interest. Like all general linear models, factorial ANOVA is a specific case of multiple regression. However, we may choose to use an ANOVA framework for the sake of interpretability.

Write the `lm()` code that performs the factorial ANOVA analysis.
```{r}
## your code here
```

## Results in regression framework

```{r}
## your code here
```

>**Question:** Write the full model with parameter estimates filled in below:

$$Liking_i = $$
>**Question:** What do each of the regression coefficient estimates mean? 

It may be helpful to refer back to the table of means.
```{r}
table_of_means
```




### Recoding the Factors

Obviously, the means being compared by b1 and b2 do not represent our main effects. The main effect of RS would be a comparison of the marginal means for the low and high conditions. The main effect of partner would be a comparison of the marginal means for the neutral and happy conditions.

With the way we have rs and partner coded in the model, the regression coefficient estimates don't correspond simply to the mean difference between conditions. We can recode the factors so that the parameter estimates, b1 and b2, correspond to a test of the main effects.

>**Question:** How should we code the levels of rejection sensitivity and partner emotion so that the parameter estimates, b1 and b2, correspond to a test of their main effects?

Assign new codes to the levels of each factor.
```{r}
## your code here
```

Re-run the analysis with newly coded factors.
```{r}
## your code here
```


## Using Traditional ANOVA: `anova()`

You can also perform a traditional ANOVA and get straightforward results by passing the model to the `anova()` function.

```{r}
model <- lm(liking ~ rs*partner, data = reject)
anova(model)
```

>**Question:** How would you interpret these results? 



# Effect Sizes 

Effect size is really important to report. It provides an idea of the *size* of the effect of a predictor variable on an outcome variable. Reporting effect size in publications is also very helpful for other researchers wanting to conduct a priori power analyses.

We can get effect sizes by passing our model to the `etaSquared()` function.
```{r}
etaSquared(model)
```

>**Question:** What is the difference between eta-squared and partial eta-squared? 




# Visualization

Visualizing the results is helpful for understanding your results and interpreting significant (or non-significant) effects. Especially when there is a significant interaction effect, a plot can help illustrate what the interaction effect is.  

## Main effects

* Remember that main effects correspond to differences in marginal means. To plot main effects, we can use the plot_model() function from the sjPlot package. This function takes three arguments: 1) the fitted model object, 2) because we want to plot marginal means, we specify `type = emm`, and 3) Specify in the `terms = ` argument which variable(s) you want the marginal means for.

### Plotting the Main Effect of Rejection Sensitivity
```{r}
plot_model(model, type = "emm", terms = "rs") # bars correspond to 95% CIs
```

### Plotting the Main Effect of Partner's Emotional Expression
```{r}
plot_model(model, type = "emm", terms = "partner") 
```

### Plotting the Interaction Effect
```{r}
plot_model(model, type = "emm", terms = c("rs", "partner")) 
```

Switch how the interaction is visualized by switching the order of `terms`.

>**Question:** How would you describe in a paper what the significant interaction effect means? 


# Simple effects

Simple effects are the effect of some factor (e.g., interaction partner's expression) at each level of another factor (e.g., at high and low RS separately).

We can get these simple effects using the `emmeans` function and calculating the differences between every pair of cell means.

```{r}
# Comparing cell means
emmeans(model, pairwise~rs*partner, adjust = "bonferroni")

# for info on types of p-value adjustments you can choose from:
?summary.emmGrid # scroll to p-value adjustments
```

>**Question:** Does partner emotion have an effect on perceived liking for people low on rejection sensitivity? For people high on rejection sensitivty? 


>**Question:** Does rejection sensitivity have an effect on perceived liking when partner shows neutral emotions? When partner shows happy emotions? 



You could also use the `emmeans` function to compare marginal means for each factor individually.
```{r}
# Comparing marginal means for partner emotion
emmeans(model, pairwise~partner, adjust = "bonferroni")

# Comparing marginal means for partner emotion
emmeans(model, pairwise~rs, adjust = "bonferroni")
```

