Bastão de Asclépio & Distribuição Normal

Bastão de Asclépio & Distribuição Normal 

invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
options(warn=-1)
suppressMessages(library(bruceR, warn.conflicts=FALSE))
suppressMessages(library(car, warn.conflicts=FALSE))
suppressMessages(library(effectsize, warn.conflicts=FALSE))
suppressMessages(library(emmeans, warn.conflicts=FALSE))
suppressMessages(library(estimatr, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(gplots, warn.conflicts=FALSE))
suppressMessages(library(lmerTest, warn.conflicts=FALSE))
suppressMessages(library(multcomp, warn.conflicts=FALSE))
suppressMessages(library(MuMIn, warn.conflicts=FALSE))
suppressMessages(library(patchwork, warn.conflicts=FALSE))
suppressMessages(library(phia, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts=FALSE))
suppressMessages(library(pwr2, warn.conflicts=FALSE))
suppressMessages(library(RcmdrMisc, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
suppressMessages(library(twowaytests, warn.conflicts=FALSE))
suppressMessages(library(WebPower, warn.conflicts=FALSE))
source("summarySEwithin2.R")

Material

Objetivos

  • Discorrer sobre uma extensão da ANOVA que inclui dois ou mais fatores.
  • Descrever e indicar três delineamentos diferentes da ANOVA com dois fatores:
    • independente: entre participantes,
    • relacionada: intraparticipantes,
    • mista: um fator entre e outro fator intraparticipantes.
  • Proceder com a análise descritiva numérica e gráfica dos dados.
  • Formular e implementar os modelos de ANOVA adequados a cada situação.
  • Determinar as significâncias estatística e prática dos efeitos de interação e principais da ANOVA.
  • Executar uma análise de efeito simples.
  • Executar e interpretar os testes post hoc, quando adequados.

Introdução

ANOVA bifatorial permite testar efeitos principais e de interação numa única análise.

As variáveis envolvidas no teste estatístico são:

  • duas ou mais VI (variável independente) com duas ou mais categorias (fatores nominais dicotômicos ou politômicos): multifatorial
  • VD (variável dependente) intervalar com distribuição normal por condição: univariada

ANOVA multifatorial pode ser:

  • independente (fator entre participantes)
  • relacionada (fator intraparticipantes)
  • mista (fatores entre e intraparticipantes)

ANOVA multifatorial é uma extensão do ANOVA unifatorial.

As suposições de normalidade homocedasticidade são condições suficientes e a suposição de independência é necessária para ANOVA multifatorial.

ANOVA multifatorial em R

  1. Planejamento de ANOVA
  1. Bifatorial Independente de Fisher
  • pwr2::pwr.2way
  1. Multifatorial Independente de Fisher
  • WebPower::wp.kanova
  1. ANOVA multifatorial independente
  1. Fisher (homocedástico)
  • teste omnibus: lm, car::Anova, effectsize::eta_squared
    • teste de efeito principal simples: phia::testInteractions, emmeans::emmip
      • teste post hoc: emmeans::emmeans, multcomp::cld, emmeans::emmip
  1. Fisher-White (heterocedástico)
  • teste omnibus: lm, car::Anova(white.adjust="hc2", ...)), estimatr::lm_robust(se_type="HC2", ...)
    • teste post hoc: ?
  1. Parametric Bootstrap based Generalized Test (robusta a hetrocedasticidade e não normalidade)
  • teste omnibus: twowaytests::gpTwoWay
    • teste post hoc: twowaytests::paircompTwoWay
  1. ANOVA multifatorial relacionada ou para medidas repetidas
  • teste omnibus: lmerTest::lmer e car::Anova(test.statistic="F", ...), effectsize::eta_squared
    • teste de efeito principal simples: phia::testInteractions, emmeans::emmip
      • teste post hoc: emmeans::emmeans, multcomp::cld, emmeans::emmip
  1. ANOVA multifatorial mista
  • teste omnibus: lmerTest::lmer e car::Anova(test.statistic="F", ...), effectsize::eta_squared
    • teste de efeito principal simples: phia::testInteractions, emmeans::emmip
      • teste post hoc: emmeans::emmeans, multcomp::cld, emmeans::emmip

ANOVA independente bifatorial

Suponha que conduzimos um estudo com delineamento totalmente entre participantes, balanceado e não-randomizado com 48 participantes para investigar o efeito da cafeína na habilidade de dirigir sob o efeito do álcool.

Encontramos um artigo com este tema em Howland et al. (2010), disponível em https://pubmed.ncbi.nlm.nih.gov/21134017/; lamentavelmente, não conseguimos acesso aos dados publicados para que pudéssemos replicar a análise. Usaremos, aqui, dados hipotéticos que trarão conclusão similar.

“Vale a pena observar que a MANOVA é ainda um teste válido, mesmo com modestas violações na suposição de normalidade, particularmente quando os tamanhos dos grupos são iguais e existe um número razoável de participantes em cada grupo; por “razoável” entendemos que, em um delineamento completamente entre participantes, deve haver pelo menos 12 participantes por grupo[…]”

Dancey & Reidy, 2019, p. 472

O estudo analisa os efeitos de presença e ausência de álcool e de cafeína na quantidade de infrações (erros) cometidos ao dirigir um simulador de direção veicular durante 20 minutos.

Dada a antiga crença de que o café auxilia a mantermo-nos alertas, podemos fazer algumas previsões para esse experimento:

  • Altos níveis de álcool diminuem a capacidade de dirigir.
  • Altos níveis de cafeína podem melhorar a habilidade de dirigir devido ao efeito estimulante.
  • Um aumento de cafeína reduz a influência do álcool na habilidade de dirigir.

Se beber, tomar café para dirigir?

Os dados estão no arquivo CafEntre_AlcEntre.rds. Foram obtidos e adaptados de Dancey & Reidy (2006).

Dados <- data.frame(readxl::read_excel("CafEntre_AlcEntre.xlsx"))
Dados$UE <- factor(Dados$UE)
Dados$Cafeina <- factor(Dados$Cafeina, 
                        levels=c("SemC","ComC"))
Dados$Alcool <- factor(Dados$Alcool, 
                       levels=c("SemA","ComA"))
saveRDS(Dados, "CafEntre_AlcEntre.rds")
Dados <- readRDS("CafEntre_AlcEntre.rds")
str(Dados)
'data.frame':   48 obs. of  4 variables:
 $ UE      : Factor w/ 48 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ Cafeina : Factor w/ 2 levels "SemC","ComC": 1 1 1 1 1 1 1 1 1 1 ...
 $ Alcool  : Factor w/ 2 levels "SemA","ComA": 1 1 1 1 1 1 1 1 1 1 ...
 $ NumErros: num  4 9 10 8 6 11 2 11 11 10 ...
print.data.frame(Dados)
   UE Cafeina Alcool NumErros
1   1    SemC   SemA        4
2   2    SemC   SemA        9
3   3    SemC   SemA       10
4   4    SemC   SemA        8
5   5    SemC   SemA        6
6   6    SemC   SemA       11
7   7    SemC   SemA        2
8   8    SemC   SemA       11
9   9    SemC   SemA       11
10 10    SemC   SemA       10
11 11    SemC   SemA        3
12 12    SemC   SemA       10
13 13    SemC   ComA       28
14 14    SemC   ComA       22
15 15    SemC   ComA       21
16 16    SemC   ComA       27
17 17    SemC   ComA       21
18 18    SemC   ComA       20
19 19    SemC   ComA       19
20 20    SemC   ComA       16
21 21    SemC   ComA       25
22 22    SemC   ComA       17
23 23    SemC   ComA       19
24 24    SemC   ComA       20
25 25    ComC   SemA       15
26 26    ComC   SemA        8
27 27    ComC   SemA       17
28 28    ComC   SemA        0
29 29    ComC   SemA       15
30 30    ComC   SemA       11
31 31    ComC   SemA       11
32 32    ComC   SemA        6
33 33    ComC   SemA        0
34 34    ComC   SemA       15
35 35    ComC   SemA       17
36 36    ComC   SemA       15
37 37    ComC   ComA       10
38 38    ComC   ComA       11
39 39    ComC   ComA       27
40 40    ComC   ComA       15
41 41    ComC   ComA       27
42 42    ComC   ComA       10
43 43    ComC   ComA       21
44 44    ComC   ComA       15
45 45    ComC   ComA       19
46 46    ComC   ComA       21
47 47    ComC   ComA       15
48 48    ComC   ComA       15
print(summary(Dados[-1]))
 Cafeina    Alcool      NumErros    
 SemC:24   SemA:24   Min.   : 0.00  
 ComC:24   ComA:24   1st Qu.:10.00  
                     Median :15.00  
                     Mean   :14.29  
                     3rd Qu.:19.25  
                     Max.   :28.00  
print(xtabs(~Cafeina+Alcool, data=Dados))
       Alcool
Cafeina SemA ComA
   SemC   12   12
   ComC   12   12
print(ftable(Dados[-4]))
           Alcool SemA ComA
UE Cafeina                 
1  SemC              1    0
   ComC              0    0
2  SemC              1    0
   ComC              0    0
3  SemC              1    0
   ComC              0    0
4  SemC              1    0
   ComC              0    0
5  SemC              1    0
   ComC              0    0
6  SemC              1    0
   ComC              0    0
7  SemC              1    0
   ComC              0    0
8  SemC              1    0
   ComC              0    0
9  SemC              1    0
   ComC              0    0
10 SemC              1    0
   ComC              0    0
11 SemC              1    0
   ComC              0    0
12 SemC              1    0
   ComC              0    0
13 SemC              0    1
   ComC              0    0
14 SemC              0    1
   ComC              0    0
15 SemC              0    1
   ComC              0    0
16 SemC              0    1
   ComC              0    0
17 SemC              0    1
   ComC              0    0
18 SemC              0    1
   ComC              0    0
19 SemC              0    1
   ComC              0    0
20 SemC              0    1
   ComC              0    0
21 SemC              0    1
   ComC              0    0
22 SemC              0    1
   ComC              0    0
23 SemC              0    1
   ComC              0    0
24 SemC              0    1
   ComC              0    0
25 SemC              0    0
   ComC              1    0
26 SemC              0    0
   ComC              1    0
27 SemC              0    0
   ComC              1    0
28 SemC              0    0
   ComC              1    0
29 SemC              0    0
   ComC              1    0
30 SemC              0    0
   ComC              1    0
31 SemC              0    0
   ComC              1    0
32 SemC              0    0
   ComC              1    0
33 SemC              0    0
   ComC              1    0
34 SemC              0    0
   ComC              1    0
35 SemC              0    0
   ComC              1    0
36 SemC              0    0
   ComC              1    0
37 SemC              0    0
   ComC              0    1
38 SemC              0    0
   ComC              0    1
39 SemC              0    0
   ComC              0    1
40 SemC              0    0
   ComC              0    1
41 SemC              0    0
   ComC              0    1
42 SemC              0    0
   ComC              0    1
43 SemC              0    0
   ComC              0    1
44 SemC              0    0
   ComC              0    1
45 SemC              0    0
   ComC              0    1
46 SemC              0    0
   ComC              0    1
47 SemC              0    0
   ComC              0    1
48 SemC              0    0
   ComC              0    1

Este é um formato de arquivo wide, no qual cada linha contêm todos os dados de cada unidade experimental (UE). A identificação de cada unidade experimental, que no caso é um participante, é o elemento-chave que distingue os delineamentos. Neste caso, há 48 participantes e 48 linhas no arquivo.

Modelagem da ANOVA

Para abordar este tipo de situação, precisamos testar todos os efeitos principais e de interação numa única análise. Utilizaremos ANOVA bifatorial (two-way ANOVA), uma extensão da ANOVA unifatorial na qual:

  • a variável dependente (VD) é intervalar.
  • as variáveis independentes (VI) são fatores nominais com pelo menos dois níveis.

Dependendo do delineamento do estudo, três formas de executar a ANOVA são necessárias (independente, relacionada e mista).

Neste exemplo aplica-se a ANOVA bifatorial independente, na qual cada unidade observacional/experimental é alocada em apenas um tratamento. Cada tratamento é uma combinação dos níveis dos fatores. No exemplo acima, há quatro tratamentos:

  • com cafeína e com álcool,
  • com cafeína e sem álcool,
  • sem cafeína e com álcool,
  • sem cafeína e sem álcool.

Quando existem dois fatores há efeitos de interação e efeitos principais (de cada um dos fatores em separado). Portanto, o propósito da ANOVA é testar todos os efeitos dos fatores, que são fontes de variação para identificar quanto da variação total nos escores da VD pode ser atribuída a cada um destes efeitos.

Neste exemplo de fatores entre participantes, há quatro efeitos fixos:

  • interação entre Álcool e Cafeína,
  • principal de Álcool,
  • principal de Cafeína,
  • outros fatores não considerados no modelo, chamado de termo de erro.

Note que a forma mais correta é dizer que testamos efeitos dos fatores e não, como coloquialmente costuma-se dizer, que testam-se os fatores.

O caso deste exemplo é um modelo completo (full model) porque os quatro efeitos são considerados.

Em capítulos subsequentes, abordaremos os modelos parciais (custom models) nos quais podemos testar apenas o efeito de interação, os efeitos principais, parte dos efeitos de interação (quando há três ou mais fatores) ou parte dos efeitos principais de interesse do pesquisador. No modelo parcial, o efeito do erro sempre aparece, contendo tudo o que não é incorporado ao modelo.

Planejamento de ANOVA bifatorial independente de Fisher

Nas condições deste exemplo (serve apenas para ANOVA independente bifatorial), podemos saber o poder a priori usando a função pwr2::pwr.2way:

print(pwr2::pwr.2way(a=2, 
                     b=2, 
                     alpha=0.05,
                     size.A=24, 
                     size.B=24,
                     f.A=0.25, 
                     f.B=0.25))

     Balanced two-way analysis of variance power calculation 

              a = 2
              b = 2
            n.A = 24
            n.B = 24
      sig.level = 0.05
        power.A = 0.6787416
        power.B = 0.6787416
          power = 0.6787416

NOTE: power is the minimum power among two factors

ANOVA bifatorial independente de Fisher

Testes:

  • Omnibus: lm, car::Anova, effectsize::eta_squared
    • Efeito principal simples: phia::testInteractions, emmeans::emmip
      • Post hoc: emmeans::emmeans, multcomp::cld, emmeans::emmip

O código completo está implementado em demo_ANOVA2f_Indep_Fisher_AlcCaf.R.

alfa <- 0.05

Dados <- readRDS("CafEntre_AlcEntre.rds")
cat("\nDescriptive Statistics")

Descriptive Statistics
cat("\nBy Alcool")

By Alcool
print(psych::describeBy(NumErros~Alcool,
                        digits=2,
                        mat=2,
                        data=Dados))
          item group1 vars  n  mean   sd median trimmed  mad min max range
NumErros1    1   SemA    1 24  9.38 5.04   10.0    9.55 5.93   0  17    17
NumErros2    2   ComA    1 24 19.21 5.27   19.5   19.30 5.93  10  28    18
           skew kurtosis   se
NumErros1 -0.30    -0.93 1.03
NumErros2 -0.03    -0.91 1.08
cat("\nBy Cafeina")

By Cafeina
print(psych::describeBy(NumErros~Cafeina,
                        digits=2,
                        mat=2,
                        data=Dados))
          item group1 vars  n  mean   sd median trimmed  mad min max range
NumErros1    1   SemC    1 24 14.58 7.63   13.5    14.5 8.90   2  28    26
NumErros2    2   ComC    1 24 14.00 6.72   15.0    14.1 5.93   0  27    27
           skew kurtosis   se
NumErros1  0.06    -1.25 1.56
NumErros2 -0.16    -0.02 1.37
cat("\nBy Alcool and Cafeina")

By Alcool and Cafeina
print(psych::describeBy(NumErros~Alcool+Cafeina,
                        digits=2,
                        mat=2,
                        data=Dados))
          item group1 group2 vars  n  mean   sd median trimmed  mad min max
NumErros1    1   SemA   SemC    1 12  7.92 3.32    9.5     8.2 2.22   2  11
NumErros2    2   ComA   SemC    1 12 21.25 3.72   20.5    21.1 2.22  16  28
NumErros3    3   SemA   ComC    1 12 10.83 6.12   13.0    11.3 4.45   0  17
NumErros4    4   ComA   ComC    1 12 17.17 5.92   15.0    16.9 6.67  10  27
          range  skew kurtosis   se
NumErros1     9 -0.63    -1.35 0.96
NumErros2    12  0.49    -1.04 1.07
NumErros3    17 -0.70    -1.08 1.77
NumErros4    17  0.41    -1.24 1.71
boxplot(NumErros~Cafeina+Alcool, 
        data=Dados,
        ylab="NumErros")

alfaBonf <- alfa/(length(unique(Dados$Cafeina))*
                  length(unique(Dados$Alcool)))
gplots::plotmeans(NumErros~interaction(Cafeina, Alcool),
                  data=Dados,
                  p=1-alfaBonf,
                  main=paste("CI",round((1-alfaBonf)*100,4),"%",sep=""),
                  barcol="black",
                  connect=FALSE,)

cat("\nAssumptions")

Assumptions
cat("\nNormality")

Normality
norm.test <- by(data=Dados$NumErros, 
                INDICES=list("Alcool" = Dados$Alcool, 
                     "Cafeina" = Dados$Cafeina), 
                FUN=shapiro.test)
print(norm.test)
Alcool: SemA
Cafeina: SemC

    Shapiro-Wilk normality test

data:  dd[x, ]
W = 0.84029, p-value = 0.02792

------------------------------------------------------------ 
Alcool: ComA
Cafeina: SemC

    Shapiro-Wilk normality test

data:  dd[x, ]
W = 0.93036, p-value = 0.384

------------------------------------------------------------ 
Alcool: SemA
Cafeina: ComC

    Shapiro-Wilk normality test

data:  dd[x, ]
W = 0.84711, p-value = 0.03382

------------------------------------------------------------ 
Alcool: ComA
Cafeina: ComC

    Shapiro-Wilk normality test

data:  dd[x, ]
W = 0.89979, p-value = 0.1576
cat("\nHomoscedasticity")

Homoscedasticity
homoc.test <- car::leveneTest(NumErros~Alcool*Cafeina,
                              data=Dados)
print(homoc.test)
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  3  1.3621 0.2668
      44               
# ANOVA bifatorial independente de Fisher
cat("\nFisher's two-factor between-subjects ANOVA")

Fisher's two-factor between-subjects ANOVA
cat("\nGLM")

GLM
modelo <- lm(NumErros ~ Alcool*Cafeina,  
             data=Dados)

cat("\nANOVA table")

ANOVA table
print(anv <- car::Anova(modelo))
Anova Table (Type II tests)

Response: NumErros
                Sum Sq Df F value   Pr(>F)    
Alcool         1160.33  1 47.6924 1.57e-08 ***
Cafeina           4.08  1  0.1678  0.68403    
Alcool:Cafeina  147.00  1  6.0420  0.01798 *  
Residuals      1070.50 44                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(reg <- summary(modelo))

Call:
lm(formula = NumErros ~ Alcool * Cafeina, data = Dados)

Residuals:
    Min      1Q  Median      3Q     Max 
-10.833  -2.396   0.125   3.771   9.833 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)               7.917      1.424   5.560 1.49e-06 ***
AlcoolComA               13.333      2.014   6.621 4.11e-08 ***
CafeinaComC               2.917      2.014   1.448    0.155    
AlcoolComA:CafeinaComC   -7.000      2.848  -2.458    0.018 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.932 on 44 degrees of freedom
Multiple R-squared:  0.5506,    Adjusted R-squared:  0.5199 
F-statistic: 17.97 on 3 and 44 DF,  p-value: 9.277e-08
# cat("\nOmnibus F test by regression")
print(reg <- summary(modelo))

Call:
lm(formula = NumErros ~ Alcool * Cafeina, data = Dados)

Residuals:
    Min      1Q  Median      3Q     Max 
-10.833  -2.396   0.125   3.771   9.833 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)               7.917      1.424   5.560 1.49e-06 ***
AlcoolComA               13.333      2.014   6.621 4.11e-08 ***
CafeinaComC               2.917      2.014   1.448    0.155    
AlcoolComA:CafeinaComC   -7.000      2.848  -2.458    0.018 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.932 on 44 degrees of freedom
Multiple R-squared:  0.5506,    Adjusted R-squared:  0.5199 
F-statistic: 17.97 on 3 and 44 DF,  p-value: 9.277e-08
fobs <- reg$fstatistic[1]
dfn <- reg$fstatistic[2]
dfd <- reg$fstatistic[3]
fc <- qf(1-alfa, dfn, dfd)
p <- 1-pf(fobs, dfn, dfd)
if (p < 1e-4)
{
  p <- sprintf("%.2e",p)
} else
{
  p <- sprintf("%.4f",p)
}
f <- seq(0,1.4*max(fc,fobs),length.out=300)
densf <- df(f, dfn, dfd)
plot(f, densf, 
     main="Fisher Two-way ANOVA: Omnibus test",
     xlab="F", ylab="Density", 
     lwd=2, type="l")
abline(v=fc, lty=3)
abline(v=fobs, lty=4)
legend("topright",
       c("H0: there is not model",
         paste("F(",dfn,",",dfd,",",1-alfa,") = ",round(fc,3),sep=""),
         paste("F(",dfn,",",dfd,") = ",round(fobs,3),"\n",
               "p = ",p,sep="")
       ),
       lwd=c(2,1,1), lty=c(1,3,4),
       cex=0.8)

cat("\nEffect size analysis\n")

Effect size analysis
cat("R^2 = eta^2 omnibus =", round(reg$r.squared, 4))
R^2 = eta^2 omnibus = 0.5506
print(effectsize::interpret_eta_squared(reg$r.squared,
                                        rules="cohen1992"))
[1] "large"
(Rules: cohen1992)
eta2 <- effectsize::eta_squared(anv,
                                partial=FALSE,
                                generalized=FALSE,
                                ci=1-alfa/3,
                                alternative="two.sided",
                                verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=4)
# Effect Size for ANOVA (Type II)

Parameter      |   Eta2 |      98.3333% CI |  interpret
-------------------------------------------------------
Alcool         | 0.4871 | [0.2219, 0.6653] |      large
Cafeina        | 0.0017 | [0.0000, 0.0035] | very small
Alcool:Cafeina | 0.0617 | [0.0000, 0.2789] |     medium
cat("\nGráfico de perfil de médias")

Gráfico de perfil de médias
g1 <- emmeans::emmip(object=modelo,
                     formula=~Cafeina ~Alcool,
                     level=1-alfa,
                     adjust="holm",
                     CIs=TRUE,
                     dodge=.3) + 
      ggplot2::theme_bw() 
g2 <- emmeans::emmip(object=modelo,
                     formula=~Alcool ~Cafeina,
                     level=1-alfa,
                     adjust="holm",
                     CIs=TRUE,
                     dodge=.3) + 
      ggplot2::theme_bw()
g3 <- emmeans::emmip(modelo, 
                     ~Cafeina,
                     level=1-alfa,
                     adjust="holm",
                     CIs=TRUE)
NOTE: Results may be misleading due to involvement in interactions
g4 <- emmeans::emmip(modelo, 
                     ~Alcool, 
                     level=1-alfa,
                     adjust="holm",
                     CIs=TRUE)
NOTE: Results may be misleading due to involvement in interactions
p <- (g1 | g2) / (g3 | g4)   # 2 x 2
print(p)       

cat("\nPost hoc test se efeito de interação Alcool:Cafeina não é significante")

Post hoc test se efeito de interação Alcool:Cafeina não é significante
cat("\nAlcool")

Alcool
EMM.A <- emmeans::emmeans(modelo, 
                          specs=pairwise~"Alcool", 
                          level=1-alfa,
                          adjust="holm")
NOTE: Results may be misleading due to involvement in interactions
print(summary(EMM.A$emmeans))
 Alcool emmean   SE df lower.CL upper.CL
 SemA     9.38 1.01 44     7.35     11.4
 ComA    19.21 1.01 44    17.18     21.2

Results are averaged over the levels of: Cafeina 
Confidence level used: 0.95 
print(summary(EMM.A$contrasts, 
              infer=TRUE))
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemA - ComA    -9.83 1.42 44    -12.7    -6.96  -6.906  <.0001

Results are averaged over the levels of: Cafeina 
Confidence level used: 0.95 
print(plot(EMM.A$emmeans, 
           colors="black"))

print(plot(EMM.A$contrasts, 
           colors="black"))

print(multcomp::cld(object=EMM.A$emmeans,
                    level=1-alfa,
                    alpha=alfa,
                    adjust="holm",
                    Letters=letters))
 Alcool emmean   SE df lower.CL upper.CL .group
 SemA     9.38 1.01 44     7.04     11.7  a    
 ComA    19.21 1.01 44    16.87     21.5   b   

Results are averaged over the levels of: Cafeina 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 2 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
cat("\nCafeina")

Cafeina
EMM.C <- emmeans::emmeans(modelo, 
                          specs=pairwise~"Cafeina", 
                          level=1-alfa,
                          adjust="holm")
NOTE: Results may be misleading due to involvement in interactions
print(summary(EMM.C$emmeans))
 Cafeina emmean   SE df lower.CL upper.CL
 SemC      14.6 1.01 44     12.6     16.6
 ComC      14.0 1.01 44     12.0     16.0

Results are averaged over the levels of: Alcool 
Confidence level used: 0.95 
print(summary(EMM.C$contrasts, 
              infer=TRUE))
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemC - ComC    0.583 1.42 44    -2.29     3.45   0.410  0.6840

Results are averaged over the levels of: Alcool 
Confidence level used: 0.95 
print(plot(EMM.C$emmeans, 
           colors="black"))

print(plot(EMM.C$contrasts, 
           colors="black"))

print(multcomp::cld(object=EMM.C$emmeans,
                    level=1-alfa,
                    alpha=alfa,
                    adjust="holm",
                    Letters=letters))
 Cafeina emmean   SE df lower.CL upper.CL .group
 ComC      14.0 1.01 44     11.7     16.3  a    
 SemC      14.6 1.01 44     12.2     16.9  a    

Results are averaged over the levels of: Alcool 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 2 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
cat("\nPost hoc teste se efeito de interação Alcool:Cafeina é significante")

Post hoc teste se efeito de interação Alcool:Cafeina é significante
cat("\nPost hoc test: interaction Cafeina:Alcool")

Post hoc test: interaction Cafeina:Alcool
plot(emmeans::emmip(modelo, 
                    ~Cafeina:Alcool,
                    level=1-alfa,
                    adjust="holm",
                    CIs=TRUE))

EMM.AC <- emmeans::emmeans(modelo, 
                           specs=pairwise~"Alcool:Cafeina",
                           level=1-alfa,
                           adjust="holm")
print(summary(EMM.AC$emmeans))
 Alcool Cafeina emmean   SE df lower.CL upper.CL
 SemA   SemC      7.92 1.42 44     5.05     10.8
 ComA   SemC     21.25 1.42 44    18.38     24.1
 SemA   ComC     10.83 1.42 44     7.96     13.7
 ComA   ComC     17.17 1.42 44    14.30     20.0

Confidence level used: 0.95 
print(summary(EMM.AC$contrasts, infer=TRUE))
 contrast              estimate   SE df lower.CL upper.CL t.ratio p.value
 SemA SemC - ComA SemC   -13.33 2.01 44   -18.90    -7.77  -6.621  <.0001
 SemA SemC - SemA ComC    -2.92 2.01 44    -8.48     2.65  -1.448  0.1546
 SemA SemC - ComA ComC    -9.25 2.01 44   -14.81    -3.69  -4.594  0.0001
 ComA SemC - SemA ComC    10.42 2.01 44     4.85    15.98   5.173  <.0001
 ComA SemC - ComA ComC     4.08 2.01 44    -1.48     9.65   2.028  0.0973
 SemA ComC - ComA ComC    -6.33 2.01 44   -11.90    -0.77  -3.145  0.0089

Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 6 estimates 
P value adjustment: holm method for 6 tests 
print(plot(EMM.AC$emmeans, 
           colors="black"))

print(plot(EMM.AC$contrasts, 
           colors="black"))

print(multcomp::cld(object=EMM.AC$emmeans,
                    level=1-alfa,
                    alpha=alfa,
                    adjust="holm",
                    Letters=letters))
 Alcool Cafeina emmean   SE df lower.CL upper.CL .group
 SemA   SemC      7.92 1.42 44     4.21     11.6  a    
 SemA   ComC     10.83 1.42 44     7.12     14.5  a    
 ComA   ComC     17.17 1.42 44    13.46     20.9   b   
 ComA   SemC     21.25 1.42 44    17.54     25.0   b   

Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 4 estimates 
P value adjustment: holm method for 6 tests 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
cat("\nSimple main effect test: Alcool effect by Cafeina levels\n")

Simple main effect test: Alcool effect by Cafeina levels
print(phia::testInteractions(modelo, 
                             across="Alcool",
                             fixed="Cafeina",
                             adjustment="holm")) 
F Test: 
P-value adjustment method: holm
             Value     SE Df Sum of Sq      F    Pr(>F)    
SemC      -13.3333 2.0137  1   1066.67 43.842 8.225e-08 ***
ComC       -6.3333 2.0137  1    240.67  9.892  0.002973 ** 
Residuals                 44   1070.50                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(emmeans::emmip(modelo,
                    ~Alcool|Cafeina,
                    level=1-alfa,
                    adjust="holm",
                    CIs=TRUE)) 

cat("\nPost hoc test of the simple main effect of Alcool")

Post hoc test of the simple main effect of Alcool
model_means <- emmeans::emmeans(object=modelo,
                                specs=pairwise~Alcool|Cafeina:Alcool,
                                level=1-alfa,
                                adjust="holm")
print(summary(model_means$emmeans, infer=TRUE))
Cafeina = SemC:
 Alcool emmean   SE df lower.CL upper.CL t.ratio p.value
 SemA     7.92 1.42 44     5.05     10.8   5.560  <.0001
 ComA    21.25 1.42 44    18.38     24.1  14.924  <.0001

Cafeina = ComC:
 Alcool emmean   SE df lower.CL upper.CL t.ratio p.value
 SemA    10.83 1.42 44     7.96     13.7   7.608  <.0001
 ComA    17.17 1.42 44    14.30     20.0  12.056  <.0001

Confidence level used: 0.95 
print(summary(model_means$contrasts, infer=TRUE))
Cafeina = SemC:
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemA - ComA   -13.33 2.01 44    -17.4    -9.28  -6.621  <.0001

Cafeina = ComC:
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemA - ComA    -6.33 2.01 44    -10.4    -2.28  -3.145  0.0030

Confidence level used: 0.95 
plot(model_means$emmeans,col="black")

plot(model_means$contrasts,col="black")

print(multcomp::cld(object=model_means$emmeans,
                    adjust="holm",
                    Letters=letters,
                    alpha=alfa))
Cafeina = SemC:
 Alcool emmean   SE df lower.CL upper.CL .group
 SemA     7.92 1.42 44     4.61     11.2  a    
 ComA    21.25 1.42 44    17.95     24.6   b   

Cafeina = ComC:
 Alcool emmean   SE df lower.CL upper.CL .group
 SemA    10.83 1.42 44     7.53     14.1  a    
 ComA    17.17 1.42 44    13.86     20.5   b   

Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 2 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
cat("\nSimple main effect test: Cafeina effect by Alcool levels\n")

Simple main effect test: Cafeina effect by Alcool levels
print(phia::testInteractions(modelo, 
                             across="Cafeina",
                             fixed="Alcool",
                             adjustment="holm"))
F Test: 
P-value adjustment method: holm
            Value     SE Df Sum of Sq      F  Pr(>F)  
SemA      -2.9167 2.0137  1     51.04 2.0979 0.15459  
ComA       4.0833 2.0137  1    100.04 4.1119 0.09733 .
Residuals                44   1070.50                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(emmeans::emmip(modelo, 
                    ~Cafeina|Alcool,
                    level=1-alfa,
                    adjust="holm",
                    CIs=TRUE))

cat("\nPost hoc test of the simple main effect of Cafeina")

Post hoc test of the simple main effect of Cafeina
model_means <- emmeans::emmeans(object=modelo,
                                specs=pairwise~Cafeina|Cafeina:Alcool,
                                level=1-alfa,
                                adjust="holm")
print(summary(model_means$emmeans, infer=TRUE))
Alcool = SemA:
 Cafeina emmean   SE df lower.CL upper.CL t.ratio p.value
 SemC      7.92 1.42 44     5.05     10.8   5.560  <.0001
 ComC     10.83 1.42 44     7.96     13.7   7.608  <.0001

Alcool = ComA:
 Cafeina emmean   SE df lower.CL upper.CL t.ratio p.value
 SemC     21.25 1.42 44    18.38     24.1  14.924  <.0001
 ComC     17.17 1.42 44    14.30     20.0  12.056  <.0001

Confidence level used: 0.95 
print(summary(model_means$contrasts, infer=TRUE))
Alcool = SemA:
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemC - ComC    -2.92 2.01 44   -6.975     1.14  -1.448  0.1546

Alcool = ComA:
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemC - ComC     4.08 2.01 44    0.025     8.14   2.028  0.0487

Confidence level used: 0.95 
plot(model_means$emmeans,col="black")

plot(model_means$contrasts,col="black")

print(multcomp::cld(object=model_means$emmeans,
                    level=1-alfa,
                    alpha=alfa,
                    adjust="holm",
                    Letters=letters))
Alcool = SemA:
 Cafeina emmean   SE df lower.CL upper.CL .group
 SemC      7.92 1.42 44     4.61     11.2  a    
 ComC     10.83 1.42 44     7.53     14.1  a    

Alcool = ComA:
 Cafeina emmean   SE df lower.CL upper.CL .group
 ComC     17.17 1.42 44    13.86     20.5  a    
 SemC     21.25 1.42 44    17.95     24.6   b   

Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 2 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same.

Adiante veremos ANOVA relacionada (com dois fatores intraparticipantes) e mista (com um fator intra e outro entre participantes):

  • As suposições de normalidade e homocedasticidade não são testáveis para a ANOVA relacionada.
  • Elas podem ser testadas apenas para o fator entre participantes na ANOVA mista.

A função lm implementa o GLM Neste caso é um General Linear Model (GLM) porque testa efeitos fixos de fatores, tratando NumErros em função (~) de Alcool interagindo (*) com Cafeina. NumErros é a VD, intervalar. Alcool e Cafeina são os fatores nominais.

Esta forma de escrever a interação é equivalente à sintaxe na qual os efeitos são explicitamente exibidos (exceto o termo de erro, que é inerente ao modelo):

modelo <- lm(NumErros ~ Alcool+Cafeina+Alcool:Cafeina, data=Dados)
modelo <- lm(NumErros ~ Alcool*Cafeina, data=Dados)
modelo <- lm(NumErros ~ (Alcool+Cafeina)^2, data=Dados)

Portanto, NumErros é tratado em função dos efeitos principais do Alcool e da Cafeina, e também em função do efeito de interação Alcool:Cafeina.

Este é o teste omnibus, e utilizamos apenas o valor p associado à estatística \(F\) da saída da função lm (neste caso, rejeita-se a hipótese nula de ausência de modelo).

O valor de \(R^2\) que aparece nesta saída corresponde à estimativa pontual do tamanho de efeito omnibus a posteriori.

Observe que obtivemos esta saída estabelecendo o modelo com lm e exibindo no formato de regressão com summary. Vamos denominar este procedimento de “Regressão da ANOVA”.

A regressão funciona como um teste omnibus, cuja hipótese nula é a não existência de um modelo.

Para o nível de significância adotado (\(\alpha=0.05\)), o valor p rejeita a hipótese nula (a hipótese alternativa é a violação de qualquer igualdade), portanto existe modelo, significando que a estrutura proposta detecta a existência de associação entre os efeitos fixos dos fatores e a VD.

\(R^2\) (Multiple R-squared) é a porcentagem da variância da VD explicada pelos efeitos do modelo e, portanto, é um tamanho de efeito que equivale ao \(\eta^2\) global.

Tabela ANOVA da ANOVA Bifatorial

Como existe modelo, podemos verificar os efeitos na saída ANOVA, para explicitamente podermos definir a rejeição ou não rejeição de cada uma das hipóteses. Utilizamos o resultado de lm no formato de ANOVA exibido com a função car::Anova, que denominaremos, aqui, de “ANOVA da ANOVA”.

Há parâmetros para a função car::Anova. O parâmetro type=2 é o default. Estes tipos são maneiras diferentes de particionar a variância total da VD entre os efeitos elencados no modelo.

Existe também o parâmetro white.adjust=TRUE serve para executar método mais robusto, com correção para heterocedasticidade (a mesma que utilizamos anteriormente com ANOVA unifatorial). Não o utilizamos aqui por coerência, pois os testes post hoc que são realizados na ANOVA bifatorial independente de Fisher não usam o mesmo tipo de correção.

Caso prefira, existe a função anova que está disponível no RStudio sem que se chame seu pacote stats, mais limitada que a car::Anova (por não oferecer correção para heterocedasticidade) e que também utiliza tipo 2 por default (portanto, neste exemplo, mostraria resultados idênticos). Cuidado com a função anova, pois os valores p podem variar em função da ordem dos fatores na fórmula do modelo. A função car::Anova faz os valores p serem invariantes com relação à ordem dos fatores na fórmula do modelo.

Efeito de interação

O teste do modelo como um todo implica em testar implicitamente todas as hipóteses nulas envolvidas na ANOVA:

\[ \begin{cases} H_{0}^\text{Alcool}: &\mu_{\text{SemA}} = \mu_{\text{ComA}}\\ H_{0}^\text{Cafeina}: &\mu_{\text{SemC}} = \mu_{\text{ComC}}\\ H_{0}^\text{Alcool:Cafeina}: &\text{ausência de efeito de interação populacional} \end{cases} \]

A primeira hipótese a ser verificada é o efeito de interação dos fatores em Alcool:Cafeina. O seu valor p é 0.01798.

Se o nível de significância adotado é 0.01, a hipótese nula do efeito de interação não é rejeitada. Nesta situação, os valores p dos efeitos principais dos fatore Alcool e Cafeína pode ser comparados com o nível de significância adotado. Se o fator tem três ou mais níveis, podem ser realizados os testes post hoc.

Se o nível de significância adotado é 0.05, rejeita-se a hipótese nula de ausência de efeito. Na presença de efeito de interação entre Cafeína e Álcool, não podemos considerar a significância estatística dos efeitos isolados que aparecem nas linhas Alcool e Cafeina nesta saída.

Podemos, apenas, considerar seus tamanhos de efeito, os \(\eta^2\) parciais, que são as porcentagens da variância da VD que são explicadas pelos respectivos efeitos. Cada \(\eta\) parcial é uma correlação de Pearson parcial absoluta implícita entre o respectivo efeito e a VD. Como há interação de efeitos, existem redundâncias que fazem com que a soma dos \(\eta^2\) parciais não seja necessariamente igual ao \(\eta^2\) global.

Efeito principal simples

Como houve efeito de interação entre Cafeína e Álcool, precisamos estudar os efeitos principais simples, que é a verificação de cada um efeitos de um fator em todos os níveis do outro fator.

Neste exemplo, em que ocorreu interação usando o nível de significância de 0.05, a verificação dos efeitos principais simples de Álcool e Cafeína é feito com funções R que fornecem os valores p já ajustados, portanto podemos compará-los diretamente com \(\alpha=0.05\). Neste exemplo, verificamos que o efeito do Álcool é significante nos dois níveis de cafeína, portanto podemos dizer que o efeito principal do Álcool é genuíno. Ao contrário, o efeito da Cafeína não foi significante nas condições com ou sem Álcool, portanto não se trata de um efeito principal genuíno. Isto equivale a dizer Cafeína não modifica a VD de acordo com o nível de Álcool.

O gráfico foi produzido pela função nativa plot que interpreta o objeto gerado por phia::interactionMeans. Os intervalos de confiança 95% neste gráfico são pós-modelo (observe que estão equalizados). Na diagonal principal vemos os efeitos simples de álcool e cafeína. Na diagonal secundária observamos as interações; já definimos que, para \(\alpha=0.05\), supomos interação de álcool e cafeína, o que é visto graficamente pelos segmentos de reta não paralelos.

Os testes de efeito principal simples com a função phia::testInteractions e emmeans::emmeans.

Portanto, a pergunta que foi feita inicialmente: Se beber, tomar café para dirigir? é respondida, de acordo com esta análise:

Segundo este estudo, não há evidência amostral de que a ingestão de Cafeína contrabalança o efeito do Álcool para que alguém possa dirigir alcoolizado.

Em geral, os trabalhos publicados não têm mais do que três fatores. Imagine um estudo com Álcool, Cafeína e Glicose, cada um com dois níveis. O total de efeitos é \(2×2×2=8\):

  • efeitos principais:
    • Álcool (com ou sem)
    • Cafeína (com ou sem)
    • Glicose (com ou sem)
  • efeitos de interação de ordem 2:
    • Álcool:Cafeína
    • Álcool:Glicose
    • Cafeína:Glicose
  • efeito de interação de ordem 3:
    • Álcool:Cafeína:Glicose
  • efeito do termo de erro.

Genericamente, uma ANOVA \(a_1×a_2×\cdots×a_k\) tem \(2^k\) efeitos, sendo que \(a_i\) é o número de níveis do fator \(i\). Assim, \(a_1×a_2×\cdots×a_k\) é número de tratamentos/condições experimentais, divididos em:

  • \(k\) efeitos principais, i.e., \(C(k,1)\) ou choose(k,1)
  • \(k\) efeitos de interação de ordem 2, i.e., \(C(k,2)\) ou choose(k,2)

…

  • \(k\) efeitos de interação de ordem \(k-1\), i.e., \(C(k,k-1)\) ou choose(k,k-1)
  • \(1\) efeito de interação de ordem \(k\), i.e., \(C(k,k)\) ou choose(k,k)
  • \(1\) efeito de erro, i.e., \(C(k,0)\) ou choose(k,0)

Como regra heurística em ANOVA com mais de três fatores analisam-se os efeitos de interação até ordem 3, sendo que os efeitos de interação podem ser escolhidos para análise.

Como foi visto, os testes de efeitos de ordem inferior só podem ser considerados se não há efeitos de ordem superior estatisticamente significantes mas podem ser considerados diretamente se o efeito de interação não é estatisticamente significante.

> Howell & Lacroix, 2012

ANOVA bifatorial independente de Fisher-White

Testes:

  • teste omnibus: lm, car::Anova(white.adjust="hc2", ...)), estimatr::lm_robust(se_type="HC2", ...)

ANOVA bifatorial independente de Fisher-White é um teste robusto à heterocedasticidade. Não há testes de efeito principal simples e post hoc em R.

O código completo está implementado em demo_ANOVA2f_Indep_Fisher-White_CafAlc.R.


Fisher-White two-factor between-subjects ANOVA
Statistical analysis: omnibus test
ANOVA
Coefficient covariances computed by hccm()
Analysis of Deviance Table (Type II tests)

Response: NumErros
               Df      F    Pr(>F)    
Alcool          1 86.504 5.951e-12 ***
Cafeina         1  0.158   0.69297    
Alcool:Cafeina  1  6.042   0.01798 *  
Residuals      44                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Call:
estimatr::lm_robust(formula = NumErros ~ Alcool * Cafeina, data = Dados, 
    se_type = "HC2")

Standard error type:  HC2 

Coefficients:
                       Estimate Std. Error t value  Pr(>|t|) CI Lower CI Upper
(Intercept)               7.917     0.9571   8.272 1.649e-10    5.988    9.846
AlcoolComA               13.333     1.4386   9.269 6.591e-12   10.434   16.233
CafeinaComC               2.917     2.0087   1.452 1.536e-01   -1.132    6.965
AlcoolComA:CafeinaComC   -7.000     2.8478  -2.458 1.798e-02  -12.739   -1.261
                       DF
(Intercept)            44
AlcoolComA             44
CafeinaComC            44
AlcoolComA:CafeinaComC 44

Multiple R-squared:  0.5506 ,   Adjusted R-squared:  0.5199 
F-statistic: 30.86 on 3 and 44 DF,  p-value: 6.729e-11
R^2 = eta^2 omnibus = 0.550572[1] "substantial"
(Rules: cohen1988)
# Effect Size for ANOVA (Type I)

Parameter      |     Eta2 |          98.3333% CI
------------------------------------------------
Alcool         | 0.487143 | [0.221913, 0.665274]
Cafeina        | 0.001714 | [0.000000, 0.003486]
Alcool:Cafeina | 0.061715 | [0.000000, 0.278859]

ANOVA bifatorial independente por Parametric Bootstrap based Generalized Test

Testes:

  • teste omnibus: twowaytests::gpTwoWay
    • testes post hoc: twowaytests::paircompTwoWay

ANOVA bifatorial independente por Parametric Bootstrap based Generalized Test é um teste robusto à normalidade e heterocedasticadade. Não há teste post hoc em R.

O código completo está implementado em demo_ANOVA2f_Indep_PBoot_CafAlc.R.

alfa <- 0.05

Dados <- readRDS("CafEntre_AlcEntre.rds")

# ANOVA multifatorial independente por _Parametric Bootstrap based Generalized Test_:   
omni_Caf <- twowaytests::gpTwoWay(NumErros ~ Cafeina*Alcool,
                                  alpha = alfa,
                                  data=Dados)

  Generalized p-value by PB method (alpha = 0.05)
--------------------------------------------------
         Factor P.value     Result
        Cafeina  0.0681 Not reject
         Alcool  0.0000     Reject
 Cafeina:Alcool  0.0213     Reject
--------------------------------------------------
# Teste de efeito simples de Cafeina
twowaytests::paircompTwoWay(omni_Caf,
                            alpha = alfa,
                            adjust.method = "holm")

  Holm correction (alpha = 0.05) for subgroups of each Cafeina level 
-------------------------------------------------------------------------------- 
  Cafeina Alcool(a) Alcool(b)      P.value   No difference
1    SemC      SemA      ComA 5.316791e-09          Reject
2    ComC      SemA      ComA 1.721035e-02          Reject
-------------------------------------------------------------------------------- 
omni_Alc <- twowaytests::gpTwoWay(NumErros ~ Alcool*Cafeina, 
                                  alpha = alfa,
                                  data=Dados)

  Generalized p-value by PB method (alpha = 0.05)
--------------------------------------------------
         Factor P.value     Result
         Alcool  0.0000     Reject
        Cafeina  0.0684 Not reject
 Alcool:Cafeina  0.0205     Reject
--------------------------------------------------
# Teste de efeito simples de Alcool
twowaytests::paircompTwoWay(omni_Alc,
                            alpha = alfa,
                            adjust.method = "holm")

  Holm correction (alpha = 0.05) for subgroups of each Alcool level 
-------------------------------------------------------------------------------- 
  Alcool Cafeina(a) Cafeina(b)    P.value   No difference
1   SemA       SemC       ComC 0.16475000      Not reject
2   ComA       SemC       ComC 0.05778413      Not reject
--------------------------------------------------------------------------------

ANOVA bifatorial relacionada

Testes:

  • Omnibus: lmerTest::lmer e car::Anova(test.statistic="F", ...), effectsize::eta_squared
    • Efeito principal simples: phia::testInteractions, emmeans::emmip
      • Post hoc: emmeans::emmeans, multcomp::cld, emmeans::emmip

Usando os mesmos valores do exemplo anterior (para comparação dos resultados), imagine que tenhamos apenas 12 participantes, cada um deles sendo submetido às quatro condições experimentais.

Neste caso, é usado General Linear Mixed Model (GLMM) por meio da função lmerTest::lmer.

“Vale a pena observar que a MANOVA é ainda um teste válido, mesmo com modestas violações na suposição de normalidade, particularmente quando os tamanhos dos grupos são iguais e existe um número razoável de participantes em cada grupo; por “razoável” entendemos que [para um delineamento] completamente intraparticipantes, [deve haver] pelo menos 22 participantes ao todo.”

Dancey & Reidy, 2019, p. 472

Neste exemplo aplica-se a ANOVA relacionada bifatorial. Os tratamentos são os mesmos, com a diferença importante de que todos os participantes são submetidos às quatro condições:

  • com cafeína e com álcool,
  • com cafeína e sem álcool,
  • sem cafeína e com álcool,
  • sem cafeína e sem álcool.

Dizemos, portanto, que ambos os fatores são intra participantes. Nesta modelagem da ANOVA é preciso controlar pelo participante, incluindo a variável de identificação como um efeito aleatório.

Os dados estão em CafIntra_AlcIntra.rds.

Dados <- data.frame(readxl::read_excel("CafIntra_AlcIntra.xlsx"))
Dados$UE <- factor(Dados$UE)
Dados$Cafeina <- factor(Dados$Cafeina, 
                        levels=c("SemC","ComC"))
Dados$Alcool <- factor(Dados$Alcool, 
                       levels=c("SemA","ComA"))
saveRDS(Dados, "CafIntra_AlcIntra.rds")
Dados <- readRDS("CafIntra_AlcIntra.rds")
str(Dados)
'data.frame':   48 obs. of  4 variables:
 $ UE      : Factor w/ 12 levels "1","2","3","4",..: 1 1 1 1 2 2 2 2 3 3 ...
 $ Alcool  : Factor w/ 2 levels "SemA","ComA": 1 1 2 2 1 1 2 2 1 1 ...
 $ Cafeina : Factor w/ 2 levels "SemC","ComC": 1 2 1 2 1 2 1 2 1 2 ...
 $ NumErros: num  4 15 28 10 9 8 22 11 10 17 ...
print.data.frame(Dados)
   UE Alcool Cafeina NumErros
1   1   SemA    SemC        4
2   1   SemA    ComC       15
3   1   ComA    SemC       28
4   1   ComA    ComC       10
5   2   SemA    SemC        9
6   2   SemA    ComC        8
7   2   ComA    SemC       22
8   2   ComA    ComC       11
9   3   SemA    SemC       10
10  3   SemA    ComC       17
11  3   ComA    SemC       21
12  3   ComA    ComC       27
13  4   SemA    SemC        8
14  4   SemA    ComC        0
15  4   ComA    SemC       27
16  4   ComA    ComC       15
17  5   SemA    SemC        6
18  5   SemA    ComC       15
19  5   ComA    SemC       21
20  5   ComA    ComC       27
21  6   SemA    SemC       11
22  6   SemA    ComC       11
23  6   ComA    SemC       20
24  6   ComA    ComC       10
25  7   SemA    SemC        2
26  7   SemA    ComC       11
27  7   ComA    SemC       19
28  7   ComA    ComC       21
29  8   SemA    SemC       11
30  8   SemA    ComC        6
31  8   ComA    SemC       16
32  8   ComA    ComC       15
33  9   SemA    SemC       11
34  9   SemA    ComC        0
35  9   ComA    SemC       25
36  9   ComA    ComC       19
37 10   SemA    SemC       10
38 10   SemA    ComC       15
39 10   ComA    SemC       17
40 10   ComA    ComC       21
41 11   SemA    SemC        3
42 11   SemA    ComC       17
43 11   ComA    SemC       19
44 11   ComA    ComC       15
45 12   SemA    SemC       10
46 12   SemA    ComC       15
47 12   ComA    SemC       20
48 12   ComA    ComC       15
print(summary(Dados[-1]))
  Alcool   Cafeina      NumErros    
 SemA:24   SemC:24   Min.   : 0.00  
 ComA:24   ComC:24   1st Qu.:10.00  
                     Median :15.00  
                     Mean   :14.29  
                     3rd Qu.:19.25  
                     Max.   :28.00  
print(xtabs(~Cafeina+Alcool, data=Dados))
       Alcool
Cafeina SemA ComA
   SemC   12   12
   ComC   12   12
print(ftable(Dados[-4]))
          Cafeina SemC ComC
UE Alcool                  
1  SemA              1    1
   ComA              1    1
2  SemA              1    1
   ComA              1    1
3  SemA              1    1
   ComA              1    1
4  SemA              1    1
   ComA              1    1
5  SemA              1    1
   ComA              1    1
6  SemA              1    1
   ComA              1    1
7  SemA              1    1
   ComA              1    1
8  SemA              1    1
   ComA              1    1
9  SemA              1    1
   ComA              1    1
10 SemA              1    1
   ComA              1    1
11 SemA              1    1
   ComA              1    1
12 SemA              1    1
   ComA              1    1

Da mesma forma que fizemos no exemplo anterior, o código completo está em demo_ANOVA2f_CafIntra_AlcIntra.R.

Estatística descritiva

Os dados estão no formato long. Observe a coluna das unidades experimentais (UE). Cada observação aparece em uma linha, mas a identificação das unidades experimentais é repetida, de forma que há apenas 12 participantes, cada um deles submetido aos quatro tratamentos.

Os gráficos que foram utilizados para a condição com dois fatores independentes não devem ser usados neste contexto, pois os intervalos de confiança dependem dos desvios-padrão e estes, por sua vez, pressupõem observações independentes. O gráfico de intervalo de confiança apresentado é o mais adequado para medidas repetidas.

ANOVA bifatorial relacionada e tamanho de efeito

A linha que define este modelo é

modelo <- lmerTest::lmer(NumErros ~ Alcool*Cafeina + (1|UE), 
                         data=Dados)

Note o termo (1|UE). Em delineamento com medidas repetidas, esta é forma de incluir a identificação do participante como um controle por meio de efeito aleatório.

A ANOVA é computada por

print(anv <- car::Anova(modelo,
                        test.statistic="F"))

O tamanho de efeito global (\(\eta^2\)) é dado por

eta2 <- as.numeric(MuMIn::r.squaredGLMM(modelo)[1])

Estas funções são adequadas para o participante modelado como efeito aleatório.

ANOVA bifatorial relacionada mostra que há efeito de interação entre Álcool e Cafeína para \(\alpha=5\%\).

Análise dos efeitos principais simples

A análise, neste caso, utiliza as mesmas funções do caso independente. Porém, o objeto veio de lmerTest::lmer.

Por ser um objeto diferente, a estatística de teste mudou, bem como os respectivos valores p. A conclusão ainda é a mesma: há efeito de álcool nos dois níveis de cafeína, mas não há efeito de cafeína nos dois níveis de álcool. O gráfico das médias marginais estimadas é praticamente o mesmo.

Testes post hoc

Os graus de liberdade são diferentes do caso com fatores independentes. Os intervalos, neste exemplo, ficaram praticamente iguais ao anterior.

source("summarySEwithin2.R")
Dados <- readRDS("CafIntra_AlcIntra.rds")
alfa <- 0.05
# str(Dados)
# print.data.frame(Dados)
# print(summary(Dados[-1]))
# print(xtabs(~Cafeina+Alcool, data=Dados))
# print(ftable(Dados[-4])) 
print(psych::describe(NumErros ~ Alcool*Cafeina,
                      data=Dados), 
      digits=2)

 Descriptive statistics by group 
Alcool: SemA
Cafeina: SemC
         vars  n mean   sd median trimmed  mad min max range  skew kurtosis
NumErros    1 12 7.92 3.32    9.5     8.2 2.22   2  11     9 -0.63    -1.35
           se
NumErros 0.96
------------------------------------------------------------ 
Alcool: ComA
Cafeina: SemC
         vars  n  mean   sd median trimmed  mad min max range skew kurtosis
NumErros    1 12 21.25 3.72   20.5    21.1 2.22  16  28    12 0.49    -1.04
           se
NumErros 1.07
------------------------------------------------------------ 
Alcool: SemA
Cafeina: ComC
         vars  n  mean   sd median trimmed  mad min max range skew kurtosis
NumErros    1 12 10.83 6.12     13    11.3 4.45   0  17    17 -0.7    -1.08
           se
NumErros 1.77
------------------------------------------------------------ 
Alcool: ComA
Cafeina: ComC
         vars  n  mean   sd median trimmed  mad min max range skew kurtosis
NumErros    1 12 17.17 5.92     15    16.9 6.67  10  27    17 0.41    -1.24
           se
NumErros 1.71
alfaBonf <- alfa/(length(unique(Dados$Alcool))*
                  length(unique(Dados$Cafeina)))
ic <- summarySEwithin2(Dados,
                       measurevar="NumErros",
                       withinvars=c("Alcool","Cafeina"),
                       idvar="UE",
                       na.rm=TRUE,
                       conf.interval=1-alfaBonf)
print(ic, digits=3)
  Alcool Cafeina  N NumErros NumErrosNormed   sd   se   ci
1   ComA    ComC 12    17.17          17.17 5.23 1.51 4.50
2   ComA    SemC 12    21.25          21.25 5.03 1.45 4.33
3   SemA    ComC 12    10.83          10.83 5.91 1.71 5.09
4   SemA    SemC 12     7.92           7.92 4.41 1.27 3.79
grf <- ggplot2::ggplot(ic,
                       ggplot2::aes(x=Alcool,
                                    y=NumErros,
                                    colour=Cafeina)) +
  ggplot2::geom_errorbar(position=ggplot2::position_dodge(.9),
                         width=.1,
                         ggplot2::aes(ymin=NumErros-ci,
                                      ymax=NumErros+ci)) +
  ggplot2::geom_point(shape=21,
                      size=3,
                      fill="white",
                      position=ggplot2::position_dodge(.9)) +
  ggplot2::ylab("NumErros") +
  ggplot2::ggtitle("Álcool & Cafeína: Número de Erros\nWithin-subject CI95% Bonferroni") +
  ggplot2::theme_bw()
print(grf)

# ANOVA bifatorial relacionada
cat("\nTwo-factor within-subjects ANOVA")

Two-factor within-subjects ANOVA
cat("\nGLMM")

GLMM
modelo <- lmerTest::lmer(NumErros ~ Alcool*Cafeina + (1|UE), 
                         data=Dados)
boundary (singular) fit: see help('isSingular')
cat("\nANOVA table")

ANOVA table
print(anv <- car::Anova(modelo,test.statistic="F"))
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: NumErros
                     F Df Df.res    Pr(>F)    
Alcool         47.6924  1     33 6.888e-08 ***
Cafeina         0.1678  1     33   0.68469    
Alcool:Cafeina  6.0420  1     33   0.01939 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(summary(modelo, correl=FALSE))
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: NumErros ~ Alcool * Cafeina + (1 | UE)
   Data: Dados

REML criterion at convergence: 275.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.19632 -0.48572  0.02534  0.76449  1.99358 

Random effects:
 Groups   Name        Variance Std.Dev.
 UE       (Intercept)  0.00    0.000   
 Residual             24.33    4.932   
Number of obs: 48, groups:  UE, 12

Fixed effects:
                       Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)               7.917      1.424 44.000   5.560 1.49e-06 ***
AlcoolComA               13.333      2.014 44.000   6.621 4.11e-08 ***
CafeinaComC               2.917      2.014 44.000   1.448    0.155    
AlcoolComA:CafeinaComC   -7.000      2.848 44.000  -2.458    0.018 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
cat("\nEffect size analysis")

Effect size analysis
eta2g <- as.numeric(MuMIn::r.squaredGLMM(modelo)[1])
cat("\nTamanho de efeito: eta^2 omnibus =", round(eta2g,4))

Tamanho de efeito: eta^2 omnibus = 0.5342
es <- effectsize::interpret_eta_squared(eta2g,
                                        rules="cohen1992")
names(es) <- c("Tamanho de efeito omnibus: estimativa pontual")
print(es)
Tamanho de efeito omnibus: estimativa pontual 
                                      "large" 
(Rules: cohen1992)
eta2 <- effectsize::eta_squared(anv,
                                partial=FALSE,
                                generalized=FALSE,
                                ci=1-alfa/3,
                                alternative="two.sided",
                                verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=4)
# Effect Size for ANOVA (Type II)

Parameter      | Eta2 (partial) |      98.3333% CI |  interpret
---------------------------------------------------------------
Alcool         |         0.5910 | [0.2991, 0.7528] |      large
Cafeina        |         0.0051 | [0.0000, 0.0103] | very small
Alcool:Cafeina |         0.1548 | [0.0000, 0.4254] |      large
cat("\nGráfico de perfil de médias")

Gráfico de perfil de médias
g1 <- emmeans::emmip(object=modelo,
                     formula=~Cafeina ~Alcool,
                     level=1-alfa,
                     adjust="holm",
                     CIs=TRUE,
                     dodge=.3,
                     lmer.df="satterthwaite",
                     lmerTest.limit=nrow(Dados)) + 
      ggplot2::theme_bw() 
g2 <- emmeans::emmip(object=modelo,
                     formula=~Alcool ~Cafeina,
                     level=1-alfa,
                     adjust="holm",
                     CIs=TRUE,
                     dodge=.3,
                     lmer.df="satterthwaite",
                     lmerTest.limit=nrow(Dados)) + 
      ggplot2::theme_bw() 
g3 <- emmeans::emmip(modelo, 
                     ~Cafeina,
                     level=1-alfa,
                     adjust="holm",
                     CIs=TRUE,
                     lmer.df="satterthwaite",
                     lmerTest.limit=nrow(Dados))
NOTE: Results may be misleading due to involvement in interactions
g4 <- emmeans::emmip(modelo, 
                     ~Alcool, 
                     level=1-alfa,
                     adjust="holm",
                     CIs=TRUE,
                     lmer.df="satterthwaite",
                     lmerTest.limit=nrow(Dados))
NOTE: Results may be misleading due to involvement in interactions
p <- (g1 | g2) / (g3 | g4)   # 2 x 2
print(p)

cat("\nPost hoc test se efeito de interação Alcool:Cafeina não é significante")

Post hoc test se efeito de interação Alcool:Cafeina não é significante
cat("\n\tAlcool")

    Alcool
EMM.A <- emmeans::emmeans(modelo, 
                          specs=pairwise~"Alcool", 
                          level=1-alfa,
                          adjust="holm",
                          lmer.df="satterthwaite",
                          lmerTest.limit=nrow(Dados))
NOTE: Results may be misleading due to involvement in interactions
print(summary(EMM.A$emmeans))
 Alcool emmean   SE df lower.CL upper.CL
 SemA     9.38 1.01 44     7.35     11.4
 ComA    19.21 1.01 44    17.18     21.2

Results are averaged over the levels of: Cafeina 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(summary(EMM.A$contrasts, 
              infer=TRUE))
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemA - ComA    -9.83 1.42 44    -12.7    -6.96  -6.906  <.0001

Results are averaged over the levels of: Cafeina 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(plot(EMM.A$emmeans, 
           colors="black"))

print(plot(EMM.A$contrasts, 
           colors="black"))

print(multcomp::cld(object=EMM.A$emmeans,
                    level=1-alfa,
                    alpha=alfa,
                    adjust="holm",
                    Letters=letters))
 Alcool emmean   SE df lower.CL upper.CL .group
 SemA     9.38 1.01 44     7.04     11.7  a    
 ComA    19.21 1.01 44    16.87     21.5   b   

Results are averaged over the levels of: Cafeina 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 2 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
cat("\n\tCafeina")

    Cafeina
EMM.C <- emmeans::emmeans(modelo, 
                          specs=pairwise~"Cafeina", 
                          level=1-alfa,
                          adjust="holm",
                          lmer.df="satterthwaite",
                          lmerTest.limit=nrow(Dados))
NOTE: Results may be misleading due to involvement in interactions
print(summary(EMM.C$emmeans))
 Cafeina emmean   SE df lower.CL upper.CL
 SemC      14.6 1.01 44     12.6     16.6
 ComC      14.0 1.01 44     12.0     16.0

Results are averaged over the levels of: Alcool 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(summary(EMM.C$contrasts, 
              infer=TRUE))
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemC - ComC    0.583 1.42 44    -2.29     3.45   0.410  0.6840

Results are averaged over the levels of: Alcool 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(plot(EMM.C$emmeans, 
           colors="black"))

print(plot(EMM.C$contrasts, 
           colors="black"))

print(multcomp::cld(object=EMM.C$emmeans,
                    level=1-alfa,
                    alpha=alfa,
                    adjust="holm",
                    Letters=letters))
 Cafeina emmean   SE df lower.CL upper.CL .group
 ComC      14.0 1.01 44     11.7     16.3  a    
 SemC      14.6 1.01 44     12.2     16.9  a    

Results are averaged over the levels of: Alcool 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 2 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
cat("\nPost hoc teste se efeito de interação Alcool:Cafeina é significante")

Post hoc teste se efeito de interação Alcool:Cafeina é significante
cat("\nPost hoc test: interaction Cafeina:Alcool")

Post hoc test: interaction Cafeina:Alcool
plot(emmeans::emmip(modelo, 
                    ~Cafeina:Alcool,
                    level=1-alfa,
                    adjust="holm",
                    CIs=TRUE))

EMM.AC <- emmeans::emmeans(modelo, 
                           specs=pairwise~"Alcool:Cafeina",
                           level=1-alfa,
                           adjust="holm",
                           lmer.df="satterthwaite",
                           lmerTest.limit=nrow(Dados))
print(summary(EMM.AC$emmeans))
 Alcool Cafeina emmean   SE df lower.CL upper.CL
 SemA   SemC      7.92 1.42 44     5.05     10.8
 ComA   SemC     21.25 1.42 44    18.38     24.1
 SemA   ComC     10.83 1.42 44     7.96     13.7
 ComA   ComC     17.17 1.42 44    14.30     20.0

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(summary(EMM.AC$contrasts, infer=TRUE))
 contrast              estimate   SE df lower.CL upper.CL t.ratio p.value
 SemA SemC - ComA SemC   -13.33 2.01 44   -18.90    -7.77  -6.621  <.0001
 SemA SemC - SemA ComC    -2.92 2.01 44    -8.48     2.65  -1.448  0.1546
 SemA SemC - ComA ComC    -9.25 2.01 44   -14.81    -3.69  -4.594  0.0001
 ComA SemC - SemA ComC    10.42 2.01 44     4.85    15.98   5.173  <.0001
 ComA SemC - ComA ComC     4.08 2.01 44    -1.48     9.65   2.028  0.0973
 SemA ComC - ComA ComC    -6.33 2.01 44   -11.90    -0.77  -3.145  0.0089

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 6 estimates 
P value adjustment: holm method for 6 tests 
print(plot(EMM.AC$emmeans, 
           colors="black"))

print(plot(EMM.AC$contrasts, 
           colors="black"))

print(multcomp::cld(object=EMM.AC$emmeans,
                    level=1-alfa,
                    alpha=alfa,
                    adjust="holm",
                    Letters=letters))
 Alcool Cafeina emmean   SE df lower.CL upper.CL .group
 SemA   SemC      7.92 1.42 44     4.21     11.6  a    
 SemA   ComC     10.83 1.42 44     7.12     14.5  a    
 ComA   ComC     17.17 1.42 44    13.46     20.9   b   
 ComA   SemC     21.25 1.42 44    17.54     25.0   b   

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 4 estimates 
P value adjustment: holm method for 6 tests 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
cat("\nSimple main effect test: Alcool effect by Cafeina levels\n")

Simple main effect test: Alcool effect by Cafeina levels
print(phia::testInteractions(modelo, 
                             across="Alcool",
                             fixed="Cafeina",
                             adjustment="holm")) 
Chisq Test: 
P-value adjustment method: holm
        Value     SE Df  Chisq Pr(>Chisq)    
SemC -13.3333 2.0137  1 43.842  7.118e-11 ***
ComC  -6.3333 2.0137  1  9.892    0.00166 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(emmeans::emmip(modelo,
                    ~Alcool|Cafeina,
                    level=1-alfa,
                    adjust="holm",
                    CIs=TRUE,
                    lmer.df="satterthwaite",
                    lmerTest.limit=nrow(Dados))) 

cat("\nPost hoc test of the simple main effect of Alcool")

Post hoc test of the simple main effect of Alcool
model_means <- emmeans::emmeans(object=modelo,
                                specs=pairwise~Alcool|Cafeina:Alcool,
                                level=1-alfa,
                                adjust="holm",
                                lmer.df="satterthwaite",
                                lmerTest.limit=nrow(Dados))
print(summary(model_means$emmeans, infer=TRUE))
Cafeina = SemC:
 Alcool emmean   SE df lower.CL upper.CL t.ratio p.value
 SemA     7.92 1.42 44     5.05     10.8   5.560  <.0001
 ComA    21.25 1.42 44    18.38     24.1  14.924  <.0001

Cafeina = ComC:
 Alcool emmean   SE df lower.CL upper.CL t.ratio p.value
 SemA    10.83 1.42 44     7.96     13.7   7.608  <.0001
 ComA    17.17 1.42 44    14.30     20.0  12.056  <.0001

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(summary(model_means$contrasts, infer=TRUE))
Cafeina = SemC:
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemA - ComA   -13.33 2.01 44    -17.4    -9.28  -6.621  <.0001

Cafeina = ComC:
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemA - ComA    -6.33 2.01 44    -10.4    -2.28  -3.145  0.0030

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
plot(model_means$emmeans,col="black")

plot(model_means$contrasts,col="black")

print(multcomp::cld(object=model_means$emmeans,
                    level=1-alfa,
                    alpha=alfa,
                    adjust="holm",
                    Letters=letters))
Cafeina = SemC:
 Alcool emmean   SE df lower.CL upper.CL .group
 SemA     7.92 1.42 44     4.61     11.2  a    
 ComA    21.25 1.42 44    17.95     24.6   b   

Cafeina = ComC:
 Alcool emmean   SE df lower.CL upper.CL .group
 SemA    10.83 1.42 44     7.53     14.1  a    
 ComA    17.17 1.42 44    13.86     20.5   b   

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 2 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
cat("\nSimple main effect test: Cafeina effect by Alcool levels\n")

Simple main effect test: Cafeina effect by Alcool levels
print(phia::testInteractions(modelo, 
                             across="Cafeina",
                             fixed="Alcool",
                             adjustment="holm"))
Chisq Test: 
P-value adjustment method: holm
       Value     SE Df  Chisq Pr(>Chisq)  
SemA -2.9167 2.0137  1 2.0979    0.14750  
ComA  4.0833 2.0137  1 4.1119    0.08516 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(emmeans::emmip(modelo, 
                    ~Cafeina|Alcool,
                    level=1-alfa,
                    adjust="holm",
                    CIs=TRUE,
                    lmer.df="satterthwaite",
                    lmerTest.limit=nrow(Dados))) 

cat("\nPost hoc test of the simple main effect of Cafeina")

Post hoc test of the simple main effect of Cafeina
model_means <- emmeans::emmeans(object=modelo,
                                specs=pairwise~Cafeina|Cafeina:Alcool,
                                level=1-alfa,
                                adjust="holm",
                                lmer.df="satterthwaite",
                                lmerTest.limit=nrow(Dados))
print(summary(model_means$emmeans, infer=TRUE))
Alcool = SemA:
 Cafeina emmean   SE df lower.CL upper.CL t.ratio p.value
 SemC      7.92 1.42 44     5.05     10.8   5.560  <.0001
 ComC     10.83 1.42 44     7.96     13.7   7.608  <.0001

Alcool = ComA:
 Cafeina emmean   SE df lower.CL upper.CL t.ratio p.value
 SemC     21.25 1.42 44    18.38     24.1  14.924  <.0001
 ComC     17.17 1.42 44    14.30     20.0  12.056  <.0001

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(summary(model_means$contrasts, infer=TRUE))
Alcool = SemA:
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemC - ComC    -2.92 2.01 44   -6.975     1.14  -1.448  0.1546

Alcool = ComA:
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemC - ComC     4.08 2.01 44    0.025     8.14   2.028  0.0487

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
plot(model_means$emmeans,col="black")

plot(model_means$contrasts,col="black")

print(multcomp::cld(object=model_means$emmeans,
                    level=1-alfa,
                    alpha=alfa,
                    adjust="holm",
                    Letters=letters))
Alcool = SemA:
 Cafeina emmean   SE df lower.CL upper.CL .group
 SemC      7.92 1.42 44     4.61     11.2  a    
 ComC     10.83 1.42 44     7.53     14.1  a    

Alcool = ComA:
 Cafeina emmean   SE df lower.CL upper.CL .group
 ComC     17.17 1.42 44    13.86     20.5  a    
 SemC     21.25 1.42 44    17.95     24.6   b   

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 2 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same.

ANOVA bifatorial mista

Testes:

  • Omnibus: lmerTest::lmer e car::Anova(test.statistic="F", ...), effectsize::eta_squared
    • Efeito principal simples: phia::testInteractions, emmeans::emmip
      • Post hoc: emmeans::emmeans, multcomp::cld, emmeans::emmip

Na ANOVA bifatorial mista, a identificação de cada um dos participantes (UE) aparece duas vezes, de tal forma que há apenas 24 participantes no total, alocados nos quatro tratamentos da seguinte maneira:

  • 12 participantes sem Álcool, submetidos às duas condições:
    • com Cafeína,
    • sem Cafeína
  • 12 participantes com Álcool, submetidos às duas condições:
    • com Cafeína,
    • sem Cafeína.

Por isso, o Álcool é o fator entre participantes, mas a Cafeína é intra participantes.

Sinonímia:

  • Two-factor mixed-design ANOVA
  • Two-way mixed-design ANOVA
  • Two-factor repeated measures ANOVA with split-plot structure
  • Two-way ANOVA with split-plot design
  • Two-way split-plot ANOVA

Neste caso, é usado General Linear Mixed Model (GLMM) por meio da função lmerTest::lmer.

Os dados estão em CafIntra_AlcEntre.rds.

Dados <- data.frame(readxl::read_excel("CafIntra_AlcEntre.xlsx"))
Dados$UE <- factor(Dados$UE)
Dados$Cafeina <- factor(Dados$Cafeina, 
                        levels=c("SemC","ComC"))
Dados$Alcool <- factor(Dados$Alcool, 
                       levels=c("SemA","ComA"))
saveRDS(Dados, "CafIntra_AlcEntre.rds")
Dados <- readRDS("CafIntra_AlcEntre.rds")
str(Dados)
'data.frame':   48 obs. of  4 variables:
 $ UE      : Factor w/ 24 levels "1","2","3","4",..: 1 1 2 2 3 3 4 4 5 5 ...
 $ NumErros: num  4 15 9 8 10 17 8 0 6 15 ...
 $ Alcool  : Factor w/ 2 levels "SemA","ComA": 1 1 1 1 1 1 1 1 1 1 ...
 $ Cafeina : Factor w/ 2 levels "SemC","ComC": 1 2 1 2 1 2 1 2 1 2 ...
print.data.frame(Dados)
   UE NumErros Alcool Cafeina
1   1        4   SemA    SemC
2   1       15   SemA    ComC
3   2        9   SemA    SemC
4   2        8   SemA    ComC
5   3       10   SemA    SemC
6   3       17   SemA    ComC
7   4        8   SemA    SemC
8   4        0   SemA    ComC
9   5        6   SemA    SemC
10  5       15   SemA    ComC
11  6       11   SemA    SemC
12  6       11   SemA    ComC
13  7        2   SemA    SemC
14  7       11   SemA    ComC
15  8       11   SemA    SemC
16  8        6   SemA    ComC
17  9       11   SemA    SemC
18  9        0   SemA    ComC
19 10       10   SemA    SemC
20 10       15   SemA    ComC
21 11        3   SemA    SemC
22 11       17   SemA    ComC
23 12       10   SemA    SemC
24 12       15   SemA    ComC
25 13       28   ComA    SemC
26 13       10   ComA    ComC
27 14       22   ComA    SemC
28 14       11   ComA    ComC
29 15       21   ComA    SemC
30 15       27   ComA    ComC
31 16       27   ComA    SemC
32 16       15   ComA    ComC
33 17       21   ComA    SemC
34 17       27   ComA    ComC
35 18       20   ComA    SemC
36 18       10   ComA    ComC
37 19       19   ComA    SemC
38 19       21   ComA    ComC
39 20       16   ComA    SemC
40 20       15   ComA    ComC
41 21       25   ComA    SemC
42 21       19   ComA    ComC
43 22       17   ComA    SemC
44 22       21   ComA    ComC
45 23       19   ComA    SemC
46 23       15   ComA    ComC
47 24       20   ComA    SemC
48 24       15   ComA    ComC
print(summary(Dados[-1]))
    NumErros      Alcool   Cafeina  
 Min.   : 0.00   SemA:24   SemC:24  
 1st Qu.:10.00   ComA:24   ComC:24  
 Median :15.00                      
 Mean   :14.29                      
 3rd Qu.:19.25                      
 Max.   :28.00                      
print(xtabs(~Cafeina+Alcool, data=Dados))
       Alcool
Cafeina SemA ComA
   SemC   12   12
   ComC   12   12
print(ftable(Dados[-2]))
          Cafeina SemC ComC
UE Alcool                  
1  SemA              1    1
   ComA              0    0
2  SemA              1    1
   ComA              0    0
3  SemA              1    1
   ComA              0    0
4  SemA              1    1
   ComA              0    0
5  SemA              1    1
   ComA              0    0
6  SemA              1    1
   ComA              0    0
7  SemA              1    1
   ComA              0    0
8  SemA              1    1
   ComA              0    0
9  SemA              1    1
   ComA              0    0
10 SemA              1    1
   ComA              0    0
11 SemA              1    1
   ComA              0    0
12 SemA              1    1
   ComA              0    0
13 SemA              0    0
   ComA              1    1
14 SemA              0    0
   ComA              1    1
15 SemA              0    0
   ComA              1    1
16 SemA              0    0
   ComA              1    1
17 SemA              0    0
   ComA              1    1
18 SemA              0    0
   ComA              1    1
19 SemA              0    0
   ComA              1    1
20 SemA              0    0
   ComA              1    1
21 SemA              0    0
   ComA              1    1
22 SemA              0    0
   ComA              1    1
23 SemA              0    0
   ComA              1    1
24 SemA              0    0
   ComA              1    1

Implementado em demo_ANOVA2f_CafIntra_AlcEntre.R. Este é o mesmo Rscript usado para a ANOVA relacionada bifatorial. Somente a estrutura da planilha (a coluna com a variável UE é a diferença) basta para que o comportamento da análise seja adaptado por lmerTest::lmer:

source("summarySEwithin2.R")
Dados <- readRDS("CafIntra_AlcEntre.rds")
alfa <- 0.05
# str(Dados)
# print.data.frame(Dados)
# print(summary(Dados[-1]))
# print(xtabs(~Cafeina+Alcool, data=Dados))
# print(ftable(Dados[-2]))
print(psych::describe(NumErros ~ Alcool*Cafeina, 
                      data=Dados), digits=2) 

 Descriptive statistics by group 
Alcool: SemA
Cafeina: SemC
         vars  n mean   sd median trimmed  mad min max range  skew kurtosis
NumErros    1 12 7.92 3.32    9.5     8.2 2.22   2  11     9 -0.63    -1.35
           se
NumErros 0.96
------------------------------------------------------------ 
Alcool: ComA
Cafeina: SemC
         vars  n  mean   sd median trimmed  mad min max range skew kurtosis
NumErros    1 12 21.25 3.72   20.5    21.1 2.22  16  28    12 0.49    -1.04
           se
NumErros 1.07
------------------------------------------------------------ 
Alcool: SemA
Cafeina: ComC
         vars  n  mean   sd median trimmed  mad min max range skew kurtosis
NumErros    1 12 10.83 6.12     13    11.3 4.45   0  17    17 -0.7    -1.08
           se
NumErros 1.77
------------------------------------------------------------ 
Alcool: ComA
Cafeina: ComC
         vars  n  mean   sd median trimmed  mad min max range skew kurtosis
NumErros    1 12 17.17 5.92     15    16.9 6.67  10  27    17 0.41    -1.24
           se
NumErros 1.71
alfaBonf <- alfa/(length(unique(Dados$Alcool))*
                  length(unique(Dados$Cafeina)))
ic <- summarySEwithin2(Dados,
                       measurevar="NumErros",
                       withinvars=c("Alcool","Cafeina"),
                       idvar="UE",
                       na.rm=TRUE,
                       conf.interval=1-alfaBonf)
print(ic)
  Alcool Cafeina  N  NumErros NumErrosNormed       sd       se       ci
1   ComA    ComC 12 17.166667       12.25000 4.468464 1.289934 3.845130
2   ComA    SemC 12 21.250000       16.33333 4.468464 1.289934 3.845130
3   SemA    ComC 12 10.833333       15.75000 4.562418 1.317057 3.925978
4   SemA    SemC 12  7.916667       12.83333 4.562418 1.317057 3.925978
grf <- ggplot2::ggplot(ic,
                       ggplot2::aes(x=Alcool,
                                    y=NumErros,
                                    colour=Cafeina)) +
  ggplot2::geom_errorbar(position=ggplot2::position_dodge(.9),
                         width=.1,
                         ggplot2::aes(ymin=NumErros-ci,
                                      ymax=NumErros+ci)) +
  ggplot2::geom_point(shape=21,
                      size=3,
                      fill="white",
                      position=ggplot2::position_dodge(.9)) +
  ggplot2::ylab("NumErros") +
  ggplot2::ggtitle("Álcool & Cafeína: Número de Erros\nWithin-subject CI95% Bonferroni") +
  ggplot2::theme_bw()
print(grf)

# ANOVA bifatorial mista
cat("\nTwo-factor mixed-desing ANOVA")

Two-factor mixed-desing ANOVA
cat("\nGLMM")

GLMM
modelo <- lmerTest::lmer(NumErros ~ Alcool*Cafeina + (1|UE), 
                         data=Dados)
boundary (singular) fit: see help('isSingular')
cat("\nANOVA table")

ANOVA table
print(anv <- car::Anova(modelo,test.statistic="F"))
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: NumErros
                     F Df Df.res    Pr(>F)    
Alcool         47.6924  1     22 6.191e-07 ***
Cafeina         0.1678  1     22   0.68601    
Alcool:Cafeina  6.0420  1     22   0.02231 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(summary(modelo, correl=FALSE))
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: NumErros ~ Alcool * Cafeina + (1 | UE)
   Data: Dados

REML criterion at convergence: 275.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.19632 -0.48572  0.02534  0.76449  1.99358 

Random effects:
 Groups   Name        Variance Std.Dev.
 UE       (Intercept)  0.00    0.000   
 Residual             24.33    4.932   
Number of obs: 48, groups:  UE, 24

Fixed effects:
                       Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)               7.917      1.424 44.000   5.560 1.49e-06 ***
AlcoolComA               13.333      2.014 44.000   6.621 4.11e-08 ***
CafeinaComC               2.917      2.014 44.000   1.448    0.155    
AlcoolComA:CafeinaComC   -7.000      2.848 44.000  -2.458    0.018 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
cat("\nEffect size analysis")

Effect size analysis
eta2g <- as.numeric(MuMIn::r.squaredGLMM(modelo)[1])
cat("\nTamanho de efeito: eta^2 omnibus =", round(eta2g,4))

Tamanho de efeito: eta^2 omnibus = 0.5342
es <- effectsize::interpret_eta_squared(eta2g,
                                        rules="cohen1992")
names(es) <- c("Tamanho de efeito omnibus: estimativa pontual")
print(es)
Tamanho de efeito omnibus: estimativa pontual 
                                      "large" 
(Rules: cohen1992)
eta2 <- effectsize::eta_squared(anv,
                                partial=FALSE,
                                generalized=FALSE,
                                ci=1-alfa/3,
                                alternative="two.sided",
                                verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=4)
# Effect Size for ANOVA (Type II)

Parameter      | Eta2 (partial) |      98.3333% CI |  interpret
---------------------------------------------------------------
Alcool         |         0.6843 | [0.3511, 0.8293] |      large
Cafeina        |         0.0076 | [0.0000, 0.0153] | very small
Alcool:Cafeina |         0.2155 | [0.0000, 0.5306] |      large
cat("\nGráfico de perfil de médias")

Gráfico de perfil de médias
g1 <- emmeans::emmip(object=modelo,
                     formula=~Cafeina ~Alcool,
                     level=1-alfa,
                     adjust="holm",
                     CIs=TRUE,
                     dodge=.3,
                     lmer.df="satterthwaite",
                     lmerTest.limit=nrow(Dados)) + 
      ggplot2::theme_bw() 
g2 <- emmeans::emmip(object=modelo,
                     formula=~Alcool ~Cafeina,
                     level=1-alfa,
                     adjust="holm",
                     CIs=TRUE,
                     dodge=.3,
                     lmer.df="satterthwaite",
                     lmerTest.limit=nrow(Dados)) + 
      ggplot2::theme_bw() 
g3 <- emmeans::emmip(modelo, 
                     ~Cafeina,
                     level=1-alfa,
                     adjust="holm",
                     CIs=TRUE,
                     lmer.df="satterthwaite",
                     lmerTest.limit=nrow(Dados))
NOTE: Results may be misleading due to involvement in interactions
g4 <- emmeans::emmip(modelo, 
                     ~Alcool, 
                     level=1-alfa,
                     adjust="holm",
                     CIs=TRUE,
                     lmer.df="satterthwaite",
                     lmerTest.limit=nrow(Dados))
NOTE: Results may be misleading due to involvement in interactions
p <- (g1 | g2) / (g3 | g4)   # 2 x 2
print(p)

cat("\nPost hoc test se efeito de interação Alcool:Cafeina não é significante")

Post hoc test se efeito de interação Alcool:Cafeina não é significante
cat("\n\tAlcool")

    Alcool
EMM.A <- emmeans::emmeans(modelo, 
                          specs=pairwise~"Alcool", 
                          level=1-alfa,
                          adjust="holm",
                          lmer.df="satterthwaite",
                          lmerTest.limit=nrow(Dados))
NOTE: Results may be misleading due to involvement in interactions
print(summary(EMM.A$emmeans))
 Alcool emmean   SE df lower.CL upper.CL
 SemA     9.38 1.01 44     7.35     11.4
 ComA    19.21 1.01 44    17.18     21.2

Results are averaged over the levels of: Cafeina 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(summary(EMM.A$contrasts, 
              infer=TRUE))
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemA - ComA    -9.83 1.42 44    -12.7    -6.96  -6.906  <.0001

Results are averaged over the levels of: Cafeina 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(plot(EMM.A$emmeans, 
           colors="black"))

print(plot(EMM.A$contrasts, 
           colors="black"))

print(multcomp::cld(object=EMM.A$emmeans,
                    level=1-alfa,
                    alpha=alfa,
                    adjust="holm",
                    Letters=letters))
 Alcool emmean   SE df lower.CL upper.CL .group
 SemA     9.38 1.01 44     7.04     11.7  a    
 ComA    19.21 1.01 44    16.87     21.5   b   

Results are averaged over the levels of: Cafeina 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 2 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
cat("\n\tCafeina")

    Cafeina
EMM.C <- emmeans::emmeans(modelo, 
                          specs=pairwise~"Cafeina", 
                          level=1-alfa,
                          adjust="holm",
                          lmer.df="satterthwaite",
                          lmerTest.limit=nrow(Dados))
NOTE: Results may be misleading due to involvement in interactions
print(summary(EMM.C$emmeans))
 Cafeina emmean   SE df lower.CL upper.CL
 SemC      14.6 1.01 44     12.6     16.6
 ComC      14.0 1.01 44     12.0     16.0

Results are averaged over the levels of: Alcool 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(summary(EMM.C$contrasts, 
              infer=TRUE))
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemC - ComC    0.583 1.42 44    -2.29     3.45   0.410  0.6840

Results are averaged over the levels of: Alcool 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(plot(EMM.C$emmeans, 
           colors="black"))

print(plot(EMM.C$contrasts, 
           colors="black"))

print(multcomp::cld(object=EMM.C$emmeans,
                    level=1-alfa,
                    alpha=alfa,
                    adjust="holm",
                    Letters=letters))
 Cafeina emmean   SE df lower.CL upper.CL .group
 ComC      14.0 1.01 44     11.7     16.3  a    
 SemC      14.6 1.01 44     12.2     16.9  a    

Results are averaged over the levels of: Alcool 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 2 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
cat("\nPost hoc teste se efeito de interação Alcool:Cafeina é significante")

Post hoc teste se efeito de interação Alcool:Cafeina é significante
cat("\nPost hoc test: interaction Cafeina:Alcool")

Post hoc test: interaction Cafeina:Alcool
plot(emmeans::emmip(modelo, 
                    ~Cafeina:Alcool,
                    level=1-alfa,
                    adjust="holm",
                    CIs=TRUE))

EMM.AC <- emmeans::emmeans(modelo, 
                           specs=pairwise~"Alcool:Cafeina",
                           level=1-alfa,
                           adjust="holm",
                           lmer.df="satterthwaite",
                           lmerTest.limit=nrow(Dados))
print(summary(EMM.AC$emmeans))
 Alcool Cafeina emmean   SE df lower.CL upper.CL
 SemA   SemC      7.92 1.42 44     5.05     10.8
 ComA   SemC     21.25 1.42 44    18.38     24.1
 SemA   ComC     10.83 1.42 44     7.96     13.7
 ComA   ComC     17.17 1.42 44    14.30     20.0

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(summary(EMM.AC$contrasts, infer=TRUE))
 contrast              estimate   SE df lower.CL upper.CL t.ratio p.value
 SemA SemC - ComA SemC   -13.33 2.01 44   -18.90    -7.77  -6.621  <.0001
 SemA SemC - SemA ComC    -2.92 2.01 44    -8.48     2.65  -1.448  0.1546
 SemA SemC - ComA ComC    -9.25 2.01 44   -14.81    -3.69  -4.594  0.0001
 ComA SemC - SemA ComC    10.42 2.01 44     4.85    15.98   5.173  <.0001
 ComA SemC - ComA ComC     4.08 2.01 44    -1.48     9.65   2.028  0.0973
 SemA ComC - ComA ComC    -6.33 2.01 44   -11.90    -0.77  -3.145  0.0089

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 6 estimates 
P value adjustment: holm method for 6 tests 
print(plot(EMM.AC$emmeans, 
           colors="black"))

print(plot(EMM.AC$contrasts, 
           colors="black"))

print(multcomp::cld(object=EMM.AC$emmeans,
                    level=1-alfa,
                    alpha=alfa,
                    adjust="holm",
                    Letters=letters))
 Alcool Cafeina emmean   SE df lower.CL upper.CL .group
 SemA   SemC      7.92 1.42 44     4.21     11.6  a    
 SemA   ComC     10.83 1.42 44     7.12     14.5  a    
 ComA   ComC     17.17 1.42 44    13.46     20.9   b   
 ComA   SemC     21.25 1.42 44    17.54     25.0   b   

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 4 estimates 
P value adjustment: holm method for 6 tests 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
cat("\nSimple main effect test: Alcool effect by Cafeina levels\n")

Simple main effect test: Alcool effect by Cafeina levels
print(phia::testInteractions(modelo, 
                             across="Alcool",
                             fixed="Cafeina",
                             adjustment="holm")) 
Chisq Test: 
P-value adjustment method: holm
        Value     SE Df  Chisq Pr(>Chisq)    
SemC -13.3333 2.0137  1 43.842  7.118e-11 ***
ComC  -6.3333 2.0137  1  9.892    0.00166 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(emmeans::emmip(modelo,
                    ~Alcool|Cafeina,
                    level=1-alfa,
                    adjust="holm",
                    CIs=TRUE,
                    lmer.df="satterthwaite",
                    lmerTest.limit=nrow(Dados))) 

cat("\nPost hoc test of the simple main effect of Alcool")

Post hoc test of the simple main effect of Alcool
model_means <- emmeans::emmeans(object=modelo,
                                specs=pairwise~Alcool|Cafeina:Alcool,
                                level=1-alfa,
                                adjust="holm",
                                lmer.df="satterthwaite",
                                lmerTest.limit=nrow(Dados))
print(summary(model_means$emmeans, infer=TRUE))
Cafeina = SemC:
 Alcool emmean   SE df lower.CL upper.CL t.ratio p.value
 SemA     7.92 1.42 44     5.05     10.8   5.560  <.0001
 ComA    21.25 1.42 44    18.38     24.1  14.924  <.0001

Cafeina = ComC:
 Alcool emmean   SE df lower.CL upper.CL t.ratio p.value
 SemA    10.83 1.42 44     7.96     13.7   7.608  <.0001
 ComA    17.17 1.42 44    14.30     20.0  12.056  <.0001

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(summary(model_means$contrasts, infer=TRUE))
Cafeina = SemC:
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemA - ComA   -13.33 2.01 44    -17.4    -9.28  -6.621  <.0001

Cafeina = ComC:
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemA - ComA    -6.33 2.01 44    -10.4    -2.28  -3.145  0.0030

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
plot(model_means$emmeans,col="black")

plot(model_means$contrasts,col="black")

print(multcomp::cld(object=model_means$emmeans,
                    level=1-alfa,
                    alpha=alfa,
                    adjust="holm",
                    Letters=letters))
Cafeina = SemC:
 Alcool emmean   SE df lower.CL upper.CL .group
 SemA     7.92 1.42 44     4.61     11.2  a    
 ComA    21.25 1.42 44    17.95     24.6   b   

Cafeina = ComC:
 Alcool emmean   SE df lower.CL upper.CL .group
 SemA    10.83 1.42 44     7.53     14.1  a    
 ComA    17.17 1.42 44    13.86     20.5   b   

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 2 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
cat("\nSimple main effect test: Cafeina effect by Alcool levels\n")

Simple main effect test: Cafeina effect by Alcool levels
print(phia::testInteractions(modelo, 
                             across="Cafeina",
                             fixed="Alcool",
                             adjustment="holm"))
Chisq Test: 
P-value adjustment method: holm
       Value     SE Df  Chisq Pr(>Chisq)  
SemA -2.9167 2.0137  1 2.0979    0.14750  
ComA  4.0833 2.0137  1 4.1119    0.08516 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(emmeans::emmip(modelo, 
                    ~Cafeina|Alcool,
                    level=1-alfa,
                    adjust="holm",
                    CIs=TRUE,
                    lmer.df="satterthwaite",
                    lmerTest.limit=nrow(Dados))) 

cat("\nPost hoc test of the simple main effect of Cafeina")

Post hoc test of the simple main effect of Cafeina
model_means <- emmeans::emmeans(object=modelo,
                                specs=pairwise~Cafeina|Cafeina:Alcool,
                                level=1-alfa,
                                adjust="holm",
                                lmer.df="satterthwaite",
                                lmerTest.limit=nrow(Dados))
print(summary(model_means$emmeans, infer=TRUE))
Alcool = SemA:
 Cafeina emmean   SE df lower.CL upper.CL t.ratio p.value
 SemC      7.92 1.42 44     5.05     10.8   5.560  <.0001
 ComC     10.83 1.42 44     7.96     13.7   7.608  <.0001

Alcool = ComA:
 Cafeina emmean   SE df lower.CL upper.CL t.ratio p.value
 SemC     21.25 1.42 44    18.38     24.1  14.924  <.0001
 ComC     17.17 1.42 44    14.30     20.0  12.056  <.0001

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(summary(model_means$contrasts, infer=TRUE))
Alcool = SemA:
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemC - ComC    -2.92 2.01 44   -6.975     1.14  -1.448  0.1546

Alcool = ComA:
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemC - ComC     4.08 2.01 44    0.025     8.14   2.028  0.0487

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
plot(model_means$emmeans,col="black")

plot(model_means$contrasts,col="black")

print(multcomp::cld(object=model_means$emmeans,
                    level=1-alfa,
                    alpha=alfa,
                    adjust="holm",
                    Letters=letters))
Alcool = SemA:
 Cafeina emmean   SE df lower.CL upper.CL .group
 SemC      7.92 1.42 44     4.61     11.2  a    
 ComC     10.83 1.42 44     7.53     14.1  a    

Alcool = ComA:
 Cafeina emmean   SE df lower.CL upper.CL .group
 ComC     17.17 1.42 44    13.86     20.5  a    
 SemC     21.25 1.42 44    17.95     24.6   b   

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 2 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same.

Na ANOVA independente é requerido pelo menos duas observações em cada um dos tratamentos para testar efeito de interação.

Para ANOVA relacionada ou mista é requerido pelo menos uma observação de cada tratamento para testar efeito de interação.

Quando um modelo completo como

modelo <- lm(NumErros ~ Alcool+Cafeina+Alcool:Cafeina,
             data=Dados)

não pode ser testado ou não se deseja testar, uma alternativa é usar um modelo de efeitos principais:

modelo <- lm(NumErros ~ Alcool+Cafeina, data=Dados)

ANOVA bifatorial para medidas repetidas ou mista com dados faltantes

Testes:

  • Omnibus: lmerTest::lmer e car::Anova(test.statistic="F", ...), effectsize::eta_squared
    • Efeito principal simples: phia::testInteractions, emmeans::emmip
      • Post hoc: emmeans::emmeans, multcomp::cld, emmeans::emmip
Para todos códigos R de ANOVA bifatorial fornecidos neste capítulo, dados desbalanceados podem ser utilizados.

Neste caso, é usado General Linear Mixed Model (GLMM) por meio da função lmerTest::lmer.

Por exemplo, trocando o arquivo de dados no caso de ANOVA bifatorial mista para CafIntra_AlcEntre_Missing.rds (com valor faltante), a saída gerada por demo_ANOVA2f_AlcEntre_CafeIntra_Missing.R resulta em:

Dados <- data.frame(readxl::read_excel("CafIntra_AlcEntre_Missing.xlsx"))
Dados$UE <- factor(Dados$UE)
Dados$Cafeina <- factor(Dados$Cafeina, 
                        levels=c("SemC","ComC"))
Dados$Alcool <- factor(Dados$Alcool, 
                       levels=c("SemA","ComA"))
saveRDS(Dados, "CafIntra_AlcEntre_Missing.rds")
Dados <- readRDS("CafIntra_AlcEntre_Missing.rds")
str(Dados)
'data.frame':   47 obs. of  4 variables:
 $ UE      : Factor w/ 24 levels "1","2","3","4",..: 1 2 2 3 3 4 4 5 5 6 ...
 $ NumErros: num  15 9 8 10 17 8 0 6 15 11 ...
 $ Alcool  : Factor w/ 2 levels "SemA","ComA": 1 1 1 1 1 1 1 1 1 1 ...
 $ Cafeina : Factor w/ 2 levels "SemC","ComC": 2 1 2 1 2 1 2 1 2 1 ...
print.data.frame(Dados)
   UE NumErros Alcool Cafeina
1   1       15   SemA    ComC
2   2        9   SemA    SemC
3   2        8   SemA    ComC
4   3       10   SemA    SemC
5   3       17   SemA    ComC
6   4        8   SemA    SemC
7   4        0   SemA    ComC
8   5        6   SemA    SemC
9   5       15   SemA    ComC
10  6       11   SemA    SemC
11  6       11   SemA    ComC
12  7        2   SemA    SemC
13  7       11   SemA    ComC
14  8       11   SemA    SemC
15  8        6   SemA    ComC
16  9       11   SemA    SemC
17  9        0   SemA    ComC
18 10       10   SemA    SemC
19 10       15   SemA    ComC
20 11        3   SemA    SemC
21 11       17   SemA    ComC
22 12       10   SemA    SemC
23 12       15   SemA    ComC
24 13       28   ComA    SemC
25 13       10   ComA    ComC
26 14       22   ComA    SemC
27 14       11   ComA    ComC
28 15       21   ComA    SemC
29 15       27   ComA    ComC
30 16       27   ComA    SemC
31 16       15   ComA    ComC
32 17       21   ComA    SemC
33 17       27   ComA    ComC
34 18       20   ComA    SemC
35 18       10   ComA    ComC
36 19       19   ComA    SemC
37 19       21   ComA    ComC
38 20       16   ComA    SemC
39 20       15   ComA    ComC
40 21       25   ComA    SemC
41 21       19   ComA    ComC
42 22       17   ComA    SemC
43 22       21   ComA    ComC
44 23       19   ComA    SemC
45 23       15   ComA    ComC
46 24       20   ComA    SemC
47 24       15   ComA    ComC
print(summary(Dados[-1]))
    NumErros      Alcool   Cafeina  
 Min.   : 0.00   SemA:23   SemC:23  
 1st Qu.:10.00   ComA:24   ComC:24  
 Median :15.00                      
 Mean   :14.51                      
 3rd Qu.:19.50                      
 Max.   :28.00                      
print(xtabs(~Cafeina+Alcool, data=Dados))
       Alcool
Cafeina SemA ComA
   SemC   11   12
   ComC   12   12
print(ftable(Dados[,-2]))
          Cafeina SemC ComC
UE Alcool                  
1  SemA              0    1
   ComA              0    0
2  SemA              1    1
   ComA              0    0
3  SemA              1    1
   ComA              0    0
4  SemA              1    1
   ComA              0    0
5  SemA              1    1
   ComA              0    0
6  SemA              1    1
   ComA              0    0
7  SemA              1    1
   ComA              0    0
8  SemA              1    1
   ComA              0    0
9  SemA              1    1
   ComA              0    0
10 SemA              1    1
   ComA              0    0
11 SemA              1    1
   ComA              0    0
12 SemA              1    1
   ComA              0    0
13 SemA              0    0
   ComA              1    1
14 SemA              0    0
   ComA              1    1
15 SemA              0    0
   ComA              1    1
16 SemA              0    0
   ComA              1    1
17 SemA              0    0
   ComA              1    1
18 SemA              0    0
   ComA              1    1
19 SemA              0    0
   ComA              1    1
20 SemA              0    0
   ComA              1    1
21 SemA              0    0
   ComA              1    1
22 SemA              0    0
   ComA              1    1
23 SemA              0    0
   ComA              1    1
24 SemA              0    0
   ComA              1    1
source("summarySEwithin2.R")
Dados <- readRDS("CafIntra_AlcEntre_Missing.rds")
alfa <- 0.05
# str(Dados)
# print.data.frame(Dados)
# print(summary(Dados[-1]))
# print(xtabs(~Cafeina+Alcool, data=Dados))
# print(ftable(Dados[-2]))
print(psych::describe(NumErros ~ Alcool*Cafeina, 
                      data=Dados), digits=2) 

 Descriptive statistics by group 
Alcool: SemA
Cafeina: SemC
         vars  n mean   sd median trimmed  mad min max range  skew kurtosis
NumErros    1 11 8.27 3.23     10    8.67 1.48   2  11     9 -0.89    -0.87
           se
NumErros 0.97
------------------------------------------------------------ 
Alcool: ComA
Cafeina: SemC
         vars  n  mean   sd median trimmed  mad min max range skew kurtosis
NumErros    1 12 21.25 3.72   20.5    21.1 2.22  16  28    12 0.49    -1.04
           se
NumErros 1.07
------------------------------------------------------------ 
Alcool: SemA
Cafeina: ComC
         vars  n  mean   sd median trimmed  mad min max range skew kurtosis
NumErros    1 12 10.83 6.12     13    11.3 4.45   0  17    17 -0.7    -1.08
           se
NumErros 1.77
------------------------------------------------------------ 
Alcool: ComA
Cafeina: ComC
         vars  n  mean   sd median trimmed  mad min max range skew kurtosis
NumErros    1 12 17.17 5.92     15    16.9 6.67  10  27    17 0.41    -1.24
           se
NumErros 1.71
alfaBonf <- alfa/(length(unique(Dados$Alcool))*
                  length(unique(Dados$Cafeina)))
ic <- summarySEwithin2(Dados,
                       measurevar="NumErros",
                       withinvars=c("Alcool","Cafeina"),
                       idvar="UE",
                       na.rm=TRUE,
                       conf.interval=1-alfaBonf)
print(ic)
  Alcool Cafeina  N  NumErros NumErrosNormed       sd       se       ci
1   ComA    ComC 12 17.166667       12.46897 4.468464 1.289934 3.845130
2   ComA    SemC 12 21.250000       16.55230 4.468464 1.289934 3.845130
3   SemA    ComC 12 10.833333       15.51064 4.334499 1.251262 3.729852
4   SemA    SemC 11  8.272727       13.41973 4.530034 1.365857 4.149805
grf <- ggplot2::ggplot(ic,
                       ggplot2::aes(x=Alcool,
                                    y=NumErros,
                                    colour=Cafeina)) +
  ggplot2::geom_errorbar(position=ggplot2::position_dodge(.9),
                         width=.1,
                         ggplot2::aes(ymin=NumErros-ci,
                                      ymax=NumErros+ci)) +
  ggplot2::geom_point(shape=21,
                      size=3,
                      fill="white",
                      position=ggplot2::position_dodge(.9)) +
  ggplot2::ylab("NumErros") +
  ggplot2::ggtitle("Álcool & Cafeína: Número de Erros\nWithin-subject CI95% Bonferroni") +
  ggplot2::theme_bw()
print(grf)

# ANOVA bifatorial mista
cat("\nTwo-factor mixed-desing ANOVA")

Two-factor mixed-desing ANOVA
cat("\nGLMM")

GLMM
modelo <- lmerTest::lmer(NumErros ~ Alcool*Cafeina + (1|UE), 
                         data=Dados)
boundary (singular) fit: see help('isSingular')
cat("\nANOVA table")

ANOVA table
print(anv <- car::Anova(modelo,test.statistic="F"))
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: NumErros
                     F Df Df.res    Pr(>F)    
Alcool         43.8255  1 21.732 1.249e-06 ***
Cafeina         0.3390  1 21.707   0.56643    
Alcool:Cafeina  5.2725  1 21.730   0.03169 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(summary(modelo, correl=FALSE))
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: NumErros ~ Alcool * Cafeina + (1 | UE)
   Data: Dados

REML criterion at convergence: 269.4

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.18839 -0.45681  0.03367  0.76594  1.98638 

Random effects:
 Groups   Name        Variance Std.Dev.
 UE       (Intercept)  0.00    0.00    
 Residual             24.51    4.95    
Number of obs: 47, groups:  UE, 24

Fixed effects:
                       Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)               8.273      1.493 43.000   5.543 1.69e-06 ***
AlcoolComA               12.977      2.066 43.000   6.280 1.43e-07 ***
CafeinaComC               2.561      2.066 43.000   1.239   0.2220    
AlcoolComA:CafeinaComC   -6.644      2.890 43.000  -2.299   0.0264 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
cat("\nEffect size analysis")

Effect size analysis
eta2g <- as.numeric(MuMIn::r.squaredGLMM(modelo)[1])
cat("\nTamanho de efeito: eta^2 omnibus =", round(eta2g,4))

Tamanho de efeito: eta^2 omnibus = 0.5197
es <- effectsize::interpret_eta_squared(eta2g,
                                        rules="cohen1992")
names(es) <- c("Tamanho de efeito omnibus: estimativa pontual")
print(es)
Tamanho de efeito omnibus: estimativa pontual 
                                      "large" 
(Rules: cohen1992)
eta2 <- effectsize::eta_squared(anv,
                                partial=FALSE,
                                generalized=FALSE,
                                ci=1-alfa/3,
                                alternative="two.sided",
                                verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=4)
# Effect Size for ANOVA (Type II)

Parameter      | Eta2 (partial) |      98.3333% CI | interpret
--------------------------------------------------------------
Alcool         |         0.6685 | [0.3237, 0.8211] |     large
Cafeina        |         0.0154 | [0.0000, 0.2889] |     small
Alcool:Cafeina |         0.1953 | [0.0000, 0.5155] |     large
cat("\nGráfico de perfil de médias")

Gráfico de perfil de médias
g1 <- emmeans::emmip(object=modelo,
                     formula=~Cafeina ~Alcool,
                     level=1-alfa,
                     adjust="holm",
                     CIs=TRUE,
                     dodge=.3,
                     lmer.df="satterthwaite",
                     lmerTest.limit=nrow(Dados)) + 
      ggplot2::theme_bw() 
g2 <- emmeans::emmip(object=modelo,
                     formula=~Alcool ~Cafeina,
                     level=1-alfa,
                     adjust="holm",
                     CIs=TRUE,
                     dodge=.3,
                     lmer.df="satterthwaite",
                     lmerTest.limit=nrow(Dados)) + 
      ggplot2::theme_bw() 
g3 <- emmeans::emmip(modelo, 
                     ~Cafeina,
                     level=1-alfa,
                     adjust="holm",
                     CIs=TRUE,
                     lmer.df="satterthwaite",
                     lmerTest.limit=nrow(Dados))
NOTE: Results may be misleading due to involvement in interactions
g4 <- emmeans::emmip(modelo, 
                     ~Alcool, 
                     level=1-alfa,
                     adjust="holm",
                     CIs=TRUE,
                     lmer.df="satterthwaite",
                     lmerTest.limit=nrow(Dados))
NOTE: Results may be misleading due to involvement in interactions
p <- (g1 | g2) / (g3 | g4)   # 2 x 2
print(p)

cat("\nPost hoc test se efeito de interação Alcool:Cafeina não é significante")

Post hoc test se efeito de interação Alcool:Cafeina não é significante
cat("\n\tAlcool")

    Alcool
EMM.A <- emmeans::emmeans(modelo, 
                          specs=pairwise~"Alcool", 
                          level=1-alfa,
                          adjust="holm",
                          lmer.df="satterthwaite",
                          lmerTest.limit=nrow(Dados))
NOTE: Results may be misleading due to involvement in interactions
print(summary(EMM.A$emmeans))
 Alcool emmean   SE df lower.CL upper.CL
 SemA     9.55 1.03 43     7.47     11.6
 ComA    19.21 1.01 43    17.17     21.2

Results are averaged over the levels of: Cafeina 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(summary(EMM.A$contrasts, 
              infer=TRUE))
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemA - ComA    -9.66 1.45 43    -12.6    -6.74  -6.681  <.0001

Results are averaged over the levels of: Cafeina 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(plot(EMM.A$emmeans, 
           colors="black"))

print(plot(EMM.A$contrasts, 
           colors="black"))

print(multcomp::cld(object=EMM.A$emmeans,
                    level=1-alfa,
                    alpha=alfa,
                    adjust="holm",
                    Letters=letters))
 Alcool emmean   SE df lower.CL upper.CL .group
 SemA     9.55 1.03 43     7.15     12.0  a    
 ComA    19.21 1.01 43    16.86     21.6   b   

Results are averaged over the levels of: Cafeina 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 2 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
cat("\n\tCafeina")

    Cafeina
EMM.C <- emmeans::emmeans(modelo, 
                          specs=pairwise~"Cafeina", 
                          level=1-alfa,
                          adjust="holm",
                          lmer.df="satterthwaite",
                          lmerTest.limit=nrow(Dados))
NOTE: Results may be misleading due to involvement in interactions
print(summary(EMM.C$emmeans))
 Cafeina emmean   SE df lower.CL upper.CL
 SemC      14.8 1.03 43     12.7     16.8
 ComC      14.0 1.01 43     12.0     16.0

Results are averaged over the levels of: Alcool 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(summary(EMM.C$contrasts, 
              infer=TRUE))
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemC - ComC    0.761 1.45 43    -2.15     3.68   0.527  0.6010

Results are averaged over the levels of: Alcool 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(plot(EMM.C$emmeans, 
           colors="black"))

print(plot(EMM.C$contrasts, 
           colors="black"))

print(multcomp::cld(object=EMM.C$emmeans,
                    level=1-alfa,
                    alpha=alfa,
                    adjust="holm",
                    Letters=letters))
 Cafeina emmean   SE df lower.CL upper.CL .group
 ComC      14.0 1.01 43     11.7     16.3  a    
 SemC      14.8 1.03 43     12.4     17.2  a    

Results are averaged over the levels of: Alcool 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 2 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
cat("\nPost hoc teste se efeito de interação Alcool:Cafeina é significante")

Post hoc teste se efeito de interação Alcool:Cafeina é significante
cat("\nPost hoc test: interaction Cafeina:Alcool")

Post hoc test: interaction Cafeina:Alcool
plot(emmeans::emmip(modelo, 
                    ~Cafeina:Alcool,
                    level=1-alfa,
                    adjust="holm",
                    CIs=TRUE))

EMM.AC <- emmeans::emmeans(modelo, 
                           specs=pairwise~"Alcool:Cafeina",
                           level=1-alfa,
                           adjust="holm",
                           lmer.df="satterthwaite",
                           lmerTest.limit=nrow(Dados))
print(summary(EMM.AC$emmeans))
 Alcool Cafeina emmean   SE df lower.CL upper.CL
 SemA   SemC      8.27 1.49 43     5.26     11.3
 ComA   SemC     21.25 1.43 43    18.37     24.1
 SemA   ComC     10.83 1.43 43     7.95     13.7
 ComA   ComC     17.17 1.43 43    14.28     20.0

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(summary(EMM.AC$contrasts, infer=TRUE))
 contrast              estimate   SE df lower.CL upper.CL t.ratio p.value
 SemA SemC - ComA SemC   -12.98 2.07 43   -18.69   -7.262  -6.280  <.0001
 SemA SemC - SemA ComC    -2.56 2.07 43    -8.28    3.155  -1.239  0.2220
 SemA SemC - ComA ComC    -8.89 2.07 43   -14.61   -3.179  -4.304  0.0004
 ComA SemC - SemA ComC    10.42 2.02 43     4.83   16.006   5.154  <.0001
 ComA SemC - ComA ComC     4.08 2.02 43    -1.51    9.673   2.020  0.0992
 SemA ComC - ComA ComC    -6.33 2.02 43   -11.92   -0.744  -3.134  0.0093

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 6 estimates 
P value adjustment: holm method for 6 tests 
print(plot(EMM.AC$emmeans, 
           colors="black"))

print(plot(EMM.AC$contrasts, 
           colors="black"))

print(multcomp::cld(object=EMM.AC$emmeans,
                    level=1-alfa,
                    alpha=alfa,
                    adjust="holm",
                    Letters=letters))
 Alcool Cafeina emmean   SE df lower.CL upper.CL .group
 SemA   SemC      8.27 1.49 43     4.38     12.2  a    
 SemA   ComC     10.83 1.43 43     7.11     14.6  a    
 ComA   ComC     17.17 1.43 43    13.44     20.9   b   
 ComA   SemC     21.25 1.43 43    17.52     25.0   b   

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 4 estimates 
P value adjustment: holm method for 6 tests 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
cat("\nSimple main effect test: Alcool effect by Cafeina levels\n")

Simple main effect test: Alcool effect by Cafeina levels
print(phia::testInteractions(modelo, 
                             across="Alcool",
                             fixed="Cafeina",
                             adjustment="holm")) 
Chisq Test: 
P-value adjustment method: holm
        Value     SE Df   Chisq Pr(>Chisq)    
SemC -12.9773 2.0664  1 39.4401  6.766e-10 ***
ComC  -6.3333 2.0210  1  9.8207   0.001726 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(emmeans::emmip(modelo,
                    ~Alcool|Cafeina,
                    level=1-alfa,
                    adjust="holm",
                    CIs=TRUE,
                    lmer.df="satterthwaite",
                    lmerTest.limit=nrow(Dados))) 

cat("\nPost hoc test of the simple main effect of Alcool")

Post hoc test of the simple main effect of Alcool
model_means <- emmeans::emmeans(object=modelo,
                                specs=pairwise~Alcool|Cafeina:Alcool,
                                level=1-alfa,
                                adjust="holm",
                                lmer.df="satterthwaite",
                                lmerTest.limit=nrow(Dados))
print(summary(model_means$emmeans, infer=TRUE))
Cafeina = SemC:
 Alcool emmean   SE df lower.CL upper.CL t.ratio p.value
 SemA     8.27 1.49 43     5.26     11.3   5.543  <.0001
 ComA    21.25 1.43 43    18.37     24.1  14.870  <.0001

Cafeina = ComC:
 Alcool emmean   SE df lower.CL upper.CL t.ratio p.value
 SemA    10.83 1.43 43     7.95     13.7   7.581  <.0001
 ComA    17.17 1.43 43    14.28     20.0  12.013  <.0001

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(summary(model_means$contrasts, infer=TRUE))
Cafeina = SemC:
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemA - ComA   -12.98 2.07 43    -17.1    -8.81  -6.280  <.0001

Cafeina = ComC:
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemA - ComA    -6.33 2.02 43    -10.4    -2.26  -3.134  0.0031

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
plot(model_means$emmeans,col="black")

plot(model_means$contrasts,col="black")

print(multcomp::cld(object=model_means$emmeans,
                    level=1-alfa,
                    alpha=alfa,
                    adjust="holm",
                    Letters=letters))
Cafeina = SemC:
 Alcool emmean   SE df lower.CL upper.CL .group
 SemA     8.27 1.49 43     4.81     11.7  a    
 ComA    21.25 1.43 43    17.93     24.6   b   

Cafeina = ComC:
 Alcool emmean   SE df lower.CL upper.CL .group
 SemA    10.83 1.43 43     7.51     14.2  a    
 ComA    17.17 1.43 43    13.85     20.5   b   

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 2 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same. 
cat("\nSimple main effect test: Cafeina effect by Alcool levels\n")

Simple main effect test: Cafeina effect by Alcool levels
print(phia::testInteractions(modelo, 
                             across="Cafeina",
                             fixed="Alcool",
                             adjustment="holm"))
Chisq Test: 
P-value adjustment method: holm
       Value     SE Df  Chisq Pr(>Chisq)  
SemA -2.5606 2.0664  1 1.5355    0.21529  
ComA  4.0833 2.0210  1 4.0823    0.08667 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(emmeans::emmip(modelo, 
                    ~Cafeina|Alcool,
                    level=1-alfa,
                    adjust="holm",
                    CIs=TRUE,
                    lmer.df="satterthwaite",
                    lmerTest.limit=nrow(Dados))) 

cat("\nPost hoc test of the simple main effect of Cafeina")

Post hoc test of the simple main effect of Cafeina
model_means <- emmeans::emmeans(object=modelo,
                                specs=pairwise~Cafeina|Cafeina:Alcool,
                                level=1-alfa,
                                adjust="holm",
                                lmer.df="satterthwaite",
                                lmerTest.limit=nrow(Dados))
print(summary(model_means$emmeans, infer=TRUE))
Alcool = SemA:
 Cafeina emmean   SE df lower.CL upper.CL t.ratio p.value
 SemC      8.27 1.49 43     5.26     11.3   5.543  <.0001
 ComC     10.83 1.43 43     7.95     13.7   7.581  <.0001

Alcool = ComA:
 Cafeina emmean   SE df lower.CL upper.CL t.ratio p.value
 SemC     21.25 1.43 43    18.37     24.1  14.870  <.0001
 ComC     17.17 1.43 43    14.28     20.0  12.013  <.0001

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
print(summary(model_means$contrasts, infer=TRUE))
Alcool = SemA:
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemC - ComC    -2.56 2.07 43 -6.72790     1.61  -1.239  0.2220

Alcool = ComA:
 contrast    estimate   SE df lower.CL upper.CL t.ratio p.value
 SemC - ComC     4.08 2.02 43  0.00764     8.16   2.020  0.0496

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
plot(model_means$emmeans,col="black")

plot(model_means$contrasts,col="black")

print(multcomp::cld(object=model_means$emmeans,
                    level=1-alfa,
                    alpha=alfa,
                    adjust="holm",
                    Letters=letters))
Alcool = SemA:
 Cafeina emmean   SE df lower.CL upper.CL .group
 SemC      8.27 1.49 43     4.81     11.7  a    
 ComC     10.83 1.43 43     7.51     14.2  a    

Alcool = ComA:
 Cafeina emmean   SE df lower.CL upper.CL .group
 ComC     17.17 1.43 43    13.85     20.5  a    
 SemC     21.25 1.43 43    17.93     24.6   b   

Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 2 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same.

Exemplo de ANOVA multifatorial e GLMM em Takayanagi, Siqueira, Silveira, Valentova, 2024

Título: O que Diferentes Pessoas Buscam em um Parceiro? Efeitos do Sexo, da Orientação Sexual e das Estratégias de Acasalamento nas Preferências de Parceiro

A artigo apresenta a análise estatística GLMM com alguns fatores entre e intraparticipantes com três ou mais níveis.

Os dados e código R estão disponíveis em OSF Home.

Resumo: As preferências de parceiro são um diferencial importante na formação de relacionamentos e na aptidão evolutiva, variando de acordo com fatores individuais, ecológicos e sociais. Neste estudo, avaliamos a variação na preferência por inteligência, bondade, atratividade física, saúde e nível socioeconômico entre indivíduos de diferentes sexos e orientações sexuais em uma amostra brasileira. Analisamos as pontuações de preferência de 778 homens e mulheres heterossexuais, bissexuais e homossexuais em três tarefas de orçamento de parceiro (baixo vs. médio vs. alto orçamento) e sua associação com sociosexualidade, estilos de apego, homogamia e disposição para engajar-se em relacionamentos de curto e longo prazo. Os resultados indicaram uma ordem global de preferência por traços, com a inteligência em primeiro lugar, seguida por bondade, atratividade física, saúde e, por último, status socioeconômico. Diferenças típicas de sexo foram observadas principalmente no grupo heterossexual, e combinações específicas de sexo e orientação sexual estiveram associadas a variações na preferência por atratividade física, bondade e status socioeconômico. Também encontramos associações únicas de outras variáveis com preferências de parceiro e com a disposição para engajar-se em relacionamentos de curto ou longo prazo. Ao explorar as preferências de parceiro de indivíduos não heterossexuais de um país latino-americano, um grupo sub-representado na pesquisa em psicologia evolucionista, nossos resultados ajudam a compreender os fatores universais e específicos que orientam as preferências de parceiro e o comportamento sexual humano.

Análises Estatísticas: Os dados foram analisados no software R, versão 4.2.1 (Core Team, 2022). Computamos modelos lineares mistos gerais (GLMM) utilizando a função lmer do pacote lmerTest (Kuznetsova et al., 2017).

Na análise principal, utilizamos um GLMM para avaliar como o nível de preferência por cada traço variou em função das características individuais. Usamos a pontuação de preferência pelo traço como variável dependente; o identificador do participante como efeito aleatório; traço do parceiro, sexo e orientação sexual como variáveis fixas nominais; e tamanho do orçamento, sociosexualidade, apego ansioso, apego evitativo e pontuação de autoavaliação do traço como variáveis fixas intervalares. Idade, estado civil, renda mensal domiciliar per capita e nível educacional foram incluídos como variáveis de controle.

Nas análises adicionais, comparamos as pontuações de disposição para engajar-se em relacionamentos de curto e longo prazo entre condições de orçamento, sexo e orientações sexuais, utilizando uma análise de variância (ANOVA) de medidas repetidas com quatro fatores. Também usamos GLMM para explorar como características individuais e diferenças nas preferências de parceiro influenciaram a inclinação para formar relacionamentos de curto ou longo prazo com o parceiro hipotético. Usamos o nível de disposição para engajar-se em relacionamentos de curto ou longo prazo como variável dependente de cada modelo; o identificador do participante como efeito aleatório; sexo e orientação sexual como variáveis fixas nominais; e tamanho do orçamento, sociosexualidade, apego ansioso, apego evitativo e pontuações de preferência por traços como variáveis fixas intervalares. Mais uma vez, idade, estado civil, renda mensal domiciliar per capita e nível educacional foram incluídos como variáveis de controle.

Referências

  • Dancey CP & Reidy J (2006) Estatística sem Matemática para Psicologia. 3a ed. Porto Alegre: Bookman.
  • Dancey CP & Reidy J (2019) Estatística sem Matemática para Psicologia. 7a ed. Porto Alegre: PENSO.
  • Howell GT & Lacroix GL (2012) Decomposing interactions using GLM in combination with the COMPARE, LMATRIX and MMATRIX subcommands in SPSS. Tutorials in Quantitative Methods for Psychology 8(1):1-22. DOI:10.20982/tqmp.08.1.p001
  • Howland J et al. (2011) The acute effects of caffeinated versus non-caffeinated alcoholic beverage on driving performance and attention/reaction time. Addiction. 106(2):335-41. doi:10.1111/j.1360-0443.2010.03219.x
  • Takayanagi, JFGB, Siqueira, JO, Silveira, PSP, Valentova, JV (2024). What Do Different People Look for in a Partner? Effects of Sex, Sexual Orientation, and Mating Strategies on Partner Preferences. Arch Sex Behav 53, 981–1000 (2024). https://doi.org/10.1007/s10508-023-02767-4. Dados e código R: https://osf.io/pc6rn/?view%2520only=38097a82dd274d34bbddc8dcc2258bb9