EDUC 640

class: center, middle, inverse, title-slide

# EDUC 640
## Two-Way ANOVA
### Chris Ives

---

# Contents

* Foundations
  + Crosstabs
  + Clustered Boxplots

* ANOVA
  + Two-way ANOVA
  + Marginal Means
  + Simple Effects
  + Interaction Plots

---

## PSA - Scientific Notation

Sometimes I find scientific notation to be a bit annoying to work with in R. If you want to turn it off for your session run the following code:

```r
options(scipen = 999)
```

For the sake of table formatting, I'll keep mine on for the slides.
---
## Crosstabs

The following code shows how to make crosstabs using `dplyr` functions. `group_by` and `summarise` are good to know, so I'd suggest running these lines one by one so you can get a sense of what they're doing. `summarise` creates a row for each group we've specified, then we specify the column contents.

```r
data %>% 
  group_by(inst, meth) %>% 
  summarise(n = n()) %>% #new column "n" = row count of each factor grouping
  spread(meth, n) 
```

```
## `summarise()` has grouped output by 'inst'. You can override using the `.groups` argument.
```

```
## # A tibble: 3 x 3
## # Groups:   inst [3]
##   inst             computer standard
##   <fct>               <int>    <int>
## 1 physical science        6        6
## 2 social science          6        6
## 3 history                 6        6
```

---
## Check Descriptives

```r
rmarkdown::paged_table(
  describe(data)
  )
```

## Grouped Descriptives

Setting `mat = TRUE` gives you the output for each level in one data frame.

```r
rmarkdown::paged_table(
  describeBy(x = data$vocab, group = data$inst, 
             mat = TRUE, data = data)
  )
```

---
## Separate Boxplots

Since we already covered normal boxplots, I'll just demonstrate a way to plot multiple boxplots if the need every arises. You'll need the `gridExtra` package.

```r
inst_plot <- ggplot(data, aes(x = inst, y = vocab)) +
  geom_boxplot()
meth_plot <- ggplot(data, aes(x = meth, y = vocab)) +
  geom_boxplot()
gridExtra::grid.arrange(inst_plot, meth_plot, nrow = 1) ## nrow = 1 so they are arranged side by side
```

![](wk4_files/figure-html/unnamed-chunk-5-1.png)

---
## Clustered Boxplots
Same process but you'll need to specify the other IV using `color =`, or `fill =`.

```r
inst_color <- ggplot(data, aes(x = inst, y = vocab, color = meth)) +
  geom_boxplot()
meth_fill <- ggplot(data, aes(x = meth, y = vocab, fill = inst)) +
  geom_boxplot()
gridExtra::grid.arrange(inst_color, meth_fill, nrow = 1)
```

![](wk4_files/figure-html/unnamed-chunk-6-1.png)

---
class: inverse-blue middle
# Two-Way ANOVA

```r
needs(rstatix, emmeans)
```
---

## Two-Way ANOVA
For a two-way ANOVA we just add our additional IV to the right-side of the formula. We specify the interaction using `inst*meth`.

```r
m1 <- anova_test(data, formula = vocab ~ inst + meth + inst:meth, 
                 detailed = TRUE, type = 3, effect.size = "pes")
```

```
## Coefficient covariances computed by hccm()
```

```r
m1
```

```
## ANOVA Table (type III tests)
## 
##        Effect   SSn  SSd DFn DFd       F        p p<.05   pes
## 1 (Intercept) 40401 1668   1  30 726.637 1.39e-22     * 0.960
## 2        inst  1194 1668   2  30  10.737 3.04e-04     * 0.417
## 3        meth  2209 1668   1  30  39.730 5.99e-07     * 0.570
## 4   inst:meth   722 1668   2  30   6.493 5.00e-03     * 0.302
```

---
## Levene's Test

```r
car::leveneTest(vocab ~ inst*meth, data = data, center = "mean")
```

```
## Levene's Test for Homogeneity of Variance (center = "mean")
##       Df F value  Pr(>F)  
## group  5  2.0579 0.09877 .
##       30                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---
## Marginal Means

To get estimated marginal means for an interaction, we need to group by one of the IVs and then run our `emmeans_test` with the other. This just gives us the mean differences and significance tests.

```r
means <- data %>%
  group_by(meth) %>% 
 emmeans_test(vocab ~ inst, p.adjust.method = "holm", detailed = TRUE)
paged_table(means)
```

<div data-pagedtable="false">
  <script data-pagedtable-source type="application/json">
{"columns":[{"label":[""],"name":["_rn_"],"type":[""],"align":["left"]},{"label":["meth"],"name":[1],"type":["chr"],"align":["left"]},{"label":["term"],"name":[2],"type":["chr"],"align":["left"]},{"label":[".y."],"name":[3],"type":["chr"],"align":["left"]},{"label":["group1"],"name":[4],"type":["chr"],"align":["left"]},{"label":["group2"],"name":[5],"type":["chr"],"align":["left"]},{"label":["null.value"],"name":[6],"type":["dbl"],"align":["right"]},{"label":["estimate"],"name":[7],"type":["dbl"],"align":["right"]},{"label":["se"],"name":[8],"type":["dbl"],"align":["right"]},{"label":["df"],"name":[9],"type":["dbl"],"align":["right"]},{"label":["conf.low"],"name":[10],"type":["dbl"],"align":["right"]},{"label":["conf.high"],"name":[11],"type":["dbl"],"align":["right"]},{"label":["statistic"],"name":[12],"type":["dbl"],"align":["right"]},{"label":["p"],"name":[13],"type":["dbl"],"align":["right"]},{"label":["p.adj"],"name":[14],"type":["dbl"],"align":["right"]},{"label":["p.adj.signif"],"name":[15],"type":["chr"],"align":["left"]}],"data":[{"1":"computer","2":"inst","3":"vocab","4":"physical science","5":"social science","6":"0","7":"6","8":"4.305036","9":"30","10":"-2.7920561","11":"14.79206","12":"1.3937166","13":"1.736427e-01","14":"3.472853e-01","15":"ns","_rn_":"1"},{"1":"computer","2":"inst","3":"vocab","4":"physical science","5":"history","6":"0","7":"8","8":"4.305036","9":"30","10":"-0.7920561","11":"16.79206","12":"1.8582888","13":"7.296550e-02","14":"2.188965e-01","15":"ns","_rn_":"2"},{"1":"computer","2":"inst","3":"vocab","4":"social science","5":"history","6":"0","7":"2","8":"4.305036","9":"30","10":"-6.7920561","11":"10.79206","12":"0.4645722","13":"6.455920e-01","14":"6.455920e-01","15":"ns","_rn_":"3"},{"1":"standard","2":"inst","3":"vocab","4":"physical science","5":"social science","6":"0","7":"22","8":"4.305036","9":"30","10":"13.2079439","11":"30.79206","12":"5.1102943","13":"1.706157e-05","14":"5.118471e-05","15":"****","_rn_":"4"},{"1":"standard","2":"inst","3":"vocab","4":"physical science","5":"history","6":"0","7":"3","8":"4.305036","9":"30","10":"-5.7920561","11":"11.79206","12":"0.6968583","13":"4.912565e-01","14":"4.912565e-01","15":"ns","_rn_":"5"},{"1":"standard","2":"inst","3":"vocab","4":"social science","5":"history","6":"0","7":"-19","8":"4.305036","9":"30","10":"-27.7920561","11":"-10.20794","12":"-4.4134360","13":"1.212903e-04","14":"2.425807e-04","15":"***","_rn_":"6"}],"options":{"columns":{"min":{},"max":[10]},"rows":{"min":[10],"max":[10]},"pages":{}}}
  </script>
</div>
---
## Marginal Means

To get the actual marginal means, we just run `get_emmeans()` on our `means` object.

```r
paged_table(
  get_emmeans(means)
)
```

<div data-pagedtable="false">
  <script data-pagedtable-source type="application/json">
{"columns":[{"label":["meth"],"name":[1],"type":["fct"],"align":["left"]},{"label":["inst"],"name":[2],"type":["fct"],"align":["left"]},{"label":["emmean"],"name":[3],"type":["dbl"],"align":["right"]},{"label":["se"],"name":[4],"type":["dbl"],"align":["right"]},{"label":["df"],"name":[5],"type":["dbl"],"align":["right"]},{"label":["conf.low"],"name":[6],"type":["dbl"],"align":["right"]},{"label":["conf.high"],"name":[7],"type":["dbl"],"align":["right"]},{"label":["method"],"name":[8],"type":["chr"],"align":["left"]}],"data":[{"1":"computer","2":"physical science","3":"46","4":"3.04412","5":"30","6":"39.783078","7":"52.21692","8":"Emmeans test"},{"1":"computer","2":"social science","3":"40","4":"3.04412","5":"30","6":"33.783078","7":"46.21692","8":"Emmeans test"},{"1":"computer","2":"history","3":"38","4":"3.04412","5":"30","6":"31.783078","7":"44.21692","8":"Emmeans test"},{"1":"standard","2":"physical science","3":"34","4":"3.04412","5":"30","6":"27.783078","7":"40.21692","8":"Emmeans test"},{"1":"standard","2":"social science","3":"12","4":"3.04412","5":"30","6":"5.783078","7":"18.21692","8":"Emmeans test"},{"1":"standard","2":"history","3":"31","4":"3.04412","5":"30","6":"24.783078","7":"37.21692","8":"Emmeans test"}],"options":{"columns":{"min":{},"max":[10]},"rows":{"min":[10],"max":[10]},"pages":{}}}
  </script>
</div>

---
## Univariate Comparisons

Since we have an interaction, we are going to check the significance of method (`meth`) on each level of instruction (`inst`). We are telling it to use SS and df from our full model with `error = model`.

```r
model <- lm(vocab ~ inst + meth + inst:meth, data = data)

data %>%
  group_by(inst) %>%
  anova_test(vocab ~ meth, error = model)
```

```
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
```

We can instead test the effect of *instruction* on each level of *method* by switching the two IVs.

```r
data %>%
  group_by(meth) %>%
  anova_test(vocab ~ inst, error = model)
```

```
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
```

```
## # A tibble: 2 x 8
##   meth     Effect   DFn   DFd     F         p `p<.05`   ges
## * <fct>    <chr>  <dbl> <dbl> <dbl>     <dbl> <chr>   <dbl>
## 1 computer inst       2    30  1.87 0.172     ""      0.111
## 2 standard inst       2    30 15.4  0.0000255 "*"     0.506
```
---
## Interaction Plots

This the most straightforward way to produce an interaction plot. For those interested, I'll add a slide on how to extract means for use in `ggplot` by early next week.

```r
with(data, {
interaction.plot(x.factor = inst,
                 trace.factor = meth,
                 response = vocab,
                 type = "l")
})
```

<img src="wk4_files/figure-html/unnamed-chunk-15-1.png" style="display: block; margin: auto;" />
---

## Final Notes
The tough thing about running this in R is having to ask for everything and then having lots of separate output to review. For now it'll be good to know how to ask for something specific and not rely on something too automated. However, I'll look into finding or writing some custom functions that help tidy up all this output.