In mjskay/ARTool: Aligned Rank Transform

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) 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: The test data we will use to compare a linear model against ART
2. Partial eta-squared: Calculation of partial eta-squared (effect size for F tests)
3. Cohen's d: Calculation of standardized mean differences (Cohen's d; effect size for t tests), including confidence intervals

Libraries needed for this

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

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 vignette("art-contrasts"):

data(InteractionTestData, package = "ARTool")
df = InteractionTestData    #save some typing

Let's fit a linear model:

#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:

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:

m.art.anova = anova(m.art)
print(m.art.anova, verbose=TRUE)

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

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):

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:

x2.contrasts = summary(pairs(emmeans(m.linear, ~ X2)))

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

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:

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:

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:

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):

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()

mjskay/ARTool documentation built on Oct. 15, 2021, 1:27 p.m.