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))RPubsConforme 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
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 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.
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.
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:
Na sequência há vários gráficos para que são adequados para o propósito da sua pesquisa.
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.
Além das suposições da regressão linear simples e da ANOVA, a ANCOVA supõe:
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:
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:
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.
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.
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.
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.
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
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:
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()
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.
|
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.
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.”
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:
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 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
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ávelSegundo 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.
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
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
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
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.
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.
# 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
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
R^2 = eta^2 omnibus = 0.55
[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
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
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
R^2 = eta^2 omnibus = 0.54
[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
Superpower. https://aaroncaldwell.us/SuperpowerBook/