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.
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
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)
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.
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).
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.
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.
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()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.