Análise de Variância Multivariada (MANOVA) em R
MSP4112: Estatística Aplicada à Medicina
28 novembro 2025 00:05h
Bastão de Asclépio & Distribuição Normal
suppressMessages(library(car, warn.conflicts=FALSE))
suppressMessages(library(DescTools, 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(heplots, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(MASS, warn.conflicts=FALSE))
suppressMessages(library(multcomp, warn.conflicts=FALSE))
suppressMessages(library(mvdalab, warn.conflicts=FALSE))
suppressMessages(library(MVN, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
suppressMessages(library(rrcov, warn.conflicts=FALSE))
suppressMessages(library(SHT, warn.conflicts=FALSE))
suppressMessages(library(tidyr, warn.conflicts=FALSE))source("eiras.bartitle.R")
source("eiras.bivCI.R")
source("eiras.correg.R")
source("eiras.friendlycolor.R")
source("eiras.numeric.summary.R")
source("eiras.showdataframe.R")Material
- HTML de Rmarkdown em
RPubs
Incluir
bruceRbruceR::CorrbruceR::DescribebruceR::RECODEbruceR::model_summarybruceR::regressbruceR::EMMEANSbruceR::GLM_summarybruceR::HLM_summarybruceR::TTESTbruceR::MANOVA: https://psychbruce.github.io/bruceR/reference/MANOVA.htmlbruceR::lavaan_summarybruceR::PCAbruceR::EFAbruceR::cfa
Objetivos
- Conceituar MANOVA, distinguindo-a de ANOVA e ANCOVA.
- Conceituar e testar as suposições da MANOVA.
- Formular e implementar os estes multivariados:
- \(T^2\) de Hotelling para uma condição
- \(T^2\) de Hotelling independente para duas condições: homocedástico e heterocedástico
- MANOVA unifatorial independente
- MANCOVA
- Determinar as significâncias estatística e prática do efeito do fator fixo.
Disciplina de pós-graduação do PBE-UEM
- DES4033: Análise Multivariada em Biologia
O que é MANOVA
Técnicas estatísticas multivariadas são extensões das técnicas univariadas e bivariadas para situações nas quais temos mais do que uma variável dependente ou medida:
- \(T^2\) de Hotelling: testar centróide(s) de medidas intervalares multivariadas observadas em uma ou duas condições independentes.
- MANOVA: Análise de Variância Multivariada Independente
- MANCOVA: Análise de Covariância Multivariada
- rmMANOVA: Análise de Variância Multivariada para Medidas Repetidas
- MANOVA mista (um fator entre participantes e um fator intraparticipantes)
O teste \(T^2\) de Hotelling para uma condição é o primeiro e mais básico dos testes multivariados.
O teste \(T^2\) de Hotelling para duas condições independentes testa a igualdade de dois centróides multivariados populacionais.
Os testes \(T^2\) de Hotelling foram concebidos por Hotelling (1931, 1951).
Em MANOVA, concebida por Wilks (1932), o efeito de um fator com três ou mais condições pode ser testado e, em MANCOVA, incluiremos covariáveis.
Com estes testes, ampliamos o uso do modelo linear geral (GLM).
Por que utilizar MANOVA?
Um construto pode ser medido por meio de alguns itens intervalares, por exemplo, compleição que poderia ser decomposta com medidas de massa corporal total (MCT) e estatura. Em geral, as medidas são correlacionadas. Além disto, para indicar o uso de MANOVA independente, deve existir um fator como, por exemplo, sexo (masculino e feminino). Queremos determinar se existe diferença de compleição masculina e feminina.
Diferentemente dos testes feitos anteriomente, a variável de desfecho (VD), aqui, não é composta de uma medida única: neste exemplo a VD é a combinação das medidas de MCT e estatura. O fator nominal que explicaria esta diferença é sexo. Outros fatores (nominais) ou covariáveis (intervalares) também podem ser adicionados a estes modelos como variáveis explicativas.
Da mesma forma que ocorre com testes t em relação à ANOVA, tentar substituir uma MANOVA pela realização de um conjunto de ANOVAs aumenta a probabilidade do erro do tipo I (falso-positivo). MANOVA permite tratar o conjunto em uma única análise, aproveitando as correlações entre as medidas.
Os testes básicos e o GLM
Os testes básicos podem ser vistos como manifestações do GLM com graus crescentes de complexidade:
- Teste t
1 medida intervalar & 1 fator entre ou intraparticipantes com duas condições - ANOVA unifatorial
1 medida intervalar & 1 fator com três ou mais condições - Regressão univariada simples
1 medida intervalar & 1 VE intervalar - ANCOVA unifatorial simples
1 medida intervalar & 1 fator e 1 VE intervalar (covariável) - Regressão univariada múltipla
1 medida intervalar & 2 ou mais VEs intervalares - Regressão multivariada múltipla
2 ou mais medidas intervalares & 2 ou mais VEs intervalares - \(T^2\) de Hotelling
2 ou mais medidas intervalares & 1 fator entre ou intraparticipantes (duplamente multivariada) com duas condições - MANOVA unifatorial
2 ou mais medidas intervalares & 1 fator entre ou intraparticipantes (duplamente multivariada) com três ou mais condições - MANCOVA unifatorial simples
2 ou mais medidas intervalares & 1 fator entre ou intraparticipantes (duplamente multivariada) e 1 VE intervalar (covariável) - …
Modelos podem evoluir daqui em diante. Se as medidas fossem estatura e MCT, e as variáveis explicativas fossem sexo e raça (ambas nominais), teremos uma MANOVA independente bifatorial. Se adicionarmos idade e tempo de prática esportiva (intervalares), é uma MANCOVA independente bifatorial múltipla. Caso as medidas sejam repetidas, uma MANCOVA bifatorial múltipla relacionada.
A lógica da MANOVA
Na ANOVA, particionamos a variância da medida que pode ser atribuída aos efeitos principais e de interação das VIs. Na MANOVA, a situação é mais complexa porque há mais de uma medida envolvida e, portanto, há uma matriz de covariância compondo a variabilidade total das medidas. Em nosso exemplo, a variabilidade total depende da variância da estatura, da variância da massa corporal e da covariância entre elas. Além disto, se há um fator com dois níveis (e.g., sexos masculino e feminino), então há dois destes conjuntos de variância a serem levados em conta.
No processo de cálculo da MANOVA, as diversas medidas são internamente convertidas em um escore intervalar similar ao que acontece com um modelo de análise fatorial unidimensional e utiliza este escore, chamado de combinação linear ótima das medidas, criando uma única VD para a análise. A análise informa se existe(m) efeito(s) das VIs na combinação linear das VDs.
As suposições subjacentes ao uso da MANOVA
As duas suposições são de normalidade multivariada e homocedasticidade multivariada das medidas nos níveis do fator.
Teste \(T^2\) de Hotelling para uma condição
O exemplo de estatura e MCT tem um fator entre participantes (sexo) com dois níveis e duas medidas (estatura e MCT).
Vamos supor, no entanto, que estejamos interessados apenas nos participantes do sexo masculino. Adiante verificaremos se estes estudantes são brasileiros típicos de acordo com a medida de Moreira et al. (2013) que indica, em média, estatura de 172 cm e MCT de 76.68 kg.
O arquivo de dados usado neste caso é Adm2008.rds.
Dados <- data.frame(readxl::read_excel("Adm2008.xlsx"))
Dados$Nome <- factor(Dados$Nome)
Dados$Sexo <- factor(Dados$Sexo)
Dados$Estatura <- Dados$Estatura*100
saveRDS(Dados, "Adm2008.rds")Iniciando pela estatística descritiva, sugere-se o que está
implementado em demo_Adm2008Masc_descritiva.R:
source("eiras.friendlycolor.R")
source("eiras.correg.R")
source("eiras.bivCI.R")
source("eiras.showdataframe.R")
source("eiras.numeric.summary.R")
source("eiras.bartitle.R")
H0.M <- c(172,76.68)
alfa <- 0.05
B <- 1e3
colM <- friendlycolor(7) # azul
pchM <- 24
alpha <- 0.05
cat(bartitle("Data"))
Dados <- readRDS("Adm2008.rds")
Masc <- subset(Dados, Sexo=="Masculino", select=c(Estatura, MCT))
showdataframe(Masc, head=5, tail=3)
print(summary(Masc))
cat(bartitle("Summary"))
print(psych::describe(Masc[,c("Estatura","MCT")], data=Masc))
cat(bartitle("Correlation matrix",3))
print(cor(Masc), digits=2)
cat(bartitle("Covariance matrix",3))
print(cov(Masc), digits=2)
correg(Masc$Estatura,
Masc$MCT,
method="raw", conf.band = FALSE, lowess=FALSE,
alpha=0.05, B=0,
xlab="Estatura (cm)", ylab="MCT (kg)",
main="Estudante de Administração Masculino",
col=colM, bg="transparent", pch=pchM,
show.equation = FALSE, suppress.text = TRUE)
# elipses de confianca do centroide populacional
lines(bivCI(s = var(Masc),
xbar = colMeans(Masc),
n = nrow(Masc),
alpha = alpha, m = 1e4),
type = "l", col = colM, lwd = 2)
cx <- mean(Dados$Estatura[Dados$Sexo=="Masculino"], na.rm=TRUE)
cy <- mean(Dados$MCT[Dados$Sexo=="Masculino"], na.rm=TRUE)
points(cx,cy,pch=21,col=colM,bg=colM)
points(H0.M[1],H0.M[2],pch=10)
text(H0.M[1],H0.M[2],expression(H[0]),pos=2)
text(cx,cy, pos=3, "M", col=colM)
----
Data
----
Estatura MCT
188 82
191 83
172 58
172 82
182 77
... ...
178 92
168 71
176 78
Estatura MCT
Min. :156.0 Min. : 48.00
1st Qu.:172.0 1st Qu.: 66.00
Median :176.0 Median : 73.00
Mean :176.4 Mean : 74.47
3rd Qu.:182.0 3rd Qu.: 82.00
Max. :193.0 Max. :105.00
-------
Summary
-------
vars n mean sd median trimmed mad min max range skew kurtosis
Estatura 1 51 176.35 8.08 176 176.49 8.90 156 193 37 -0.17 -0.36
MCT 2 51 74.47 12.10 73 73.93 11.86 48 105 57 0.39 -0.30
se
Estatura 1.13
MCT 1.69
------------------
Correlation matrix
------------------
Estatura MCT
Estatura 1.00 0.65
MCT 0.65 1.00
-----------------
Covariance matrix
-----------------
Estatura MCT
Estatura 65 63
MCT 63 146
A novidade é a região elíptica de confiança de 95% do centróide populacional. Como a hipótese nula está fora desta região elíptica, rejeita-se \(H_0\) para qualquer nível de significância convencional. O estudante do sexo masculino do curso noturno de Aministração da USP de 2008 têm compleição média diferente da população brasileira geral de 2013: tem estatura média maior e mesma massa média.
Suposições
MANOVA supõe normalidade multivariada das medidas (neste caso bivariada porque só temos duas medidas: estatura e MCT, na condição sexo masculino).
Os testes estão implementados em demo_Adm2008Masc_normalidade.R:
source("eiras.bartitle.R")
Dados <- readRDS("Adm2008.rds")
Masc <- subset(Dados, Sexo=="Masculino", select=c(Estatura, MCT))
result <- MVN::mvn(data = Masc,
# R=5e4,
# mvnTest = "energy",
mvn_test = "hz",
univariate_test = "SW",)
cat(bartitle("Binormality test"))
print(result$multivariate_normality)
----------------
Binormality test
----------------
Test Statistic p.value Method MVN
1 Henze-Zirkler 1.126 0.01 asymptotic ✗ Not normalEstes testes foram realizados com MVN::mvn e mostram os testes de uninormalidade e
binormalidade para o nível de significância de 5%.
Univariadamente, as distribuições podem ser consideradas normais. Houve rejeição da binormalidade para o grupo Masculino (mas não seria rejeitada com \(\alpha=1\%\)). No entanto, a amostra tem 51 homens e podemos invocar o teorema central do limite.
Teste da hipótese nula multivariada
Admitindo que as medidas relatadas em Moreira et al. (2013) são estimativas da população brasileira, as hipóteses nula e alternativa são:
\[
\begin{cases}
H_0:
\begin{bmatrix}
\mu_{\text{Est}}\\
\mu_{\text{MCT}}
\end{bmatrix} =
\begin{bmatrix}
172\\
76.68
\end{bmatrix}\\
H_1:
\begin{bmatrix}
\mu_{\text{Est}}\\
\mu_{\text{MCT}}
\end{bmatrix} \ne
\begin{bmatrix}
172\\
76.68
\end{bmatrix}
\end{cases}\\
\alpha=5\%
\] Implementado em demo_T2Hotelling_1cond.R:
source("eiras.bartitle.R")
source("eiras.bivCI.R")
H0.M <- c(172, 76.68)
alfa <- 0.05
B <- 1e4
Dados <- readRDS("Adm2008.rds")
Masc <- subset(Dados, Sexo=="Masculino", select=c(Estatura, MCT))
rm(Dados)
p <- dim(Masc)[2]
n.M <- dim(Masc)[1]
cov.M <- cov(Masc)
mean.M <- colMeans(Masc)
Sinv.M <- solve(cov.M)
F.M <- (n.M*(n.M-p)/((n.M-1)*p))*
(mean.M-H0.M)%*%Sinv.M%*%(mean.M-H0.M)
p.M <- 1-pf(q=F.M,df1=p,df2=n.M-p)
cat("\nMANOVA para uma condição:\n", sep="")
cat("\tF(",p,",",n.M-p,") = ",F.M,", p = ",p.M,"\n", sep="")
T2One <- MVTests::OneSampleHT2(Masc,
mu=H0.M,
alpha=alfa)
print(summary(T2One), digits=2)
T2crit <- ((n.M-1)*p/(n.M-p))*qf(1-alfa, p, n.M-p)
c <- sqrt(T2crit/n.M)
car::ellipse(center=mean.M,
shape=cov.M,
radius=c,
fill=TRUE,
fill.alpha=0.1,
grid=FALSE,
col="black",
add=FALSE,
xlim=c(170, 180),
ylim=c(70, 80),
xlab=expression(mu[Est]),
ylab=expression(mu[MCT]),
main="Região elíptica de confiança de 95%\nEstudante Masculino")
points(mean.M[1], mean.M[2], col = "black", pch = 19)
points(H0.M[1], H0.M[2], pch=3, col="black")
text(H0.M[1], H0.M[2], pos=1, expression(H[0]), col="black")
text(mean.M[1],mean.M[2], pos=2, expression(bar(x)[M]), col="black")
mvdalab::MVcis(Masc,
level=1-alfa,
Vars2Plot=c(1, 2),
include.zero=FALSE)
MANOVA para uma condição:
F(2,49) = 19.34467, p = 6.421062e-07
One Sample Hotelling T Square Test
Hotelling T Sqaure Statistic = 39.47892
F value = 19.345 , df1 = 2 , df2 = 49 , p-value: 6.42e-07
Descriptive Statistics
Estatura MCT
N 51.000000 51.00000
Means 176.352941 74.47059
Sd 8.076691 12.09852
Detection important variable(s)
Lower Upper Mu0 Important Variables?
Estatura 173.46882 179.23706 172.00 *TRUE*
MCT 70.15031 78.79087 76.68 FALSE
NULL
[,1] [,2]
Estatura 173.46882 179.23706
MCT 70.15031 78.79087Rejeita-se \(H_0\) para qualquer nível de significância convencional. O estudante do sexo masculino do curso noturno de Aministração da USP de 2008 têm compleição média diferente da população brasileira geral de 2013: tem estatura média maior e mesma massa média.
Teste \(T^2\) de Hotelling para duas condições independentes
Teste \(T^2\) de Hotelling
independente tem um fator dicotômico entre participantes (sexo) e duas
medidas (estatura e MCT). A estatística descritiva está implementada em
demo_Adm2008_descritiva.R:
suppressMessages(suppressWarnings(
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
))
source("eiras.friendlycolor.R")
source("eiras.correg.R")
source("eiras.bivCI.R")
source("eiras.showdataframe.R")
source("eiras.numeric.summary.R")
source("eiras.bartitle.R")
colM <- friendlycolor(7) # azul
colF <- friendlycolor(30) # vermelho
pchM <- 24
pchF <- 21
alpha <- 0.05
cat(bartitle("Data"))
Dados <- readRDS("Adm2008.rds")
Dados <- Dados[, c("Sexo","Estatura","MCT")]
showdataframe(Dados, head=5, tail=3)
cat(bartitle("Summary"))
cat(bartitle("Masculino",2))
cat(bartitle("Estatura",3))
print(numeric.summary(Dados$Estatura[Dados$Sexo=="Masculino"]))
cat(bartitle("MCT",3))
print(numeric.summary(Dados$MCT[Dados$Sexo=="Masculino"]))
dt_tmp <- data.frame(Dados$Estatura[Dados$Sexo=="Masculino"],
Dados$MCT[Dados$Sexo=="Masculino"])
names(dt_tmp) <- c("Estatura","MCT")
cat(bartitle("Correlation matrix",3))
print(cor(dt_tmp), digits=2)
cat(bartitle("Covariance matrix",3))
print(cov(dt_tmp), digits=2)
cat(bartitle("Feminino",2))
cat(bartitle("Estatura",3))
print(numeric.summary(Dados$Estatura[Dados$Sexo=="Feminino"]))
cat(bartitle("MCT",3))
print(numeric.summary(Dados$MCT[Dados$Sexo=="Feminino"]))
dt_tmp <- data.frame(Dados$Estatura[Dados$Sexo=="Feminino"],
Dados$MCT[Dados$Sexo=="Feminino"])
names(dt_tmp) <- c("Estatura","MCT")
cat(bartitle("Correlation matrix",3))
print(cor(dt_tmp), digits=2)
cat(bartitle("Covariance matrix",3))
print(cov(dt_tmp), digits=2)
print(psych::describeBy(Dados[,2:3], g = Dados$Sexo,
mat=1, digits = 2))
Dados.long <- reshape2::melt(Dados,
id=c("Sexo"),
measured=c("Estatura","MCT"))
names(Dados.long) <- c('Group', 'Measure', 'Value')
Dados.long.Boxplot <- ggplot2::ggplot(Dados.long,
ggplot2::aes(Group,
Value,
color = Measure)) +
ggplot2::geom_boxplot() +
ggplot2::labs(x="Sexo",
y="Valor",
color="Medida") +
ggplot2::theme_bw()+
ggplot2::theme(panel.grid = ggplot2::element_blank())
plot(Dados.long.Boxplot)
rm(Dados.long)
suppressMessages(suppressWarnings(
invisible(print(GGally::ggpairs(data = Dados,
mapping = ggplot2::aes(color = Sexo,
alpha = .5),
progress = FALSE) +
ggplot2::theme_bw()+
ggplot2::theme(panel.grid = ggplot2::element_blank())))))
correg(Dados$Estatura[Dados$Sexo=="Masculino"],
Dados$MCT[Dados$Sexo=="Masculino"],
method="raw", conf.band = FALSE, lowess=FALSE,
alpha=0.05, B=0,
xlab="Estatura (cm)", ylab="MCT (kg)",
main="Elipses de confiança de 95% Bonferroni\nEstudante de Administração",
col=colM, bg="transparent", pch=pchM,
show.equation = FALSE, suppress.text = TRUE)
correg(Dados$Estatura[Dados$Sexo=="Feminino"],
Dados$MCT[Dados$Sexo=="Feminino"],
method="raw", conf.band = FALSE, lowess=FALSE,
alpha=0.05, B=0,
xlab="Estatura (cm)", ylab="MCT (kg)",
col=colF, bg="transparent", pch=pchF,
show.equation = FALSE, suppress.text = TRUE, add=TRUE)
# elipses de confianca do centroide populacional com correçao de Bonferroni
df_tmp <- data.frame(Dados$Estatura[Dados$Sexo=="Masculino"],
Dados$MCT[Dados$Sexo=="Masculino"])
names(df_tmp) <- c("Estatura","MCT")
df_tmp <- df_tmp[!is.na(df_tmp$Estatura),]
df_tmp <- df_tmp[!is.na(df_tmp$MCT),]
lines(bivCI(s = var(df_tmp),
xbar = colMeans(df_tmp),
n = nrow(df_tmp),
alpha = alpha/2, m = 10000),
type = "l", col = colM, lwd = 2)
cx <- mean(Dados$Estatura[Dados$Sexo=="Masculino"], na.rm=TRUE)
cy <- mean(Dados$MCT[Dados$Sexo=="Masculino"], na.rm=TRUE)
points(cx,cy,pch=21,col=colM,bg=colM)
text(cx,cy, pos=2, "M", col=colM)
df_tmp <- data.frame(Dados$Estatura[Dados$Sexo=="Feminino"],
Dados$MCT[Dados$Sexo=="Feminino"])
names(df_tmp) <- c("Estatura","MCT")
df_tmp <- df_tmp[!is.na(df_tmp$Estatura),]
df_tmp <- df_tmp[!is.na(df_tmp$MCT),]
lines(bivCI(s = var(df_tmp),
xbar = colMeans(df_tmp),
n = nrow(df_tmp),
alpha = alpha/2, m = 10000),
type = "l", col = colF, lwd = 2)
cx <- mean(Dados$Estatura[Dados$Sexo=="Feminino"], na.rm=TRUE)
cy <- mean(Dados$MCT[Dados$Sexo=="Feminino"], na.rm=TRUE)
points(cx,cy,pch=21,col=colF,bg=colF)
text(cx,cy, pos=2, "F", col=colF)
----
Data
----
Sexo Estatura MCT
Feminino 161 53
Feminino 156 50
Feminino 172 60
Feminino 157 44
Feminino 168 57
... ... ...
Masculino 178 92
Masculino 168 71
Masculino 176 78
-------
Summary
-------
---------
Masculino
---------
--------
Estatura
--------
Min. 1st Qu. Median 3rd Qu. Max. Mean St.Dev. n NA
156 172 176 182 193 176.3529 8.076691 51 0
---
MCT
---
Min. 1st Qu. Median 3rd Qu. Max. Mean St.Dev. n NA
48 66 73 82 105 74.47059 12.09852 51 0
------------------
Correlation matrix
------------------
Estatura MCT
Estatura 1.00 0.65
MCT 0.65 1.00
-----------------
Covariance matrix
-----------------
Estatura MCT
Estatura 65 63
MCT 63 146
--------
Feminino
--------
--------
Estatura
--------
Min. 1st Qu. Median 3rd Qu. Max. Mean St.Dev. n NA
150 160 164 169.75 178 164.1316 6.798931 38 0
---
MCT
---
Min. 1st Qu. Median 3rd Qu. Max. Mean St.Dev. n NA
43 50.25 53.5 59.75 70 54.63158 6.188136 38 0
------------------
Correlation matrix
------------------
Estatura MCT
Estatura 1.00 0.64
MCT 0.64 1.00
-----------------
Covariance matrix
-----------------
Estatura MCT
Estatura 46 27
MCT 27 38
item group1 vars n mean sd median trimmed mad min max
Estatura1 1 Feminino 1 38 164.13 6.80 164.0 164.22 8.15 150 178
Estatura2 2 Masculino 1 51 176.35 8.08 176.0 176.49 8.90 156 193
MCT1 3 Feminino 2 38 54.63 6.19 53.5 54.47 5.93 43 70
MCT2 4 Masculino 2 51 74.47 12.10 73.0 73.93 11.86 48 105
range skew kurtosis se
Estatura1 28 -0.18 -0.81 1.10
Estatura2 37 -0.17 -0.36 1.13
MCT1 27 0.34 -0.48 1.00
MCT2 57 0.39 -0.30 1.69
Este gráfico mostra o plano das medidas brutas. As observações para os dois grupos foram representados separadamente. Assinalamos o centróide (um par ordenado) para cada grupo.
As regiões elípticas de confiança de 95% (com correção de Bonferroni) dos centróides populacionais podem ser vistos como uma extensão bivariada do intervalo de confiança de 95% das médias populacionais de estatura e MCT combinadas. Podemos perceber que as elipses dos dois grupos estão bem distantes: quando não há sobreposição, podemos esperar que a hipótese nula de igualdade dos centróides populacionais será rejeitada.
Suposição de normalidade multivariada
Implementado em demo_Adm2008_normalidade.R:
source("eiras.bartitle.R")
Dados <- readRDS("Adm2008.rds")
Dados <- Dados[, c("Sexo","Estatura","MCT")]
result <- MVN::mvn(data = Dados,
subset="Sexo",
# R=5e4,
# mvnTest = "energy",
univariate_test = "SW")
cat(bartitle("Binormality test"))
print(result$multivariate_normality)
----------------
Binormality test
----------------
Group Test Statistic p.value MVN
1 Feminino Henze-Zirkler 0.287 0.877 ✓ Normal
2 Masculino Henze-Zirkler 1.126 0.010 ✗ Not normalEstes testes foram feitos com MVN::mvn e mostram os testes de uninormalidade e
binormalidade em cada nível do fator. Para o nível de significância de
5%. Isoladamente, as distribuições são consideradas normais. Houve
rejeição apenas da binormalidade para o grupo Masculino (mas não seria
rejeitada com \(\alpha=1\%\)). No
entanto, a amostra tem 51 homens e 38 mulheres, e podemos apelar para o
teorema central do limite (TCL).
“O teste da suposição de multinormalidade, em geral, é necessário para amostra pequena. Se a amostra é suficientemente grande (maior que \(20 - p\) unidade experimentais, sendo que \(p\) é o número de medidas), o TCL entra em ação e a suposição de multinormalidade não é mais necessária.”
Johnson & Wichern, 2007, p. 294
Suposição de homocedasticidade multivariada
Avaliamos a homogeneidade das matrizes de covariância. Implementado
em demo_Adm2008_homocedasticidade.R:
source("eiras.bartitle.R")
Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8")
Sys.setlocale("LC_ALL", "pt_BR.UTF-8")
options(warn=-1) # disable warnings
Dados <- readRDS("Adm2008.rds")
print(res <- heplots::boxM(Dados[,c("Estatura","MCT")], Dados$Sexo))
plot(res, main="Estudante de Administração\nEstatura e MCT")
options(warn=0) # enable warnings
Box's M-test for Homogeneity of Covariance Matrices
data: Dados[, c("Estatura", "MCT")]
Chi-Sq (approx.) = 19.918, df = 3, p-value = 0.0001765
Há heterocedasticidade multivariada entre os níveis do fator, porém…
“… com amostras razoavelmente grandes, os testes de médias ou de efeito de tratamentos por MANOVA são bastante robustos à não normalidade. Assim, o teste M de Box pode rejeitar a hipótese nula em alguns casos sem normalidade sem prejuízo para MANOVA. Além disso, com tamanhos de amostra iguais, alguma diferença nas matrizes de covariância têm pouco efeito na MANOVA. Para resumir, podemos decidir continuar com os testes de MANOVA habituais, mesmo que o teste M de Box leve à rejeição da hipótese nula.”
Johnson & Wichern, 2007, p. 291-6
Teste da hipótese nula multivariada
A hipótese nula a seguir é a igualdade de dois centróides populacionais.
Formalmente é:
\[ \begin{cases} H_0: \left[ \mu^{F}_{Est}~ \;\mu^{F}_{MCT} \right] = \left[ \mu^{M}_{Est}~\;\mu^{M}_{MCT} \right]\\H_1: \left[ \mu^{F}_{Est}~\;\mu^{F}_{MCT} \right] \ne \left[ \mu^{M}_{Est}~\;\mu^{M}_{MCT} \right] \end{cases}\\ \alpha=5\% \]
A matriz que conjuga Estatura e MCT para cada um dos sexos é a representação da compleição. Portanto, as hipóteses nula e alternativa podem ser reescritas como:
\[ \begin{cases} H_0: \text{Compleição média feminina} = \text{Compleição média masculina}\\H_1: \text{Compleição média feminina} \ne \text{Compleição média masculina} \end{cases}\\ \alpha=5\% \]
Implementado em demo_T2Hotelling_2cond.R:
source("eiras.bartitle.R")
cat(bartitle("Data"))
Dados <- readRDS("Adm2008.rds")
Dados <- Dados[, c("Sexo","Estatura","MCT")]
print(summary(Dados))
print(FSA::Summarize(Estatura~Sexo,
digits=2,
data=Dados))
print(FSA::Summarize(MCT~Sexo,
digits=2,
data=Dados))
cat(bartitle("T^2 de Hotelling (homocedástico)"))
# Solução 1: DescTools::HotellingsT2Test
# T.2 é F
print(DescTools::HotellingsT2Test(cbind(Estatura,MCT)~Sexo,
data = Dados))
Dados.Mm <- as.matrix(subset(Dados, Sexo=="Masculino", select=c(Estatura, MCT)))
Dados.Fm <- as.matrix(subset(Dados, Sexo=="Feminino", select=c(Estatura, MCT)))
print(SHT::mean2.1931Hotelling(Dados.Mm, Dados.Fm))
cat(bartitle("T^2 de Hotelling (homocedástico)"))
# Solução 2: manova
fit <- lm(cbind(Estatura,MCT)~Sexo,
data = Dados)
print(anv <- car::Anova(fit,
test="Pillai"),
digits=5)
alpha <- 0.05
eta2 <- effectsize::eta_squared(anv,
partial=TRUE,
generalized=FALSE,
ci=1-alpha,
alternative="two.sided",
verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=2)
----
Data
----
Sexo Estatura MCT
Feminino :38 Min. :150.0 Min. : 43
Masculino:51 1st Qu.:164.0 1st Qu.: 55
Median :171.0 Median : 65
Mean :171.1 Mean : 66
3rd Qu.:178.0 3rd Qu.: 76
Max. :193.0 Max. :105 Registered S3 methods overwritten by 'FSA':
method from
confint.boot car
hist.boot car Sexo n mean sd min Q1 median Q3 max
1 Feminino 38 164.13 6.80 150 160 164 169.75 178
2 Masculino 51 176.35 8.08 156 172 176 182.00 193
Sexo n mean sd min Q1 median Q3 max
1 Feminino 38 54.63 6.19 43 50.25 53.5 59.75 70
2 Masculino 51 74.47 12.10 48 66.00 73.0 82.00 105
--------------------------------
T^2 de Hotelling (homocedástico)
--------------------------------
Hotelling's two sample T2-test
data: cbind(Estatura, MCT) by Sexo
T.2 = 44.633, df1 = 2, df2 = 86, p-value = 5.063e-14
alternative hypothesis: true location difference is not equal to c(0,0)
Hotelling's T-squared Test for Independent Samples with Equal
Covariance Assumption.
data: Dados.Mm and Dados.Fm
T2 = 90.304, p-value = 5.064e-14
alternative hypothesis: true means are different.
--------------------------------
T^2 de Hotelling (homocedástico)
--------------------------------
Type II MANOVA Tests: Pillai test statistic
Df test stat approx F num Df den Df Pr(>F)
Sexo 1 0.50932 44.633 2 86 5.064e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Effect Size for ANOVA (Type II)
Parameter | Eta2 (partial) | 95% CI | interpret
-----------------------------------------------------
Sexo | 0.51 | [0.36, 0.62] | largeUtilizamos, como em capítulos anteriores, a função lm
(ou equivalentemente manova) nativa do R, chamada com a
sintaxe:
A função
A hipótese nula de igualdade das compleições de homem e mulher foi
rejeitada pelos testes \(T^2\) de
Hotelling (DescTools::HotellingsT2Test) e MANOVA
(lm). Verifica-se que o teste \(T^2\) de Hotelling é um caso particular de
MANOVA.
O código, por razões históricas, aplica primeiro o \(T^2\) de Hotelling ainda usado na
literatura e implementado com a função
DescTools::HotellingsT2Test. Este teste não é robusto
heterocedasticidade multivariada. Formas robustas à heterocedasticidade
multivariada são descritas e comparadas em Lin (1992). O pacote
SHT implementa dois testes robustas à heterocedasticidade
multivariada.
Como foi feito em ANCOVA, a saída produz uma “Tabela ANOVA da MANOVA”
que o teste de Pillai (default), sendo que o mais utilizado na
literatura é o teste de Wilks. Neste caso particular, todos os métodos
apresentam o mesmo resultado. No teste de Pillai, a coluna
test stat coincide com o \(\eta^2_{\text{parcial}}\) obtido com
effectsize::eta_squared. No teste de Wilks, este valor é
chamado de \(\lambda\) de Wilks, que é
igual a \(1-\eta^2_{\text{parcial}}\)
(estas estimativas dos dois métodos, portanto, somam um). Para qualquer
nível de significância habitual, \(H_0\) foi rejeitada.
Teste \(T^2\) de Welch por Johansen: robusto à heterocedasticidade multivariada
As funções SHT::mean2.1980Johansen e
SHT::mean2.1965Yao executam a MANOVA robusta à
heterocedasticidade multivariada. Estes testes rejeitam a hipótese nula
de igualdade de compleição entre homem e mulher deste exemplo.
source("eiras.bartitle.R")
cat(bartitle("Data"))
Dados <- readRDS("Adm2008.rds")
Dados <- Dados[, c("Sexo","Estatura","MCT")]
result <- try(MVN::mvn(data = Dados,
subset="Sexo",
# R=5e4,
# mvnTest = "energy",
mvn_test = "hz",
univariate_test = "SW"))
print(result$multivariate_normality)
print(res <- heplots::boxM(Dados[,2:3], Dados[,1]))
plot(res, main="Adm2008")
# Robust One-way MANOVA (heterocedástica)
Dados.Mm <- as.matrix(subset(Dados, Sexo=="Masculino", select=c(Estatura, MCT)))
Dados.Fm <- as.matrix(subset(Dados, Sexo=="Feminino", select=c(Estatura, MCT)))
colMeans(Dados.Mm)
colMeans(Dados.Fm)
print(SHT::mean2.1980Johansen(Dados.Mm, Dados.Fm))
print(SHT::mean2.1965Yao(Dados.Mm, Dados.Fm))
print(SHT::mean2.1931Hotelling(Dados.Mm, Dados.Fm))
# Na saída, T.2 é F
print(DescTools::HotellingsT2Test(cbind(Estatura,
MCT)~Sexo,
test="f",
data=Dados))
----
Data
----
Group Test Statistic p.value MVN
1 Feminino Henze-Zirkler 0.287 0.877 ✓ Normal
2 Masculino Henze-Zirkler 1.126 0.010 ✗ Not normal
Box's M-test for Homogeneity of Covariance Matrices
data: Dados[, 2:3]
Chi-Sq (approx.) = 19.918, df = 3, p-value = 0.0001765
Two-sample Test for Multivariate Means by Johansen (1980)
data: Dados.Mm and Dados.Fm
T2 = 104.91, p-value = 7.291e-15
alternative hypothesis: true means are different.
Two-sample Test for Multivariate Means by Yao (1965).
data: Dados.Mm and Dados.Fm
T2 = 104.91, p-value = 3.442e-15
alternative hypothesis: true means are different.
Hotelling's T-squared Test for Independent Samples with Equal
Covariance Assumption.
data: Dados.Mm and Dados.Fm
T2 = 90.304, p-value = 5.064e-14
alternative hypothesis: true means are different.
Hotelling's two sample T2-test
data: cbind(Estatura, MCT) by Sexo
T.2 = 44.633, df1 = 2, df2 = 86, p-value = 5.063e-14
alternative hypothesis: true location difference is not equal to c(0,0)MANOVA com três ou mais condições independentes
O \(T^2\) de Hotelling permite duas ou mais medidas, mas o fator é limitado em dois grupos.
MANOVA (assim como t em relação à one-way ANOVA)
permite mais de dois grupos. Para exemplificar, usaremos o famoso
arquivo iris, usando código similar ao
anterior.
Estes dados estão sempre disponíveis em R (veja a documentação com
?iris) e contém as medidas de comprimento e largura de
sépalas e pétalas (quatro medidas) de flores de três espécies de
Iris sp (fator com três níveis: setosa,
versicolor e virginica).
A estatística descritiva está implementada em demo_iris_descritiva.R:
source("eiras.friendlycolor.R")
source("eiras.correg.R")
source("eiras.bivCI.R")
source("eiras.showdataframe.R")
source("eiras.numeric.summary.R")
source("eiras.bartitle.R")
colM <- friendlycolor(7) # azul
colF <- friendlycolor(30) # vermelho
pchM <- 24
pchF <- 21
alpha <- 0.05
cat(bartitle("Data"))
showdataframe(iris, head=5, tail=3)
cat(bartitle("Summary"))
print(
psych::describeBy(iris[,c("Sepal.Length","Sepal.Width",
"Petal.Length","Petal.Width")],
group=iris$Species,
mat=1,
digits=2)
)
cat(bartitle("Correlation & covariance"))
species <- unique(iris$Species)
colors <- c("#666666","#888888","#cccccc")
for (s in species)
{
cat(bartitle(s,2))
cat(bartitle("Correlation matrix",3))
print(cor(subset(iris,Species=s,select=-Species)), digits=2)
cat(bartitle("Covariance matrix",3))
print(cov(subset(iris,Species=s,select=-Species)), digits=2)
}
Dados.long <- reshape2::melt(iris,
id=c("Species"),
measured=c("Sepal.Length",
"Sepal.Width",
"Petal.Length",
"Petal.Width"))
names(Dados.long) <- c('Group', 'Measure', 'Value')
Dados.long.Boxplot <- ggplot2::ggplot(Dados.long,
ggplot2::aes(Group,
Value,
color = Measure)) +
ggplot2::geom_boxplot() +
ggplot2::labs(x="Species",
y="Value (cm)",
color="Medida") +
ggplot2::theme_bw()+
ggplot2::theme(panel.grid =
ggplot2::element_blank())
plot(Dados.long.Boxplot)
rm(Dados.long)
suppressMessages(suppressWarnings(
invisible(print(GGally::ggpairs(data = iris,
mapping = ggplot2::aes(color = Species,
#fill = Species,
alpha = .5),
progress = FALSE) +
ggplot2::theme_bw()+
ggplot2::theme(panel.grid =
ggplot2::element_blank())))))
# grafico
car::scatterplotMatrix(iris[,c("Sepal.Length",
"Sepal.Width",
"Petal.Length",
"Petal.Width")],
groups=iris$Species,
regLine=FALSE,
smooth=FALSE,
boxplots=TRUE,
by.groups=TRUE,
ellipse=list(levels=c(0.95),
robust=TRUE,
fill=FALSE),
grid=FALSE,
col=colors,
cex=0.5,
cex.lab=1)
# elipses de confianca do centroide populacional
for (i1 in 1:3)
{
for (i2 in (i1+1):4)
{
dt_tmp <- iris[,c(i1,i2,5)]
dt_tmp <- na.omit(dt_tmp)
n1 <- names(dt_tmp)[1]
n2 <- names(dt_tmp)[2]
car::scatterplot(dt_tmp[,n1], dt_tmp[,n2],
xlab=n1, ylab=n2,
group=dt_tmp[,3],
regLine=list(lwd=0.5),
ellipse=FALSE,
smooth=FALSE,
grid=FALSE,
col=colors,
legend=TRUE,
data=dt_tmp)
for (s in 1:length(species))
{
dt_tmp2 <- dt_tmp[dt_tmp[3]==as.character(species[s]),c(1,2)]
lines(bivCI(s = var(dt_tmp2),
xbar = colMeans(dt_tmp2),
n = nrow(dt_tmp2),
alpha = alpha/3, m = 10000),
type = "l", col = colors[s], lwd = 2)
cx <- mean(dt_tmp2[,1], na.rm=TRUE)
cy <- mean(dt_tmp2[,2], na.rm=TRUE)
points(cx,cy,pch=21,col=colors[s],bg=colors[s])
}
}
}
----
Data
----
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
5.1 3.5 1.4 0.2 setosa
4.9 3 1.4 0.2 setosa
4.7 3.2 1.3 0.2 setosa
4.6 3.1 1.5 0.2 setosa
5 3.6 1.4 0.2 setosa
... ... ... ... ...
6.5 3 5.2 2 virginica
6.2 3.4 5.4 2.3 virginica
5.9 3 5.1 1.8 virginica
-------
Summary
-------
item group1 vars n mean sd median trimmed mad min max
Sepal.Length1 1 setosa 1 50 5.01 0.35 5.00 5.00 0.30 4.3 5.8
Sepal.Length2 2 versicolor 1 50 5.94 0.52 5.90 5.94 0.52 4.9 7.0
Sepal.Length3 3 virginica 1 50 6.59 0.64 6.50 6.57 0.59 4.9 7.9
Sepal.Width1 4 setosa 2 50 3.43 0.38 3.40 3.42 0.37 2.3 4.4
Sepal.Width2 5 versicolor 2 50 2.77 0.31 2.80 2.78 0.30 2.0 3.4
Sepal.Width3 6 virginica 2 50 2.97 0.32 3.00 2.96 0.30 2.2 3.8
Petal.Length1 7 setosa 3 50 1.46 0.17 1.50 1.46 0.15 1.0 1.9
Petal.Length2 8 versicolor 3 50 4.26 0.47 4.35 4.29 0.52 3.0 5.1
Petal.Length3 9 virginica 3 50 5.55 0.55 5.55 5.51 0.67 4.5 6.9
Petal.Width1 10 setosa 4 50 0.25 0.11 0.20 0.24 0.00 0.1 0.6
Petal.Width2 11 versicolor 4 50 1.33 0.20 1.30 1.32 0.22 1.0 1.8
Petal.Width3 12 virginica 4 50 2.03 0.27 2.00 2.03 0.30 1.4 2.5
range skew kurtosis se
Sepal.Length1 1.5 0.11 -0.45 0.05
Sepal.Length2 2.1 0.10 -0.69 0.07
Sepal.Length3 3.0 0.11 -0.20 0.09
Sepal.Width1 2.1 0.04 0.60 0.05
Sepal.Width2 1.4 -0.34 -0.55 0.04
Sepal.Width3 1.6 0.34 0.38 0.05
Petal.Length1 0.9 0.10 0.65 0.02
Petal.Length2 2.1 -0.57 -0.19 0.07
Petal.Length3 2.4 0.52 -0.37 0.08
Petal.Width1 0.5 1.18 1.26 0.01
Petal.Width2 0.8 -0.03 -0.59 0.03
Petal.Width3 1.1 -0.12 -0.75 0.04
------------------------
Correlation & covariance
------------------------
------
setosa
------
------------------
Correlation matrix
------------------ Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length 1.00 -0.12 0.87 0.82
Sepal.Width -0.12 1.00 -0.43 -0.37
Petal.Length 0.87 -0.43 1.00 0.96
Petal.Width 0.82 -0.37 0.96 1.00
-----------------
Covariance matrix
----------------- Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length 0.686 -0.042 1.27 0.52
Sepal.Width -0.042 0.190 -0.33 -0.12
Petal.Length 1.274 -0.330 3.12 1.30
Petal.Width 0.516 -0.122 1.30 0.58
----------
versicolor
----------
------------------
Correlation matrix
------------------ Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length 1.00 -0.12 0.87 0.82
Sepal.Width -0.12 1.00 -0.43 -0.37
Petal.Length 0.87 -0.43 1.00 0.96
Petal.Width 0.82 -0.37 0.96 1.00
-----------------
Covariance matrix
----------------- Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length 0.686 -0.042 1.27 0.52
Sepal.Width -0.042 0.190 -0.33 -0.12
Petal.Length 1.274 -0.330 3.12 1.30
Petal.Width 0.516 -0.122 1.30 0.58
---------
virginica
---------
------------------
Correlation matrix
------------------ Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length 1.00 -0.12 0.87 0.82
Sepal.Width -0.12 1.00 -0.43 -0.37
Petal.Length 0.87 -0.43 1.00 0.96
Petal.Width 0.82 -0.37 0.96 1.00
-----------------
Covariance matrix
----------------- Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length 0.686 -0.042 1.27 0.52
Sepal.Width -0.042 0.190 -0.33 -0.12
Petal.Length 1.274 -0.330 3.12 1.30
Petal.Width 0.516 -0.122 1.30 0.58
Uma das novidades que aparece nesta saída é o gráfico produzido pela
função car::scatterplotMatrix que exibe density
plots na diagonal principal e scatterplots para as
combinações dois a dois das quatro medidas. Estes dados estão associados
com elipses de predição de 95% dos pontos observados.
Suposição de normalidade multivariada
Implementado em demo_iris_normalidade.R:
source("eiras.bartitle.R")
data(iris)
result <- try(MVN::mvn(data = iris,
subset="Species",
#B=1e4,
#mvn_test = "energy",
mvn_test = "hz",
univariate_test = "SW",
multivariate_outlier_method ="none"))
print(result$multivariate_normality)
result <- try(MVN::mvn(data = iris,
subset="Species",
#B=1e4,
#mvn_test = "energy",
mvn_test = "doornik_hansen",
univariate_test = "SW",
multivariate_outlier_method ="none"))
print(result$multivariate_normality) Group Test Statistic p.value MVN
1 setosa Henze-Zirkler 0.949 0.050 ✗ Not normal
2 versicolor Henze-Zirkler 0.839 0.226 ✓ Normal
3 virginica Henze-Zirkler 0.757 0.497 ✓ Normal
Group Test Statistic p.value MVN
1 setosa Doornik-Hansen 126.558 <0.001 ✗ Not normal
2 versicolor Doornik-Hansen 20.262 0.009 ✗ Not normal
3 virginica Doornik-Hansen 24.393 0.002 ✗ Not normalNo teste de multinormalidade, a hipótese nula de tetranormalidade não foi rejeitada para as três espécies se considerarmos correção de Bonferroni (três testes simultâneos de multinormalidade, \(\alpha/3 \approx 0.167\)). Adicionalmente, podemos invocar o Teorema Central do Limite e prosseguirmos com a análise.
Suposição de homocedasticidade multivariada
Avaliamos a homogeneidade das matrizes de covariância. Implementado
em demo_iris_homocedasticidade.R:
source("eiras.bartitle.R")
print(res <- heplots::boxM(iris[,1:4], iris$Species))
plot(res, main="Iris")
Box's M-test for Homogeneity of Covariance Matrices
data: iris[, 1:4]
Chi-Sq (approx.) = 140.94, df = 20, p-value < 2.2e-16
Apesar de ser heterocedástico, os grupos são perfeitamente balanceados e, por isso, podemos prosseguir com a MANOVA.
Teste MANOVA omnibus
Implementado em demo_MANOVA_iris.R:
source("eiras.bartitle.R")
cat(bartitle("MANOVA"))
cat(bartitle("Tabela ANOVA da MANOVA",2))
fit <- lm(cbind(Sepal.Length,Sepal.Width,
Petal.Length,Petal.Width) ~ Species,
data=iris)
print(anv <- car::Anova(fit,
test="Pillai"),
digits=3)
alpha <- 0.05
eta2 <- effectsize::eta_squared(anv,
partial=TRUE,
generalized=FALSE,
ci=1-alpha,
alternative="two.sided",
verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=2)
------
MANOVA
------
----------------------
Tabela ANOVA da MANOVA
----------------------
Type II MANOVA Tests: Pillai test statistic
Df test stat approx F num Df den Df Pr(>F)
Species 2 1.19 53.5 8 290 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Effect Size for ANOVA (Type II)
Parameter | Eta2 (partial) | 95% CI | interpret
-----------------------------------------------------
Species | 0.60 | [0.52, 0.65] | largePelo teste de Pillai, rejeita-se \(H_0\).
Os supostos testes post hoc estão implementados aqui por ilustração, exibindo a saída por duas formas diferentes (para o leitor escolher sua preferida), mas (como já dissemos acima) são testes inadequados porque desprezam as correlações entre as medidas.
Teste MANOVA post hoc
A forma adequada de fazer o teste post hoc é comparar as
condições independentes (grupos) dois a dois mas considerando
conjuntamente todas as medidas usadas para o teste omnibus.
Implementado em demo_MANOVA_iris_posthoc.R:
# Testes post hoc
source("eiras.bartitle.R")
alpha <- 0.05
# setosa & versicolor
cat(bartitle("setosa & versicolor"))
Dados1 <- subset(iris, Species!="virginica")
fit <- lm(cbind(Sepal.Length, Sepal.Width,
Petal.Length, Petal.Width) ~ Species,
data=Dados1)
print(anv <- car::Anova(fit,
test="Pillai"),
digits=3)
eta2 <- effectsize::eta_squared(anv,
partial=TRUE,
generalized=FALSE,
ci=1-alpha,
alternative="two.sided",
verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=2)
# setosa & virginica
cat(bartitle("setosa & virginica"))
Dados2 <- subset(iris, Species!="versicolor")
fit <- lm(cbind(Sepal.Length, Sepal.Width,
Petal.Length, Petal.Width) ~ Species,
data=Dados2)
print(anv <- car::Anova(fit,
test="Pillai"),
digits=3)
eta2 <- effectsize::eta_squared(anv,
partial=TRUE,
generalized=FALSE,
ci=1-alpha,
alternative="two.sided",
verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=2)
# versicolor & virginica
cat(bartitle("versicolor & virginica"))
Dados3 <- subset(iris, Species!="setosa")
fit <- lm(cbind(Sepal.Length, Sepal.Width,
Petal.Length, Petal.Width) ~ Species,
data=Dados3)
print(anv <- car::Anova(fit,
test="Pillai"),
digits=3)
eta2 <- effectsize::eta_squared(anv,
partial=TRUE,
generalized=FALSE,
ci=1-alpha,
alternative="two.sided",
verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=2)
-------------------
setosa & versicolor
-------------------
Type II MANOVA Tests: Pillai test statistic
Df test stat approx F num Df den Df Pr(>F)
Species 1 0.963 625 4 95 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Effect Size for ANOVA (Type II)
Parameter | Eta2 (partial) | 95% CI | interpret
-----------------------------------------------------
Species | 0.96 | [0.95, 0.97] | large
------------------
setosa & virginica
------------------
Type II MANOVA Tests: Pillai test statistic
Df test stat approx F num Df den Df Pr(>F)
Species 1 0.98 1183 4 95 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Effect Size for ANOVA (Type II)
Parameter | Eta2 (partial) | 95% CI | interpret
-----------------------------------------------------
Species | 0.98 | [0.97, 0.98] | large
----------------------
versicolor & virginica
----------------------
Type II MANOVA Tests: Pillai test statistic
Df test stat approx F num Df den Df Pr(>F)
Species 1 0.784 86.1 4 95 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Effect Size for ANOVA (Type II)
Parameter | Eta2 (partial) | 95% CI | interpret
-----------------------------------------------------
Species | 0.78 | [0.71, 0.83] | largeUtilizando correção de Bonferroni para o nível de significância de 5% (i.e., \(\alpha \approx 0.0167\)), conclui-se que os tamanhos totais das três espécies são diferentes.
Podemos imaginar como indicador o produto do comprimento e largura das pétalas e sépalas (que são uma aproximação das respectivas áreas) para colocar as flores das três espécies em ordem de tamanho:
Conforme Grice & Iwasaki (2007), a maneira incorreta de fazer o teste post hoc é realizar uma ANOVA para cada medida.
“Resumo: Muitos pesquisadores realizam uma Análise Multivariada de Variância (MANOVA) em seus dados, mas não reconhecem plenamente a natureza verdadeiramente multivariada de seus efeitos. O erro mais comum é seguir a MANOVA com análises univariadas das variáveis dependentes. Uma das causas desse erro é a falta de materiais pedagógicos claros para identificar e testar os efeitos multivariados obtidos na análise. Este artigo revisa as diferenças fundamentais entre MANOVA e ANOVA univariada e apresenta um conjunto coerente de métodos para explorar a natureza multivariada de um conjunto de dados. Um exemplo completo, com dados reais, é fornecido, incluindo estimativas de tamanho de efeito e intervalos de confiança. Também é apresentada uma seção de resultados de exemplo, redigida segundo o estilo técnico da American Psychological Association. Diversas questões relativas aos métodos atuais também são discutidas.”
Dados <- iris
fit <- lm(cbind(Sepal.Length, Sepal.Width,
Petal.Length, Petal.Width) ~ Species,
data=Dados)
print(car::Anova(fit))
Type II MANOVA Tests: Pillai test statistic
Df test stat approx F num Df den Df Pr(>F)
Species 2 1.1919 53.466 8 290 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Response Sepal.Length :
Df Sum Sq Mean Sq F value Pr(>F)
Species 2 63.212 31.606 119.26 < 2.2e-16 ***
Residuals 147 38.956 0.265
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Response Sepal.Width :
Df Sum Sq Mean Sq F value Pr(>F)
Species 2 11.345 5.6725 49.16 < 2.2e-16 ***
Residuals 147 16.962 0.1154
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Response Petal.Length :
Df Sum Sq Mean Sq F value Pr(>F)
Species 2 437.10 218.551 1180.2 < 2.2e-16 ***
Residuals 147 27.22 0.185
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Response Petal.Width :
Df Sum Sq Mean Sq F value Pr(>F)
Species 2 80.413 40.207 960.01 < 2.2e-16 ***
Residuals 147 6.157 0.042
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1MANCOVA
Retomando o exemplo da estatura e MCT dos estudantes de
administração, podemos controlar por uma covariável intervalar que,
neste exemplo, é a idade. Implementamos em demo_MANCOVA.R:
source("eiras.bartitle.R")
Dados <- readRDS("Adm2008.rds")
cat(bartitle("MANCOVA"))
fit <- lm(cbind(Estatura,MCT)~Sexo+Idade,
data=Dados)
print(anv <- car::Anova(fit,
univariate=FALSE,
multivariate=TRUE,
test="Pillai"),
digits=3)
alpha <- 0.05
eta2 <- effectsize::eta_squared(anv,
partial=TRUE,
generalized=FALSE,
ci=1-alpha,
alternative="two.sided",
verbose=TRUE)
eta2$interpret <- effectsize::interpret_eta_squared(eta2$Eta2)
print(eta2, digits=2)
-------
MANCOVA
-------
Type II MANOVA Tests: Pillai test statistic
Df test stat approx F num Df den Df Pr(>F)
Sexo 1 0.478 38.9 2 85 1e-12 ***
Idade 1 0.135 6.7 2 85 0.0021 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Effect Size for ANOVA (Type II)
Parameter | Eta2 (partial) | 95% CI | interpret
-----------------------------------------------------
Sexo | 0.48 | [0.32, 0.59] | large
Idade | 0.14 | [0.02, 0.27] | mediumO valor p que interessa é o que está associado com
Sexo e, portanto, rejeita-se \(H_0\). A MANCOVA é, simplesmente, feita
com
lm(cbind(Estatura, MCT) ~ Sexo + Idade, data=Dados)
A VD é a composição cbind(Estatura, MCT) e a
Idade é adicionada como covariável, da mesma forma que foi
feito na ANCOVA.
Referências
- Borm, GF et al. (2007) A simple sample size formula for analysis of covariance in randomized clinical trials. Journal of Clinical Epidemiology 60: 1234-8.
- Crawley, MJ (2012) The R Book. 2nd ed.
- Grice, JW & Iwasaki, M (2007) A truly multivariate approach to MANOVA. Applied Multivariate Research 12(30): 199-226.
- Hotelling, H (1931), The generalization of Student’s ratio. Annals of Mathematical Statistics 2(3): 360-378.
- Hotelling, H (1951) A generalized T test and measure of multivariate dispersion. Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability, 23–41. University of California Press.
- Johansen, S (1980) The Welch-James approximation to the distribution of the residual sum of squares in a weighted linear regression. Biometrika 67(1). doi:10.1093/biomet/67.1.85.
- Johnson, R & Wichern, D (2007) Applied Multivariate Statistical Analysis. 6ª ed. NJ: Prentice-Hall.
- Lin, W-Y (1992) An Overview of the Performance of Four Alternatives to Hotelling’s T Square. Educational Research Journal 7:110-114.
- Moreira, NF et al. (2018) Self-reported weight and height are valid measures to determine weight status: Results from the Brazilian National Health Survey (PNS 2013). Cadernos de Saúde Pública 34(5). https://doi.org/10.1590/0102-311x00063917
- Sakae, TM et al. (2016) Efeitos da anestesia geral da cognição do idoso. Arquivos Catarinenses de Medicina 45(3): 107-16.
- Wilks, SS (1932) Certain generalizations in the analysis of variance. Biometrika 24(3/4): 471-494. https://doi.org/10.1093/biomet/24.3-4.471
- Xiaofeng, L (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