Método no paramétrico ART en diseños factoriales en bloques y su ajuste post hoc
Autor: Aldo Richard Meza Rodríguez
15/2/2022
En esta investigación se desarrolla una prueba no paramétrica para diseños factoriales cuando no se cumplen los supuestos básico de normalidad o varianza constante. Esta presentación contiene funciones propias (creadas bajo el enfoque de la literatura), también funciones por defecto que se encuentran en el R. La descripción de las funciones y métodos se encuentra en el siguiente artículo: DOI: 10.21704/ac.v82i2.1795
CONTENIDO
CARGANDO LIBRERÍAS
library(ART)
library(ARTool)## Warning: package 'ARTool' was built under R version 4.0.4## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod carlibrary(ggpubr)## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 4.0.5library(ggplot2)
library(rstatix)##
## Attaching package: 'rstatix'## The following object is masked from 'package:stats':
##
## filterlibrary(agricolae)## Warning: package 'agricolae' was built under R version 4.0.5library(splitstackshape)
library(emmeans)
library(nortest)
library(car)## Loading required package: carDatalibrary(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(reshape)
library(tibble)## Warning: package 'tibble' was built under R version 4.0.5library(stringr)
library(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:rstatix':
##
## selectlibrary(DescTools)## Warning: package 'DescTools' was built under R version 4.0.4##
## Attaching package: 'DescTools'## The following object is masked from 'package:car':
##
## Recodelibrary(cowplot)##
## Attaching package: 'cowplot'## The following object is masked from 'package:reshape':
##
## stamp## The following object is masked from 'package:ggpubr':
##
## get_legendELABORACIÓN DEL CONJUNTO DE DATOS DEL DISEÑO FACTORIAL
Factor A (VARIEDAD): Varidades de quinua (Salcedo, Altiplano y Quillahuamán)
Factor B (FERTILIZANTE): Tipos de fertilizante (DF1=250N-150P-120K kg/ha, DF2=180N-120P-100K kg/ha, DF3=200N-130P-90K kg/ha)
Bloque: Densidades de semilla (Block I, Block II, Block III, Block IV)
Variable de respuesta: Altura de la planta
Height <- c(166.68, 177.99, 170.74, 188.72, 153.60, 159.40, 154.18, 178.54, 145.98, 160.60, 162.11,
169.79, 163.85, 169.94, 164.72, 187.63, 188.80, 197.21, 201.27, 209.97, 151.68, 159.51,
163.86, 176.62, 196.19, 205.18, 197.06, 214.17, 209.24, 215.33, 210.11, 230.12, 234.69,
221.42, 211.99, 233.74)
Variety <- c("Salcedo","Salcedo","Salcedo","Salcedo","Salcedo","Salcedo","Salcedo","Salcedo",
"Salcedo","Salcedo","Salcedo","Salcedo","Altiplano","Altiplano","Altiplano",
"Altiplano","Altiplano","Altiplano","Altiplano","Altiplano","Altiplano","Altiplano",
"Altiplano","Altiplano","Quillahuaman","Quillahuaman","Quillahuaman","Quillahuaman",
"Quillahuaman","Quillahuaman","Quillahuaman","Quillahuaman","Quillahuaman",
"Quillahuaman","Quillahuaman","Quillahuaman")
Fertilization <- c("DF1","DF1","DF1","DF1","DF2","DF2","DF2","DF2","DF3","DF3","DF3","DF3","DF1",
"DF1","DF1","DF1","DF2","DF2","DF2","DF2","DF3","DF3","DF3","DF3","DF1","DF1",
"DF1","DF1","DF2","DF2","DF2","DF2","DF3","DF3","DF3","DF3")
Block <- c("I", "II", "III", "IV","I", "II", "III", "IV","I", "II", "III", "IV","I", "II", "III",
"IV", "I", "II", "III", "IV","I", "II", "III", "IV","I", "II", "III", "IV","I", "II",
"III", "IV", "I", "II", "III", "IV")
dat <- data.frame(Block,Variety,Fertilization,Height)
dat$Variety=as.factor(dat$Variety)
dat$Fertilization=as.factor(dat$Fertilization)
dat$Block=as.factor(dat$Block)
dat## Block Variety Fertilization Height
## 1 I Salcedo DF1 166.68
## 2 II Salcedo DF1 177.99
## 3 III Salcedo DF1 170.74
## 4 IV Salcedo DF1 188.72
## 5 I Salcedo DF2 153.60
## 6 II Salcedo DF2 159.40
## 7 III Salcedo DF2 154.18
## 8 IV Salcedo DF2 178.54
## 9 I Salcedo DF3 145.98
## 10 II Salcedo DF3 160.60
## 11 III Salcedo DF3 162.11
## 12 IV Salcedo DF3 169.79
## 13 I Altiplano DF1 163.85
## 14 II Altiplano DF1 169.94
## 15 III Altiplano DF1 164.72
## 16 IV Altiplano DF1 187.63
## 17 I Altiplano DF2 188.80
## 18 II Altiplano DF2 197.21
## 19 III Altiplano DF2 201.27
## 20 IV Altiplano DF2 209.97
## 21 I Altiplano DF3 151.68
## 22 II Altiplano DF3 159.51
## 23 III Altiplano DF3 163.86
## 24 IV Altiplano DF3 176.62
## 25 I Quillahuaman DF1 196.19
## 26 II Quillahuaman DF1 205.18
## 27 III Quillahuaman DF1 197.06
## 28 IV Quillahuaman DF1 214.17
## 29 I Quillahuaman DF2 209.24
## 30 II Quillahuaman DF2 215.33
## 31 III Quillahuaman DF2 210.11
## 32 IV Quillahuaman DF2 230.12
## 33 I Quillahuaman DF3 234.69
## 34 II Quillahuaman DF3 221.42
## 35 III Quillahuaman DF3 211.99
## 36 IV Quillahuaman DF3 233.741. ANÁLISIS EXPLORATORIO
pp1=aggregate(Height ~ Variety + Fertilization,dat, mean )
pp2=aggregate(Height ~ Variety + Fertilization,dat, sd)
dat2=data.frame(pp1,pp2[,3])
names(dat2)=c("Variety","Fertilization","media","desv" )
Plot.1<-ggplot(dat2, aes(x=Variety, y=media, group=Fertilization))+
geom_line(size=1.2, aes(color=Fertilization))+
geom_point(aes(colour=Fertilization), size=3)+
geom_ribbon(aes(ymin=media-desv, ymax=media+desv,fill=Fertilization),alpha=.2)+
#geom_errorbar(aes(ymin = media-desv, ymax=media+desv,color=Fertilization), width = 0.2)+
theme_classic()+
ggpubr::color_palette("jco")+
scale_fill_manual(values = c("#0073C2FF", "#EFC000FF", "#868686FF", "#FC4E07"))+
ylab("Height (cm)")+
xlab("Variety")+
#ggtitle("Tutors and Gender as GPA Predictors")+
theme_bw()+
theme(text = element_text(size=11),
legend.text = element_text(size=11),
legend.direction = "horizontal",
#panel.grid.major = element_blank(),
#panel.grid.minor = element_blank(),
legend.position="top")
Plot.2<-ggplot(dat2, aes(x=Fertilization, y=media, group=Variety))+
geom_line(size=1.2, aes(color=Variety))+
geom_point(aes(colour=Variety), size=3)+
geom_ribbon(aes(ymin=media-desv, ymax=media+desv,fill=Variety),alpha=.2)+
#geom_errorbar(aes(ymin = media-desv, ymax=media+desv,color=Variety), width = 0.2)+
theme_bw()+
ggpubr::color_palette("jco")+
scale_fill_manual(values = c("#0073C2FF", "#EFC000FF", "#868686FF", "#FC4E07"))+
ylab("Height (cm)")+
xlab("Fertilization")+
#ggtitle("Tutors and Gender as GPA Predictors")+
theme(text = element_text(size=11),
legend.text = element_text(size=11),
legend.direction = "horizontal",
#panel.grid.major = element_blank(),
#panel.grid.minor = element_blank(),
legend.position="top")
figure2=ggarrange(Plot.1,Plot.2,nrow = 1,common.legend = F)
final2<-annotate_figure(figure2, top = text_grob("Interaction graph point-line with standard deviation",
color = "black", face = "bold", size = 11),
right = text_grob(bquote(""), rot = 90,size = 6))
final2
bxp <- ggline(
dat, x = "Fertilization", y = "Height",
color = "Variety", palette = "jco",size = 0.8,shape = 19,point.size = 1.5,
facet.by = "Block", short.panel.labs = FALSE,ylab = "Height (cm)",y.size =10,
)
bxp
2. ANOVA CLÁSICO
md1<-lm(Height~Block+Variety*Fertilization,data=dat)
summary(aov(md1))## Df Sum Sq Mean Sq F value Pr(>F)
## Block 3 2087 696 24.95 1.49e-07 ***
## Variety 2 15778 7889 282.96 < 2e-16 ***
## Fertilization 2 681 341 12.22 0.000219 ***
## Variety:Fertilization 4 3865 966 34.66 1.20e-09 ***
## Residuals 24 669 28
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13. PRUEBAS DE SUPUESTOS DEL MODELO
3.1 Gráfico de residuales
i= 1:length(md1$residuals)
fi <- (i - 0.5) / length(md1$residuals)
x.norm <- sort(qnorm(fi))
da=data.frame(y=sort(md1$residuals),yvar=md1$residuals,x=x.norm, xvar=md1$fitted.values)
qqnorm=ggplot(data = da, aes(x = x, y = y,color="blue")) +
geom_point(size = 4, alpha = 0.5)+
scale_color_manual("Origin",values = c("#0073C2FF"))+
#geom_smooth(method = "lm",se = F, size = 0.2, linetype = 2)+
geom_abline(intercept = -1.496e-15, slope = 2.8, color="#0073C2FF", size=0.5,linetype = 2)+
#lm(da$y~da$x)
#lims( y = c(-3,3))+
theme_light()+
theme(legend.position="none")+
labs(x="Theoritical Quantiles",y="Standarized residuals")+
ggtitle("Normal Q-Q") +
theme(plot.title = element_text(hjust = 0.5,size=11))
homog=ggplot(data = da, aes(x = xvar, y = yvar,color="blue")) +
geom_point(size = 4, alpha = 0.5)+
scale_color_manual("Origin",values = c("#0073C2FF"))+
geom_smooth(method = "lm",formula=y ~ poly(x,4), se = F, size = 0.5, color= "red",linetype = 2)+
geom_abline(intercept = -6.702e-15, slope =3.560e-17, color="black", size=0.5,linetype = 2)+
#lm(y ~ x)
#lims( y = c(-3,3))+
theme_light()+
theme(legend.position="none")+
labs(x="Fitted values",y="Residuals")+
ggtitle("Residual vs Fitted") +
theme(plot.title = element_text(hjust = 0.5,size=11))
fig<-ggarrange(homog, qqnorm,ncol = 2)
fig1<-annotate_figure(fig,
top = text_grob("Basic assumptions of normality and homogeneity of variances"
, color = "black", face = "bold", size = 11))
fig1
3.2 Pruebas de normalidad
shapiro.test(rstandard(md1))##
## Shapiro-Wilk normality test
##
## data: rstandard(md1)
## W = 0.89049, p-value = 0.001894ad.test(rstandard(md1))##
## Anderson-Darling normality test
##
## data: rstandard(md1)
## A = 1.1186, p-value = 0.0054493.3 Prueba de varianza constate
ncvTest(md1)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 4.625421, Df = 1, p = 0.031502gqtest(md1)##
## Goldfeld-Quandt test
##
## data: md1
## GQ = 3.3472, df1 = 6, df2 = 6, p-value = 0.08359
## alternative hypothesis: variance increases from segment 1 to 24. Modelo ART (Rango Transformado Alineado)
4.1 Modelo con la función art
model <- anova(art(Height~Variety*Fertilization+Error(Block) ,data=dat))
model## Analysis of Variance of Aligned Rank Transformed Data
##
## Table Type: Repeated Measures Analysis of Variance Table (Type I)
## Model: Repeated Measures (aov)
## Response: art(Height)
##
## Error Df Df.res F value Pr(>F)
## 1 Variety Withn 2 24 126.829 1.7394e-13 ***
## 2 Fertilization Withn 2 24 8.241 0.0018854 **
## 3 Variety:Fertilization Withn 4 24 30.064 4.9486e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14.2 Modelo completo (con todas sus fuente de variabilidad)
art.fact.block <- function(formula,data){
if (mode(formula)!="call") stop("invalid formula")
resp <- rstatix:::get_formula_left_hand_side(formula)
elem <- rstatix:::get_formula_right_hand_side(formula)
fact0 <- as.vector(levels(as.factor(unlist(strsplit(elem, "[:]")))))
int <- gsub(":", "*", elem[length(elem)])
fact <- c(as.vector(unlist(strsplit(int, "[*]"))))
nfac <-length(fact)
block <-fact0[(fact0 %in% fact)==F]
form <- as.formula(paste0(resp,"~",int,"+ Error(",block,")"))
art. <- art(form,data=data) # Model art factorial in Block
anv <- anova(art.) # Model anova in art factorial in Block
df <- data.frame(data,art.$aligned.ranks)
rsl <- list()
for(i in 1:dim(art.$aligned.ranks)[2]){
rsl[[i]] <- summary(aov(as.formula(paste0(names(df[,
(dim(data)[2]+1):(dim(df)[2])][i]),
"~",int,"+ Error(",block,")")),data=df))}
rsl1 <- matrix(0,2^nfac-1,6)
for(i in 1:dim(art.$aligned.ranks)[2]){
Blq <- c(unlist(rsl[[i]][[1]])[1:3],unlist(rsl[[i]][[2]])[(2^nfac)*3])
Fval <- round(Blq[3]/Blq[4],3)
Pval <- 1-pf(as.numeric(Fval),as.numeric(Blq[1]),anv$Df.res[1])
rsl1[i,] <- as.numeric(c(round(Blq,3),Fval,Pval=Pval))}
anv.blk <- data.frame(Df.res=anv$Df.res,rsl1)
fac<-cbind(round(cbind(anv$Df,anv$`Sum Sq`,anv$`Mean Sq`,
anv$`F value`),3),anv$`Pr(>F)`)
art.block<-cbind(round(cbind(anv.blk$X1,anv.blk$X2,anv.blk$X3,
anv.blk$X5),3),anv.blk$X6)
res<-cbind(anv$Df.res,round(anv$`Sum Sq.res`,3), anv.blk$X4,NA,NA)
aov.art=rbind(fac,art.block,res)
Fv=c(anv$Term, paste0(block,"(",anv$Term,")"),
paste0("residual","(",anv$Term,")"))
aov.fct.bl=data.frame(Fv,aov.art)
names(aov.fct.bl) <- c("Fv","Df","Sum.Sq", "Mean.Sq", "F.value", "Pvalue")
aov.fct.bl$Sig <- with(aov.fct.bl, ifelse(Pvalue <= 0.001,"***",
ifelse(Pvalue <= 0.01,"**", ifelse(Pvalue <= 0.05,"*",
ifelse(Pvalue <= 0.1,".","ns")))))
aov.fct.bl$Sig[is.na(aov.fct.bl$Sig)]<-""
return(c(art.=list(art.),anova=list(anv),Block=list(aov.fct.bl)))
}
art.fact.block(Height~Variety*Fertilization+Block,dat)## $art.
## Aligned Rank Transform of Factorial Model
##
## Call:
## art(formula = form, data = data)
##
## Column sums of aligned responses (should all be ~0):
## Variety Fertilization Variety:Fertilization
## 0 0 0
##
## F values of ANOVAs on aligned responses not of interest (should all be ~0):
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
##
## $anova
## Analysis of Variance of Aligned Rank Transformed Data
##
## Table Type: Repeated Measures Analysis of Variance Table (Type I)
## Model: Repeated Measures (aov)
## Response: art(Height)
##
## Error Df Df.res F value Pr(>F)
## 1 Variety Withn 2 24 126.829 1.7394e-13 ***
## 2 Fertilization Withn 2 24 8.241 0.0018854 **
## 3 Variety:Fertilization Withn 4 24 30.064 4.9486e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Block
## Fv Df Sum.Sq Mean.Sq F.value Pvalue Sig
## 1 Variety 2 2944.667 1472.333 126.829 1.739436e-13 ***
## 2 Fertilization 2 630.167 315.083 8.241 1.885388e-03 **
## 3 Variety:Fertilization 4 2222.500 555.625 30.064 4.948556e-09 ***
## 4 Block(Variety) 3 659.889 219.963 18.948 1.608522e-06 ***
## 5 Block(Fertilization) 3 2317.889 772.630 20.208 9.373761e-07 ***
## 6 Block(Variety:Fertilization) 3 1207.444 402.481 21.778 4.939273e-07 ***
## 7 residual(Variety) 24 278.611 11.609 NA NA
## 8 residual(Fertilization) 24 917.611 38.234 NA NA
## 9 residual(Variety:Fertilization) 24 443.556 18.481 NA NA5. Efecto del diseño del modelo ART
Eff.size.art<-function(formula,data){
rs <- art.fact.block(formula,data)
nter<-length(rs$anova$Term)
scb <- rs$Block[(nter+1):(2*nter),3]
eta.sq.part.art <- rs$anova$`Sum Sq`/(rs$anova$`Sum Sq` +
rs$anova$`Sum Sq.res`)
eta.sq.gn.art <- rs$anova$`Sum Sq`/(1*rs$anova$`Sum Sq` +
scb + rs$anova$`Sum Sq.res`)
eta.sq.gn.art.block <- scb /(0*rs$anova$`Sum Sq`+
scb + rs$anova$`Sum Sq.res`)
eta2.lm <- EtaSq(aov(rs$art.$formula,data = dat), type=1)[,2:3]
eta2.fn<-data.frame(eta2.lm,eta.sq.part.art,eta.sq.gn.art,
eta.sq.gn.art.block)
return(eta2.fn)
}
Eff.size.art(Height~Variety*Fertilization + Block,dat)## eta.sq.part eta.sq.gen eta.sq.part.art eta.sq.gn.art
## Variety 0.9593167 0.8512844 0.9135628 0.7583158
## Fertilization 0.5045177 0.1981939 0.4071429 0.1630163
## Variety:Fertilization 0.8524428 0.5837529 0.8336285 0.5737706
## eta.sq.gn.art.block
## Variety 0.7031315
## Fertilization 0.7163928
## Variety:Fertilization 0.73134126. Eficiencia relativa del modelo ART
ER.art<-function(formula,data){
art.bl<-art.fact.block(formula,data)
rs <- art.fact.block(formula,data)
nter<-length(art.bl$anova$Term)
scb <- rs$Block[(nter+1),2]
nb<-art.bl$Block[(nter+1),2]+1
elem <- rstatix:::get_formula_right_hand_side(formula)
fact <- data[,c(as.vector(unlist(strsplit(gsub(":", "*",
elem[length(elem)]), "[*]")))),]
n<-c()
f<-list()
for (i in 1:dim(fact)[2]){
for (j in 1:dim(fact)[2]){
n[i]<-nlevels(fact[,i])
f[[j]]<-apply(t(combn(n,j)), 1, prod)}}
r<-unlist(f)
CM.Error<-art.bl$Block[(2*nter+1):(3*nter),4]
CM.Bloq<-art.bl$Block[(nter+1):(2*nter),4]
ER <- (((nb - 1)*CM.Bloq) + (nb*(r - 1)*CM.Error))/(((nb*r) - 1)*CM.Error)
gl1<-(r-1)*(nb-1)
gl2<-r*nb-r
ER.adj<-((gl1+1)*(gl2+3))/((gl1+3)*(gl2+1))*ER
res<-data.frame(FV=art.bl$anova$Term,ER,ER.adj)
res}
ER.art(Height~Variety*Fertilization + Block,dat)## FV ER ER.adj
## 1 Variety 5.894807 5.501820
## 2 Fertilization 6.238526 5.822625
## 3 Variety:Fertilization 2.780980 2.7589087. Post hoc (comparaciones múltiples)
7.1 Efectos principales
m.art<-art(Height~Variety*Fertilization+Error(Block),data=dat)
art.con(m.art,"Variety", adjust="tukey")## NOTE: Results may be misleading due to involvement in interactions## contrast estimate SE df t.ratio p.value
## Altiplano - Quillahuaman -14.17 1.39 24 -10.185 <.0001
## Altiplano - Salcedo 7.67 1.39 24 5.512 <.0001
## Quillahuaman - Salcedo 21.83 1.39 24 15.696 <.0001
##
## Results are averaged over the levels of: Fertilization
## P value adjustment: tukey method for comparing a family of 3 estimatesart.con(m.art,"Fertilization", adjust="tukey")## NOTE: Results may be misleading due to involvement in interactions## contrast estimate SE df t.ratio p.value
## DF1 - DF2 -8.8333 2.52 24 -3.499 0.0051
## DF1 - DF3 0.0833 2.52 24 0.033 0.9994
## DF2 - DF3 8.9167 2.52 24 3.532 0.0047
##
## Results are averaged over the levels of: Variety
## P value adjustment: tukey method for comparing a family of 3 estimates7.2 Efectos de interacción
art.con(m.art,"Variety:Fertilization", adjust="tukey")## contrast estimate SE df t.ratio p.value
## Altiplano,DF1 - Altiplano,DF2 -11.25 1.78 24 -6.307 <.0001
## Altiplano,DF1 - Altiplano,DF3 4.75 1.78 24 2.663 0.2133
## Altiplano,DF1 - Quillahuaman,DF1 -12.25 1.78 24 -6.868 <.0001
## Altiplano,DF1 - Quillahuaman,DF2 -17.25 1.78 24 -9.671 <.0001
## Altiplano,DF1 - Quillahuaman,DF3 -20.25 1.78 24 -11.353 <.0001
## Altiplano,DF1 - Salcedo,DF1 -2.75 1.78 24 -1.542 0.8251
## Altiplano,DF1 - Salcedo,DF2 5.75 1.78 24 3.224 0.0725
## Altiplano,DF1 - Salcedo,DF3 6.00 1.78 24 3.364 0.0539
## Altiplano,DF2 - Altiplano,DF3 16.00 1.78 24 8.970 <.0001
## Altiplano,DF2 - Quillahuaman,DF1 -1.00 1.78 24 -0.561 0.9997
## Altiplano,DF2 - Quillahuaman,DF2 -6.00 1.78 24 -3.364 0.0539
## Altiplano,DF2 - Quillahuaman,DF3 -9.00 1.78 24 -5.046 0.0010
## Altiplano,DF2 - Salcedo,DF1 8.50 1.78 24 4.765 0.0021
## Altiplano,DF2 - Salcedo,DF2 17.00 1.78 24 9.531 <.0001
## Altiplano,DF2 - Salcedo,DF3 17.25 1.78 24 9.671 <.0001
## Altiplano,DF3 - Quillahuaman,DF1 -17.00 1.78 24 -9.531 <.0001
## Altiplano,DF3 - Quillahuaman,DF2 -22.00 1.78 24 -12.334 <.0001
## Altiplano,DF3 - Quillahuaman,DF3 -25.00 1.78 24 -14.016 <.0001
## Altiplano,DF3 - Salcedo,DF1 -7.50 1.78 24 -4.205 0.0079
## Altiplano,DF3 - Salcedo,DF2 1.00 1.78 24 0.561 0.9997
## Altiplano,DF3 - Salcedo,DF3 1.25 1.78 24 0.701 0.9983
## Quillahuaman,DF1 - Quillahuaman,DF2 -5.00 1.78 24 -2.803 0.1659
## Quillahuaman,DF1 - Quillahuaman,DF3 -8.00 1.78 24 -4.485 0.0041
## Quillahuaman,DF1 - Salcedo,DF1 9.50 1.78 24 5.326 0.0005
## Quillahuaman,DF1 - Salcedo,DF2 18.00 1.78 24 10.091 <.0001
## Quillahuaman,DF1 - Salcedo,DF3 18.25 1.78 24 10.231 <.0001
## Quillahuaman,DF2 - Quillahuaman,DF3 -3.00 1.78 24 -1.682 0.7512
## Quillahuaman,DF2 - Salcedo,DF1 14.50 1.78 24 8.129 <.0001
## Quillahuaman,DF2 - Salcedo,DF2 23.00 1.78 24 12.894 <.0001
## Quillahuaman,DF2 - Salcedo,DF3 23.25 1.78 24 13.034 <.0001
## Quillahuaman,DF3 - Salcedo,DF1 17.50 1.78 24 9.811 <.0001
## Quillahuaman,DF3 - Salcedo,DF2 26.00 1.78 24 14.576 <.0001
## Quillahuaman,DF3 - Salcedo,DF3 26.25 1.78 24 14.716 <.0001
## Salcedo,DF1 - Salcedo,DF2 8.50 1.78 24 4.765 0.0021
## Salcedo,DF1 - Salcedo,DF3 8.75 1.78 24 4.905 0.0015
## Salcedo,DF2 - Salcedo,DF3 0.25 1.78 24 0.140 1.0000
##
## P value adjustment: tukey method for comparing a family of 9 estimates7.3 Efectos simples
formula=Height~Variety*Fertilization
dat$Yart<-artlm.con(art(formula,data=dat), "Variety")$model$.y
model=lm(Yart~Variety*Fertilization+Block,data=dat)
simple<- stats::as.formula(paste(" ~ ", "Variety"," | ","Fertilization"))
contrast(emmeans(model, simple), method="pairwise", adjust="tukey")## Fertilization = DF1:
## contrast estimate SE df t.ratio p.value
## Altiplano - Quillahuaman -14.00 2.41 24 -5.811 <.0001
## Altiplano - Salcedo 7.75 2.41 24 3.217 0.0099
## Quillahuaman - Salcedo 21.75 2.41 24 9.028 <.0001
##
## Fertilization = DF2:
## contrast estimate SE df t.ratio p.value
## Altiplano - Quillahuaman -14.00 2.41 24 -5.811 <.0001
## Altiplano - Salcedo 8.00 2.41 24 3.321 0.0078
## Quillahuaman - Salcedo 22.00 2.41 24 9.132 <.0001
##
## Fertilization = DF3:
## contrast estimate SE df t.ratio p.value
## Altiplano - Quillahuaman -14.50 2.41 24 -6.019 <.0001
## Altiplano - Salcedo 7.25 2.41 24 3.009 0.0161
## Quillahuaman - Salcedo 21.75 2.41 24 9.028 <.0001
##
## Results are averaged over the levels of: Block
## P value adjustment: tukey method for comparing a family of 3 estimatessimple<- stats::as.formula(paste(" ~ ", "Fertilization"," | ","Variety"))
contrast(emmeans(model, simple), method="pairwise", adjust="tukey")## Variety = Altiplano:
## contrast estimate SE df t.ratio p.value
## DF1 - DF2 0.00 2.41 24 0.000 1.0000
## DF1 - DF3 0.50 2.41 24 0.208 0.9766
## DF2 - DF3 0.50 2.41 24 0.208 0.9766
##
## Variety = Quillahuaman:
## contrast estimate SE df t.ratio p.value
## DF1 - DF2 0.00 2.41 24 0.000 1.0000
## DF1 - DF3 0.00 2.41 24 0.000 1.0000
## DF2 - DF3 0.00 2.41 24 0.000 1.0000
##
## Variety = Salcedo:
## contrast estimate SE df t.ratio p.value
## DF1 - DF2 0.25 2.41 24 0.104 0.9941
## DF1 - DF3 0.00 2.41 24 0.000 1.0000
## DF2 - DF3 -0.25 2.41 24 -0.104 0.9941
##
## Results are averaged over the levels of: Block
## P value adjustment: tukey method for comparing a family of 3 estimates7.4 Gráfica de Efectos principales y efectos de interacción
Graph.EFF_Main<-function(formula,data,main,graph,p.adjust.method,resul.var)
{
df=data
main<-gsub("\\*",":",main)
main1=unlist(strsplit(main, "[:]"))
vars.dep <- rstatix:::get_formula_left_hand_side(formula)
vars.ind <- rstatix:::get_formula_right_hand_side(formula)
inter<-vars.ind[length(vars.ind)]
inter.<-gsub("\\:", "*",inter)
main0=ifelse(main==inter,"inter",main)
main01=ifelse(main==inter|main==inter.,inter.,main)
inter..=unlist(strsplit(inter., "[*]"))
df$inter =interaction(df[,inter..[1]], df[,inter..[2]], sep = ":")
data$inter =interaction(df[,inter..[1]], df[,inter..[2]], sep = ":")
formula2<-as.formula(paste0(vars.dep," ~ ",inter.))
df[,vars.dep]<-artlm.con(art(formula2,data=df), main)$model$.y
model2=lm(formula,data=df)
formula0<-as.formula(paste0(" ~ ",main01))
comp=contrast(emmeans(model2, formula0), method = "pairwise",adjust=p.adjust.method)
mes <- tapply.stat(df[,vars.dep], df[,main0], stat = "mean")
comp1<-data.frame(comp)
ntr=nlevels(df[,main0])
alpha=0.05
Q <- matrix(1, ncol = ntr, nrow = ntr)
p <- comp1$p.value
k <- 0
for (i in 1:(ntr - 1)) {
for (j in (i + 1):ntr) {
k <- k + 1
Q[i, j] <- p[k]
Q[j, i] <- p[k]
}
}
Comp.fin <- data.frame(orderPvalue(mes[, 1], mes[, 2], alpha, Q))
Dep=data.frame(t=rownames(Comp.fin),Comp.fin,
id=rep(1:nlevels(as.factor(df[,main0]))),row.names=NULL)
mean. <- tapply(data[,vars.dep], data[,main0],"mean")
sd. <- tapply(data[,vars.dep], data[,main0],"sd")
max. <- tapply(data[,vars.dep], data[,main0],"max")
measure<-cbind(mean.,sd.,max.)
Dep1=data.frame(t=rownames(measure),measure,
row.names=NULL)
mean.. <- tapply(df[,vars.dep], df[,main0],"mean")
sd.. <- tapply(df[,vars.dep], df[,main0],"sd")
max.. <- tapply(df[,vars.dep], df[,main0],"max")
measure1<-cbind(mean..,sd..,max..)
Dep2=data.frame(t=rownames(measure1),measure1,
row.names=NULL)
if (resul.var == "Y.orig"){
hh=merge(Dep1, Dep, by.x = 't', by.y='t')
medida<-paste(vars.dep,"(cm)")}
if (resul.var == "Y.art"){
hh=merge(Dep2, Dep, by.x = 't', by.y='t')
medida<-paste(vars.dep," (aligned ranks)")}
ff=data.frame(cSplit(hh, "t", "."))[,c(6,7,1:5)][order(hh$id),]
mm=rep(vars.dep,each=nlevels(as.factor(data[,main0])))
df0=data.frame(ff,f=mm)[,-c(1,6)]
names(df0)=c(main0,"mean","desv","max","groups","Medidas")
if (resul.var == "Y.orig"){
df1=data.frame(merge(data[,c(main0,vars.dep)], df0, by =main0))
names(df1)=c("trat","response",names(df1)[-c(1,2)])}
if (resul.var == "Y.art"){
df1=data.frame(merge(df[,c(main0,vars.dep)], df0, by =main0))
names(df1)=c("trat","response",names(df1)[-c(1,2)])}
p1=(ggplot (df1, aes(x=reorder(trat,-response), y=response,color=groups))+
geom_boxplot() +
geom_point(size=2) +
#stat_summary(fun=mean, geom="point", size=4) +
ggpubr::color_palette("jco")+
theme_classic()+
theme(legend.position="none") +
geom_text(aes(x=trat, y=max+desv,label=groups,color=groups),
position=position_dodge(width=0.9),alpha=1, size=4)+
#ggtitle(paste0("Post hoc",method," interaction Variety:Fertilization")) +
geom_hline(yintercept = median(df1$response), linetype = 2,color="light blue")+
ylab(medida)+ xlab("")+
theme(plot.title = element_text(hjust = 0.5,color = "black",size=11,face="bold"),
axis.title=element_text(size=10))+
rotate_x_text())
p2=(ggplot (df1, aes(x=reorder(trat,-response), y=response,color=groups))+
geom_violin(trim=F) +
#geom_boxplot(width=0.15) +
geom_point(pch="-", size=4) +
stat_summary(fun=mean, geom="point", size=4) +
ggpubr::color_palette("jco")+
theme_classic()+
#theme_bw()+
theme(legend.position="none",
panel.border = element_rect(colour = "#868686FF", fill=NA, size=0.5),
axis.line=element_line()) +
geom_text(aes(x=trat, y=1.07*max+desv,label=groups,color=groups),
position=position_dodge(width=0.9),alpha=1, size=4)+
#ggtitle(paste0("Post hoc",method," interaction Variety:Fertilization")) +
geom_hline(yintercept = median(df1$response), linetype = 2,color="blue")+
ylab(medida) + xlab("") +
ggtitle(paste0("")) +
theme(plot.title = element_text(hjust = 0.5,color = "black",size=11,face="bold"),
axis.title=element_text(size=10)) + rotate_x_text()
)
if (graph == "violin") {
return(p2)
}
if (graph == "boxplot") {
return(p1)
}
if (graph == "ambos"){
p<-ggarrange(p1,p2,ncol = 2, nrow = 1)
return(p)
}
}
Graph.EFF_Main(Height~Variety*Fertilization+Block,dat,main="Variety","violin","tukey","Y.orig")## NOTE: Results may be misleading due to involvement in interactions
Graph.EFF_Main(Height~Variety*Fertilization+Block,dat,main="Fertilization","violin","tukey","Y.orig")## NOTE: Results may be misleading due to involvement in interactions
Graph.EFF_Main(Height~Variety*Fertilization+Block,dat,main="Variety:Fertilization","violin","tukey","Y.orig")
7.5 Gráfica de Efectos simples
Artool_EFFS2<-function(data,formula,Block,grouping.vars,graph,p.adjust.method,resul.var)
{
df0<-dat
vars1 = names(attr(aov(formula,data)$terms,"dataClasses"))
vars=vars1[vars1!= paste0(grouping.vars)]
group <- vars[2]
data <- data %>% rstatix:::.as_factor(group)
data <- data %>% rstatix:::.as_factor(grouping.vars)
intr<- rstatix:::get_formula_right_hand_side(formula)[3]
inter<-paste0(vars1[2],"*",vars1[3])
resp<-rstatix:::get_formula_left_hand_side(formula)
form0 <- as.formula(paste0(resp,"~",inter,"+ Error(",Block,")"))
data$vartool <- artlm.con(art(formula,data=data), intr)$model$.y
data=data[ , !(names(data) %in% vars[1])]
outcome<-"vartool"
formula <- stats::as.formula(paste(outcome, " ~ ", inter,"+",Block))
date=data.frame(data)
names(date)[names(date) == group] <- "groupp"
names(df0)[names(df0) == group] <- "groupp"
names(df0)[names(df0) == resp] <- "vartool"
date$t=as.factor(paste(date[,grouping.vars] , date$groupp, sep="."))
df0$t=as.factor(paste(df0[,grouping.vars] , df0$groupp, sep="."))
model.fact.bloq <- lm(formula, data=data)
simple<- stats::as.formula(paste(" ~ ", group," | ",grouping.vars))
pwc=data.frame(melt(contrast(emmeans(model.fact.bloq, simple), method="pairwise", adjust=p.adjust.method)))
names(pwc)=gsub("\\value.","",names(pwc))
dd1=melt(do.call(cbind,lapply(split(date,date[,grouping.vars]),
function(x) with(x,by(vartool,groupp,mean)))))
dd2=melt(do.call(cbind,lapply(split(date,date[,grouping.vars]),
function(x) with(x,by(vartool,groupp,sd)))))
dd3=melt(do.call(cbind,lapply(split(date,date[,grouping.vars]),
function(x) with(x,by(vartool,groupp,max)))))
dg=cbind(dd1,dd2[,3],dd3[,3],vars[1])
dg$t=paste0(dg$X2,".",dg$X1)
df1=split(dg,dg[,"X2"])
groupp<-"groupp"
fd=function(p.value){
ntr=nlevels(date[,groupp])
Q <- matrix(1, ncol = ntr, nrow = ntr)
p <- p.value
k <- 0
for (i in 1:(ntr - 1)) {
for (j in (i + 1):ntr) {
k <- k + 1
Q[i, j] <- p[k]
Q[j, i] <- p[k]
}
}
return(Q)
}
ddd=lapply(split(pwc,pwc[,grouping.vars]),
function(x) with(x
,fd(p.value)))
alpha=0.05
re=list()
for(i in levels(date[,grouping.vars])){
re[[i]]=data.frame(orderPvalue(df1[[i]]$X1,df1[[i]]$value,alpha,ddd[[i]]),
groupp=rownames(orderPvalue(df1[[i]]$X1,df1[[i]]$value,alpha,ddd[[i]])),
row.names=NULL)
}
df2=melt(re)
df2$t <- paste(df2$L1 , df2$groupp, sep=".")
df2=cbind(df2,id=rep(1:dim(df2)[1]))
d1=melt(do.call(cbind,lapply(split(df0,df0[,grouping.vars]),
function(x) with(x,by(vartool,groupp,mean)))))
d2=melt(do.call(cbind,lapply(split(df0,df0[,grouping.vars]),
function(x) with(x,by(vartool,groupp,sd)))))
d3=melt(do.call(cbind,lapply(split(df0,df0[,grouping.vars]),
function(x) with(x,by(vartool,groupp,max)))))
dg1=cbind(d1,d2[,3],d3[,3],vars[1])
dg1$t=paste0(dg1$X2,".",dg1$X1)
if (resul.var == "Y.art") {
hh=merge(dg,df2, by.x = 't', by.y='t')
ff=data.frame(hh[order(hh$id),][,c(1,3,2,4:6,8,7)])
ff$t=as.factor(ff$t)
names(ff)=c("t","grouping.vars","group1","mean","desv","max","groups","Medidas")
df3=data.frame(merge(date[,c("t","vartool")], ff, by ="t"))
medida<-paste(resp," (aligned ranks)")
}
if (resul.var == "Y.orig") {
hh=merge(dg1,df2, by.x = 't', by.y='t')
ff=data.frame(hh[order(hh$id),][,c(1,3,2,4:6,8,7)])
ff$t=as.factor(ff$t)
names(ff)=c("t","grouping.vars","group1","mean","desv","max","groups","Medidas")
df3=data.frame(merge(df0[,c("t","vartool")], ff, by ="t"))
medida<-paste(resp,"(cm)")
}
p1 <- list()
for (j in levels(df3$grouping.vars)) {
p1[[j]] <-
(ggplot (data <- subset(df3,grouping.vars==j), aes(x=reorder(group1,-mean), y=mean,color=groups)) +
geom_pointrange(aes(ymin=mean-desv, ymax=mean+desv,color=groups),
size=.5,position=position_dodge(.5)) +
ggpubr::color_palette("jco")+
rotate_x_text()+ theme_bw()+
geom_point()+
geom_hline(yintercept = mean(df3$vartool), linetype = 2,color="light blue")+
geom_text(aes(x=group1, y=(mean+desv)+mean(df3$mean)*0.15,label=groups,color=groups),
position=position_dodge(width=0.9),alpha=1, size=4)+
theme(plot.margin=margin(unit(c(8,4,1,-8),"cm")))+ # 1 top, 2 rigth, 3 left, 4 down # cambiar a mm si es en decimales
facet_grid(. ~ grouping.vars)+
theme(legend.position="none")+
lims( y = c(min(aggregate(mean ~ grouping.vars+group1,
df3, mean )$mean-aggregate(desv ~ grouping.vars+group1, df3, mean )$desv),
max(df3$vartool)+0.24*(mean(df3$vartool))))+
labs(x="",y="")+
theme(plot.title = element_text(hjust = 0.5)))}
figure1<-ggarrange(plotlist = p1,ncol = nlevels(df3$grouping.vars))
final1<-annotate_figure(figure1,
top = text_grob(paste("Post hoc",p.adjust.method,"Simple Effect:",
group,"|",grouping.vars)
, color = "black", face = "bold", size = 11),
left = text_grob(medida, color = "black",size = 10, rot = 90),
right = text_grob(bquote(""), rot = 90,size = 6))
p2 <- list()
for (j in levels(df3$grouping.vars)) {
p2[[j]] <-
(ggplot (data <- subset(df3,grouping.vars==j), aes(x=reorder(group1,-vartool), y=vartool,color=groups)) +
geom_boxplot() +
ggpubr::color_palette("jco")+
rotate_x_text()+ theme_bw()+
geom_point()+
geom_hline(yintercept = median(df3$vartool), linetype = 2,color="light blue")+
geom_text(aes(x=group1, y=max+mean(df3$max)*0.15,label=groups,color=groups),
position=position_dodge(width=0.9),alpha=1, size=4)+
theme(plot.margin=margin(unit(c(8,4,1,-8),"cm")))+ # 1 top, 2 rigth, 3 left, 4 down # cambiar a mm si es en decimales
facet_grid(. ~ grouping.vars)+
theme(legend.position="none")+
lims( y = c(min(df3$vartool)-0.08*(mean(df3$vartool)), max(df3$vartool)+0.24*(mean(df3$vartool))))+
labs(x="",y="")+
theme(plot.title = element_text(hjust = 0.5)))}
figure2<-ggarrange(plotlist = p2,ncol = nlevels(df3$grouping.vars))
final2<-annotate_figure(figure2,
top = text_grob(paste("Post hoc",p.adjust.method,"Simple Effect:",
group,"|",grouping.vars)
, color = "black", face = "bold",size = 11), #
left = text_grob(medida, color = "black",size = 10, rot = 90))
p3 <- list()
for (j in levels(df3$grouping.vars)) {
p3[[j]] <-
(ggplot (data <- subset(df3,grouping.vars==j), aes(x=reorder(group1,-vartool), y=vartool,color=groups)) +
geom_violin(trim=F,adjust=1.3) +
geom_point(pch="-", size=4) +
stat_summary(fun=mean, geom="point", size=4) +
ggpubr::color_palette("jco")+
#rotate_x_text()+
theme_bw()+
#theme_classic()+
geom_hline(yintercept = median(df3$vartool), linetype = 2,color="blue")+
geom_text(aes(x=group1, y=max+3.5*desv,label=groups,color=groups),
position=position_dodge(width=0.9),alpha=1, size=4)+
theme(plot.margin=margin(unit(c(8,4,1,-8),"cm")))+ # 1 top, 2 rigth, 3 left, 4 down # cambiar a mm si es en decimales
facet_grid(. ~ grouping.vars)+
theme(legend.position="none")+
lims( y = c(min(df3$vartool)-0.25*(mean(df3$vartool)), max(df3$vartool)+1.6*(sd(df3$vartool))))+
labs(x="",y="")+
theme(plot.title = element_text(hjust = 0.5)))}
figure3<-ggarrange(plotlist = p3,ncol = nlevels(df3$grouping.vars))
final3<-annotate_figure(figure3,
top = text_grob(paste("Post hoc",p.adjust.method,"Simple Effect:",
group,"|",grouping.vars)
, color = "black", face = "bold",size = 11), #
left = text_grob(medida, color = "black",size = 10, rot = 90))
#, right = text_grob(bquote(""), rot = 90,size = 1)
if (graph == "violin") {
return(final3)
}
if (graph == "boxplot") {
return(final2)
}
if (graph == "pointrange") {
return(final1)
}
if (graph == "ambos"){
p<-ggarrange(final1,final2,ncol = 1, nrow = 2)
return(p)
}
}
Artool_EFFS2(dat,Height~Variety*Fertilization,Block="Block", "Variety","violin","tukey","Y.orig")## Using groups, groupp as id variables
## Using groups, groupp as id variables
## Using groups, groupp as id variables
Artool_EFFS2(dat,Height~Variety*Fertilization,Block="Block","Fertilization","violin","tukey","Y.orig")## Using groups, groupp as id variables
## Using groups, groupp as id variables
## Using groups, groupp as id variables
8. REFERENCIAS BIBLIOGRÁFICAS
Beasley T.M. & Zumbo B.D. (2009). Aligned Rank Tests for Interactions in Split- Plot Designs: Distributional Assumptions and Stochastic Heterogeneity. Journal of Modern Applied Statistical Methods, 8(1), 16-50.
Conover W. & Iman R.L. (1981). Rank transformations as a bridge between parametric and nonparametric statistics. . American Statistician, 35. (3), 124-129.
De Neve J. & Thas O. (2017). A Mann–Whitney type effect measure of interaction for factorial designs. Communications in Statistics - Theory and Methods, 46(22), 11243-11260.
Elkin L.A, Kay. M., Higgins. J.J. & Wobbrock, J.O. (2021). An Aligned Rank Transform Procedure for Multifactor Contrast Tests. Stat.ME Cornell University. https://arxiv.org/abs/2102.11824
Feys J. (2016). Nonparametric Tests for the Interaction in Two-way Factorial Designs Using R. The R Journal, 8(1), 367-378.
Friedrich S., Konietschke F. & Pauly M. (2017). GFD: An R Package for the Analysis of General Factorial Designs. Journal of Statistical Software, 79(1), 1-18.
Funder D.C. & Ozer D.J. (2019). Evaluating Effect Size in Psychological Research: Sense and Nonsense. Advances in Methods and Practices in Psychological Science, 2 ,156-168.
Gao X. & Alvo M. (2005). A nonparametric test for interaction in two‐way layouts. Canadian Journal of Statistics, 33(4), 529-543.
Higgins J. & Tashtoush S. (1994). An aligned rank transform test for interaction. Nonlinear World, 1(2), 201-211
Hinkelmann K. & Kempthorne O. (2007). Designs with Repeated Measures. En Design and Analysis of Experiments, Introduction to Experimental Design(2da ed., pp. 278-291). Wiley & Sons.
Kay M., Elkin L.A. & Wobbrock J.O. (2021). Contrast tests with ART.
Kay M., Elkin, L.A., Higgins. J.J. & Wobbrock J. (2021). ARTool: Aligned Rank Transform. R package version 0.11.0.
Kutner M., Nachtsheim C., Neter J. & Li W. (2005). Randomized Complete Block Designs. En Applied Linear Statistical Models (5th ed., pp. 911-912). McGraw-Hill Irwin.
Lakens D. (2013). Calculating and reporting effect sizes to facilitate cumulative science: a practical primer for t-tests and ANOVAs. Frontiers in Psychology, 4, 863.
Luepsen, H. (2018). Comparison of nonparametric analysis of variance methods: A vote for van der Waerden. Journal Communications in Statistics - Simulation and Computation, 47(9), 2547-2576.
Mansouri H. (1998). Multifactor analysis of variance based on the aligned rank transform technique. Computational Statistics & Data Analysis, 29(2), 177–189.
Melo O., López L. & Melo S. (2020). Diseño de experimentos. Métodos y aplicaciones. Segunda edición. Universidad Nacional de Colombia, Colombia 292-293 p.
Meza, A. (2021). Métodos alternativos ante la violación de supuestos en diseños de experimentos factoriales. Anales Científicos, 81(2), 318-335.
Norouzian R. & Plonsky L. (2018). Eta- and partial eta-squared in L2 research: A cautionary review and guide to more appropriate usage. Second Language Research, 34(2), 257-271. Olejnik S. & Algina J. (2003). Generalized eta and omega squared statistics: Measures of effect size for some common research designs. Psychological Methods, 8(4), 434-447.
Saste S.V., Sananse S.L. & Sonar C. (2016). On parametric and nonparametric analysis of two factor factorial experiment. International Journal of Applied Research, 2(7), 653–656.
Sullivan L.M. (2008). Repeated measures. Circulation, 117(9), 1238-1243.
Thongteeraparp A. (2019). The comparison of nonparametric statistical tests for interaction effects in factorial design. Decision Science Letters, 8(3), 309-316.
Trigo M.E. & Martínez R.J. (2016). Generalized eta squared for multiple comparisons on between-groups designs. Psicothema, 28(3), 340-345.
Uyanto S.S. (2019). Monte Carlo power comparison of seven most commonly used heteroscedasticity tests. Communications in Statistics - Simulation and Computation.
Wijekularathna D.K., Manage A.B. & Scariano S.M. (2019). Power analysis of several normality tests: A Monte Carlo simulation study. Communications in Statistics - Simulation and Computation.
Wobbrock J.O, Findlater L., Gergle D. & Higgins J.J. (2011). The Aligned Rank Transform for Nonparametric Factorial Analyses Using Only Anova Procedures. In Proceedings of the SIGCHI Conference on Human Factors in Computing Systems (CHI ’11). ACM, New York, NY, USA, 143–146.
Wobbrock J.O., Elkin L.A., Higgins J.J; Findlater L., Gergle D. & Kay M. (2021). Artool: Align-and-rank data for a nonparametric ANOVA.