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

options(warn=-1)
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
suppressMessages(library(psych, warn.conflicts=FALSE))
suppressMessages(library(emmeans, warn.conflicts=FALSE))
suppressMessages(library(multcomp, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(GGally, warn.conflicts=FALSE))
suppressMessages(library(effects, warn.conflicts=FALSE))
suppressMessages(library(effectsize, warn.conflicts=FALSE))
suppressMessages(library(tidyr, warn.conflicts=FALSE))
suppressMessages(library(MASS, warn.conflicts=FALSE))
suppressMessages(library(car, warn.conflicts=FALSE))
suppressMessages(library(MVN, warn.conflicts=FALSE))
suppressMessages(library(heplots, warn.conflicts=FALSE))
suppressMessages(library(rrcov, warn.conflicts=FALSE))
suppressMessages(library(DescTools, warn.conflicts=FALSE))
suppressMessages(library(SHT, warn.conflicts=FALSE))
suppressMessages(library(mvdalab, warn.conflicts=FALSE))
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")

Material

HTML de R Markdown em RPubs

Incluir

out <- bruceR::MANOVA(data = Data,
                      dv=c("Assessment"),
                      between=c("Severity", "Complexity", "Experience"))
bruceR::EMMEANS(out, "Experience")

Objetivos

Conceituar MANOVA, distinguindo-a de ANOVA e ANCOVA.
Conceituar e testar as suposições da MANOVA.
Proceder com a análise descritiva numérica e gráfica dos dados.
Formular e implementar:
- \(T^2\) de Hotelling
- MANOVA unifatorial independente
- MANCOVA
Determinar as significâncias estatística e prática do efeito do fator.

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óides de medidas intervalares multivariadas observadas em duas condições independentes.
MANOVA: Análise de Variância Multivariada Independente
MANCOVA: Análise de Covariância Multivariada
rmMANOVA unifatorial: Análise de Variância Multivariada para Medidas Repetidas (um fator intraparticipantes)
MANOVA mista (um fator entre participantes e um fator intraparticipantes)

O teste \(T^2\) de Hotelling é o primeiro e mais básico dos testes multivariados. Originalmente, era possível testar a igualdade de dois centróides multivariados populacionais.

Em MANOVA, o efeito de um fator com três ou mais condições pode ser testado e, em MANCOVA, incluiremos covariáveis. Com estes testes, completamos 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 grandes suposições são de normalidade multivariada e homocedasticidade multivariada das medidas nos níveis do fator.

MANOVA unifatorial independente

MANOVA 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:

options(warn=-1) # disable warnings

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)

options(warn=0) # enable warnings


----
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 centroíde populacional. Como a hipótese nula está fora desta região elíptica, podemos esperar que a hipótese nula será rejeitada adiante.

Suposições

A MANOVA supõe normalidade multivariada das medidas em cada tratamento (neste caso bivariada porque só temos duas medidas: estatura e MCT, nas condições masculino e feminino). Os testes estão implementados em demo_Adm2008Masc_normalidade.R:

source("eiras.bartitle.R")

options(warn=-1) # disable warnings

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",
                   multivariate_outlier_method="quan")
cat(bartitle("Binormality test"))
print(result$multivariate_normality)
cat(bartitle("Normality test"))
print(result$univariate_normality)
cat(bartitle("Multivariate Outlier"))
print(result$multivariate_outliers)
options(warn=0) # enable warnings


----------------
Binormality test
----------------
           Test Statistic p.value     Method          MVN
1 Henze-Zirkler     1.126    0.01 asymptotic ✗ Not normal

--------------
Normality test
--------------
          Test Variable Statistic p.value Normality
1 Shapiro-Wilk Estatura     0.986   0.795  ✓ Normal
2 Shapiro-Wilk      MCT     0.978   0.465  ✓ Normal

--------------------
Multivariate Outlier
--------------------
  Observation Mahalanobis.Distance
1          64               26.749
2          57               23.862
3          73               22.523
4          71               22.104
5          87               11.860

Estes testes foram realizados com MVN::mvn e mostram os testes de uninormalidade e binormalidade para o nível de significância de 5%.

Isoladamente, as distribuições são 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.

Aplicação do teste

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:

\[H_0: \left[ \mu_{Est}~~~\mu_{MCT} \right] = \left[ 172~cm~~~76.68~kg \right]\] \[H_1: \left[ \mu_{Est}~~~\mu_{MCT} \right] \ne \left[ 172~cm~~~76.68~kg \right]\]

Implementado em demo_MANOVA_umacondicao.R:

source("eiras.bartitle.R")
source("eiras.bivCI.R")

options(warn=-1) # disable warnings

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) 
summary(T2One)

# 'mcd' for minimum covariance determinant estimator
print(rrcov::T2.test(Masc, 
               mu = H0.M, 
               conf.level = 1-alfa,
               method="mcd"))

T2crit <- ((n.M-1)*p/(n.M-p))*qf(1-alfa, p, n.M-p)
eigen_result <- eigen(cov.M)
eigvl <- eigen_result$values
a <- sqrt(T2crit*eigen_result$values[1]/n.M) # (5-19)
b <- sqrt(T2crit*eigen_result$values[2]/n.M) # (5-19)
theta <- atan2(eigen_result$vectors[2, 1], eigen_result$vectors[1, 1])
theta_seq <- seq(0, 2 * pi, length.out = 100)
ellipse_x <- mean.M[1] + a * cos(theta_seq) * cos(theta) - 
  b * sin(theta_seq) * sin(theta)
ellipse_y <- mean.M[2] + a * cos(theta_seq) * sin(theta) + 
  b * sin(theta_seq) * cos(theta)
plot(ellipse_x, ellipse_y, type = "l", xlim=c(170,180),
     xlab="Estatura (m)", 
     ylab="MCT (kg)",
     main="Região elíptica de confiança de 95%\nEstudante de Administração Masculino",
     col="black")
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]))
text(mean.M[1],mean.M[2], pos=2, "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

    One-sample Hotelling test (Reweighted MCD Location)

data:  Masc
T2 = 113.557, F = 45.465, df1 = 2.000, df2 = 25.687, p-value =
3.641e-09
alternative hypothesis: true mean vector is not equal to (172, 76.68)' 

sample estimates:
             Estatura      MCT
MCD x-vector 176.2609 71.97826

              [,1]      [,2]
Estatura 173.46882 179.23706
MCT       70.15031  78.79087

Rejeita-se \(H_0\) para qualquer nível de significância convencional. Os estudantes 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.

MANOVA com duas condições independentes

O exemplo de estatura e MCT tem um fator entre participantes (sexo) e duas medidas (estatura e MCT). Iniciando pela estatística descritiva, sugere-se o que está implementado em demo_Adm2008_descritiva.R:

options(warn=-1) # disable warnings

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") 
plot(Dados.long.Boxplot)
rm(Dados.long)

print(GGally::ggpairs(data = Dados, mapping = ggplot2::aes(color = Sexo,
                                                                fill = Sexo,
                                                                alpha = .5),
                      progress = FALSE))

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)

options(warn=0) # enable warnings


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

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

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% 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 distantes: quando não há sobreposição, podemos esperar que a hipótese nula de igualdade dos centróides populacionais será rejeitada.

Suposições

Normalidade multivariada e univariada

Implementado em demo_Adm2008_normalidade.R:

source("eiras.bartitle.R")

options(warn=-1) # disable warnings

Dados <- readRDS("Adm2008.rds")
Dados <- Dados[, c("Sexo","Estatura","MCT")]

result <- MVN::mvn(data = Dados,
                   subset="Sexo",
                   # R=5e4,
                   # mvnTest = "energy",
                   univariate_test  = "SW",
                   multivariate_outlier_method="quan")
cat(bartitle("Binormality test"))
print(result$multivariate_normality)
cat(bartitle("Normality test"))
print(result$univariate_normality)
cat(bartitle("Multivariate Outlier"))
print(result$multivariate_outliers)
options(warn=0) # enable warnings


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

--------------
Normality test
--------------
      Group         Test Variable Statistic p.value Normality
1  Feminino Shapiro-Wilk Estatura     0.977   0.626  ✓ Normal
2  Feminino Shapiro-Wilk      MCT     0.979   0.695  ✓ Normal
3 Masculino Shapiro-Wilk Estatura     0.986   0.795  ✓ Normal
4 Masculino Shapiro-Wilk      MCT     0.978   0.465  ✓ Normal

--------------------
Multivariate Outlier
--------------------
      Group Observation Mahalanobis.Distance
1  Feminino          26               14.506
2 Masculino          64               26.749
3 Masculino          57               23.862
4 Masculino          73               22.523
5 Masculino          71               22.104
6 Masculino          87               11.860

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

“O teste da suposição de multinormalidade, em geral, é necessário para amostra pequena. Se a amostra é suficientemente grande (maior que ‘20 + número de medidas’ para cada condição entre participantes), o TCL entra em ação e a suposição de multinormalidade não é mais necessária.”

Johnson & Wichern, 2007, p. 294

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

Aplicação do teste

A hipótese nula a seguir é a igualdade dos centróides populacionais. Formalmente é:

\[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]\]

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:

\[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}\]

Implementado em demo_MANOVA_independente.R:

source("eiras.bartitle.R")

options(warn=-1) # disable warnings

cat(bartitle("Data"))
Dados <- readRDS("Adm2008.rds")
Dados <- Dados[, c("Sexo","Estatura","MCT")]

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

cat(bartitle("MANOVA (homocedástica)"))
# Solução 2: manova
fit <- manova(cbind(Estatura,MCT)~Sexo,
              data = Dados)
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)


----
Data
----

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


----------------------
MANOVA (homocedástica)
----------------------

Type II MANOVA Tests: Pillai test statistic
     Df test stat approx F num Df den Df  Pr(>F)    
Sexo  1     0.509     44.6      2     86 5.1e-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] |     large

Utilizamos, como em capítulos anteriores, a função lm nativa do R, chamada com a sintaxe:

manova(cbind(Estatura, MCT) ~ Sexo,  
       data=Dados)

O modelo utilizado emprega uma matriz das medidas combinadas (cbind(Estatura, MCT)) em função do sexo.

Observe o que faz a função cbind:

nomes <- c("Juvêncio", "Palmério", "Sofrêncio", "Coríntio")
idades <- c(34,23,56,10)
juntos <- cbind(nomes,idades)
cat(paste("juntos é data frame? ",is.data.frame(juntos),"\n",sep=""))
cat(paste("juntos é matrix? ",is.matrix(juntos),"\n",sep=""))
cat("juntos contém:\n")
print(juntos)

juntos é data frame? FALSE

juntos é matrix? TRUE

juntos contém:

     nomes       idades
[1,] "Juvêncio"  "34"  
[2,] "Palmério"  "23"  
[3,] "Sofrêncio" "56"  
[4,] "Coríntio"  "10"

Observe que cbind cria uma matriz e, como tal, todos os elementos precisam ser do mesmo tipo. Como o conteúdo do vetor a é alfanumérico, apesar de b ser numérico, todos os elementos foram “demovidos” a character. No entanto, é fácil trabalhar com esta matriz se quiser tratar a coluna de idades como numérica:

juntos <- as.data.frame(juntos) # converte em data frame
juntos$idades <- as.numeric(juntos$idades) # converte a coluna em numérica
cat(paste("juntos é data frame? ",is.data.frame(juntos),"\n",sep=""))
cat(paste("juntos é matrix? ",is.matrix(juntos),"\n",sep=""))
cat("juntos contém:\n")
print(juntos)

juntos é data frame? TRUE

juntos é matrix? FALSE

juntos contém:

      nomes idades
1  Juvêncio     34
2  Palmério     23
3 Sofrêncio     56
4  Coríntio     10

A hipótese nula de igualdade das compleições de homens e mulheres foi rejeitada pelos testes \(T^2\) de Hotelling e por 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 à falta de normalidade e homocedasticidade multivariadas em amostras pequenas e desbalanceadas. Formas robustas são descritas e comparadas em Lin (1992). O pacote SHT implementa as formas robustas do artigo e outras que apareceram posteriormente na literatura.

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.

Adicionalmente, podem ser computadas ANOVAs univariadas, que são ANOVAS unifatoriais para cada medida. A primeira parte é uma forma sintética para obter a ANOVA de cada medida. A segunda parte foi feita com os comandos habituais que vimos no capítulo de ANOVA unifatorial. Note que os valores p desta segunda forma coincidem com os valores p não ajustados da primeira forma.

GRICE & IWASAKI, 2007

MANOVA robusta

Existe função para executar MANOVA utilizando bootstrapping, implementada em demo_MANOVA_robusta.R:

source("eiras.bartitle.R")

options(warn=-1) # disable warnings

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",
                       multivariate_outlier_method ="quan"))
print(result$multivariate_normality)
print(result$univariate_normality)
print(result$multivariate_outliers)

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

# Robust One-way MANOVA (Bartlett Chi2) (heterocedástica e não-multinormal)
mnv <- rrcov::Wilks.test(x=Dados[,c("Estatura","MCT")],
                         grouping=Dados[,"Sexo"],
                         method="mcd",
                         nrep=1e3)
print(mnv)


----
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
      Group         Test Variable Statistic p.value Normality
1  Feminino Shapiro-Wilk Estatura     0.977   0.626  ✓ Normal
2  Feminino Shapiro-Wilk      MCT     0.979   0.695  ✓ Normal
3 Masculino Shapiro-Wilk Estatura     0.986   0.795  ✓ Normal
4 Masculino Shapiro-Wilk      MCT     0.978   0.465  ✓ Normal
      Group Observation Mahalanobis.Distance
1  Feminino          26               14.506
2 Masculino          64               26.749
3 Masculino          57               23.862
4 Masculino          73               22.523
5 Masculino          71               22.104
6 Masculino          87               11.860

    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.


    Robust One-way MANOVA (Bartlett Chi2)

data:  x
Wilks' Lambda = 0.44685, Chi2-Value = 36.6821, DF = 1.6254, p-value =
5.406e-09
sample estimates:
          Estatura      MCT
Feminino  163.7568 54.72973
Masculino 176.8182 72.29545

A função rrcov::Wilks.test computa a MANOVA reportando apenas o \(\lambda\) de Wilks, que corresponde a \(1-\eta^2_{\text{parcial}}\) e rejeita-se a hipótese nula de igualdade de compleição entre homens e mulheres deste exemplo.

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 no RStudio (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:

options(warn=-1) # disable warnings

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="Valor (cm)", 
                color="Medida") 
plot(Dados.long.Boxplot)
rm(Dados.long)

print(GGally::ggpairs(data = iris, mapping = ggplot2::aes(color = Species,
                                                           fill = Species,
                                                           alpha = .5),
                      progress = FALSE))

# 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])
    }
  }
}

options(warn=0) # enable warnings


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

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

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ções

Normalidade univariada e 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)
print(result$univariate_normality)
print(result$multivariate_outliers)

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)
print(result$univariate_normality)
print(result$multivariate_outliers)

       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     Variable Statistic p.value    Normality
1      setosa Shapiro-Wilk Sepal.Length     0.978    0.46     ✓ Normal
2      setosa Shapiro-Wilk  Sepal.Width     0.972   0.272     ✓ Normal
3      setosa Shapiro-Wilk Petal.Length     0.955   0.055     ✓ Normal
4      setosa Shapiro-Wilk  Petal.Width     0.800  <0.001 ✗ Not normal
5  versicolor Shapiro-Wilk Sepal.Length     0.978   0.465     ✓ Normal
6  versicolor Shapiro-Wilk  Sepal.Width     0.974   0.338     ✓ Normal
7  versicolor Shapiro-Wilk Petal.Length     0.966   0.158     ✓ Normal
8  versicolor Shapiro-Wilk  Petal.Width     0.948   0.027 ✗ Not normal
9   virginica Shapiro-Wilk Sepal.Length     0.971   0.258     ✓ Normal
10  virginica Shapiro-Wilk  Sepal.Width     0.967   0.181     ✓ Normal
11  virginica Shapiro-Wilk Petal.Length     0.962    0.11     ✓ Normal
12  virginica Shapiro-Wilk  Petal.Width     0.960   0.087     ✓ Normal
NULL
       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 normal
        Group         Test     Variable Statistic p.value    Normality
1      setosa Shapiro-Wilk Sepal.Length     0.978    0.46     ✓ Normal
2      setosa Shapiro-Wilk  Sepal.Width     0.972   0.272     ✓ Normal
3      setosa Shapiro-Wilk Petal.Length     0.955   0.055     ✓ Normal
4      setosa Shapiro-Wilk  Petal.Width     0.800  <0.001 ✗ Not normal
5  versicolor Shapiro-Wilk Sepal.Length     0.978   0.465     ✓ Normal
6  versicolor Shapiro-Wilk  Sepal.Width     0.974   0.338     ✓ Normal
7  versicolor Shapiro-Wilk Petal.Length     0.966   0.158     ✓ Normal
8  versicolor Shapiro-Wilk  Petal.Width     0.948   0.027 ✗ Not normal
9   virginica Shapiro-Wilk Sepal.Length     0.971   0.258     ✓ Normal
10  virginica Shapiro-Wilk  Sepal.Width     0.967   0.181     ✓ Normal
11  virginica Shapiro-Wilk Petal.Length     0.962    0.11     ✓ Normal
12  virginica Shapiro-Wilk  Petal.Width     0.960   0.087     ✓ Normal
NULL

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

Homocedasticidade multivariada

Avaliamos a homogeneidade das matrizes de covariância. Implementado em demo_iris_homocedasticidade.R:

source("eiras.bartitle.R")

options(warn=-1) # disable warnings

print(res <- heplots::boxM(iris[,1:4], iris$Species))
plot(res, main="Iris")

options(warn=0) # enable warnings


    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.

Aplicação do teste

Implementado em demo_MANOVA_iris.R:

source("eiras.bartitle.R")
options(warn=-1) # disable warnings

cat(bartitle("MANOVA"))
cat(bartitle("Tabela ANOVA da MANOVA",2))
fit <- manova(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] |     large

Pelo 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 fracos porque desprezam as correlações entre as medidas.

Testes post hoc

A forma adequada de fazer o teste é comparar os 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 <- manova(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 <- manova(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 <- manova(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] |     large

Utilizando 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:

MANOVA robusta

A versão robusta, utilizando bootstrapping, implementada em demo_MANOVA_robusta_iris.R:

source("eiras.bartitle.R")

options(warn=-1) # disable warnings

cat(bartitle("MANOVA robusta"))
mnv <- rrcov::Wilks.test(x=iris[,1:4],
                          grouping=iris[,"Species"],
                          method="mcd",nrep=1e3)
print(mnv)
data(iris)
result <- try(MVN::mvn(data = iris,
                       subset="Species",
                       # R=5e4,
                       # mvnTest = "energy",
                       mvn_test = "doornik_hansen",
                       univariate_test = "SW",
                       multivariate_outlier_method ="none"))
print(result$multivariate_normality)
print(result$univariate_normality)
print(result$multivariate_outliers)

print(res <- heplots::boxM(iris[,1:4], iris[,"Species"]))
plot(res, main="Iris")

# Robust One-way MANOVA (Bartlett Chi2) (heterocedástica e não-multinormal)
mnv <- rrcov::Wilks.test(x=iris[,1:4],
                         grouping=iris[,"Species"],
                         method="mcd",nrep=1e3)
print(mnv)


--------------
MANOVA robusta
--------------

    Robust One-way MANOVA (Bartlett Chi2)

data:  x
Wilks' Lambda = 0.013115, Chi2-Value = 393.5119, DF = 6.6648, p-value <
2.2e-16
sample estimates:
           Sepal.Length Sepal.Width Petal.Length Petal.Width
setosa         4.985106    3.419149     1.470213   0.2425532
versicolor     5.946809    2.789362     4.263830   1.3212766
virginica      6.415152    2.984848     5.381818   2.0606061

       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 normal
        Group         Test     Variable Statistic p.value    Normality
1      setosa Shapiro-Wilk Sepal.Length     0.978    0.46     ✓ Normal
2      setosa Shapiro-Wilk  Sepal.Width     0.972   0.272     ✓ Normal
3      setosa Shapiro-Wilk Petal.Length     0.955   0.055     ✓ Normal
4      setosa Shapiro-Wilk  Petal.Width     0.800  <0.001 ✗ Not normal
5  versicolor Shapiro-Wilk Sepal.Length     0.978   0.465     ✓ Normal
6  versicolor Shapiro-Wilk  Sepal.Width     0.974   0.338     ✓ Normal
7  versicolor Shapiro-Wilk Petal.Length     0.966   0.158     ✓ Normal
8  versicolor Shapiro-Wilk  Petal.Width     0.948   0.027 ✗ Not normal
9   virginica Shapiro-Wilk Sepal.Length     0.971   0.258     ✓ Normal
10  virginica Shapiro-Wilk  Sepal.Width     0.967   0.181     ✓ Normal
11  virginica Shapiro-Wilk Petal.Length     0.962    0.11     ✓ Normal
12  virginica Shapiro-Wilk  Petal.Width     0.960   0.087     ✓ Normal
NULL

    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


    Robust One-way MANOVA (Bartlett Chi2)

data:  x
Wilks' Lambda = 0.013115, Chi2-Value = 476.2760, DF = 7.9217, p-value <
2.2e-16
sample estimates:
           Sepal.Length Sepal.Width Petal.Length Petal.Width
setosa         4.985106    3.419149     1.470213   0.2425532
versicolor     5.946809    2.789362     4.263830   1.3212766
virginica      6.415152    2.984848     5.381818   2.0606061

A conclusão é a mesma neste exemplo.

MANCOVA

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

options(warn=-1) # disable warnings

Dados <- readRDS("Adm2008.rds")

cat(bartitle("MANCOVA"))

fit <- manova(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] |    medium

O valor p que interessa é o que está associado com Sexo e, portanto, rejeita-se \(H_0\). A MANCOVA é, simplesmente, feita com

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

Análise de Variância Multivariada (MANOVA)

em R

José Siqueira (siqueira@usp.br)
Paulo Silveira (silveira@usp.br)

07 julho 2025 18:12h

Material

Incluir

Objetivos

Disciplina de pós-graduação do PBE-UEM

O que é MANOVA

Por que utilizar MANOVA?

Os testes básicos e o GLM

A lógica da MANOVA

As suposições subjacentes ao uso da MANOVA

MANOVA unifatorial independente

MANOVA para uma condição

Suposições

Aplicação do teste

MANOVA com duas condições independentes

Suposições

Normalidade multivariada e univariada

Homocedasticidade multivariada

Aplicação do teste

MANOVA robusta

MANOVA com três ou mais condições independentes

Suposições

Normalidade univariada e multivariada

Homocedasticidade multivariada

Aplicação do teste

Testes post hoc

MANOVA robusta

MANCOVA

Referências

Análise de Variância Multivariada (MANOVA)

em R

José Siqueira (siqueira@usp.br) Paulo Silveira (silveira@usp.br)

07 julho 2025 18:12h

Material

Incluir

Objetivos

Disciplina de pós-graduação do PBE-UEM

O que é MANOVA

Por que utilizar MANOVA?

Os testes básicos e o GLM

A lógica da MANOVA

As suposições subjacentes ao uso da MANOVA

MANOVA unifatorial independente

MANOVA para uma condição

Suposições

Aplicação do teste

MANOVA com duas condições independentes

Suposições

Normalidade multivariada e univariada

Homocedasticidade multivariada

Aplicação do teste

MANOVA robusta

MANOVA com três ou mais condições independentes

Suposições

Normalidade univariada e multivariada

Homocedasticidade multivariada

Aplicação do teste

Testes post hoc

MANOVA robusta

MANCOVA

Referências

José Siqueira (siqueira@usp.br)
Paulo Silveira (silveira@usp.br)