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

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

suppressMessages(library(car, warn.conflicts=FALSE))
suppressMessages(library(DiagrammeR, warn.conflicts=FALSE))
suppressMessages(library(dplyr, warn.conflicts=FALSE))
suppressMessages(library(effects, warn.conflicts=FALSE))
suppressMessages(library(effectsize, warn.conflicts=FALSE))
suppressMessages(library(emmeans, warn.conflicts=FALSE))
suppressMessages(library(GGally, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(MASS, warn.conflicts=FALSE))
suppressMessages(library(multcomp, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
suppressMessages(library(rgl, warn.conflicts=FALSE))
suppressMessages(library(tidyr, warn.conflicts=FALSE))
source("eiras.bartitle.R")

Adicionar

  • ANCOVA_and_T2.R

Material

  • HTML de Rmarkdown em RPubs

Objetivos

  • Conceituar ANCOVA, distinguindo-a da ANOVA e da Regressão Linear Simples.
  • Discorrer sobre as suposições e condições nas quais ANCOVA pode ser utilizada.
  • Descrever e indicar dois delineamentos diferentes da ANCOVA:
    • com um fator (nominal) entre participantes e uma covariável (intervalar),
    • em delineamento pré-pós.
  • Proceder com a análise descritiva numérica e gráfica dos dados.
  • Formular e implementar os modelos de ANCOVA adequados a cada situação.
  • Determinar as significâncias estatística e prática.

ANCOVA independente simples

Conforme Cochran (1957),

“Como Fisher (1934) expressou, a análise de covariância ‘combina as vantagens e reconcilia as exigências de dois procedimentos muito amplamente aplicáveis, conhecidos como regressão e análise de variância’.”

ANCOVA independente simples é uma combinação de regressão linear simples com ANOVA independente unifatorial. Este é o primeiro passo para o modelo linear geral univariado (GLM: General Linear Model).

ANOVA unifatorial simples

ANOVA unifatorial simples

Existem duas razões principais para usar ANCOVA:

1. Reduzir a variância do resíduo:

A ANCOVA introduz uma covariável (intervalar) que também explica parcialmente a variância da VD. Supõe-se que esta covariável não esteja correlacionada com o fator. Com isso, a variância da VD é particionada em dois tamanhos de efeito parciais, uma parte atribuída ao fator, \(\eta^2_{\text{Fator}}\) e outra parte associada com a covariável, \(\eta^2_{\text{Cov}}\). São estes tamanhos de efeito parciais que aparecerão na análise, adiante.

ANCOVA unifatorial simples

ANCOVA unifatorial simples

A introdução desta covariável permite o controlar o efeito do fator: a parcela da variância da VD atribuída à covariável é subtraída pelo modelo.

ANCOVA unifatorial simples

ANCOVA unifatorial simples

Como os tamanhos de efeitos são porcentagens explicadas da variância da VD pela covariável ou pelo fator, subtraindo-se a parte atribuída à covariável, altera-se o denominador. Com isto, a porcentagem remanescente explicada pelo fator é proporcionalmente maior do que aquela computada pela ANOVA. O modelo ganha em poder explicativo, e este é o tamanho de efeito global, denotado por \(\eta^2_{\text{Global}}\) ou simplesmente \(\eta^2\) que aparece nas análises adiante.

Pela ilustração também é fácil de perceber que a soma dos tamanhos de efeito parciais não correspondem necessariamente ao tamanho de efeito global.

2. Ajustar as médias da covariável, de modo que o valor médio da covariável seja o mesmo em todos os grupos e depois ajustar as médias e erros-padrão da VD nos grupos (equalização estatística de condições entre participantes):

O uso de uma covariável é uma forma de obter controle estatístico quando o controle experimental, obtido pelo delineamento do estudo, não for possível. Testes mais simples, como teste t ou ANOVA, muitas vezes podem ser usados em situações laboratoriais. Em um experimento com camundongos, havendo recursos financeiros, o pesquisador pode ter, por exemplo, animais com a mesma idade ou até mesmo com a mesma genética.

Suponha, no entanto, que os animais disponíveis têm idades variadas. Controlar o efeito do experimento pela idade pode ser uma boa ideia. Críticos da ANCOVA dirão que esta dificuldade poderia ser resolvida com randomização dos animais entre os grupos ou condições experimentais (na linguagem do GLM: fator). Porém, randomização mitiga a discrepância entre as condições experimentais em média (com a intenção, neste exemplo, de equalizar as idades), mas não garante esta homogeneidade completamente em um determinado experimento, especialmente se a amostra for pequena.

Enquanto a ANCOVA pode ser preterida em relação à ANOVA em delineamentos experimentais nos quais o pesquisador assume que a randomização é suficiente, em estudos quase-experimentais ou observacionais, nos quais randomização não é realizada, ANCOVA é uma estratégia mais defensável.

GLM é uma das formas mais simples para se operacionalizar este controle estatístico na análise para testar o efeito fixo do fator entre participantes.

Por “ajustar as médias da covariável, de modo que o valor médio da covariável seja o mesmo em todos os grupos”, GLM computa a média da covariável em cada condição do fator, e então utiliza a média destas médias como a média representativa para a covariável.

O segundo passo é “ajustar as médias e erros-padrão da VD nos grupos (equalização estatística de condições entre participantes)”, que é alterar os valores das médias da VD para cada uma das condições do fator, o que é conhecido como médias marginais ajustadas, que também aparecerão nas análises da ANCOVA. Da mesma forma, os erros-padrão são homogeneizados. Veremos, adiante, como estes valores de média e erro-padrão são computados na aplicação do modelo.

O centramento da covariável não muda o resultado estatístico da ANCOVA (mesmo valor p), mas melhora a interpretação dos parâmetros do modelo.

Exemplo: Álcool na direção com experiência

Em Dancey & Reidy (2019, p. 424), um experimento foi projetado para investigar se o álcool (nas condições placebo, ingestão em pequena e em grande quantidade) afeta o desempenho na condução veicular, medido pelo número de erros cometidos em simulador de direção veicular durante um determinado período de tempo.

No entanto, uma suposição perfeitamente razoável é que a experiência de direção também está relacionada aos erros cometidos ao dirigir, mesmo num simulador. Quanto mais experiente for o motorista, menor será a quantidade de erros. Assim, a experiência do motorista está negativamente associada aos erros ao dirigir.

O delineamento é entre participantes. O pesquisador randomizou participantes entre os três grupos mas, ainda assim, gostaria de mitigar o efeito da experiência dos motoristas para estudar o efeito do álcool. O dados estão disponíveis em AlcoholExperienceDrivingErrors.rds.

Dados <- data.frame(readxl::read_excel("AlcoholExperienceDrivingErrors.xlsx"))
Dados$Grupo <- factor(Dados$Grupo,
                      levels=c("Placebo","PoucoAlcool","MuitoAlcool"))
saveRDS(Dados,"AlcoholExperienceDrivingErrors.rds")

Estatística descritiva, testes das suposições e testes omnibus e post hoc da ANCOVA estão implementados em demo_ANCOVA.R.

Estatística descritiva

Implementado em demo_ANCOVA_01descritiva.R:

suppressMessages(suppressWarnings(
  invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
source("eiras.bartitle.R")
Dados <- readRDS("AlcoholExperienceDrivingErrors.rds")

cat(bartitle("Dados"))
print.data.frame(Dados, row.names=FALSE)
cat(bartitle("Análise descritiva"))
print(xtabs(~Grupo, data=Dados))

cat("\n\nSumário:\n")
print(psych::describeBy(ErrosDirecao~Grupo,
                        mat=1,
                        digits=2,
                        data=Dados))
print(GGally::ggpairs(Dados[-1], 
                      ggplot2::aes(colour=Grupo, alpha=0.5))+ 
        ggplot2::theme_bw() + 
        ggplot2::theme(panel.grid = 
                         ggplot2::element_blank()))

-----
Dados
-----
 ID       Grupo ErrosDirecao Experiencia
  1     Placebo            5          12
  2     Placebo           10           5
  3     Placebo            7           9
  4     Placebo            3          24
  5     Placebo            5          15
  6     Placebo            7           6
  7     Placebo           11           3
  8     Placebo            2          30
  9     Placebo            3          20
 10     Placebo            5          10
 11     Placebo            6           7
 12     Placebo            6           8
 13 PoucoAlcool            5          21
 14 PoucoAlcool            7          16
 15 PoucoAlcool            9           7
 16 PoucoAlcool            8          15
 17 PoucoAlcool            2          30
 18 PoucoAlcool            5          21
 19 PoucoAlcool            6          12
 20 PoucoAlcool            6          13
 21 PoucoAlcool            4          26
 22 PoucoAlcool            4          24
 23 PoucoAlcool            8           9
 24 PoucoAlcool           10           6
 25 MuitoAlcool            8          29
 26 MuitoAlcool           10           8
 27 MuitoAlcool            8          26
 28 MuitoAlcool            9          20
 29 MuitoAlcool           11          18
 30 MuitoAlcool           15           6
 31 MuitoAlcool            7          12
 32 MuitoAlcool           11           7
 33 MuitoAlcool            8          15
 34 MuitoAlcool            8           9
 35 MuitoAlcool           17           3
 36 MuitoAlcool           11           7

------------------
Análise descritiva
------------------
Grupo
    Placebo PoucoAlcool MuitoAlcool 
         12          12          12 


Sumário:
              item      group1 vars  n  mean   sd median trimmed  mad min max
ErrosDirecao1    1     Placebo    1 12  5.83 2.69    5.5     5.7 2.22   2  11
ErrosDirecao2    2 PoucoAlcool    1 12  6.17 2.33    6.0     6.2 2.97   2  10
ErrosDirecao3    3 MuitoAlcool    1 12 10.25 3.05    9.5     9.9 2.22   7  17
              range  skew kurtosis   se
ErrosDirecao1     9  0.47    -0.83 0.78
ErrosDirecao2     8 -0.03    -1.15 0.67
ErrosDirecao3    10  0.98    -0.30 0.88
Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

O arquivo de dados contém as informações de cada voluntário (identificados pelo ID) em uma linha. O delineamento é completamente entreparticipantes (participantes diferentes foram submetidos à ingestão de Placebo, PoucoAlcool e MuitoAlcool). Este placebo é uma bebida idealizada (neste exemplo fictício) com características idênticas em aspecto e sabor à bebida alcóolica, mas que não contém álcool.

Três tabelas mostram, respectivamente, em cada uma das condições do fator:

  • distribuição dos erros cometidos,
  • distribuição da experiência dos participantes,
  • distribuição simultânea dos erro e experiência simultaneamente observados,
  • sumário dos erros de direção cometidos.

Na sequência há vários gráficos para que são adequados para o propósito da sua pesquisa.

Regressões lineares simples

Implementado em demo_ANCOVA_02lm.R:

suppressMessages(suppressWarnings(
  invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
source("eiras.bartitle.R")
Dados <- readRDS("AlcoholExperienceDrivingErrors.rds")
alfa <- 0.05

cat(bartitle("Regressão Linear Simples"))

car::scatterplot(ErrosDirecao~Experiencia|Grupo, 
                 grid=FALSE,
                 data=Dados,
                 regLine=TRUE, 
                 smooth=FALSE, 
                 col=c("black","blue","red"))

alfa <- 0.05
k <- nlevels(Dados$Grupo)
alfaBonf <- alfa/k
graf <- ggplot2::ggplot(Dados,
                        ggplot2::aes(y = ErrosDirecao,
                                     x = Experiencia,
                                     group = Grupo,
                                     linetype = Grupo,
                                     shape = Grupo)) +
  ggplot2::geom_point(alpha = 1/4, size = 1) +
  ggplot2::theme_bw() +
  ggplot2::theme(panel.grid = ggplot2::element_blank()) +
  ggplot2::labs(title = "RLS por Grupo",
                subtitle = "BC95% Bonferroni",
                x = "Experiência (ano)",
                y = "#ErrosDirecao") +
  #ggplot2::geom_smooth(method = stats::lm,
  #ggplot2::geom_smooth(method = MASS::rlm,
  ggplot2::geom_smooth(method = estimatr::lm_robust,
                       formula = y ~ x, 
                       na.rm = TRUE,
                       color = "black", 
                       ggplot2::aes(fill = Grupo),
                       level = 1 - alfaBonf)
plot(graf)

for (g in unique(Dados$Grupo))
{
  cat(bartitle(g,2))
  dt_tmp <- subset(Dados, Grupo==g)
  modelo <- lm(ErrosDirecao~Experiencia, data=dt_tmp)
  print(summary(modelo))
}

------------------------
Regressão Linear Simples
------------------------


    -------
    Placebo
    -------

Call:
lm(formula = ErrosDirecao ~ Experiencia, data = dt_tmp)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.5153 -0.9831 -0.3742  0.6088  2.5093 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.33732    0.74221  12.580 1.87e-07 ***
Experiencia -0.28220    0.05034  -5.606 0.000226 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.387 on 10 degrees of freedom
Multiple R-squared:  0.7586,    Adjusted R-squared:  0.7345 
F-statistic: 31.43 on 1 and 10 DF,  p-value: 0.0002257


    -----------
    PoucoAlcool
    -----------

Call:
lm(formula = ErrosDirecao ~ Experiencia, data = dt_tmp)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.4829 -0.3483  0.0555  0.5106  1.3633 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 10.86731    0.59864  18.153 5.52e-09 ***
Experiencia -0.28204    0.03281  -8.595 6.25e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8434 on 10 degrees of freedom
Multiple R-squared:  0.8808,    Adjusted R-squared:  0.8689 
F-statistic: 73.88 on 1 and 10 DF,  p-value: 6.247e-06


    -----------
    MuitoAlcool
    -----------

Call:
lm(formula = ErrosDirecao ~ Experiencia, data = dt_tmp)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5482 -1.5513 -0.2127  1.3885  4.4392 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 13.23169    1.42182   9.306 3.06e-06 ***
Experiencia -0.22363    0.09149  -2.444   0.0346 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.53 on 10 degrees of freedom
Multiple R-squared:  0.374, Adjusted R-squared:  0.3114 
F-statistic: 5.974 on 1 and 10 DF,  p-value: 0.0346

São apresentados dois tipos de gráfico. O primeiro é gerado pela função car::scatterplot que plota retas de regressão automaticamente, usando-se o parâmetro regLine=TRUE. O segundo usa recursos do pacote ggplot2, parametrizado para exibir as linhas de regressão com bandas de confiança 95%. Observe que as retas parecem paralelas, mas quando produzimos as três bandas de confiança de 95% com correção de Bonferroni, não conseguimos visualizar um segmento de reta comum às três retas no domínio observado comum de Experiencia, sugerindo que não há uma reta populacional de associação linear entre erros de direção e experiência em direção que atenda aos três grupos (i.e., possivelmente o álcool distingue as três populações).

Temos, então, três regressões lineares simples, mostrando que há associação negativa significante da experiência em direção com o número de erros cometidos no simulador em cada uma das condições do fator. Visualmente, as três regressões parecem paralelas, o que nos leva à discussão sobre as suposições necessárias para a aplicação da ANCOVA.

Suposições

Além das suposições da regressão linear simples e da ANOVA, a ANCOVA supõe:

  • dissociação entre o fator e a covariável,
  • covariável sem erro de mensuração (não testável, assumida pelo modelo),
  • linearidade entre a VD e a covariável em todos os níveis do fator,
  • igualdade das inclinações das regressões lineares simples.

Implementado em demo_ANCOVA_03suposicoes.R:

suppressMessages(suppressWarnings(
  invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
source("eiras.bartitle.R")
Dados <- readRDS("AlcoholExperienceDrivingErrors.rds")

cat(bartitle("ANCOVA"))

cat(bartitle("Suposições",2))

## ANCOVA: Teste de dissociacao entre fator e covariavel
cat(bartitle("Dissociação entre Grupo (fator) e Experiencia (covariável)",3))
modelo <- lm(Experiencia ~ Grupo, 
             data=Dados)
suppressMessages(
  print(Anova <- car::Anova(modelo))
)

## ANCOVA: Teste de linearidade da VD e covariavel nos niveis do fator
cat(bartitle("Linearidade entre ErrosDirecao (VD) e Experiencia (covariável)",3))
levels <- as.vector(unique(Dados$Grupo))
for (l in levels)
{
  modelo <- lm(ErrosDirecao ~ Experiencia, 
               data=Dados[Dados$Grupo==l,])
  hc <- lmtest::harvtest(modelo) # Harvey-Collier test for linearity
  cat(bartitle(l,4))
  print(hc)
}

## ANCOVA: Teste de igualdade das inclinacoes das retas de regressao
cat(bartitle("Igualdade das inclinações das regressões",3))
modelo <- lm(ErrosDirecao ~ Grupo + Experiencia + Grupo:Experiencia, 
             data=Dados)
suppressMessages(
  print(Anova <- car::Anova(modelo))
)

------
ANCOVA
------

    ----------
    Suposições
    ----------

        ----------------------------------------------------------
        Dissociação entre Grupo (fator) e Experiencia (covariável)
        ----------------------------------------------------------
Anova Table (Type II tests)

Response: Experiencia
           Sum Sq Df F value Pr(>F)
Grupo      120.06  2  0.9069 0.4136
Residuals 2184.25 33               

        --------------------------------------------------------------
        Linearidade entre ErrosDirecao (VD) e Experiencia (covariável)
        --------------------------------------------------------------

            -------
            Placebo
            -------

    Harvey-Collier test

data:  modelo
HC = 0.48977, df = 9, p-value = 0.636


            -----------
            PoucoAlcool
            -----------

    Harvey-Collier test

data:  modelo
HC = 0.75424, df = 9, p-value = 0.47


            -----------
            MuitoAlcool
            -----------

    Harvey-Collier test

data:  modelo
HC = 0.00074081, df = 9, p-value = 0.9994


        ----------------------------------------
        Igualdade das inclinações das regressões
        ----------------------------------------
Anova Table (Type II tests)

Response: ErrosDirecao
                   Sum Sq Df F value    Pr(>F)    
Grupo             136.334  2 22.6338 1.018e-06 ***
Experiencia       149.531  1 49.6493 7.827e-08 ***
Grupo:Experiencia   1.701  2  0.2823     0.756    
Residuals          90.352 30                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Os dois testes das suposições apresentados nesta saída são:

  • Dissociação entre Grupo (fator) and Experiencia (covariável)

A hipótese nula é a ausência de associação entre o fator Grupo e a covariável Experiencia dos participantes. Isto é feito por uma ANOVA independente com:

modelo <- lm(Experiencia ~ Grupo, data=Dados)
print(Anova <- car::Anova(modelo))
Anova Table (Type II tests)

Response: Experiencia
           Sum Sq Df F value Pr(>F)
Grupo      120.06  2  0.9069 0.4136
Residuals 2184.25 33               

O valor \(p=0.4136\) não rejeita \(H_0\) e, portanto, assumimos dissociação entre o fator e a covariável.

  • Linearidade entre ErrosDirecao (VD) e Experiencia (covariável) em cada nível de Grupo (fator)

As hipóteses nulas para o teste de Harvey-Collier são de linearidade entre as variáveis testadas. São feitas regressões segmentando os dados para cada grupo com:

levels <- as.vector(unique(Dados$Grupo))
for (l in levels)
{
  modelo <- lm(ErrosDirecao ~ Experiencia, data=Dados[Dados$Grupo==l,])
  hc <- lmtest::harvtest(modelo) # Harvey-Collier test for linearity
  cat(bartitle(l,4))
  print(hc)
}

            -------
            Placebo
            -------

    Harvey-Collier test

data:  modelo
HC = 0.48977, df = 9, p-value = 0.636


            -----------
            PoucoAlcool
            -----------

    Harvey-Collier test

data:  modelo
HC = 0.75424, df = 9, p-value = 0.47


            -----------
            MuitoAlcool
            -----------

    Harvey-Collier test

data:  modelo
HC = 0.00074081, df = 9, p-value = 0.9994

Para os três níveis do fator, não se rejeitou a hipótese nula de ausência de linearidade.

  • Homogeneidade da inclinação das regressões lineares simples

A hipótese nula é o paralelismo das três regressões feitas, em cada nível do fator, entre a experiência (covariável) e o número de erros cometidos no simulador (VD). Isto é feito por um GLM que considera o fator, a covariável e a interação entre fator e covariável com:

modelo <- lm(ErrosDirecao ~ Grupo + Experiencia + Grupo:Experiencia, 
             data=Dados)
suppressMessages(
  print(Anova <- car::Anova(modelo))
)
Anova Table (Type II tests)

Response: ErrosDirecao
                   Sum Sq Df F value    Pr(>F)    
Grupo             136.334  2 22.6338 1.018e-06 ***
Experiencia       149.531  1 49.6493 7.827e-08 ***
Grupo:Experiencia   1.701  2  0.2823     0.756    
Residuals          90.352 30                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A inclusão da interação Grupo:Experiencia é um artifício para o teste do paralelismo das retas. O valor desta interação, \(p=0.756\) é o único relevante para a decisão estatística. Neste caso, não rejeitamos \(H_0\) e, portanto, as três regressões lineares simples são assumidas como retas populacionais paralelas.

Teste do efeito do fator

A ANCOVA propriamente dita tem hipótese nula de ausência de efeito do fator sobre a VD:

\[H_0: \mu^{\prime}_1 = \mu^{\prime}_2 = \cdots = \mu^{\prime}_k\]

em que \(k\) é o número de níveis do fator e \(\mu^{\prime}_i\) são as médias populacionais da VD ajustadas pela covariável.

Uma vez que as três regressões são assumidas paralelas, detectamos o efeito do fator (níveis de álcool) sobre a VD (número de erros de direção), controlando-se para a covariável (experiência em dirigir) caso as regressões tenham diferentes interceptos.

Implementado em demo_ANCOVA_04ANCOVA.R:

suppressMessages(suppressWarnings(
  invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
source("eiras.bartitle.R")
Dados <- readRDS("AlcoholExperienceDrivingErrors.rds")

cat(bartitle("Teste do efeito do fator fixo",2))
## ANCOVA: Teste do efeito do fator fixo: 
##          Se as declividades sao iguais, testar se os interceptos sao iguais.
cat(bartitle("Igualdade dos interceptos das regressões\n(assumindo igualdade das inclinações)\n- Teste de efeito do fator -",3))
ancova.fit <- lm(ErrosDirecao ~ Grupo + Experiencia, 
                 data=Dados)
print(Anova <- car::Anova(ancova.fit))
print(summary(ancova.fit))

    -----------------------------
    Teste do efeito do fator fixo
    -----------------------------

        ----------------------------------------
        Igualdade dos interceptos das regressões
(assumindo igualdade das inclinações)
- Teste de efeito do fator -
        ----------------------------------------
Anova Table (Type II tests)

Response: ErrosDirecao
             Sum Sq Df F value    Pr(>F)    
Grupo       136.334  2  23.697 4.851e-07 ***
Experiencia 149.531  1  51.981 3.452e-08 ***
Residuals    92.053 32                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Call:
lm(formula = ErrosDirecao ~ Grupo + Experiencia, data = Dados)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5989 -0.9540 -0.0951  0.8361  4.0463 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       9.08210    0.66541  13.649 6.85e-15 ***
GrupoPoucoAlcool  1.44533    0.70939   2.037   0.0499 *  
GrupoMuitoAlcool  4.65651    0.69322   6.717 1.38e-07 ***
Experiencia      -0.26165    0.03629  -7.210 3.45e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.696 on 32 degrees of freedom
Multiple R-squared:  0.762, Adjusted R-squared:  0.7397 
F-statistic: 34.15 on 3 and 32 DF,  p-value: 4.316e-10

Note que a sintaxe é similar à anterior, exceto pela remoção da interação Grupo:Experiencia. Computam-se os erros de direção em função do fator (Grupo) e da covariável (Experiencia); é uma tabela ANOVA da ANCOVA.

O valor \(p=4.851 \times 10^{-7}\) relevante aqui é aquele associado com o fator (Grupo), rejeitando-se a hipótese nula da ANCOVA que foi enunciado acima. Portanto, considera-se que a ingestão de álcool afeta significantemente o número de erros cometidos na direção.

Na sequência exibimos, também, a regressão linear associada à ANCOVA. O valor que aparece na coluna Estimate na linha da covariável Experiencia é a inclinação comum estimada para as três retas de regressão, que será utilizada adiante.

As últimas duas linhas da saída da regressão da ANCOVA mostram o efeito global (\(\eta^2=R^2=0.76\)) e o valor \(p=4.316 \times 10^{-10}\) que rejeita a hipótese nula da ausência de modelo ANCOVA; portanto, a ANCOVA proposta é um modelo válido para esta análise.

Tamanho de efeito

Como descrito anteriormente, computam-se o tamanho de efeito global e os tamanhos de efeito, implementadas em demo_ANCOVA_05tamefeito.R:

suppressMessages(suppressWarnings(
  invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
source("eiras.bartitle.R")
Dados <- readRDS("AlcoholExperienceDrivingErrors.rds")
alfa <- 0.05
cat(bartitle("Análise de tamanho de efeito",2))
cat(bartitle("Omnibus",2))
eta2 <- summary(ancova.fit)$r.squared
cat("R^2 = eta^2 omnibus = ", round(eta2,2), "\n", sep="")
print(effectsize::interpret_eta_squared(eta2))
eta2 <- effectsize::eta_squared(Anova,
                                partial=TRUE,
                                generalized=FALSE,
                                ci=1-alfa,
                                alternative="two.sided",
                                verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=4)

    ----------------------------
    Análise de tamanho de efeito
    ----------------------------

    -------
    Omnibus
    -------
R^2 = eta^2 omnibus = 0.76
[1] "large"
(Rules: field2013)
# Effect Size for ANOVA (Type II)

Parameter   | Eta2 (partial) |           95% CI | interpret
-----------------------------------------------------------
Grupo       |         0.5969 | [0.3520, 0.7331] |     large
Experiencia |         0.6190 | [0.3922, 0.7505] |     large

Média marginal ajustada

Os próximos itens da saída são as médias marginais ajustadas e erros-padrão ajustados. Estes são os valores estimados depois que os grupos foram estatisticamente equalizados pela covariável (i.e., os valores estimados para os erros de direção associados ao álcool sem o efeito da experiência). Por isso, são diferentes dos valores observados.

Este trecho está implementado em demo_ANCOVA_06emmeans.R:

suppressMessages(suppressWarnings(
  invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
source("eiras.bartitle.R")
Dados <- readRDS("AlcoholExperienceDrivingErrors.rds")
alfa <- 0.05

cat(bartitle("Médias marginais estimadas ajustadas",2))

EMM <- emmeans::emmeans(ancova.fit, 
                        pairwise~Grupo, 
                        adjust="holm",
                        level=1-alfa)
print(EMM$emmeans)
print(plot(EMM$emmeans,
           colors="black", 
           xlab="ErrosDirecao controlado por Experiencia",
           xlim=c(4,12),
           ylab="Grupo")+ggplot2::theme_bw())

    ------------------------------------
    Médias marginais estimadas ajustadas
    ------------------------------------
 Grupo       emmean    SE df lower.CL upper.CL
 Placebo       5.38 0.494 32     4.38     6.39
 PoucoAlcool   6.83 0.498 32     5.81     7.84
 MuitoAlcool  10.04 0.490 32     9.04    11.04

Confidence level used: 0.95 

Esta correção dos valores é realizada a partir da regressão linear associada à ANCOVA:

ancova.fit <- lm(ErrosDirecao ~ Grupo + Experiencia, 
                 data=Dados)
print(reg <- summary(ancova.fit))

Call:
lm(formula = ErrosDirecao ~ Grupo + Experiencia, data = Dados)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5989 -0.9540 -0.0951  0.8361  4.0463 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       9.08210    0.66541  13.649 6.85e-15 ***
GrupoPoucoAlcool  1.44533    0.70939   2.037   0.0499 *  
GrupoMuitoAlcool  4.65651    0.69322   6.717 1.38e-07 ***
Experiencia      -0.26165    0.03629  -7.210 3.45e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.696 on 32 degrees of freedom
Multiple R-squared:  0.762, Adjusted R-squared:  0.7397 
F-statistic: 34.15 on 3 and 32 DF,  p-value: 4.316e-10

Sendo uma regressão, a coluna Estimate contém o intercepto associado ao primeiro nível do fator (9.0821029), os interceptos dos outros níveis relativos ao primeiro nível (1.4453283 e 4.6565087), e a inclinação da reta comum a todos os níveis (-0.2616459). O objetivo é obter três retas com a mesma inclinação.

Os valores absolutos dos interceptos são dados pela soma dos valores relativos sobre o intercepto do primeiro nível:

Intercepts:
    Placebo = 9.082103
    PoucoAlcool = 9.082103 + 1.445328 = 10.52743
    MuitoAlcool = 9.082103 + 4.656509 = 13.73861

A inclinação comum às três retas é dada por

Common slope = -0.2616459

e a média comum da covariável é a média das médias da experiência dos participantes em cada um dos grupos:

        Grupo Experiencia
1     Placebo    12.41667
2 PoucoAlcool    16.66667
3 MuitoAlcool    13.33333
Common Covariate Mean = 14.13889

Finalmente, as médias marginais ajustadas correspondem aos valores do número de erros de direção para a média comum da covariável estimados pelas três regressões obtidas com os interceptos ajustados e inclinação comum:

Adjusted marginal means:
    Placebo = 9.082103 + -0.2616459 * 14.13889 = 5.382721
    PoucoAlcool = 10.52743 + -0.2616459 * 14.13889 = 6.828049
    MuitoAlcool = 13.73861 + -0.2616459 * 14.13889 = 10.03923

Uma visão gráfica destas médias marginais ajustadas é:

em que \(14.139\) é a média comum da covariável e os valores das médias marginais correspondem aos que foram obtidos com a função emmeans::emmeans()

EMM <- emmeans::emmeans(ancova.fit,
                        pairwise~Grupo,
                        adjust="holm",
                        level=1-alfa)
print(EMM$emmeans)
 Grupo       emmean    SE df lower.CL upper.CL
 Placebo       5.38 0.494 32     4.38     6.39
 PoucoAlcool   6.83 0.498 32     5.81     7.84
 MuitoAlcool  10.04 0.490 32     9.04    11.04

Confidence level used: 0.95 

No mesmo processo de cálculo das médias marginais, foram computados também erros-padrão ajustados (SE, standard errors), que costumam ter valores quase iguais para cada nível do fator. O cálculo para estes erros-padrão é um pouco mais trabalhoso (Liu, 2011; Crawley, 2012). Basta dizer aqui que estes valores são os utilizados para estimar o intervalo de confiança que aparece no gráfico correspondente.

O código completo para calcular as médias marginais ajustadas está em demo_AdjustedMarginalMeans.R.

Teste post hoc

Resta apenas localizar as diferenças dos números de erros entre os diferentes níveis do fator Grupo. Usamos os contrastes de Tukey da mesma forma que fizemos com ANOVA nos capítulos anteriores. Este trecho está implementado em demo_ANCOVA_07posthoc.R:

suppressMessages(suppressWarnings(
  invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
source("eiras.bartitle.R")
Dados <- readRDS("AlcoholExperienceDrivingErrors.rds")
alfa <- 0.05

cat(bartitle("Teste post hoc: Grupo (fator) ajustado por Experiencia (covariável)"))

print(summary(EMM$contrasts, infer=TRUE))
print(plot(EMM$contrasts, 
           colors="black")+ggplot2::theme_bw())
print(multcomp::cld(object=EMM$emmeans,
                    level=1-alfa,
                    adjust="holm",
                    Letters=letters,
                    alpha=alfa))

cat(bartitle("Teste post hoc: Grupo (fator) ajustado por Experiencia (covariável): Dunnett"))
EMM.contrast <- emmeans::contrast(EMM, 
                                  method="trt.vs.ctrl", 
                                  ref="Placebo",
                                  adjust="holm",
                                  level=1-alfa)
print(EMM.contrast)
print(plot(EMM.contrast,
           colors="black")+ggplot2::theme_bw())

cat(bartitle("Teste post hoc: Grupo (fator) ajustado por Experiencia (covariável): Consecutivo"))
EMM.contrast <- emmeans::contrast(EMM, 
                                  method="consec",
                                  adjust="holm",
                                  level=1-alfa)
print(EMM.contrast)
print(plot(EMM.contrast,
           colors="black")+ggplot2::theme_bw())

-------------------------------------------------------------------
Teste post hoc: Grupo (fator) ajustado por Experiencia (covariável)
-------------------------------------------------------------------
 contrast                  estimate    SE df lower.CL upper.CL t.ratio p.value
 Placebo - PoucoAlcool        -1.45 0.709 32    -3.24    0.347  -2.037  0.0499
 Placebo - MuitoAlcool        -4.66 0.693 32    -6.41   -2.905  -6.717  <.0001
 PoucoAlcool - MuitoAlcool    -3.21 0.703 32    -4.99   -1.435  -4.568  0.0001

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

 Grupo       emmean    SE df lower.CL upper.CL .group
 Placebo       5.38 0.494 32     4.14     6.63  a    
 PoucoAlcool   6.83 0.498 32     5.57     8.09   b   
 MuitoAlcool  10.04 0.490 32     8.80    11.28    c  

Confidence level used: 0.95 
Conf-level adjustment: bonferroni method for 3 estimates 
P value adjustment: holm method for 3 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. 

----------------------------------------------------------------------------
Teste post hoc: Grupo (fator) ajustado por Experiencia (covariável): Dunnett
----------------------------------------------------------------------------
 contrast              estimate    SE df t.ratio p.value
 PoucoAlcool - Placebo     1.45 0.709 32   2.037  0.0499
 MuitoAlcool - Placebo     4.66 0.693 32   6.717  <.0001

P value adjustment: holm method for 2 tests 


--------------------------------------------------------------------------------
Teste post hoc: Grupo (fator) ajustado por Experiencia (covariável): Consecutivo
--------------------------------------------------------------------------------
 contrast                  estimate    SE df t.ratio p.value
 PoucoAlcool - Placebo         1.45 0.709 32   2.037  0.0499
 MuitoAlcool - PoucoAlcool     3.21 0.703 32   4.568  0.0001

P value adjustment: holm method for 2 tests 

Concluímos, portanto que os níveis A (Placebo) e B (PoucoAlcool) são não estatisticamente diferentes, mas o nível C (MuitoAlcool) é o motivo da rejeição da hipótese nula e, portanto, Álcool tem efeito no número de erros cometidos no simulador de direção veicular.

ANCOVA mista

Em outra situação, os mesmos participantes são avaliados antes e depois de uma intervenção (intraparticipantes), mas são divididos em grupos que recebem tratamentos diferentes (entre participantes). O interesse é avaliar o efeito das tratamentos.

Segundo Borm et al., 2007:

“Uma vantagem do uso da ANCOVA é que ela ajusta as diferenças basais entre os grupos de tratamento. ANCOVA também possui maior poder estatístico do que o teste t, portanto requer um tamanho amostral menor.”

Exemplo: Cirurgia de grande porte com anestesia geral

Dez pacientes idosos foram submetidos a uma cirurgia de grande porte com anestesia geral. Divididos em dois grupos, cada um tratado com uma droga diferente: uma droga tradicional (padrão) e uma nova proposta terapêutica (nova). Um teste cognitivo foi aplicado alguns dias antes da cirurgia e no período pós-operatório. Espera-se que as drogas mitiguem delírio pós-operatório (DPO) e disfunção cognitiva pós-operatória (DCPO). Este estudo é duplo cego para a comparação de duas drogas, de forma que serão denominadas de droga A e droga B, sem que o analista saiba qual delas corresponde à droga padrão e qual corresponde à droga nova. Os dados obtidos estão em PrePos_2CondIndep.rds, organizados no formato long.

Dados.long <- data.frame(readxl::read_excel("PrePos_2CondIndep.xlsx"))
Dados.long$Grupo <- factor(Dados.long$Grupo,
                           levels=c("droga A","droga B"))
Dados.long$Momento <- factor(Dados.long$Momento,
                             levels=c("Pre","Pos"))
print.data.frame(Dados.long)
   Paciente   Grupo Momento Escore
1         1 droga A     Pre    150
2         2 droga A     Pre    130
3         3 droga A     Pre    125
4         4 droga A     Pre    152
5         5 droga A     Pre    160
6         6 droga B     Pre    174
7         7 droga B     Pre    110
8         8 droga B     Pre    180
9         9 droga B     Pre    145
10       10 droga B     Pre    140
11        1 droga A     Pos     51
12        2 droga A     Pos     50
13        3 droga A     Pos     40
14        4 droga A     Pos     45
15        5 droga A     Pos     60
16        6 droga B     Pos     75
17        7 droga B     Pos     41
18        8 droga B     Pos     80
19        9 droga B     Pos     60
20       10 droga B     Pos     55
# converte long -> wide
Dados.wide <- tidyr::spread(Dados.long, Momento, Escore)
print.data.frame(Dados.wide)
   Paciente   Grupo Pre Pos
1         1 droga A 150  51
2         2 droga A 130  50
3         3 droga A 125  40
4         4 droga A 152  45
5         5 droga A 160  60
6         6 droga B 174  75
7         7 droga B 110  41
8         8 droga B 180  80
9         9 droga B 145  60
10       10 droga B 140  55
lista <- list(Dados.wide, Dados.long)
names(lista) <- c("Dados.wide", "Dados.long")
saveRDS(lista, "PrePos_2CondIndep.rds")

Conforme Sakae et al. (2016), “A anestesia geral age predominantemente no sistema nervoso central, repercutindo também, em todos os aparelhos e sistemas do organismo. Seu mecanismo intrínseco de ação ainda não é completamente conhecido e por isso a possibilidade de algum prejuízo temporário ou permanente na cognição e na memória sempre foi alvo de considerações. Há uma especial preocupação quanto aos idosos, por apresentarem maior susceptibilidade às alterações da homeostasia e do meio ambiente. A cirurgia e a anestesia exercem comparativamente efeitos adversos cerebrais mais acentuados nos idosos do que nos jovens, manifestado pela maior prevalência de delírio pós-operatório e disfunção cognitiva. O delírio pós-operatório (DPO) e a disfunção cognitiva pós-operatória (DCPO) atrasam a reabilitação, e são associadas com o aumento na morbidade e na mortalidade de pacientes idosos. Embora seja difícil estabelecer metodologicamente qualquer correlação entre DPO e DCPO, um recente estudo sugere que ambas podem representar uma trajetória de insuficiência cognitiva pós-operatória, talvez como uma progressão de não reconhecida insuficiência cognitiva leve pré-operatória. […] A incidência de DCPO é elevada (40%) em pacientes com mais de 65 anos, submetidos a cirurgias de médio e grande porte e, em cirurgias cardíacas, esta porcentagem pode ser maior que 60%. […] Aproximadamente 10% dos pacientes cirúrgicos idosos desenvolvem DPO, aumentando para 30-65% após certos tipos de cirurgia, tais como cirurgia cardíaca, cirurgias de emergência e fratura de quadril.”

ANCOVA mista é similar à independente, mas o papel da covariável é feito pela medida anterior à intervenção (Pre), que funciona como controle da linha de base.

Quem faz o papel da variável de desfecho é a medida posterior à intervenção (Pos). O fator é a droga, e neste exemplo existem dois níveis (droga A e B); são dois grupos de pacientes independentes, mas cada um deles é medido duas vezes. Neste sentido, o delineamento é misto: entre participantes e intraparticipantes.

Esta análise foi implementada em demo_ANCOVA_prepos.R. O código é praticamente o mesmo que o anterior, com pequenos ajustes por causa da estatística descritiva e pelo uso dos formatos long e wide de acordo com a conveniência das funções utilizadas.

Obtém-se a seguinte saída:

suppressMessages(suppressWarnings(
  invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))))
source("eiras.bartitle.R")
Dados.long <- readRDS("PrePos_2CondIndep.rds")$Dados.long
Dados.wide <- readRDS("PrePos_2CondIndep.rds")$Dados.wide

cat(bartitle("Data"))
print.data.frame(Dados.long)

cat(bartitle("Descriptive Statistics"))
print(ftable(Dados.long[,c("Grupo","Momento","Paciente")]))
print(psych::describeBy(Escore ~ Momento*Grupo,
                        mat=1,
                        digits=2,
                        data=Dados.long))
boxplot(Escore~Momento*Grupo, data=Dados.long, ylab="Escore")
print(GGally::ggpairs(Dados.wide[-1], 
                      ggplot2::aes(colour=Grupo, alpha=0.5))+ 
        ggplot2::theme_bw() + 
        ggplot2::theme(panel.grid = 
                         ggplot2::element_blank()))

cat(bartitle("Simple linear regressions"))

# formato wide é melhor
car::scatterplot(Pos~Pre|Grupo, 
                 grid=FALSE,
                 data=Dados.wide, 
                 regLine=TRUE, 
                 smooth=FALSE, 
                 col="black")
# o mesmo com os dados no formato long
# car::scatterplot(Escore[Momento=="Pos"]~Escore[Momento=="Pre"]|Grupo[Momento=="Pre"], data=Dados.long,
#                  regLine=TRUE, smooth=FALSE, col="black")

# testes usam o formato wide
alfa <- 0.05
k <- nlevels(Dados$Grupo)
alfaBonf <- alfa/k
graf <- ggplot2::ggplot(Dados.wide, 
                        ggplot2::aes(y=Pos, 
                                     x=Pre, 
                                     group=Grupo, 
                                     linetype=Grupo,
                                     shape=Grupo))+ 
  ggplot2::geom_point(alpha=1/4,
                      size=1)+ 
  ggplot2::theme_bw()+
  ggplot2::theme(panel.grid = ggplot2::element_blank()) +
  ggplot2::labs(title="RLS por Grupo",
                subtitle="BC95% Bonferroni",
                x="Escore Pre", 
                y="Escore Pos")+
  ggplot2::geom_smooth(method=MASS::rlm,
                       na.rm=TRUE,
                       color="black",
                       level=1-alfaBonf)
plot(graf)

for (g in unique(Dados.wide$Grupo))
{
  cat(bartitle(g,3))
  dt_tmp <- subset(Dados.wide, Grupo==g)
  modelo <- lm(Pos~Pre, data=dt_tmp)
  print(summary(modelo))
}

cat(bartitle("ANCOVA"))

cat(bartitle("Assumptions",2))

## ANCOVA: Teste de dissociacao entre fator e covariavel
cat(bartitle("Dissociation between Grupo (factor) and VD Pre (covariate)",3))
modelo <- lm(Pre ~ Grupo, 
             data=Dados.wide)
print(Anova <- car::Anova(modelo))

## ANCOVA: Teste de linearidade da VD e covariavel nos niveis do fator
cat(bartitle("Linearity between Escore pos (VD) and Escore pre (covariate)",3))
levels <- as.vector(unique(Dados.wide$Grupo))
for (l in levels)
{
  modelo <- lm(Pos ~ Pre, 
               data=Dados.wide[Dados.wide$Grupo==l,])
  hc <- lmtest::harvtest(modelo) # Harvey-Collier test for linearity
  cat(bartitle(l,4))
  print(hc)
}

## ANCOVA: Teste de igualdade das inclinacoes das retas de regressao
cat(bartitle("Homogenous slope of regressions",3))
modelo <- lm(Pos ~ Grupo + Pre + Grupo:Pre, 
             data=Dados.wide)
print(Anova <- car::Anova(modelo))

## ANCOVA: Teste do efeito do fator fixo: 
##          Se as declividades sao iguais, testar se os interceptos sao iguais.
cat(bartitle("Homogenous intercepts of regressions\n(assuming homogeneous slopes)\n- Factor Effect Test -",3))
ancova.fit <- lm(Pos ~ Grupo + Pre, 
                 data=Dados.wide)
print(Anova <- car::Anova(ancova.fit))
print(summary(ancova.fit))

cat(bartitle("Effect size analysis"))
cat(bartitle("Omnibus",2))
eta2 <- summary(ancova.fit)$r.squared
cat("R^2 = eta^2 omnibus = ", round(eta2,2), "\n", sep="")
print(effectsize::interpret_eta_squared(eta2))
eta2 <- effectsize::eta_squared(Anova,
                                partial = TRUE,
                                generalized = FALSE,
                                ci = 1-alfa/2,
                                alternative = "two.sided",
                                verbose = TRUE)
print(eta2, digits = 2)
es <- effectsize::interpret_eta_squared(eta2$Eta2)
names(es) <- c("Tamanho de efeito: estimativa pontual")
print(es)

cat(bartitle("Adjusted estimated marginal means"))

EMM <- emmeans::emmeans(ancova.fit, 
                        pairwise~Grupo, 
                        adjust="holm",
                        level=1-alfa)
print(EMM$emmeans)
print(plot(EMM$emmeans,
           colors="black", 
           main="Adjusted Estimated Marginal Means",
           xlab="VD Pós controlada por VD Pré",
           xlim=c(4,12),
           ylab="Grupo")+ggplot2::theme_bw())

cat(bartitle("Post hoc tests: Grupo (factor) adjusted by VD Pre (covariate)"))

print(summary(EMM$contrasts, infer=TRUE))
print(plot(EMM$contrasts, 
           colors="black")+ggplot2::theme_bw())
print(multcomp::cld(object=EMM$emmeans,
                    level=1-alfa,
                    adjust="holm",
                    Letters=letters,
                    alpha=alfa))

----
Data
----
   Paciente   Grupo Momento Escore
1         1 droga A     Pre    150
2         2 droga A     Pre    130
3         3 droga A     Pre    125
4         4 droga A     Pre    152
5         5 droga A     Pre    160
6         6 droga B     Pre    174
7         7 droga B     Pre    110
8         8 droga B     Pre    180
9         9 droga B     Pre    145
10       10 droga B     Pre    140
11        1 droga A     Pos     51
12        2 droga A     Pos     50
13        3 droga A     Pos     40
14        4 droga A     Pos     45
15        5 droga A     Pos     60
16        6 droga B     Pos     75
17        7 droga B     Pos     41
18        8 droga B     Pos     80
19        9 droga B     Pos     60
20       10 droga B     Pos     55

----------------------
Descriptive Statistics
----------------------
                Paciente 1 2 3 4 5 6 7 8 9 10
Grupo   Momento                              
droga A Pre              1 1 1 1 1 0 0 0 0  0
        Pos              1 1 1 1 1 0 0 0 0  0
droga B Pre              0 0 0 0 0 1 1 1 1  1
        Pos              0 0 0 0 0 1 1 1 1  1
        item group1  group2 vars n  mean    sd median trimmed   mad min max
Escore1    1    Pre droga A    1 5 143.4 15.09    150   143.4 14.83 125 160
Escore2    2    Pos droga A    1 5  49.2  7.46     50    49.2  7.41  40  60
Escore3    3    Pre droga B    1 5 149.8 28.29    145   149.8 43.00 110 180
Escore4    4    Pos droga B    1 5  62.2 15.71     60    62.2 22.24  41  80
        range  skew kurtosis    se
Escore1    35 -0.18    -2.11  6.75
Escore2    20  0.20    -1.64  3.34
Escore3    70 -0.20    -1.85 12.65
Escore4    39 -0.11    -1.91  7.02

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.


-------------------------
Simple linear regressions
-------------------------

`geom_smooth()` using formula = 'y ~ x'


        -------
        droga A
        -------

Call:
lm(formula = Pos ~ Pre, data = dt_tmp)

Residuals:
      1       2       3       4       5 
-0.4715  5.4118 -2.8674 -7.1598  5.0869 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -0.1528    29.5259  -0.005    0.996
Pre           0.3442     0.2050   1.679    0.192

Residual standard error: 6.188 on 3 degrees of freedom
Multiple R-squared:  0.4844,    Adjusted R-squared:  0.3126 
F-statistic: 2.819 on 1 and 3 DF,  p-value: 0.1918


        -------
        droga B
        -------

Call:
lm(formula = Pos ~ Pre, data = dt_tmp)

Residuals:
      6       7       8       9      10 
-0.5989  0.8362  1.0790  0.4576 -1.7740 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -20.7404     3.6618  -5.664 0.010900 *  
Pre           0.5537     0.0241  22.971 0.000181 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.364 on 3 degrees of freedom
Multiple R-squared:  0.9943,    Adjusted R-squared:  0.9925 
F-statistic: 527.7 on 1 and 3 DF,  p-value: 0.0001807


------
ANCOVA
------

    -----------
    Assumptions
    -----------

        ----------------------------------------------------------
        Dissociation between Grupo (factor) and VD Pre (covariate)
        ----------------------------------------------------------
Anova Table (Type II tests)

Response: Pre
          Sum Sq Df F value Pr(>F)
Grupo      102.4  1  0.1992 0.6672
Residuals 4112.0  8               

        ------------------------------------------------------------
        Linearity between Escore pos (VD) and Escore pre (covariate)
        ------------------------------------------------------------

            -------
            droga A
            -------

    Harvey-Collier test

data:  modelo
HC = 0.30952, df = 2, p-value = 0.7862


            -------
            droga B
            -------

    Harvey-Collier test

data:  modelo
HC = 0.29829, df = 2, p-value = 0.7936


        -------------------------------
        Homogenous slope of regressions
        -------------------------------
Anova Table (Type II tests)

Response: Pos
           Sum Sq Df F value    Pr(>F)    
Grupo      232.05  1 11.5594 0.0144976 *  
Pre       1058.02  1 52.7033 0.0003473 ***
Grupo:Pre   31.13  1  1.5509 0.2594306    
Residuals  120.45  6                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

        ------------------------------------
        Homogenous intercepts of regressions
(assuming homogeneous slopes)
- Factor Effect Test -
        ------------------------------------
Anova Table (Type II tests)

Response: Pos
           Sum Sq Df F value    Pr(>F)    
Grupo      232.05  1  10.716 0.0136078 *  
Pre       1058.02  1  48.858 0.0002135 ***
Residuals  151.58  7                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Call:
lm(formula = Pos ~ Grupo + Pre, data = Dados.wide)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.5623 -1.4138  0.1841  1.9159  7.5971 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -23.53923   10.61245  -2.218 0.062053 .  
Grupodroga B   9.75362    2.97954   3.274 0.013608 *  
Pre            0.50725    0.07257   6.990 0.000213 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.653 on 7 degrees of freedom
Multiple R-squared:  0.9071,    Adjusted R-squared:  0.8806 
F-statistic: 34.18 on 2 and 7 DF,  p-value: 0.0002442


--------------------
Effect size analysis
--------------------

    -------
    Omnibus
    -------
R^2 = eta^2 omnibus = 0.91
[1] "large"
(Rules: field2013)
# Effect Size for ANOVA (Type II)

Parameter | Eta2 (partial) |     97.5% CI
-----------------------------------------
Grupo     |           0.60 | [0.00, 0.85]
Pre       |           0.87 | [0.44, 0.95]
Tamanho de efeito: estimativa pontual                                  <NA> 
                              "large"                               "large" 
(Rules: field2013)

---------------------------------
Adjusted estimated marginal means
---------------------------------
 Grupo   emmean   SE df lower.CL upper.CL
 droga A   50.8 2.09  7     45.9     55.8
 droga B   60.6 2.09  7     55.6     65.5

Confidence level used: 0.95 


-------------------------------------------------------------
Post hoc tests: Grupo (factor) adjusted by VD Pre (covariate)
-------------------------------------------------------------
 contrast          estimate   SE df lower.CL upper.CL t.ratio p.value
 droga A - droga B    -9.75 2.98  7    -16.8    -2.71  -3.274  0.0136

Confidence level used: 0.95 

 Grupo   emmean   SE df lower.CL upper.CL .group
 droga A   50.8 2.09  7     44.9     56.8  a    
 droga B   60.6 2.09  7     54.6     66.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. 

A interpretação da saída é a mesma da ANOVA independente. Neste exemplo, chamamos a atenção para algumas ocorrências:

  1. Nas regressões lineares simples, a droga \(A\) não tem inclinação significante para \(\alpha=0.05\). No entanto, o teste da homogeneidade das inclinações não rejeita a hipótese nula de paralelismo entre as regressões. Observe o gráfico correspondente: há um segmento de reta que pode ser traçado pertencendo simultaneamente às bandas de confiança de 95% com correção de Bonferroni no domínio comum da covariável.

  2. Como acontece na ANCOVA independente, existe a saída da tabela ANOVA da ANCOVA (a ANCOVA propriamente dita) e a saída da regressão da ANCOVA.

Na saída da tabela ANOVA da ANCOVA, avalia-se o valor p associado ao fator Grupo (droga), que neste caso é 0.0136. Para \(\alpha=0.05\), rejeita-se a hipótese nula de igualdade entre as médias ajustadas, concluindo que as drogas A e B têm efeitos distintos, sendo o efeito da droga B maior.

Anova Table (Type II tests)

Response: Pos
           Sum Sq Df F value    Pr(>F)    
Grupo      232.05  1  10.716 0.0136078 *  
Pre       1058.02  1  48.858 0.0002135 ***
Residuals  151.58  7                      

Na saída da regressão, o teste F global (\(F=34.18\), \(p=0.00024\)) indica que o modelo ANCOVA explica uma parcela expressiva da variabilidade da resposta. O fator Grupo tem efeito significante após o ajuste pela covariável Pre. O coeficiente de determinação (\(R^2=0.91\)) indica que 91% da variância total de Pos é explicada pelo modelo, correspondendo a um tamanho de efeito global elevado.

Multiple R-squared:  0.9071,    Adjusted R-squared:  0.8806 
F-statistic: 34.18 on 2 and 7 DF,  p-value: 0.0002442
  1. Como há somente dois níveis do fator, não haveria necessidade do teste post hoc que meramente deveria concordar com o que obtém da ANCOVA (está disponível neste código R para o caso do código ser adaptado para situações com mais níveis do fator).

Planejamento

Segundo Borm et al. (2007),

“O uso da ANCOVA pode reduzir consideravelmente o número de pacientes necessários para um ensaio. […] Uma ANCOVA com \((1 - \rho^2)n\) observações tem o mesmo poder que um teste t com \(n\) observações.”

Estes autores propõem uma fórmula simples para dimensionar uma amostra para a ANCOVA para um fator com dois níveis:

delta <- .5
rho <- .6
alfa <- .05
poder <- .80
n.porcondicao <- 2*((qnorm(1-alfa/2)+qnorm(poder))*sqrt(1-rho^2)/delta)^2
cat("n total = ",2*ceiling(n.porcondicao)," com ",ceiling(n.porcondicao)," por nível do fator\n", sep="")
n total = 82 com 41 por nível do fator

Neste código R, encontramos o tamanho da amostra para \(\alpha=0.05\) e \(\text{poder}=0.80\). Temos:

  • delta: \(d\) de Cohen populacional (tamanho de efeito populacional)
  • rho: correlação populacional de Pearson entre a VD e a covariável

Segundo Liu (2011), embora a inclusão da covariável seja geralmente aceita como capaz de reduzir o erro-padrão ajustado, consequentemente aumentando o poder da ANCOVA em comparação com a ANOVA (aumentando o valor de \(\eta^2\)), isto nem sempre se verifica. Com o aumento do tamanho da amostra, a probabilidade de haver benefício com a inclusão de uma covariável também aumenta. O autor apresenta a seguinte tabela:

Neste estudo, foram supostos dois níveis para o fator. A correlação de Pearson entre a VD e a covariável é \(0.3\) (intermediária). Há tamanho total da amostra (\(n\)) e graus de liberdade correspondentes da ANCOVA (\(\nu\)). As colunas \(p_1\) e \(p_2\) são probabilidades (calculadas por dois métodos diversos) da ANCOVA ter resultado melhor que a ANOVA nas mesmas condições: \(p_1\) é a probabilidade de reduzir o erro-padrão ajustado; \(p_2\) é probabilidade de reduzir o amplitude do intervalo de confiança ajustado. Observe que para uma probabilidade de 80% e 90% de haver benefício com a inclusão da covariável necessita-se de \(n = 16\) e \(n = 22\).

Para quem quiser modelo mais complexo, Crawley (2012) em “Chapter12 - Analysis of covariance” apresenta um delineamento com dois fatores e uma covariável.

Visualização gráfica de ANCOVA

Esse gráfico 3D ilustra de forma intuitiva o resultado de uma análise de covariância (ANCOVA) aplicada a medidas pré e pós em dois grupos, Placebo e Tratamento.

O plano azul com grade é o plano de regressão ajustado pelo modelo sob a hipótese nula, em que se assume que o efeito da variável pré sobre a variável pós é o mesmo para ambos os grupos, variando apenas o nível médio. Esse plano representa o comportamento médio previsto pelo modelo.

As duas superfícies translúcidas acima e abaixo do plano azul formam a banda de confiança simultânea de 95%. Elas indicam a região onde, com 95% de confiança, espera-se que o verdadeiro plano médio da população esteja. Quanto mais estreitas essas superfícies, maior é a precisão das estimativas do modelo.

Os pontos pretos são as observações realizadas, mostrando os valores observados de pós em função de pré e do grupo. Quanto mais próximos esses pontos estão do plano azul, melhor é o ajuste do modelo.

O plano laranja é o plano bissetriz, que representa a situação de ausência de mudança: quando os valores pós seriam iguais aos valores pré (pós = pré). Visualmente, se o plano azul estiver acima da bissetriz, há tendência de aumento; se estiver abaixo, tendência de redução.

Assim, o gráfico permite observar simultaneamente três elementos fundamentais: o comportamento médio ajustado (plano azul), a incerteza associada ao modelo (bandas de confiança) e a linha teórica de referência sem efeito (plano laranja).


Call:
lm(formula = Pos ~ Grupo.num + Pre, data = Dados)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.9397 -3.5054  0.1912  2.9852 13.6825 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.2242     3.1922   1.637 0.123996    
Grupo.num    14.5087     3.2732   4.433 0.000568 ***
Pre           0.5094     0.1595   3.194 0.006501 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.976 on 14 degrees of freedom
Multiple R-squared:  0.7883,    Adjusted R-squared:  0.7581 
F-statistic: 26.07 on 2 and 14 DF,  p-value: 1.904e-05
wgl 
  1 

ANCOVA sem dados brutos

Quando apenas estatísticas resumidas estão disponíveis (médias, tamanhos \(n_i\), desvios-padrão intragrupo de \(y\) e \(x\) por grupo de VD \(y\) e covariável \(x\) e a correlação \(r\)), ainda é possível testar a hipótese de igualdade das médias ajustadas entre grupos na ANCOVA. O teste \(F\) informa se, após ajustar \(y\) por \(x\), há evidência de diferenças entre as médias ajustadas dos grupos.

Em meta‑análises, reanálises de artigos que publicam apenas resumos, auditorias e reprodutibilidade, o teste \(F\) da ANCOVA obtido de resumos permite recuperar a inferência principal do efeito de grupo ajustado por \(x\) sem acesso aos dados brutos. Serve também para estudos de poder aproximado e checagens de consistência entre relatórios.

set.seed(123)

# --- parâmetros-alvo ---
n    <- c(20,22,18)
xbar <- c(5.0,6.0,5.5)
ybar <- c(10.0,12.0,11.0)
sxw  <- 2.0
syw  <- 3.0
rw   <- 0.60

# --- simulação coerente com momentos ---
Sigma <- matrix(c(sxw^2, rw*sxw*syw,
                  rw*sxw*syw, syw^2), 2, 2)
simula <- function(ni, mux, muy, Sigma){
  XY <- MASS::mvrnorm(ni, mu=c(mux,muy), Sigma=Sigma)
  data.frame(X=XY[,1], Y=XY[,2])
}
g <- length(n)
dados <- do.call(rbind, lapply(seq_len(g), function(i) simula(n[i], xbar[i], ybar[i], Sigma)))
dados$Grupo <- factor(rep(paste0("G",1:g), times=n))

# --- ANCOVA teórica (estatísticas-resumo) ---
ANCOVA_WORD.test <- function(n, ybar, xbar, sy, sx, r, conf = 0.95) {
  N <- sum(n); g <- length(n)
  beta <- r * (sy / sx)
  xdot <- weighted.mean(xbar, n)
  y_tilde <- ybar - beta * (xbar - xdot)
  y_tilde_dot <- weighted.mean(y_tilde, n)
  
  SSg <- sum(n * (y_tilde - y_tilde_dot)^2)
  SSw <- (N - g) * sy^2
  SSx <- r^2 * SSw
  SSE <- (1 - r^2) * SSw
  
  df_g <- g - 1
  df_x <- 1
  df_e <- N - g - 1
  MS_e <- SSE / df_e
  
  Fg <- (SSg / df_g) / MS_e
  Fx <- (SSx / df_x) / MS_e
  pg <- pf(Fg, df_g, df_e, lower.tail = FALSE)
  px <- pf(Fx, df_x, df_e, lower.tail = FALSE)
  
  # Função auxiliar para IC aproximado de η² parcial
  eta2_ci <- function(Fv, df1, df2, conf) {
    alpha <- 1 - conf
    F_lower <- qf(alpha/2, df1, df2)
    F_upper <- qf(1 - alpha/2, df1, df2)
    est <- (df1 * Fv) / (df1 * Fv + df2)
    lo  <- (df1 * F_lower) / (df1 * F_lower + df2)
    hi  <- (df1 * F_upper) / (df1 * F_upper + df2)
    c(est, lo, hi)
  }
  
  e_g <- eta2_ci(Fg, df_g, df_e, conf)
  e_x <- eta2_ci(Fx, df_x, df_e, conf)[1]  # apenas a estimativa, sem IC
  
  # R² e IC95%
  R2 <- (SSg + SSx) / (SSg + SSx + SSE)
  K  <- df_g + df_x
  alpha <- 1 - conf
  F_lower <- qf(alpha/2, K, df_e)
  F_upper <- qf(1 - alpha/2, K, df_e)
  R2_lo <- (K * F_lower) / (K * F_lower + df_e)
  R2_hi <- (K * F_upper) / (K * F_upper + df_e)
  
  tab <- data.frame(
    check.names = FALSE,
    "Sum Sq"  = c(SSg, SSx, SSE),
    "Df"      = c(df_g, df_x, df_e),
    "F value" = c(Fg, Fx, NA),
    "Pr(>F)"  = c(pg, px, NA),
    "η² parcial" = c(e_g[1], e_x, NA),
    "IC95% η² parcial (Grupo)" = c(
      sprintf("[%.3f, %.3f]", e_g[2], e_g[3]),
      NA,
      NA
    ),
    row.names = c("Group","Covariate","Residuals")
  )
  attr(tab, "R2") <- R2
  attr(tab, "R2_CI") <- c(R2_lo, R2_hi)
  class(tab) <- c("anova","data.frame")
  tab
}

tab_sum <- ANCOVA_WORD.test(n, ybar, xbar, syw, sxw, rw)
print(as.data.frame(tab_sum), digits = 4)
          Sum Sq Df F value    Pr(>F) η² parcial IC95% η² parcial (Grupo)
Group      12.68  2   1.082 3.460e-01     0.0372           [0.001, 0.123]
Covariate 184.68  1  31.500 6.432e-07     0.3600                     <NA>
Residuals 328.32 56      NA        NA         NA                     <NA>
cat(sprintf("\nR² global = %.3f; IC95%% R² = [%.3f, %.3f]\n",
            attr(tab_sum, "R2"),
            attr(tab_sum, "R2_CI")[1],
            attr(tab_sum, "R2_CI")[2]))

R² global = 0.375; IC95% R² = [0.004, 0.153]
# --- ANCOVA empírica via car::Anova ---
fit <- lm(Y ~ Grupo + X, data=dados)
print(tab_car <- car::Anova(fit))
Anova Table (Type II tests)

Response: Y
           Sum Sq Df F value    Pr(>F)    
Grupo       9.107  2  1.0373    0.3611    
X         189.742  1 43.2231 1.743e-08 ***
Residuals 245.830 56                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(summary(fit))

Call:
lm(formula = Y ~ Grupo + X, data = dados)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.6566 -1.4300 -0.1477  1.3843  4.2957 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.0514     0.9361   5.396 1.43e-06 ***
GrupoG2       0.8784     0.6636   1.324    0.191    
GrupoG3       0.1523     0.7069   0.216    0.830    
X             1.0086     0.1534   6.574 1.74e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.095 on 56 degrees of freedom
Multiple R-squared:  0.4807,    Adjusted R-squared:  0.4529 
F-statistic: 17.28 on 3 and 56 DF,  p-value: 4.586e-08
alfa <- 0.05
eta2 <- effectsize::eta_squared(tab_car,
                                partial=TRUE,
                                generalized=FALSE,
                                ci=1-alfa,
                                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) |           95% CI | interpret
---------------------------------------------------------
Grupo     |         0.0357 | [0.0000, 0.1526] |     small
X         |         0.4356 | [0.2436, 0.5836] |     large

Delineamento de Solomon

Conforme Braver & Braver (1988)

Quase quarenta anos atrás, Solomon introduziu uma nova forma de delineamento experimental, tipicamente referida hoje como o delineamento de quatro grupos de Solomon (Solomon, 1949). Campbell e Stanley (1963) discutiram esse delineamento como um dos três delineamentos experimentais com uma única condição de tratamento, sendo os outros dois o delineamento com grupo controle pré e pós-teste e o delineamento com grupo controle apenas pós-teste (ver Tabela 1). Cada um desses delineamentos é adequado para avaliar o efeito do tratamento e é imune à maioria das ameaças à validade interna.

O delineamento de quatro grupos de Solomon, entretanto, acrescenta a vantagem de ser o único dos três capaz de avaliar a presença de sensibilização ao pré-teste. Sensibilização ao pré-teste significa que “a exposição ao pré-teste aumenta a sensibilidade dos sujeitos ao tratamento experimental, impedindo assim a generalização dos resultados do grupo pré-testado para uma população não pré-testada” (Huck & Sandler, 1973, p. 54). Assim, o delineamento de quatro grupos de Solomon acrescenta um grau mais elevado de validade externa além de sua validade interna e, portanto, segundo Helmstadter (1970), é “o mais desejável de todos os delineamentos experimentais básicos” (p. 110).

Apesar de sua força, o delineamento de quatro grupos de Solomon é pouco utilizado. Quatro razões podem explicar esse subuso.

Primeiro, o delineamento requer o dobro do número de grupos em relação aos outros dois delineamentos, exigência muitas vezes mal interpretada como significando que o delineamento de Solomon requer também o dobro do número de sujeitos. Mostraremos mais adiante neste artigo que, reduzindo-se pela metade o tamanho de cada grupo, mantendo o mesmo tamanho total da amostra que os outros delineamentos, ainda é possível usufruir dos benefícios do delineamento de Solomon. Essa estratégia normalmente resulta em poder estatístico bastante adequado; de fato, o poder será maior do que o do delineamento com grupo controle apenas pós-teste.

Segundo, poucos pesquisadores têm interesse em estudar os efeitos da sensibilização ao pré-teste em si, e portanto objetam o uso de um delineamento cuja principal vantagem é justamente examinar tal sensibilização. Essa objeção perde o ponto essencial. A sensibilização ao pré-teste, assim como o efeito do experimentador, é um potencial artefato que pode limitar a generalização dos efeitos que de fato interessam ao pesquisador, e sua existência deveria ser investigada mesmo que não haja interesse direto nesse fenômeno.

Mais diretamente, pode haver a crença de que artefatos de sensibilização ao pré-teste provavelmente não existem na área de pesquisa do investigador. Essa crença é sustentada por várias revisões de literatura (Bracht & Glass, 1968; Lana, 1959; Lana, 1969; Rosnow, 1971; Willson & Putnam, 1982), mostrando que artefatos de sensibilização ao pré-teste raramente ocorrem. Ainda assim, esse artefato deve ser considerado como um efeito que pode potencialmente comprometer a validade externa de qualquer achado de pesquisa até que, e a menos que, o delineamento de Solomon seja usado para descartá-lo explicitamente.

Terceiro, como argumentaram Oliver e Berger (1980), as conclusões podem se tornar muito mais complexas ao se usar o delineamento de Solomon, devido ao número de comparações que ele permite. Essa complexidade pode desmotivar pesquisadores. Embora essa objeção seja compreensível, a complexidade de um fenômeno ou de uma análise não é razão científica válida para deixar de conduzi-la.

Quarto, e talvez o motivo mais importante para o subuso do delineamento de Solomon, é a falta de certeza sobre o tratamento estatístico apropriado para esse delineamento relativamente complexo. Campbell e Stanley (1963) fizeram algumas recomendações, mas sua apresentação foi incompleta. Huck e Sandler (1973) modificaram a análise de Campbell e Stanley, mas ainda deixaram alguns detalhes sem tratamento.

O restante deste artigo é dedicado à questão do tratamento estatístico adequado do delineamento de Solomon. A análise que recomendamos está de acordo com o trabalho de Campbell e Stanley e de Huck e Sandler, mas cobre mais contingências e apresenta maior poder estatístico (isto é, maior capacidade de detectar significância). Por uma questão de completude, toda a análise é apresentada aqui, embora as primeiras etapas já tenham sido originalmente apresentadas pelos autores anteriores.

Como apontaram Campbell & Stanley (1963), a primeira fase da análise é determinar se existe evidência de sensibilização ao pré-teste, isto é, se \(X\) afeta \(O\) apenas quando uma medida de pré-teste também é administrada. se este for o caso, \(O_2\) seria maior do que \(O_4\), mas \(O_5\) não seria maior do que \(O_6\). o teste para isso é uma análise de variância (ANOVA) 2 × 2 entre grupos sobre os quatro escores de pós-teste, conforme indicado na Tabela 2. os fatores são tratamento (sim versus não) e pré-teste (sim versus não). a evidência que demonstra sensibilização ao pré-teste é detectada pela interação (doravante denominada teste A). além disso, deve haver um efeito simples significativo de tratamento na primeira linha (pré-teste presente — teste B), mas não na segunda (pré-teste ausente — teste C). se o padrão precedente estiver presente, a análise termina com a conclusão de que há evidência de um efeito do tratamento, mas ele ocorre apenas para os grupos pré-testados; há, portanto, sensibilização ao pré-teste (um resultado improvável de agradar ao investigador).

Tabela 1. delineamentos experimentais com uma condição de tratamento

delineamento grupo pré-teste tratamento pós-teste
Solomon quatro grupos 1 \(O_1\) \(X\) \(O_2\)
2 \(O_3\) \(O_4\)
3 \(X\) \(O_5\)
4 \(O_6\)
com grupo controle & pré e pós-teste 1 \(O_1\) \(X\) \(O_2\)
2 \(O_3\) \(O_4\)
com grupo controle & apenas pós-teste 1 \(X\) \(O_5\)
2 \(O_6\)

nota. \(O\) = medida de desfecho; \(X\) = tratamento; R = randomização em todas as linhas.

Tabela 2. análise 2 × 2 dos escores de pós-teste

pré-teste tratamento (\(X\)) sim tratamento (\(X\)) não
sim \(O_2\) \(O_4\)
não \(O_5\) \(O_6\)

nota. \(O\) = medida de desfecho.

  1. Executar uma ANOVA 2×2 entre grupos nos quatro escores de pós-teste (O₂, O₄, O₅, O₆).
  • Teste A = interação entre tratamento (sim/não) e pré-teste (sim/não).
  1. Se a interação for significativa, então:
  1. Realizar o Teste B, efeito simples de tratamento nos grupos com pré-teste (Grupos 1 e 2).
  2. Realizar o Teste C, efeito simples de tratamento nos grupos sem pré-teste (Grupos 3 e 4). * Se apenas B for significativo, concluir: “O tratamento só tem efeito quando há pré-teste; existe sensibilização ao pré-teste.” * Se B e C forem significativos, concluir: “O tratamento tem efeito mesmo em grupos sem pré-teste.”
  1. Se a interação não for significativa, então realizar o Teste D, efeito principal de tratamento (média geral entre grupos com e sem pré-teste).
  • Se significativo, concluir: “O tratamento tem efeito — nenhuma qualificação necessária.”
  • Se não for significativo, prosseguir com testes mais potentes:
  1. Rodar o Teste E (ANCOVA) ou Teste F (análise de ganhos) ou Teste G (ANOVA de medidas repetidas) em Grupos 1 e 2.
  • Se significativo, concluir efeito de tratamento.
  • Se não, realizar o Teste H (teste t independente em Grupos 3 e 4).
  • Se ainda não for significativo, combinar os resultados dos Testes E/F/G e H por meta-análise (Teste I).
  1. Se a meta-análise (Teste I) for significativa, concluir efeito de tratamento; caso contrário, concluir ausência de evidência de efeito.

Exemplo: idoso com hipertensão

Pacientes idosos do sexo masculino com hipertensão estágio 2 são acompanhados em ambulatório de clínica médica. Avalia-se a eficácia de um programa supervisionado de exercícios respiratórios diafragmáticos, três sessões semanais por 12 semanas, adicionado ao cuidado padrão. O desfecho primário é a pressão arterial sistólica no pós-teste (mmHg). Em \(P=1\) mede-se \(Y_{\mathrm{pre}}\) e \(Y_{\mathrm{post}}\); em \(P=0\) mede-se apenas \(Y_{\mathrm{post}}\). Adota-se delineamento de Solomon com quatro grupos e randomização 1:1:1:1.

ANCOVA generalizada com interação tratamento:pré-teste

Modelo para testar sensibilização \(\beta_{TP}\) e efeito do tratamento \(\beta_T\):

\[ Y_{i,\mathrm{post}} = \theta_0 + \theta_T T_i + \theta_P P_i + \theta_{TP}\,T_iP_i + \gamma\,P_i\big(Y_{i,\mathrm{pre}}-\bar Y_{\mathrm{pre}}\big) + \varepsilon_i \]

A interação \(\theta_{TP}\) testa sensibilização; o efeito ajustado de tratamento é \(\theta_T\) quando \(\theta_{TP}= 0\).

Explicação da simulação:

Os parâmetros foram escolhidos para refletir um cenário clínico plausível em hipertensão estágio 2. A média basal foi fixada em \(\mu_{\mathrm{pre}}=160\) mmHg com desvio padrão \(\sigma_{\mathrm{pre}}=14\) mmHg, representando variabilidade típica ambulatorial. O erro no pós-teste foi simulado com \(\sigma_{\varepsilon}=12\) mmHg, que captura variação não explicada por tratamento e baseline. O efeito verdadeiro do tratamento foi definido como \(\tau=-20\) mmHg, refletindo uma redução clinicamente relevante com exercícios respiratórios supervisionados. Assumiu-se ausência de sensibilização do pré-teste (isto é, efeito de tratamento constante entre \(P=1\) e \(P=0\)), o que serve de base para verificar se a análise recupera o valor próximo a \(-20\) mmHg. Os tamanhos amostrais foram iguais em cada grupo (\(n_{11}=n_{01}=n_{10}=n_{00}=60\)) para simplicidade e boa potência, mantendo balanceamento típico de ensaios randomizados. As distribuições normais para \(Y_{\mathrm{pre}}\) e os erros de \(Y_{\mathrm{post}}\) são aproximações usuais em variáveis fisiológicas contínuas; ruídos independentes entre indivíduos foram adotados para focar no contraste entre tratamento e controle.

n11 <- n01 <- n10 <- n00 <- 60L
mu_pre <- 160
sd_pre <- 14
sd_eps <- 12
tau <- -20  # redução média alvo de 20 mmHg

pre_11 <- rnorm(n11, mu_pre, sd_pre)
pre_01 <- rnorm(n01, mu_pre, sd_pre)

post_11 <- pre_11 + tau + rnorm(n11, 0, sd_eps)
post_01 <- pre_01 + 0   + rnorm(n01, 0, sd_eps)
post_10 <- mu_pre + tau + rnorm(n10, 0, sd_eps)
post_00 <- mu_pre + 0   + rnorm(n00, 0, sd_eps)

df <- rbind(
  data.frame(t=1L,p=1L,y_pre=pre_11,y_post=post_11),
  data.frame(t=0L,p=1L,y_pre=pre_01,y_post=post_01),
  data.frame(t=1L,p=0L,y_pre=NA_real_,y_post=post_10),
  data.frame(t=0L,p=0L,y_pre=NA_real_,y_post=post_00)
)
row.names(df) <- NULL
head(df)
# A tibble: 6 × 4
      t     p y_pre y_post
  <int> <int> <dbl>  <dbl>
1     1     1  162.   132.
2     1     1  147.   121.
3     1     1  153.   151.
4     1     1  156.   123.
5     1     1  186.   164.
6     1     1  151.   154.
tail(df)
# A tibble: 6 × 4
      t     p y_pre y_post
  <int> <int> <dbl>  <dbl>
1     0     0    NA   157.
2     0     0    NA   162.
3     0     0    NA   181.
4     0     0    NA   156.
5     0     0    NA   164.
6     0     0    NA   157.
ybar_pre_all <- mean(df$y_pre, na.rm = TRUE)
df$y_pre_c <- ifelse(df$p==1, df$y_pre - ybar_pre_all, 0)
fit <- lm(y_post ~ t + p + t:p + I(p*y_pre_c), data=df)
print(anv <- car::Anova(fit))
Anova Table (Type II tests)

Response: y_post
               Sum Sq  Df F value Pr(>F)    
t               24708   1 167.888 <2e-16 ***
p                   9   1   0.061 0.8052    
I(p * y_pre_c)  16922   1 114.985 <2e-16 ***
t:p               509   1   3.458 0.0642 .  
Residuals       34585 235                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(summary(fit))

Call:
lm(formula = y_post ~ t + p + t:p + I(p * y_pre_c), data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-33.817  -6.589   0.139   8.612  30.321 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    161.88563    1.56615  103.37   <2e-16 ***
t              -23.20572    2.21486  -10.48   <2e-16 ***
p               -2.52580    2.21490   -1.14   0.2553    
I(p * y_pre_c)   0.86609    0.08077   10.72   <2e-16 ***
t:p              5.82493    3.13239    1.86   0.0642 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.13 on 235 degrees of freedom
Multiple R-squared:  0.5476,    Adjusted R-squared:  0.5399 
F-statistic: 71.12 on 4 and 235 DF,  p-value: < 2.2e-16
eta2 <- summary(fit)$r.squared
cat("R^2 = eta^2 omnibus = ", round(eta2,2), "\n", sep="")
R^2 = eta^2 omnibus = 0.55
print(effectsize::interpret_eta_squared(eta2))
[1] "large"
(Rules: field2013)
alfa <- 0.05
eta2 <- effectsize::eta_squared(anv,
                                partial=TRUE,
                                generalized=FALSE,
                                ci=1-alfa,
                                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) |           95% CI |  interpret
---------------------------------------------------------------
t              |         0.4167 | [0.3261, 0.4967] |      large
p              |         0.0003 | [0.0000, 0.0175] | very small
I(p * y_pre_c) |         0.3285 | [0.2363, 0.4146] |      large
t:p            |         0.0145 | [0.0000, 0.0586] |      small
# ANCOVA
fit <- lm(y_post ~ factor(t) + y_pre_c, data=df)
print(anv <- car::Anova(fit))
Anova Table (Type II tests)

Response: y_post
          Sum Sq  Df F value    Pr(>F)    
factor(t)  24708   1  166.82 < 2.2e-16 ***
y_pre_c    16971   1  114.58 < 2.2e-16 ***
Residuals  35103 237                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(summary(fit))

Call:
lm(formula = y_post ~ factor(t) + y_pre_c, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-32.554  -7.063   0.122   8.044  28.672 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 160.62282    1.11099  144.58   <2e-16 ***
factor(t)1  -20.29345    1.57121  -12.92   <2e-16 ***
y_pre_c       0.86732    0.08102   10.70   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.17 on 237 degrees of freedom
Multiple R-squared:  0.5408,    Adjusted R-squared:  0.537 
F-statistic: 139.6 on 2 and 237 DF,  p-value: < 2.2e-16
eta2 <- summary(fit)$r.squared
cat("R^2 = eta^2 omnibus = ", round(eta2,2), "\n", sep="")
R^2 = eta^2 omnibus = 0.54
print(effectsize::interpret_eta_squared(eta2))
[1] "large"
(Rules: field2013)
alfa <- 0.05
eta2 <- effectsize::eta_squared(anv,
                                partial=TRUE,
                                generalized=FALSE,
                                ci=1-alfa,
                                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) |           95% CI | interpret
---------------------------------------------------------
factor(t) |         0.4131 | [0.3227, 0.4931] |     large
y_pre_c   |         0.3259 | [0.2341, 0.4117] |     large

Referências

  • Borm et al. (2007) A simple sample size formula for analysis of covariance in randomized clinical trials. Journal of Clinical Epidemiology 60: 1234-8.
  • Braver, MCW & Braver, SL (1988) Statistical treatment of the Solomon four-group design: A meta-analytic approach. Psychological Bulletin 104(1): 150-4. https://doi.org/10.1037/0033-2909.104.1.150
  • Caldwell, AR et al. (2022) Power Analysis with Superpower. https://aaroncaldwell.us/SuperpowerBook/
  • Cochran, WG (1957) Analysis of Covariance: Its Nature and Uses. Biometrics 13(3): 261–281. https://doi.org/10.2307/2527916
  • Dancey CP & Reidy J (2019) Estatística sem Matemática para Psicologia. 7ª ed. Porto Alegre: Penso.
  • Crawley, MJ (2013) The R Book. 2nd ed. NJ: Wiley.
  • Sakae TM et al. (2016) Efeitos da anestesia geral da cognição do idoso. Arquivos Catarinenses de Medicina 45(3): 107-16.
  • Shieh, G (2020) Power Analysis and Sample Size Planning in ANCOVA Designs. Psychometrika, 85(1): 101–20. doi:10.1007/s11336-019-09692-3
  • Liu, X (2011) The Effect of a Covariate on Standard Error and Confidence Interval Width. Communications in Statistics - Theory and Methods 40(3): 449-56, DOI: 10.1080/03610920903391337