Prática R: ANOVA Bifatorial

MSP4112: Estatística Aplicada à Medicina

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(heplots, warn.conflicts=FALSE))
suppressMessages(library(lawstat, warn.conflicts=FALSE))
suppressMessages(library(lmerTest, warn.conflicts=FALSE))
suppressMessages(library(multcomp, warn.conflicts=FALSE))
suppressMessages(library(MuMIn, 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(sjPlot, warn.conflicts=FALSE))
suppressMessages(library(twowaytests, warn.conflicts=FALSE))
suppressMessages(library(WebPower, warn.conflicts=FALSE))
source("summarySEwithin2.R")

Material

• HTML de Rmarkdown em RPubs

Introdução

Esta prática é relativa aos capítulos 11 do livro-texto:

• Dancey C & Reidy J (2019) Estatística sem Matemática para Psicologia. 7^a ed. Porto Alegre: PENSO.

Diferentemente da ANOVA unifatorial, quando há dois ou mais fatores, três situações possíveis são:

• todos os fatores entre participantes
• todos os fatores intraparticipantes
• fatores intra e entre participantes

Two-way ANOVA: Dois fatores entre participantes

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.

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 na planilha CafEntre_AlcEntre.xlsx. Foram obtidos e adaptados de Dancey CP & Reidy J (2006) Estatística sem Matemática para Psicologia 3^a ed. Porto Alegre: Bookman.

O conteúdo da planilha é apresentado a seguir.

Dados.long <- readxl::read_xlsx("CafEntre_AlcEntre.xlsx")
Dados.long$UE <- factor(Dados.long$UE)
Dados.long$Cafeina <- factor(Dados.long$Cafeina, 
                             levels=c("SemC","ComC"))
Dados.long$Alcool <- factor(Dados.long$Alcool, 
                             levels=c("SemA","ComA"))
print.data.frame(Dados.long, row.names=FALSE)

 UE Cafeina Alcool NumErros
  1    SemC   SemA        4
  2    SemC   SemA        9
  3    SemC   SemA       10
  4    SemC   SemA        8
  5    SemC   SemA        6
  6    SemC   SemA       11
  7    SemC   SemA        2
  8    SemC   SemA       11
  9    SemC   SemA       11
 10    SemC   SemA       10
 11    SemC   SemA        3
 12    SemC   SemA       10
 13    SemC   ComA       28
 14    SemC   ComA       22
 15    SemC   ComA       21
 16    SemC   ComA       27
 17    SemC   ComA       21
 18    SemC   ComA       20
 19    SemC   ComA       19
 20    SemC   ComA       16
 21    SemC   ComA       25
 22    SemC   ComA       17
 23    SemC   ComA       19
 24    SemC   ComA       20
 25    ComC   SemA       15
 26    ComC   SemA        8
 27    ComC   SemA       17
 28    ComC   SemA        0
 29    ComC   SemA       15
 30    ComC   SemA       11
 31    ComC   SemA       11
 32    ComC   SemA        6
 33    ComC   SemA        0
 34    ComC   SemA       15
 35    ComC   SemA       17
 36    ComC   SemA       15
 37    ComC   ComA       10
 38    ComC   ComA       11
 39    ComC   ComA       27
 40    ComC   ComA       15
 41    ComC   ComA       27
 42    ComC   ComA       10
 43    ComC   ComA       21
 44    ComC   ComA       15
 45    ComC   ComA       19
 46    ComC   ComA       21
 47    ComC   ComA       15
 48    ComC   ComA       15

print(xtabs(~Alcool+Cafeina, data=Dados.long))

      Cafeina
Alcool SemC ComC
  SemA   12   12
  ComA   12   12

str(Dados.long)

tibble [48 × 4] (S3: tbl_df/tbl/data.frame)
 $ 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 [1:48] 4 9 10 8 6 11 2 11 11 10 ...

Este é um formato long, no qual cada observação da VD de cada unidade experimental (UE) aparece em uma linha, cada fator aparece em uma coluna e a VD em uma coluna. A identificação de cada unidade experimental, que no caso é um participante, é o elemento-chave que distingue os delineamentos. Neste caso, há 48 participantes.

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 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, há quatro efeitos fixos:

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

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

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

Testar os efeitos principais e de interação dos fatores álcool e cafeína sobre o número erros no simulador se:

Questão 1

• Os fatores Alcool e Cafeina são entre participantes;
  Os dados estão na planilha CafEntre_AlcEntre.xlsx

Solução em alcb_cafb.R

Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8")

[1] "pt_BR.UTF-8"

Sys.setlocale("LC_ALL", "pt_BR.UTF-8")

[1] "LC_COLLATE=pt_BR.UTF-8;LC_CTYPE=pt_BR.UTF-8;LC_MONETARY=pt_BR.UTF-8;LC_NUMERIC=C;LC_TIME=pt_BR.UTF-8"

alfa <- 0.05
Dados.long <- readxl::read_xlsx("CafEntre_AlcEntre.xlsx")
Dados.long$UE <- factor(Dados.long$UE)
Dados.long$Cafeina <- factor(Dados.long$Cafeina, 
                             levels=c("SemC","ComC"))
Dados.long$Alcool <- factor(Dados.long$Alcool, 
                             levels=c("SemA","ComA"))
cat("\nDados.long")

Dados.long

print.data.frame(Dados.long, row.names=FALSE)

 UE Cafeina Alcool NumErros
  1    SemC   SemA        4
  2    SemC   SemA        9
  3    SemC   SemA       10
  4    SemC   SemA        8
  5    SemC   SemA        6
  6    SemC   SemA       11
  7    SemC   SemA        2
  8    SemC   SemA       11
  9    SemC   SemA       11
 10    SemC   SemA       10
 11    SemC   SemA        3
 12    SemC   SemA       10
 13    SemC   ComA       28
 14    SemC   ComA       22
 15    SemC   ComA       21
 16    SemC   ComA       27
 17    SemC   ComA       21
 18    SemC   ComA       20
 19    SemC   ComA       19
 20    SemC   ComA       16
 21    SemC   ComA       25
 22    SemC   ComA       17
 23    SemC   ComA       19
 24    SemC   ComA       20
 25    ComC   SemA       15
 26    ComC   SemA        8
 27    ComC   SemA       17
 28    ComC   SemA        0
 29    ComC   SemA       15
 30    ComC   SemA       11
 31    ComC   SemA       11
 32    ComC   SemA        6
 33    ComC   SemA        0
 34    ComC   SemA       15
 35    ComC   SemA       17
 36    ComC   SemA       15
 37    ComC   ComA       10
 38    ComC   ComA       11
 39    ComC   ComA       27
 40    ComC   ComA       15
 41    ComC   ComA       27
 42    ComC   ComA       10
 43    ComC   ComA       21
 44    ComC   ComA       15
 45    ComC   ComA       19
 46    ComC   ComA       21
 47    ComC   ComA       15
 48    ComC   ComA       15

print(xtabs(~Alcool+Cafeina, data=Dados.long))

      Cafeina
Alcool SemC ComC
  SemA   12   12
  ComA   12   12

cat("\nStructure:\n")

Structure:

str(Dados.long)

tibble [48 × 4] (S3: tbl_df/tbl/data.frame)
 $ 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 [1:48] 4 9 10 8 6 11 2 11 11 10 ...

cat("Analisando número de erros no simulador")

Analisando número de erros no simulador

boxplot(NumErros~Cafeina*Alcool, 
        data=Dados.long,
        ylab="NumErros")

alfabonf <- alfa/(length(unique(Dados.long$Alcool))*
                    length(unique(Dados.long$Cafeina)))

print.data.frame(psych::describeBy(Dados.long$NumErros,
                                   group=list(Dados.long$Alcool,
                                              Dados.long$Cafeina), 
                                   mat=TRUE, 
                                   digits=2), 
                 row.names=FALSE)

 item group1 group2 vars  n  mean   sd median trimmed  mad min max range  skew
    1   SemA   SemC    1 12  7.92 3.32    9.5     8.2 2.22   2  11     9 -0.63
    2   ComA   SemC    1 12 21.25 3.72   20.5    21.1 2.22  16  28    12  0.49
    3   SemA   ComC    1 12 10.83 6.12   13.0    11.3 4.45   0  17    17 -0.70
    4   ComA   ComC    1 12 17.17 5.92   15.0    16.9 6.67  10  27    17  0.41
 kurtosis   se
    -1.35 0.96
    -1.04 1.07
    -1.08 1.77
    -1.24 1.71

rmiscsum <- Rmisc::summarySE(Dados.long,
                             measurevar="NumErros",
                             groupvars=c("Alcool","Cafeina"),
                             na.rm=TRUE, 
                             conf.interval=1-alfabonf)
print(rmiscsum, digits=3)

  Alcool Cafeina  N NumErros   sd    se   ci
1   SemA    SemC 12     7.92 3.32 0.957 2.85
2   SemA    ComC 12    10.83 6.12 1.766 5.26
3   ComA    SemC 12    21.25 3.72 1.074 3.20
4   ComA    ComC 12    17.17 5.92 1.709 5.10

cat("\nIntervalo de confiança ")

Intervalo de confiança 

cat("\ncom correção de Bonferroni (entre participantes)")

com correção de Bonferroni (entre participantes)

grf <- ggplot2::ggplot(rmiscsum, 
                       ggplot2::aes(x=Cafeina, 
                                    y=NumErros,
                                    colour=Alcool,
                                    group=Alcool)) +
  ggplot2::geom_line() +
  ggplot2::geom_errorbar(width=.3, 
                         ggplot2::aes(ymin=NumErros-ci, 
                                      ymax=NumErros+ci)) +
  ggplot2::geom_point(shape=21, 
                      size=3, 
                      fill="white") +
  # ggplot2::ylim(2.5,6.5) +
  ggplot2::ylab("NumErros") +
  ggplot2::ggtitle("Between-subjects CI95 Bonferroni") +
  ggplot2::theme_bw() 
print(grf)

cat("\nANOVA unifatorial independente (aov): perfeitamente balanceado: Faraway(2016)")

ANOVA unifatorial independente (aov): perfeitamente balanceado: Faraway(2016)

fit.aov <- aov(NumErros ~ Alcool + Cafeina + Alcool:Cafeina, 
               data=Dados.long)
print(summary(fit.aov))

               Df Sum Sq Mean Sq F value   Pr(>F)    
Alcool          1 1160.3  1160.3  47.692 1.57e-08 ***
Cafeina         1    4.1     4.1   0.168    0.684    
Alcool:Cafeina  1  147.0   147.0   6.042    0.018 *  
Residuals      44 1070.5    24.3                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\nANOVA bifatorial independente (lm)")

ANOVA bifatorial independente (lm)

cat("\n\tteste omnibus")

    teste omnibus

fit.lm <- lm(NumErros ~ Alcool + Cafeina + Alcool:Cafeina,
             data=Dados.long)
print(car::Anova(fit.lm))

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(summary(fit.lm))

Call:
lm(formula = NumErros ~ Alcool + Cafeina + Alcool:Cafeina, data = Dados.long)

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

print(sjPlot::tab_model(fit.lm,
                  p.adjust="holm",
                  encoding="Windows-1252",
                  df.method="satterthwaite",
                  title="p.adjust=holm",
                  show.reflvl=TRUE,
                  show.ngroups=TRUE,
                  show.stat=TRUE,
                  show.df=TRUE,
                  show.se=TRUE,
                  digits=4,
                  digits.p=4,
                  digits.rsq=6,
                  digits.re=6))

`ci_method` must be one of "residual", "wald", "normal", "profile",
  "boot", "uniroot", "betwithin" or "ml1". Using "wald" now.

options(warn=-1)
o.par <- par()
fit.means <- phia::interactionMeans(fit.lm)
plot(fit.means, 
     errorbar=paste0("ci",round((1-alfa)*100,0)))
par(o.par)

options(warn=0)

# Se ausência de ausência de efeito de interação não é rejeitado
# https://schmidtpaul.github.io/DSFAIR/compactletterdisplay.html
model_means <- emmeans::emmeans(object=fit.lm,
                                specs=pairwise~Alcool,
                                level=1-alfa,
                                adjust="holm")

NOTE: Results may be misleading due to involvement in interactions

print(summary(model_means$emmeans, infer=TRUE))

 Alcool emmean   SE df lower.CL upper.CL t.ratio p.value
 SemA     9.38 1.01 44     7.35     11.4   9.311  <.0001
 ComA    19.21 1.01 44    17.18     21.2  19.078  <.0001

Results are averaged over the levels of: Cafeina 
Confidence level used: 0.95 

print(summary(model_means$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 

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 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. 

print(emmeans::eff_size(object=model_means$emmeans, 
                        sigma=sigma(fit.lm), 
                        edf=df.residual(fit.lm)))

 contrast    effect.size    SE df lower.CL upper.CL
 SemA - ComA       -1.99 0.358 44    -2.72    -1.27

Results are averaged over the levels of: Cafeina 
sigma used for effect sizes: 4.932 
Confidence level used: 0.95 

model_means <- emmeans::emmeans(object=fit.lm,
                                specs=pairwise~Cafeina,
                                level=1-alfa,
                                adjust="holm")

NOTE: Results may be misleading due to involvement in interactions

print(summary(model_means$emmeans, infer=TRUE))

 Cafeina emmean   SE df lower.CL upper.CL t.ratio p.value
 SemC      14.6 1.01 44     12.6     16.6  14.484  <.0001
 ComC      14.0 1.01 44     12.0     16.0  13.905  <.0001

Results are averaged over the levels of: Alcool 
Confidence level used: 0.95 

print(summary(model_means$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 

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 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. 

print(emmeans::eff_size(object=model_means$emmeans, 
                        sigma=sigma(fit.lm), 
                        edf=df.residual(fit.lm)))

 contrast    effect.size    SE df lower.CL upper.CL
 SemC - ComC       0.118 0.289 44   -0.464    0.701

Results are averaged over the levels of: Alcool 
sigma used for effect sizes: 4.932 
Confidence level used: 0.95 

model_means <- emmeans::emmeans(object=fit.lm,
                                specs=pairwise~Cafeina:Alcool,
                                level=1-alfa,
                                adjust="holm")
print(summary(model_means$emmeans, infer=TRUE))

 Cafeina Alcool emmean   SE df lower.CL upper.CL t.ratio p.value
 SemC    SemA     7.92 1.42 44     5.05     10.8   5.560  <.0001
 ComC    SemA    10.83 1.42 44     7.96     13.7   7.608  <.0001
 SemC    ComA    21.25 1.42 44    18.38     24.1  14.924  <.0001
 ComC    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))

 contrast              estimate   SE df lower.CL upper.CL t.ratio p.value
 SemC SemA - ComC SemA    -2.92 2.01 44    -8.48     2.65  -1.448  0.1546
 SemC SemA - SemC ComA   -13.33 2.01 44   -18.90    -7.77  -6.621  <.0001
 SemC SemA - ComC ComA    -9.25 2.01 44   -14.81    -3.69  -4.594  0.0001
 ComC SemA - SemC ComA   -10.42 2.01 44   -15.98    -4.85  -5.173  <.0001
 ComC SemA - ComC ComA    -6.33 2.01 44   -11.90    -0.77  -3.145  0.0089
 SemC ComA - ComC ComA     4.08 2.01 44    -1.48     9.65   2.028  0.0973

Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 6 estimates 
P value adjustment: holm method for 6 tests 

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 Alcool emmean   SE df lower.CL upper.CL .group
 SemC    SemA     7.92 1.42 44     4.21     11.6  a    
 ComC    SemA    10.83 1.42 44     7.12     14.5  a    
 ComC    ComA    17.17 1.42 44    13.46     20.9   b   
 SemC    ComA    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. 

print(emmeans::eff_size(object=model_means$emmeans, 
                        sigma=sigma(fit.lm), 
                        edf=df.residual(fit.lm)))

 contrast              effect.size    SE df lower.CL upper.CL
 SemC SemA - ComC SemA      -0.591 0.413 44  -1.4238    0.241
 SemC SemA - SemC ComA      -2.703 0.500 44  -3.7102   -1.696
 SemC SemA - ComC ComA      -1.875 0.455 44  -2.7914   -0.959
 ComC SemA - SemC ComA      -2.112 0.466 44  -3.0514   -1.172
 ComC SemA - ComC ComA      -1.284 0.431 44  -2.1518   -0.416
 SemC ComA - ComC ComA       0.828 0.418 44  -0.0139    1.670

sigma used for effect sizes: 4.932 
Confidence level used: 0.95 

# Se ausência de ausência de efeito de interação é rejeitado
cat("\nTeste de efeito principal simples de Álcool")

Teste de efeito principal simples de Álcool

print(phia::testInteractions(fit.lm, 
                             fixed="Cafeina", 
                             across="Alcool",
                             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

cat("\nTeste post hoc de efeito principal simples de Álcool")

Teste post hoc de efeito principal simples de Álcool

model_means <- emmeans::emmeans(object=fit.lm,
                                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("\nTeste de efeito principal simples de Cafeína")

Teste de efeito principal simples de Cafeína

print(phia::testInteractions(fit.lm, 
                             fixed="Alcool", 
                             across="Cafeina",
                             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

cat("\nTeste post hoc de efeito principal simples de Cafeína")

Teste post hoc de efeito principal simples de Cafeína

model_means <- emmeans::emmeans(object=fit.lm,
                                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. 

Questão 2

• Os fatores Alcool e Cafeina são intraparticipantes;
  Os dados estão na planilha CafIntra_AlcIntra.xlsx

Solução em alcw_cafw.R

source("summarySEwithin2.R")
Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8")

[1] "pt_BR.UTF-8"

Sys.setlocale("LC_ALL", "pt_BR.UTF-8")

[1] "LC_COLLATE=pt_BR.UTF-8;LC_CTYPE=pt_BR.UTF-8;LC_MONETARY=pt_BR.UTF-8;LC_NUMERIC=C;LC_TIME=pt_BR.UTF-8"

alfa <- 0.05
Dados.long <- readxl::read_xlsx("CafIntra_AlcIntra.xlsx")
Dados.long$UE <- factor(Dados.long$UE)
Dados.long$Cafeina <- factor(Dados.long$Cafeina, 
                             levels=c("SemC","ComC"))
Dados.long$Alcool <- factor(Dados.long$Alcool, 
                             levels=c("SemA","ComA"))
cat("\nDados.long")

Dados.long

print.data.frame(Dados.long, row.names=FALSE)

 UE Alcool Cafeina NumErros
  1   SemA    SemC        4
  1   SemA    ComC       15
  1   ComA    SemC       28
  1   ComA    ComC       10
  2   SemA    SemC        9
  2   SemA    ComC        8
  2   ComA    SemC       22
  2   ComA    ComC       11
  3   SemA    SemC       10
  3   SemA    ComC       17
  3   ComA    SemC       21
  3   ComA    ComC       27
  4   SemA    SemC        8
  4   SemA    ComC        0
  4   ComA    SemC       27
  4   ComA    ComC       15
  5   SemA    SemC        6
  5   SemA    ComC       15
  5   ComA    SemC       21
  5   ComA    ComC       27
  6   SemA    SemC       11
  6   SemA    ComC       11
  6   ComA    SemC       20
  6   ComA    ComC       10
  7   SemA    SemC        2
  7   SemA    ComC       11
  7   ComA    SemC       19
  7   ComA    ComC       21
  8   SemA    SemC       11
  8   SemA    ComC        6
  8   ComA    SemC       16
  8   ComA    ComC       15
  9   SemA    SemC       11
  9   SemA    ComC        0
  9   ComA    SemC       25
  9   ComA    ComC       19
 10   SemA    SemC       10
 10   SemA    ComC       15
 10   ComA    SemC       17
 10   ComA    ComC       21
 11   SemA    SemC        3
 11   SemA    ComC       17
 11   ComA    SemC       19
 11   ComA    ComC       15
 12   SemA    SemC       10
 12   SemA    ComC       15
 12   ComA    SemC       20
 12   ComA    ComC       15

print(xtabs(~Alcool+Cafeina+UE, data=Dados.long))

, , UE = 1

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    1    1

, , UE = 2

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    1    1

, , UE = 3

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    1    1

, , UE = 4

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    1    1

, , UE = 5

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    1    1

, , UE = 6

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    1    1

, , UE = 7

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    1    1

, , UE = 8

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    1    1

, , UE = 9

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    1    1

, , UE = 10

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    1    1

, , UE = 11

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    1    1

, , UE = 12

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    1    1

cat("\nStructure:\n")

Structure:

str(Dados.long)

tibble [48 × 4] (S3: tbl_df/tbl/data.frame)
 $ 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 [1:48] 4 15 28 10 9 8 22 11 10 17 ...

cat("Analisando número de erros no simulador")

Analisando número de erros no simulador

alfabonf <- alfa/(length(unique(Dados.long$Alcool))*
                    length(unique(Dados.long$Cafeina)))

rmiscsum <- summarySEwithin2(Dados.long,
                             measurevar="NumErros",
                             withinvars=c("Alcool","Cafeina"),
                             idvar="UE",
                             na.rm=TRUE, 
                             conf.interval=1-alfabonf)
print(rmiscsum)

  Alcool Cafeina  N  NumErros NumErrosNormed       sd       se       ci
1   ComA    ComC 12 17.166667      17.166667 5.233198 1.510694 4.503187
2   ComA    SemC 12 21.250000      21.250000 5.034477 1.453328 4.332186
3   SemA    ComC 12 10.833333      10.833333 5.914586 1.707394 5.089523
4   SemA    SemC 12  7.916667       7.916667 4.409872 1.273020 3.794711

cat("\nIntervalo de confiança ")

Intervalo de confiança 

cat("\ncom correção de Bonferroni (intraparticipantes)")

com correção de Bonferroni (intraparticipantes)

grf <- ggplot2::ggplot(rmiscsum, 
                       ggplot2::aes(x=Cafeina, 
                                    y=NumErros,
                                    colour=Alcool,
                                    group=Alcool)) +
  ggplot2::geom_line() +
  ggplot2::geom_errorbar(width=.3, 
                         ggplot2::aes(ymin=NumErros-ci, 
                                      ymax=NumErros+ci)) +
  ggplot2::geom_point(shape=21, 
                      size=3, 
                      fill="white") +
  ggplot2::ylab("NumErros") +
  ggplot2::ggtitle("Within-subjects CI95 Bonferroni") +
  ggplot2::theme_bw() 
print(grf)

cat("\nANOVA bifatorial relacionada (aov): balanceamento perfeito: Faraway(2016)")

ANOVA bifatorial relacionada (aov): balanceamento perfeito: Faraway(2016)

fit.aov <- aov(NumErros ~ Alcool + Cafeina + Alcool:Cafeina +
               Error(UE), 
               data=Dados.long)
print(summary(fit.aov))

Error: UE
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 11  186.4   16.95               

Error: Within
               Df Sum Sq Mean Sq F value   Pr(>F)    
Alcool          1 1160.3  1160.3  43.312 1.76e-07 ***
Cafeina         1    4.1     4.1   0.152   0.6987    
Alcool:Cafeina  1  147.0   147.0   5.487   0.0253 *  
Residuals      33  884.1    26.8                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\nANOVA bifatorial relacionada (lmerTest::lmer)")

ANOVA bifatorial relacionada (lmerTest::lmer)

cat("\n\tteste omnibus")

    teste omnibus

fit.lmer <- lmerTest::lmer(NumErros ~ Alcool + Cafeina + Alcool:Cafeina +
                                      (1|UE),
                           data=Dados.long)

boundary (singular) fit: see help('isSingular')

print(car::Anova(fit.lmer))

Analysis of Deviance Table (Type II Wald chisquare tests)

Response: NumErros
                 Chisq Df Pr(>Chisq)    
Alcool         47.6924  1  4.986e-12 ***
Cafeina         0.1678  1    0.68204    
Alcool:Cafeina  6.0420  1    0.01397 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

print(summary(fit.lmer), correl = FALSE)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: NumErros ~ Alcool + Cafeina + Alcool:Cafeina + (1 | UE)
   Data: Dados.long

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

print(sjPlot::tab_model(fit.lmer,
                        p.adjust="holm",
                        encoding="Windows-1252",
                        df.method="satterthwaite",
                        title="p.adjust=holm",
                        show.reflvl=TRUE,
                        show.ngroups=TRUE,
                        show.stat=TRUE,
                        show.df=TRUE,
                        show.se=TRUE,
                        digits=4,
                        digits.p=4,
                        digits.rsq=6,
                        digits.re=6))
print(effectsize::eta_squared(fit.lmer, partial=TRUE))

# Effect Size for ANOVA (Type III)

Parameter      | Eta2 (partial) |       95% CI
----------------------------------------------
Alcool         |           0.52 | [0.35, 1.00]
Cafeina        |       3.80e-03 | [0.00, 1.00]
Alcool:Cafeina |           0.12 | [0.01, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

print(as.numeric(MuMIn::r.squaredGLMM(fit.lmer)))

[1] 0.5342025 0.5342025

options(warn=-1)
o.par <- par()
fit.means <- phia::interactionMeans(fit.lmer)
plot(fit.means, 
     errorbar=paste0("ci",round((1-alfa)*100,0)))
par(o.par)

options(warn=0)

# Se ausência de efeito de interação não é rejeitado
# https://schmidtpaul.github.io/DSFAIR/compactletterdisplay.html
model_means <- emmeans::emmeans(object=fit.lmer,
                                specs=pairwise~Alcool,
                                level=1-alfa,
                                adjust="holm",
                                lmer.df="satterthwaite",
                                lmerTest.limit=nrow(Dados.long))

NOTE: Results may be misleading due to involvement in interactions

print(summary(model_means$emmeans, infer=TRUE))

 Alcool emmean   SE df lower.CL upper.CL t.ratio p.value
 SemA     9.38 1.01 44     7.35     11.4   9.311  <.0001
 ComA    19.21 1.01 44    17.18     21.2  19.078  <.0001

Results are averaged over the levels of: Cafeina 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 

print(summary(model_means$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 

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

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

print(multcomp::cld(object=model_means$contrasts,
                    adjust="holm",
                    Letters=letters,
                    alpha=alfa))

 contrast    estimate   SE df t.ratio p.value .group
 SemA - ComA    -9.83 1.42 44  -6.906  <.0001  a    

Results are averaged over the levels of: Cafeina 
Degrees-of-freedom method: satterthwaite 
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. 

model_means <- emmeans::emmeans(object=fit.lmer,
                                specs=pairwise~Cafeina,
                                level=1-alfa,
                                adjust="holm",
                                lmer.df="satterthwaite",
                                lmerTest.limit=nrow(Dados.long))

NOTE: Results may be misleading due to involvement in interactions

print(summary(model_means$emmeans, infer=TRUE))

 Cafeina emmean   SE df lower.CL upper.CL t.ratio p.value
 SemC      14.6 1.01 44     12.6     16.6  14.484  <.0001
 ComC      14.0 1.01 44     12.0     16.0  13.905  <.0001

Results are averaged over the levels of: Alcool 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 

print(summary(model_means$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 

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

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

print(multcomp::cld(object=model_means$contrasts,
                    adjust="holm",
                    Letters=letters,
                    alpha=alfa))

 contrast    estimate   SE df t.ratio p.value .group
 SemC - ComC    0.583 1.42 44   0.410  0.6840  a    

Results are averaged over the levels of: Alcool 
Degrees-of-freedom method: satterthwaite 
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. 

model_means <- emmeans::emmeans(object=fit.lmer,
                                specs=pairwise~Cafeina:Alcool,
                                level=1-alfa,
                                adjust="holm",
                                lmer.df="satterthwaite",
                                lmerTest.limit=nrow(Dados.long))
print(summary(model_means$emmeans, infer=TRUE))

 Cafeina Alcool emmean   SE df lower.CL upper.CL t.ratio p.value
 SemC    SemA     7.92 1.42 44     5.05     10.8   5.560  <.0001
 ComC    SemA    10.83 1.42 44     7.96     13.7   7.608  <.0001
 SemC    ComA    21.25 1.42 44    18.38     24.1  14.924  <.0001
 ComC    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))

 contrast              estimate   SE df lower.CL upper.CL t.ratio p.value
 SemC SemA - ComC SemA    -2.92 2.01 44    -8.48     2.65  -1.448  0.1546
 SemC SemA - SemC ComA   -13.33 2.01 44   -18.90    -7.77  -6.621  <.0001
 SemC SemA - ComC ComA    -9.25 2.01 44   -14.81    -3.69  -4.594  0.0001
 ComC SemA - SemC ComA   -10.42 2.01 44   -15.98    -4.85  -5.173  <.0001
 ComC SemA - ComC ComA    -6.33 2.01 44   -11.90    -0.77  -3.145  0.0089
 SemC ComA - ComC ComA     4.08 2.01 44    -1.48     9.65   2.028  0.0973

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 

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

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

print(multcomp::cld(object=model_means$contrasts,
                    adjust="holm",
                    Letters=letters,
                    alpha=alfa))

 contrast              estimate   SE df t.ratio p.value .group
 SemC SemA - SemC ComA   -13.33 2.01 44  -6.621  <.0001  a    
 ComC SemA - SemC ComA   -10.42 2.01 44  -5.173  <.0001  ab   
 SemC SemA - ComC ComA    -9.25 2.01 44  -4.594  0.0001  a    
 ComC SemA - ComC ComA    -6.33 2.01 44  -3.145  0.0089  ab   
 SemC SemA - ComC SemA    -2.92 2.01 44  -1.448  0.1546   bc  
 SemC ComA - ComC ComA     4.08 2.01 44   2.028  0.0973    c  

Degrees-of-freedom method: satterthwaite 
P value adjustment: holm method for 6 tests 
P value adjustment: holm method for 15 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. 

# Se ausência de efeito de interação é rejeitado
cat("\nTeste de efeito principal simples de Álcool")

Teste de efeito principal simples de Álcool

print(phia::testInteractions(fit.lmer, 
                             fixed="Cafeina", 
                             across="Alcool",
                             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

cat("\nTeste post hoc de efeito principal simples de Álcool")

Teste post hoc de efeito principal simples de Álcool

model_means <- emmeans::emmeans(object=fit.lmer,
                                specs=pairwise~Alcool|Cafeina:Alcool,
                                level=1-alfa,
                                adjust="holm",
                                lmer.df="satterthwaite",
                                lmerTest.limit=nrow(Dados.long))
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$contrasts,
                    adjust="holm",
                    Letters=letters,
                    alpha=alfa))

Cafeina = SemC:
 contrast    estimate   SE df t.ratio p.value .group
 SemA - ComA   -13.33 2.01 44  -6.621  <.0001  a    

Cafeina = ComC:
 contrast    estimate   SE df t.ratio p.value .group
 SemA - ComA    -6.33 2.01 44  -3.145  0.0030  a    

Degrees-of-freedom method: satterthwaite 
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("\nTeste de efeito principal simples de Cafeína")

Teste de efeito principal simples de Cafeína

print(phia::testInteractions(fit.lmer, 
                             fixed="Alcool", 
                             across="Cafeina",
                             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

cat("\nTeste post hoc de efeito principal simples de Cafeína")

Teste post hoc de efeito principal simples de Cafeína

model_means <- emmeans::emmeans(object=fit.lmer,
                                specs=pairwise~Cafeina|Cafeina:Alcool,
                                level=1-alfa,
                                adjust="holm",
                                lmer.df="satterthwaite",
                                lmerTest.limit=nrow(Dados.long))
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$contrasts,
                    adjust="holm",
                    Letters=letters,
                    alpha=alfa))

Alcool = SemA:
 contrast    estimate   SE df t.ratio p.value .group
 SemC - ComC    -2.92 2.01 44  -1.448  0.1546  a    

Alcool = ComA:
 contrast    estimate   SE df t.ratio p.value .group
 SemC - ComC     4.08 2.01 44   2.028  0.0487  a    

Degrees-of-freedom method: satterthwaite 
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. 

Questão 3

• Os fatores Cafeina e Alcool e são, respectivamente, intra e entre participantes.
  Os dados estão na planilha CafIntra_AlcEntre.xlsx

Solução em alcb_cafw.R

source("summarySEwithin2.R")
Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8")

[1] "pt_BR.UTF-8"

Sys.setlocale("LC_ALL", "pt_BR.UTF-8")

[1] "LC_COLLATE=pt_BR.UTF-8;LC_CTYPE=pt_BR.UTF-8;LC_MONETARY=pt_BR.UTF-8;LC_NUMERIC=C;LC_TIME=pt_BR.UTF-8"

alfa <- 0.05
Dados.long <- readxl::read_xlsx("CafIntra_AlcEntre.xlsx")
Dados.long$UE <- factor(Dados.long$UE)
Dados.long$Cafeina <- factor(Dados.long$Cafeina, 
                             levels=c("SemC","ComC"))
Dados.long$Alcool <- factor(Dados.long$Alcool, 
                             levels=c("SemA","ComA"))
cat("\nDados.long")

Dados.long

print.data.frame(Dados.long, row.names=FALSE)

 UE NumErros Alcool Cafeina
  1        4   SemA    SemC
  1       15   SemA    ComC
  2        9   SemA    SemC
  2        8   SemA    ComC
  3       10   SemA    SemC
  3       17   SemA    ComC
  4        8   SemA    SemC
  4        0   SemA    ComC
  5        6   SemA    SemC
  5       15   SemA    ComC
  6       11   SemA    SemC
  6       11   SemA    ComC
  7        2   SemA    SemC
  7       11   SemA    ComC
  8       11   SemA    SemC
  8        6   SemA    ComC
  9       11   SemA    SemC
  9        0   SemA    ComC
 10       10   SemA    SemC
 10       15   SemA    ComC
 11        3   SemA    SemC
 11       17   SemA    ComC
 12       10   SemA    SemC
 12       15   SemA    ComC
 13       28   ComA    SemC
 13       10   ComA    ComC
 14       22   ComA    SemC
 14       11   ComA    ComC
 15       21   ComA    SemC
 15       27   ComA    ComC
 16       27   ComA    SemC
 16       15   ComA    ComC
 17       21   ComA    SemC
 17       27   ComA    ComC
 18       20   ComA    SemC
 18       10   ComA    ComC
 19       19   ComA    SemC
 19       21   ComA    ComC
 20       16   ComA    SemC
 20       15   ComA    ComC
 21       25   ComA    SemC
 21       19   ComA    ComC
 22       17   ComA    SemC
 22       21   ComA    ComC
 23       19   ComA    SemC
 23       15   ComA    ComC
 24       20   ComA    SemC
 24       15   ComA    ComC

print(xtabs(~Alcool+Cafeina+UE, data=Dados.long))

, , UE = 1

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    0    0

, , UE = 2

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    0    0

, , UE = 3

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    0    0

, , UE = 4

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    0    0

, , UE = 5

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    0    0

, , UE = 6

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    0    0

, , UE = 7

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    0    0

, , UE = 8

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    0    0

, , UE = 9

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    0    0

, , UE = 10

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    0    0

, , UE = 11

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    0    0

, , UE = 12

      Cafeina
Alcool SemC ComC
  SemA    1    1
  ComA    0    0

, , UE = 13

      Cafeina
Alcool SemC ComC
  SemA    0    0
  ComA    1    1

, , UE = 14

      Cafeina
Alcool SemC ComC
  SemA    0    0
  ComA    1    1

, , UE = 15

      Cafeina
Alcool SemC ComC
  SemA    0    0
  ComA    1    1

, , UE = 16

      Cafeina
Alcool SemC ComC
  SemA    0    0
  ComA    1    1

, , UE = 17

      Cafeina
Alcool SemC ComC
  SemA    0    0
  ComA    1    1

, , UE = 18

      Cafeina
Alcool SemC ComC
  SemA    0    0
  ComA    1    1

, , UE = 19

      Cafeina
Alcool SemC ComC
  SemA    0    0
  ComA    1    1

, , UE = 20

      Cafeina
Alcool SemC ComC
  SemA    0    0
  ComA    1    1

, , UE = 21

      Cafeina
Alcool SemC ComC
  SemA    0    0
  ComA    1    1

, , UE = 22

      Cafeina
Alcool SemC ComC
  SemA    0    0
  ComA    1    1

, , UE = 23

      Cafeina
Alcool SemC ComC
  SemA    0    0
  ComA    1    1

, , UE = 24

      Cafeina
Alcool SemC ComC
  SemA    0    0
  ComA    1    1

cat("\nStructure:\n")

Structure:

str(Dados.long)

tibble [48 × 4] (S3: tbl_df/tbl/data.frame)
 $ UE      : Factor w/ 24 levels "1","2","3","4",..: 1 1 2 2 3 3 4 4 5 5 ...
 $ NumErros: num [1:48] 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 ...

cat("Analisando número de erros no simulador")

Analisando número de erros no simulador

alfabonf <- alfa/(length(unique(Dados.long$Alcool))*
                    length(unique(Dados.long$Cafeina)))

rmiscsum <- summarySEwithin2(Dados.long,
                             measurevar="NumErros",
                             betweenvars="Alcool",
                             withinvars="Cafeina",
                             idvar="UE",
                             na.rm=TRUE, 
                             conf.interval=1-alfabonf)
print(rmiscsum)

  Alcool Cafeina  N  NumErros NumErrosNormed       sd       se       ci
1   ComA    ComC 12 17.166667       12.25000 5.472729 1.579841 4.709303
2   ComA    SemC 12 21.250000       16.33333 5.472729 1.579841 4.709303
3   SemA    ComC 12 10.833333       15.75000 5.587798 1.613058 4.808321
4   SemA    SemC 12  7.916667       12.83333 5.587798 1.613058 4.808321

cat("\nIntervalo de confiança ")

Intervalo de confiança 

cat("\ncom correção de Bonferroni (Alcool: entre, Cafeina: intra)")

com correção de Bonferroni (Alcool: entre, Cafeina: intra)

grf <- ggplot2::ggplot(rmiscsum, 
                       ggplot2::aes(x=Cafeina, 
                                    y=NumErros,
                                    colour=Alcool,
                                    group=Alcool)) +
  ggplot2::geom_line() +
  ggplot2::geom_errorbar(width=.3, 
                         ggplot2::aes(ymin=NumErros-ci, 
                                      ymax=NumErros+ci)) +
  ggplot2::geom_point(shape=21, 
                      size=3, 
                      fill="white") +
  ggplot2::ylab("NumErros") +
  ggplot2::ggtitle("Within-caffeine, between-alcohol subjects CI95 Bonferroni") +
  ggplot2::theme_bw() 
print(grf)

cat("\nANOVA bifatorial mista (aov): balanceamento perfeito: Faraway(2016)")

ANOVA bifatorial mista (aov): balanceamento perfeito: Faraway(2016)

fit.aov <- aov(NumErros ~ Alcool + Cafeina + Alcool:Cafeina + 
                          Error(UE), 
               data=Dados.long)
print(summary(fit.aov))

Error: UE
          Df Sum Sq Mean Sq F value   Pr(>F)    
Alcool     1 1160.3  1160.3   64.21 5.75e-08 ***
Residuals 22  397.6    18.1                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: Within
               Df Sum Sq Mean Sq F value Pr(>F)  
Cafeina         1    4.1    4.08   0.133 0.7183  
Alcool:Cafeina  1  147.0  147.00   4.806 0.0392 *
Residuals      22  672.9   30.59                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\nANOVA bifatorial mista (lmerTest::lmer)")

ANOVA bifatorial mista (lmerTest::lmer)

cat("\n\tteste omnibus")

    teste omnibus

fit.lmer <- lmerTest::lmer(NumErros ~ Alcool + Cafeina + Alcool:Cafeina +
                                      (1|UE),
                           data=Dados.long)

boundary (singular) fit: see help('isSingular')

print(car::Anova(fit.lmer, 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(fit.lmer, correl = FALSE))

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: NumErros ~ Alcool + Cafeina + Alcool:Cafeina + (1 | UE)
   Data: Dados.long

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

print(sjPlot::tab_model(fit.lmer,
                        p.adjust="holm",
                        encoding="Windows-1252",
                        df.method="satterthwaite",
                        title="p.adjust=holm",
                        show.reflvl=TRUE,
                        show.ngroups=TRUE,
                        show.stat=TRUE,
                        show.df=TRUE,
                        show.se=TRUE,
                        digits=4,
                        digits.p=4,
                        digits.rsq=6,
                        digits.re=6))

options(warn=-1)
o.par <- par()
fit.means <- phia::interactionMeans(fit.lmer)
plot(fit.means, 
     errorbar=paste0("ci",round((1-alfa)*100,0)))
par(o.par)

options(warn=0)

# Se ausência de efeito de interação não é rejeitado
# https://schmidtpaul.github.io/DSFAIR/compactletterdisplay.html
model_means <- emmeans::emmeans(object=fit.lmer,
                                specs=pairwise~Alcool,
                                level=1-alfa,
                                adjust="holm",
                                lmer.df="satterthwaite",
                                lmerTest.limit=nrow(Dados.long))

NOTE: Results may be misleading due to involvement in interactions

print(summary(model_means$emmeans, infer=TRUE))

 Alcool emmean   SE df lower.CL upper.CL t.ratio p.value
 SemA     9.38 1.01 44     7.35     11.4   9.311  <.0001
 ComA    19.21 1.01 44    17.18     21.2  19.078  <.0001

Results are averaged over the levels of: Cafeina 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 

print(summary(model_means$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 

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 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. 

model_means <- emmeans::emmeans(object=fit.lmer,
                                specs=pairwise~Cafeina,
                                level=1-alfa,
                                adjust="holm",
                                lmer.df="satterthwaite",
                                lmerTest.limit=nrow(Dados.long))

NOTE: Results may be misleading due to involvement in interactions

print(summary(model_means$emmeans, infer=TRUE))

 Cafeina emmean   SE df lower.CL upper.CL t.ratio p.value
 SemC      14.6 1.01 44     12.6     16.6  14.484  <.0001
 ComC      14.0 1.01 44     12.0     16.0  13.905  <.0001

Results are averaged over the levels of: Alcool 
Degrees-of-freedom method: satterthwaite 
Confidence level used: 0.95 

print(summary(model_means$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 

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 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. 

model_means <- emmeans::emmeans(object=fit.lmer,
                                specs=pairwise~Cafeina:Alcool,
                                level=1-alfa,
                                adjust="holm",
                                lmer.df="satterthwaite",
                                lmerTest.limit=nrow(Dados.long))
print(summary(model_means$emmeans, infer=TRUE))

 Cafeina Alcool emmean   SE df lower.CL upper.CL t.ratio p.value
 SemC    SemA     7.92 1.42 44     5.05     10.8   5.560  <.0001
 ComC    SemA    10.83 1.42 44     7.96     13.7   7.608  <.0001
 SemC    ComA    21.25 1.42 44    18.38     24.1  14.924  <.0001
 ComC    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))

 contrast              estimate   SE df lower.CL upper.CL t.ratio p.value
 SemC SemA - ComC SemA    -2.92 2.01 44    -8.48     2.65  -1.448  0.1546
 SemC SemA - SemC ComA   -13.33 2.01 44   -18.90    -7.77  -6.621  <.0001
 SemC SemA - ComC ComA    -9.25 2.01 44   -14.81    -3.69  -4.594  0.0001
 ComC SemA - SemC ComA   -10.42 2.01 44   -15.98    -4.85  -5.173  <.0001
 ComC SemA - ComC ComA    -6.33 2.01 44   -11.90    -0.77  -3.145  0.0089
 SemC ComA - ComC ComA     4.08 2.01 44    -1.48     9.65   2.028  0.0973

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 

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 Alcool emmean   SE df lower.CL upper.CL .group
 SemC    SemA     7.92 1.42 44     4.21     11.6  a    
 ComC    SemA    10.83 1.42 44     7.12     14.5  a    
 ComC    ComA    17.17 1.42 44    13.46     20.9   b   
 SemC    ComA    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. 

# Se ausência de efeito de interação é rejeitado
cat("\nTeste de efeito principal simples de Álcool")

Teste de efeito principal simples de Álcool

print(phia::testInteractions(fit.lmer, 
                             fixed="Cafeina", 
                             across="Alcool",
                             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

cat("\nTeste post hoc de efeito principal simples de Álcool")

Teste post hoc de efeito principal simples de Álcool

model_means <- emmeans::emmeans(object=fit.lmer,
                                specs=pairwise~Alcool|Cafeina:Alcool,
                                level=1-alfa,
                                adjust="holm",
                                lmer.df="satterthwaite",
                                lmerTest.limit=nrow(Dados.long))
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("\nTeste de efeito principal simples de Cafeína")

Teste de efeito principal simples de Cafeína

print(phia::testInteractions(fit.lmer, 
                             fixed="Alcool", 
                             across="Cafeina",
                             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

cat("\nTeste post hoc de efeito principal simples de Cafeína")

Teste post hoc de efeito principal simples de Cafeína

model_means <- emmeans::emmeans(object=fit.lmer,
                                specs=pairwise~Cafeina|Cafeina:Alcool,
                                level=1-alfa,
                                adjust="holm",
                                lmer.df="satterthwaite",
                                lmerTest.limit=nrow(Dados.long))
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.