---
title: "Effect Sizes with ART"
author: "Matthew Kay"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Effect Sizes with ART}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

## Introduction

The aligned-rank transform (ART) allows for non-parametric analyses of variance. 
But how do we derive effect sizes from ART results? 

__NOTE:__ Before embarking down the path of calculating standardized effect sizes, it is
always worth asking if that is what you really want. Carefully consider, for example, the arguments
of [Cummings (2011)](https://dx.doi.org/10.1001/archpediatrics.2011.97) against the use of
standardized effect sizes and in favor of simple (unstandardized) effect sizes. If you decide you would rather
use simple effect sizes, you may need to consider a different procedure than ART, as the ranking
procedure destroys the information necessary to calculate simple effect sizes.

## Contents

1. [Test Dataset](#test-dataset): The test data we will use to compare a linear model against ART
1. [Partial _eta_-squared](#partial-eta-squared): Calculation of partial _eta_-squared (effect size for _F_ tests)
1. [Cohen's _d_](#cohens-d): Calculation of standardized mean differences (Cohen's _d_; effect size for _t_ tests), including confidence intervals

## Libraries needed for this

```{r setup, include=FALSE}
knitr::opts_chunk$set(  #default code chunk options
    fig.width = 6,
    fig.height = 4
)           
pander::panderOptions("table.split.table", Inf)     #don't split wide tables in output
pander::panderOptions("table.style", "rmarkdown")   #table style that's supported by github
```

```{r message=FALSE}
library(dplyr)      #%>%
library(emmeans)    #emmeans
library(DescTools)  #EtaSq
library(car)        #sigmaHat
library(ARTool)     #art, artlm
library(ggplot2)    #ggplot, stat_..., geom_..., etc
```


## Test dataset

Let's load the test dataset from <code>[vignette("art-contrasts")](art-contrasts.html)</code>:

```{r}
data(InteractionTestData, package = "ARTool")
df = InteractionTestData    #save some typing
```  

Let's fit a linear model:

```{r}
#we'll be doing type 3 tests, so we want sum-to-zero contrasts
options(contrasts = c("contr.sum", "contr.poly"))
m.linear = lm(Y ~ X1*X2, data=df)
```

Now with ART:

```{r}
m.art = art(Y ~ X1*X2, data=df)
```


## Partial _eta_-squared

Note that for Fixed-effects-only models and repeated measures models
(those with `Error()` terms) ARTool also collects the sums of squares, but
does not print them by default. We can pass `verbose = TRUE` to `print()`
to print them:

```{r}
m.art.anova = anova(m.art)
print(m.art.anova, verbose=TRUE)
```

We can use the sums of squares to calculate partial _eta_-squared:

```{r}
m.art.anova$eta.sq.part = with(m.art.anova, `Sum Sq`/(`Sum Sq` + `Sum Sq.res`))
m.art.anova
```

We can compare the above results to partial _eta_-squared calculated on the
linear model (the second column below):

```{r}
EtaSq(m.linear, type=3)
```

The results are comparable.


## Cohen's _d_

We can derive Cohen's _d_ (the standardized mean difference) by dividing estimated differences by the
residual standard deviation of the model. Note that this relies somewhat on the assumption of 
constant variance across levels (aka homoscedasticity).

### in the linear model (for comparison)

As a comparison, let's first derive pairwise contrasts for
all levels of X2 in the linear model:

```{r}
x2.contrasts = summary(pairs(emmeans(m.linear, ~ X2)))
```

Then divide these estimates by the residual standard deviation to get an estimate of _d_:

```{r}
x2.contrasts$d = x2.contrasts$estimate / sigmaHat(m.linear)
x2.contrasts
```

Note that this is essentially the same as the unstandardized estimate for this model;
that is because this test dataset was generated with a residual standard deviation of 1.

### in ART

We can follow the same procedure on the ART model for factor X2:

```{r}
m.art.x2 = artlm(m.art, "X2")
x2.contrasts.art = summary(pairs(emmeans(m.art.x2, ~ X2)))
x2.contrasts.art$d = x2.contrasts.art$estimate / sigmaHat(m.art.x2)
x2.contrasts.art
```

Note how standardization is helping us now: The standardized mean differences (_d_) are
quite similar to the estimates of _d_ from the linear model above.

## Confidence intervals

We can also derive confidence intervals on these effect sizes. To do that, we'll use the `d.ci` function from the `psych` package, which also requires us to indicate how many observations were in each group for each contrast. That is easy in this case, as each group has 100 observations. Thus:

```{r}
x2.contrasts.ci = confint(pairs(emmeans(m.linear, ~ X2))) %>%
    mutate(d = estimate / sigmaHat(m.linear)) %>%
    cbind(d = plyr::ldply(.$d, psych::d.ci, n1 = 100, n2 = 100))

x2.contrasts.ci
```

And from the ART model:

```{r}
x2.contrasts.art.ci = confint(pairs(emmeans(m.art.x2, ~ X2))) %>%
    mutate(d = estimate / sigmaHat(m.art.x2)) %>%
    cbind(d = plyr::ldply(.$d, psych::d.ci, n1 = 100, n2 = 100)) 

x2.contrasts.art.ci
```

And plotting both, to compare (red dashed line is the true effect):

```{r cohens-d-comparison}
rbind(
        cbind(x2.contrasts.ci, model="linear"), 
        cbind(x2.contrasts.art.ci, model="ART")
    ) %>%
    ggplot(aes(x=model, y=d, ymin=d.lower, ymax=d.upper)) +
    geom_pointrange() +
    geom_hline(aes(yintercept = true_effect), 
      data = data.frame(true_effect = c(-2, -2, 0), contrast = c("C - D", "C - E", "D - E")), 
      linetype = "dashed", color = "red") +
    facet_grid(contrast ~ .) + 
    coord_flip()
```