ANÁLISE DE VARIÂNCIA MULTIVARIADA - MANOVA

Outubro de 2008

OBJETIVOS DA AULA

• Usar o R para realizar análises de variância univariadas (aov) e multivariadas (manova).

• Realizar comparações simultâneas no caso de rejeição da hipótese nula de ausência de efeito de tratamento.

EXEMPLO 1

Para começar vamos trabalhar com a base de dados milk.txt.

Descrição dos dados: as unidades de observação referem-se a caminhões de transporte de leite e os dados observados são custos (combustível, consertos, capital) associados ao veículo. O fator refere-se ao tipo de combustível que pode ser gasolina ou diesel.

Dados de transporte de leite

• Primeiro, é necessário verificar se as suposições básicas do modelo são plausíveis: normalidade e variância constante.

• milk=read.table(http://www.im.ufrj.br/~flavia/mad484/milk.txt,header=T)

• Para isso vamos usar as funções Shapiro.test (verifica a normalidade dos dados) e var.test (realiza um teste de comparação das variâncias nos dois tipos de combustível).

Dados de transporte de leite

• Verificadas as suposições básicas, estamos prontos para realizar a análise de variância univariada para verificar a hipótese de não haver diferença nas médias de custo de combustível.

• comb=aov(milk$x1~milk$comb)

TABELA ANOVA PARA CUSTO DE COMBUSTÍVEL

summary(comb) g.l. SQ QM F p-valor

tratamento 1 19,96 19,96 2,7874 0,1007 resíduos 55 393,80 7,16 Total 56 413,76 Portanto, não rejeitamos a hipótese nula de igualdade entre

os custos médios de combustível.

Análise de variância do custo sobre consertos

• cons=aov(milk$x2~milk$comb)• summary(cons)• g.l. SQ QM F p-valor • tratamento 1 134,34 134.34 7,1096 0.01005 *• resíduos 55 1039,26 18.90

Portanto, ao nível de significância de 5%, rejeitamos a hipótese nula de igualdade entre as médias de custo de conserto para os dois tipos de caminhão.

Análise de variância do custo sobre capital

• cap=aov(milk$x3~milk$comb)• summary(cap)• gl SQ QM F p-valor • Tratamento 1 1016,25 1016,25 39,307 5.885e-08 • Residuals 55 1421,98 25,85

• Portanto, para esse custo também rejeitamos a hipótese nula.

Análise de variância multivariada

• Agora vamos realizar a análise de variância multivariada. Observe que aqui também é necessário verificar as suposições básicas do modelo, a saber, normalidade, variância igual e independência entre as diferentes observações.

• Será necessário carregar o pacote stats do R.

ESTATÍSTICAS PARA TESTAR A HIPÓTESE DE AUSÊNCIA DE EFEITO DE TRATAMENTO

• Vimos em aula a estatística lambda de Wilks dada pela razão entre os determinantes da matriz de somas de quadrados e produtos cruzados devida aos resíduos sobre o determinante da matriz de somas de quadrados e produtos cruzados da variação total.

• Quanto menor for o valor dessa estatística, maior a evidência a favor da hipótese nula de ausência de efeito de tratamento.

BR

R

*

ESTATÍSTICAS PARA TESTAR A HIPÓTESE DE AUSÊNCIA DE EFEITO DE TRATAMENTO

• Outras estatísticas usadas para esse teste são baseadas nos auto-valores da matriz

• Sejam • Estatística de Hotelling-Lawley:• Estatística de Pillai: • Estatística de Roy:

1BRs ,...,, 21 os respectivos auto-valores

s

jjBRtrT

1

1)(

s

j j

jRBBtrV1

1

1])([

s

j j

jU1 1

ESTATÍSTICAS PARA TESTAR A HIPÓTESE DE AUSÊNCIA DE EFEITO DE TRATAMENTO

• O R calcula todas essas estatísticas.• Voltando aos dados de transporte de leite,

suponha que após análise inicial, as suposições básicas do modelo tenham sido consideradas adequadas (normalidade, variâncias iguais e independência das observações).

MANOVA

• Após carregar o pacote stats, defina o vetor-resposta Y de dimensão 3 por:

• Y=cbind(milk$x1,milk$x2,milk$x3)• Defina o fator combustível por• classe=milk$comb• Faça então:• geral=manova(Y~classe)• geral2=summary.manova(geral)

Call: manova(Y ~ classe)Terms: classe Residualsresp 1 19.9576 393.7967resp 2 134.3407 1039.2641resp 3 1016.249 1421.979Deg. of Freedom 1 55

Residual standard error: 2.675806 4.34692 5.084699

geral2$SS$classe [,1] [,2] [,3][1,] 19.95757 -51.77947 -142.4144[2,] -51.77947 134.34071 369.4910[3,] -142.41438 369.49102 1016.2490

$Residuals [,1] [,2] [,3][1,] 393.7967 186.8572 157.6213[2,] 186.8572 1039.2641 311.6113[3,] 157.6213 311.6113 1421.9791

geralW=summary.manova(geral,test="Wilks")geralP=summary.manova(geral,test="Pillai")geralR=summary.manova(geral,test="Roy")geralHL=summary.manova(geral,test="Hotelling-Lawley")

Df Wilks approx F num Df den Df Pr(>F) classe 1 0.5122 16.8262 3 53 8.358e-08 ***

Df Pillai approx F num Df den Df Pr(>F) classe 1 0.4878 16.8262 3 53 8.358e-08 ***

Df Roy approx F num Df den Df Pr(>F) classe 1 0.9524 16.8262 3 53 8.358e-08 ***

Df Hotelling-Lawley approx F num Df den Df Pr(>F) classe 1 0.9524 16.8262 3 53 8.358e-08 ***

Residuals 55

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Resultado

• Verifica-se então que os dados não trazem evidência a favor da hipótese nula, de modo que rejeita-se H0.

Comparações Múltiplas

• Quando a hipótese de ausência de efeito de tratamento é rejeitada, os efeitos que levaram à rejeição da hipótese são de interesse.

• Para comparações duas a duas, a abordagem de Bonferroni pode ser usada para construir intervalos simultâneos de confiança para as diferenças dos efeitos de tratamento tomados dois a dois.

• Esses intervalos serão mais estreitos que os intervalos simultâneos T2 obtidos para todos os contrastes.

MODELO

k

lipliliT

il

pil

lillil

H

XXXX

tesindependenN

klniX

....:

),...,,(

.),,0(~

,....,1,,...,1,

210

,,2,1

grupo. ésimo- do

componente ésima- da média a representa que

1

e componente ésima- da global média é que

1

com ˆ

que segue ˆ Como

.grupo) ésimo-l do (efeito de componente ésima-j a é

que tal,...,2,1),,...,,(Seja

1

1 1

21

l

j

xn

x

j

xn

xxxα

xx

kl

j

iijl

ljl

k

l

n

iijljjjljl

ll

ljl

plllT

l

l

p jnn

)XXVar()ααVar(

xx

jjsr

jsjrjsjr

jsjrjsjr

,...,2,1,11

ˆˆ

que Observe

tes.independen amostrais

médias duas entre diferença a é ˆˆ

Portanto,

etc.) classes,ou grupos, de (número

fator do níveis de número o ék e sobservaçõe

de totalnúmero o én R, matriz da diagonal da

elemento ésimo-j o é que em 11

:por dado é médias essas entre

diferença da variânciadaestimador Um

jjjj

sr

rkn

r

nn

mtxx

mt

pkkm

)k(kk

knjsjrjsjr

kn

21:)1,(IC

que tal2

1

obteremos e 2

)1( teremos,Bonferroni de

abordagem a usando pares, 2

1

2 são Como

B

No exemplo de transporte de leite, ahipótese nula foi rejeitada. Obtenha os intervalos de confiança deBonferroni.

Observe que como k=2 e p=3, teremos ao todo 3 contrastes a serem analisados, referindo-se às diferenças nas médias de cada uma das três componentes.

Como exercício obtenha os três intervalos e tire Suas conclusões.

Como segunda atividade vamos analisar os dados crabs sobre medidas morfológicas de duas espécies de caranguejos.

Será necessário carregar o pacote MASS para obter os dados.