Análise de Covariância (ANCOVA)
em R
02 julho 2025 21:01h
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))options(warn=-1)
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts=FALSE))
suppressMessages(library(emmeans, warn.conflicts=FALSE))
suppressMessages(library(multcomp, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(GGally, warn.conflicts=FALSE))
suppressMessages(library(effects, warn.conflicts=FALSE))
suppressMessages(library(effectsize, warn.conflicts=FALSE))
suppressMessages(library(tidyr, warn.conflicts=FALSE))
suppressMessages(library(MASS, warn.conflicts=FALSE))
suppressMessages(library(car, warn.conflicts=FALSE))
source("eiras.bartitle.R")Material
- HTML de R Markdown 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 intervenção.
- 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 das suposições e do efeito do fator.
ANCOVA independente simples
ANCOVA independente simples é uma combinação de regressão linear simples com ANOVA unifatorial independente. Este é o primeiro passo para o modelo linear geral (GLM: General Linear Model).
Na ANOVA, o tamanho de efeito global, \(R^2\), é a porcentagem da variância da variável de desfecho (VD) intervalar, explicada pelo fator nominal. Note que, como só existe um fator, o tamanho de efeito global coincide com o tamanho de efeito do fator, \(\eta^2\). A parcela da variância da VD que não é explicada pelo modelo da ANOVA (chamado de resíduo ou termo de erro) corresponde a \(1-\eta^2\).
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
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
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 é grande.
Enquanto a ANCOVA pode ser preterida em relação à ANOVA em delineamentos experimentais nas 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.
O 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”, o 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 após a aplicação do modelo.
Exemplo: Dancey & Reidy, 2019, p. 424
Um experimento foi projetado para descobrir se álcool (nas condições placebo, ingestão em pequena quantidade e ingestão em grande quantidade) afeta o desempenho de dirigir, medido por erros cometidos no simulador veicular. A expectativa é que o número de erros aumente quanto maior for a intoxicação pelo álcool.
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")Alguma estatística descritiva, os testes das suposições para a
ANCOVA, a execução da ANCOVA e testes post hoc estão
implementados em demo_ANCOVA.R.
Estatística descritiva
Implementado em demo_ANCOVA_01descritiva.R:
source("eiras.bartitle.R")
Dados <- readRDS("AlcoholExperienceDrivingErrors.rds")
cat(bartitle("Dados"))
print.data.frame(Dados)
cat(bartitle("Análise descritiva"))
print(xtabs(~Grupo, data=Dados))
ftable(Dados$ID, Dados$Grupo)
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)))
-----
Dados
-----
ID Grupo ErrosDirecao Experiencia
1 1 Placebo 5 12
2 2 Placebo 10 5
3 3 Placebo 7 9
4 4 Placebo 3 24
5 5 Placebo 5 15
6 6 Placebo 7 6
7 7 Placebo 11 3
8 8 Placebo 2 30
9 9 Placebo 3 20
10 10 Placebo 5 10
11 11 Placebo 6 7
12 12 Placebo 6 8
13 13 PoucoAlcool 5 21
14 14 PoucoAlcool 7 16
15 15 PoucoAlcool 9 7
16 16 PoucoAlcool 8 15
17 17 PoucoAlcool 2 30
18 18 PoucoAlcool 5 21
19 19 PoucoAlcool 6 12
20 20 PoucoAlcool 6 13
21 21 PoucoAlcool 4 26
22 22 PoucoAlcool 4 24
23 23 PoucoAlcool 8 9
24 24 PoucoAlcool 10 6
25 25 MuitoAlcool 8 29
26 26 MuitoAlcool 10 8
27 27 MuitoAlcool 8 26
28 28 MuitoAlcool 9 20
29 29 MuitoAlcool 11 18
30 30 MuitoAlcool 15 6
31 31 MuitoAlcool 7 12
32 32 MuitoAlcool 11 7
33 33 MuitoAlcool 8 15
34 34 MuitoAlcool 8 9
35 35 MuitoAlcool 17 3
36 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`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:
source("eiras.bartitle.R")
Dados <- readRDS("AlcoholExperienceDrivingErrors.rds")
alfa <- 0.05
cat(bartitle("Regressão Linear Simples"))
car::scatterplot(ErrosDirecao~Experiencia|Grupo,
data=Dados,
regLine=TRUE, smooth=FALSE, col=c("black","blue","red"))
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_classic()+
ggplot2::labs(title="RLS por Grupo",
subtitle=paste0("BC",round((1-alfa/
length(unique(Dados$Grupo)))*
100,2),"%"),
x="Experiência (ano)",
y="#ErrosDirecao")+
ggplot2::geom_smooth(method=MASS::rlm,
na.rm=TRUE,
color="black",
level = 1-alfa/length(unique(Dados$Grupo)))
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
------------------------
`geom_smooth()` using formula = 'y ~ x'
-------
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.0346Sã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, 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:
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 ' ' 1Os dois testes das suposições apresentados nesta saída são:
- Dissociação entre Grupo (factor) and Experiência (covariável)
A hipótese nula é a ausência de associação entre os níveis do fator e os níveis de experiência dos participantes. Isto é feito por uma ANOVA independente com:
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 todos os níveis do 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.9994Para os três níveis do fator não rejeitou-se 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 ' ' 1A 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 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:
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-10Note que a sintaxe é similar à anterior, exceto pela remoção da
interação Grupo:Experiencia. Computa-se os erros de direção
em função do fator (Grupo) e da covariável
(Experiencia); é uma 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 o número de erros
cometidos na direção.
Na sequência exibimos, também, a regressão 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.
Tamanhos de efeito
Como descrito anteriormente, computam-se o tamanho de efeito global e
os tamanhos de efeito, implementadas em demo_ANCOVA_05tamefeito.R:
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=FALSE,
generalized=FALSE,
ci=1-alfa,
alternative="two.sided",
verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=2)
----------------------------
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 | 95% CI | interpret
---------------------------------------------
Grupo | 0.36 | [0.09, 0.56] | large
Experiencia | 0.40 | [0.14, 0.59] | largeMédias marginais estimadas ajustadas
O próximo item da saída são as médias marginais 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:
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"))
------------------------------------
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 é a partir da regressão linear associada à ANCOVA:
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-10Sendo 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:
grp <- unique(Dados$Grupo)
ngrp <- length(unique(Dados$Grupo))
# interceptos adjustados
adj_intercept <- rep(NA,ngrp)
adj_intercept[1] <- reg$coefficients[1,1]
for (g in 2:ngrp)
{
adj_intercept[g] <- adj_intercept[1]+reg$coefficients[g,1]
}
cat("\n")Intercepts: Placebo = 9.082103for (g in 2:ngrp)
{
cat("\t",as.character(grp[g])," = ",adj_intercept[1]," + ",reg$coefficients[g,1]," = ",adj_intercept[g],"\n", sep="")
} PoucoAlcool = 9.082103 + 1.445328 = 10.52743
MuitoAlcool = 9.082103 + 4.656509 = 13.73861A inclinação comum às três retas é dada por
Common slope = -0.2616459e 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.33333Common Covariate Mean = 14.13889Finalmente, as médias marginais estimadas 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:
adj_mm <- rep(NA,ngrp)
for (g in 1:ngrp)
{
adj_mm[g] <- adj_intercept[g] + adj_slope*adj_cov
}
cat("Adjusted marginal means:\n")Adjusted marginal means:for (g in 1:ngrp)
{
cat("\t",as.character(grp[g])," = ",adj_intercept[g]," + ",adj_slope," * ",adj_cov," = ",adj_mm[g],"\n", sep="")
} 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.03923Uma visão gráfica destas médias marginais ajustadas é:
xreg <- rep(NA,ngrp)
yreg <- rep(NA,ngrp)
xmin <- xmax <- ymin <- ymax <- NA
for (g in 1:ngrp)
{
xmin <- min(Dados$Experiencia[Dados$Grupo==as.character(grp[g])],na.rm=TRUE)
xmax <- max(Dados$Experiencia[Dados$Grupo==as.character(grp[g])],na.rm=TRUE)
x <- seq(xmin,xmax,length.out=100)
if(is.na(xmin)) {xmin<-min(x)} else { if(xmin>min(x)){xmin<-min(x)} }
if(is.na(xmax)) {xmax<-max(x)} else { if(xmax<max(x)){xmax<-max(x)} }
y <- adj_intercept[g] + adj_slope*x
if(is.na(ymin)) {ymin<-min(y)} else { if(ymin>min(y)){ymin<-min(y)} }
if(is.na(ymax)) {ymax<-max(y)} else { if(ymax<max(y)){ymax<-max(y)} }
xreg[g] <- list(x)
yreg[g] <- list(y)
}
plot(NA,
xlim=c(xmin,xmax), ylim=c(ymin,ymax),
main="Adjusted Marginal Means",
xlab="Experiencia",
ylab="ErrosDirecao")
abline(v=adj_cov,lty=2)
text(adj_cov,ymin,round(adj_cov,3),pos=2,cex=0.6)
for (g in 1:ngrp)
{
lines(unlist(xreg[g]),unlist(yreg[g]))
x <- max(unlist(xreg[g]))
y <- adj_intercept[g] + adj_slope*x
text(x,y,as.character(grp[g]),pos=2,cex=0.8)
y <- adj_intercept[g] + adj_slope*adj_cov
lines(c(xmin,adj_cov), c(y, y), lty=2)
text(xmin,y,round(adj_mm[g],3),pos=3,cex=0.6)
}
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 estimadas ajustadas
está em demo_AdjustedMarginalMeans.R.
|
Testes post hoc
Resta apenas localizar as diferenças dos números de erros entre os
diferentes níveis do fator. 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:
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"))
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"))
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"))
-------------------------------------------------------------------
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:
“An advantage of the use of ANCOVA is that it adjusts for baseline differences between the treatment groups. ANCOVA also has more statistical power than the t-test, so sample size requirements are lower.”
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 55lista <- list(Dados.wide, Dados.long)
names(lista) <- c("Dados.wide", "Dados.long")
saveRDS(lista, "PrePos_2CondIndep.rds")“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.” (Sakae et al., 2016)
Esta ANCOVA é 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:
source("eiras.bartitle.R")
Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8")
Sys.setlocale("LC_ALL", "pt_BR.UTF-8")
alfa <- 0.05
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)))
cat(bartitle("Simple linear regressions"))
# formato wide eh melhor
car::scatterplot(Pos~Pre|Grupo, 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
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_classic()+
ggplot2::labs(title="RLS por Grupo",
subtitle=paste0("BC",
round((1-alfa/
length(unique(Dados.wide$Grupo))),4)*
100,"%"),
x="Escore Pre",
y="Escore Pos")+
ggplot2::geom_smooth(method=MASS::rlm,
na.rm=TRUE,
color="black",
level=1-alfa/length(unique(Dados.wide$Grupo)))
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,
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"))
cat(bartitle("Post hoc tests: Grupo (factor) adjusted by VD Pre (covariate)"))
print(summary(EMM$contrasts, infer=TRUE))
print(plot(EMM$contrasts,
colors="black"))
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) | 95% CI
-----------------------------------------
Grupo | 0.60 | [0.05, 0.83]
Pre | 0.87 | [0.54, 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:
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.
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 ANOVA da ANCOVA, avalia-se o valor p associado ao fator, que neste caso é 0.0136. Portanto, rejeita-se a hipótese nula para \(\alpha=0.05\) e assume-se que as drogas A e B têm efeitos diversos.
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 da ANCOVA, interessa observar as duas últimas linhas:
Multiple R-squared: 0.9071, Adjusted R-squared: 0.8806 F-statistic: 34.18 on 2 and 7 DF, p-value: 0.0002442
Aqui testa-se se existe modelo ANCOVA. O valor p=0.0002442 rejeita a hipótese nula de ausência de modelo ANCOVA. Portanto, há ANCOVA. Também é onde encontramos o tamanho de efeito global (\(\eta^2 = R^2 = 0.91\)).
- 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 Rscript para o caso do código ser adaptado para situações com mais níveis do fator).
Planejamento
Segundo BORM et al. (2007):
“The use of ANCOVA may considerably reduce the number of patients required for a trial. […] ANCOVA on \((1 - \rho^2) \cdot n\) observations has the same power as a t-test on \(n\) observations.”
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 <- .90
n.porcondicao <- 2*((qnorm(1-alfa/2)+qnorm(poder))*sqrt(1-rho^2)/delta)^2
n.total <- 2*n.porcondicao
cat("\nn total = ",n.total," com ",n.porcondicao," por nivel do fator\n")
n total = 107.596 com 53.79801 por nivel do fatorNeste Rscript encontramos o tamanho da amostra para \(\alpha=0.05\) e \(\text{poder}=0.90\). Temos:
deltaque é o \(d\) de Cohen populacional (tamanho de efeito populacional)rhoque a correlação populacional de Pearson entre a VD e a covariável.
Podemos variar o valor do poder para termos ideia do poder do estudo de nosso exemplo, com 10 idosos no total. Encontramos:
delta <- .5
rho <- .6
alfa <- .05
poder <- .17
n.porcondicao <- 2*((qnorm(1-alfa/2)+qnorm(poder))*sqrt(1-rho^2)/delta)^2
n.total <- 2*n.porcondicao
cat("\nn total = ",n.total," com ",n.porcondicao," por nivel do fator\n")
n total = 10.3591 com 5.179551 por nivel do fatorPortanto, o poder a priori de um estudo deste tamanho era bem baixo, da ordem de 17%.
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 90% de haver benefício com a inclusão da covariável necessita-se de \(n \ge 22\).
Para quem quiser modelo mais complexo, Crawley (2012) em “Chapter12 - Analysis of covariance” apresenta um delineamento com dois fatores e uma covariável.
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.
- Dancey CP & Reidy J (2019) Estatística sem Matemática para Psicologia. 7ª ed. Porto Alegre: Penso.
- Michael J. Crawley (2012) The R Book. 2nd Edition.
- Sakae TM et al. (2016) Efeitos da anestesia geral da cognição do idoso. Arquivos Catarinenses de Medicina 45(3): 107-16.
- Xiaofeng Liu (2011) The Effect of a Covariate on Standard Error and Confidence Interval Width. Communications in Statistics - Theory and Methods 40:3, 449-456, DOI: 10.1080/03610920903391337