12 ANOVA balanceada com um fator - Teste de Hipótese

Para exemplificar o teste de hipótese em modelos ANOVA balanceada com um fator, considere o seguinte caso: Três métodos (A, B, C) de armazenar alimentos congelados foram comparados. A variável resposta é a quantidade de ácido ascórbico (mg/100g) e os dados estão representados a seguir:

A	B	C
14.29	20.06	20.04
19.10	20.64	26.23
19.09	18.00	22.74
16.25	19.56	24.04
15.09	19.47	23.37
16.61	19.07	25.02
19.63	18.38	23.27

Como suposições do modelo \(\mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}\) consideraremos:

\(E(\boldsymbol{\epsilon}) = \mathbf{0}\);
\(var(\boldsymbol{\epsilon}) = \sigma^2 \mathbf{I}\);
\(cov(\mathbf{\boldsymbol{\epsilon}}) = 0\);
\(\boldsymbol{\epsilon} \sim N(\mathbf{0}, \sigma^2 \mathbf{I})\).

Novamente, utilizaremos a função ginv() do pacote MASS.

install.packages("MASS")
library(MASS)

12.1 Estimação de parâmetros

y <- as.vector(
  c(14.29, 19.10, 19.09, 16.25, 15.09, 16.61, 19.63,
    20.06, 20.64, 18.00, 19.56, 19.47, 19.07, 18.38,
    20.04, 26.23, 22.74, 24.04, 23.37, 25.02, 23.27)
)
y

 [1] 14.29 19.10 19.09 16.25 15.09 16.61 19.63 20.06 20.64 18.00 19.56 19.47
[13] 19.07 18.38 20.04 26.23 22.74 24.04 23.37 25.02 23.27

A matriz de delineamento \(\mathbf{X}\) do modelo matricial \(\mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}\) é dada por:

X <- matrix(c(
  rep(1,21),
  rep(1,7),rep(0,14),
  rep(0,7),rep(1,7),rep(0,7),
  rep(0,14),rep(1,7)
),
ncol = 4, byrow = FALSE)
X

      [,1] [,2] [,3] [,4]
 [1,]    1    1    0    0
 [2,]    1    1    0    0
 [3,]    1    1    0    0
 [4,]    1    1    0    0
 [5,]    1    1    0    0
 [6,]    1    1    0    0
 [7,]    1    1    0    0
 [8,]    1    0    1    0
 [9,]    1    0    1    0
[10,]    1    0    1    0
[11,]    1    0    1    0
[12,]    1    0    1    0
[13,]    1    0    1    0
[14,]    1    0    1    0
[15,]    1    0    0    1
[16,]    1    0    0    1
[17,]    1    0    0    1
[18,]    1    0    0    1
[19,]    1    0    0    1
[20,]    1    0    0    1
[21,]    1    0    0    1

No objeto kn, guardaremos o número de valores observados; no p, o número de variáveis; e no k, o posto de \(\mathbf{X}\).

kn <- length(y)
kn

[1] 21

p <- ncol(X)
p

[1] 4

k <- sum(diag(X %*% ginv(X)))
k

[1] 3

p - k

[1] 1

Note que a diferença entre p - k resulta no déficit de rank. Como o resultado é 1, \(\mathbf{X}\) é uma matriz de posto incompleto.

Calculando \(\mathbf{X'X}\) e \(\mathbf{X'y}\), temos:

XLX <- t(X) %*% X
XLX

     [,1] [,2] [,3] [,4]
[1,]   21    7    7    7
[2,]    7    7    0    0
[3,]    7    0    7    0
[4,]    7    0    0    7

XLy <- t(X) %*% y
XLy

       [,1]
[1,] 419.95
[2,] 120.06
[3,] 135.18
[4,] 164.71

Com os resultados anteriores, calcularemos duas estimativas de \(\boldsymbol{\beta}\) utilizando duas inversas generalizadas distintas:

\[\boldsymbol{\hat\beta} = \mathbf{(X'X)^- X' y}\]

Inversa generalizada de Moore-Penrose:

Beta <- ginv(XLX) %*% XLy
Beta

          [,1]
[1,] 14.998214
[2,]  2.153214
[3,]  4.313214
[4,]  8.531786

Inversa generalizada simples (Searle):

XLX

     [,1] [,2] [,3] [,4]
[1,]   21    7    7    7
[2,]    7    7    0    0
[3,]    7    0    7    0
[4,]    7    0    0    7

igXLX <- (1/7) * matrix(
  c(0, 0, 0, 0,
    0, 1, 0, 0,
    0, 0, 1, 0,
    0, 0, 0, 1),
  ncol = 4, byrow = TRUE
)
igXLX |> fractions()

     [,1] [,2] [,3] [,4]
[1,]   0    0    0    0 
[2,]   0  1/7    0    0 
[3,]   0    0  1/7    0 
[4,]   0    0    0  1/7

Beta2 <- igXLX %*% XLy
Beta2

         [,1]
[1,]  0.00000
[2,] 17.15143
[3,] 19.31143
[4,] 23.53000

A partir das duas soluções de \(\mathbf{\beta}'s\), verificaremos se as seguintes funções são estimáveis:

\(\beta = \alpha_1 - \alpha_2\):

L1 <- c(0, 1, -1, 0)
L1

[1]  0  1 -1  0

L1Beta <- t(L1) %*% Beta
L1Beta

      [,1]
[1,] -2.16

L1Beta2 <- t(L1) %*% Beta2
L1Beta2

      [,1]
[1,] -2.16

\(\beta = \alpha_1 + \alpha_2 + \alpha_3\):

L2 <- c(0, 1, 1, 1)
L2

[1] 0 1 1 1

L2Beta <- t(L2) %*% Beta
L2Beta

         [,1]
[1,] 14.99821

L2Beta2 <- t(L2) %*% Beta2
L2Beta2

         [,1]
[1,] 59.99286

Como se pode notar, \(\beta = \alpha_1 - \alpha_2\) é uma função estimável e \(\beta = \alpha_1 + \alpha_2 + \alpha_3\) é uma função não estimável.

options nodate nocenter ps=1000;

proc iml;
*reset print;
reset fuzz;
y = {14.29,19.10,19.09,16.25,15.09,16.61,19.63,
     20.06,20.64,18.00,19.56,19.47,19.07,18.38,
     20.04,26.23,22.74,24.04,23.37,25.02,23.27};
X = {1 1 0 0,1 1 0 0,1 1 0 0,1 1 0 0,1 1 0 0,1 1 0 0,1 1 0 0,
     1 0 1 0,1 0 1 0,1 0 1 0,1 0 1 0,1 0 1 0,1 0 1 0,1 0 1 0,
     1 0 0 1,1 0 0 1,1 0 0 1,1 0 0 1,1 0 0 1,1 0 0 1,1 0 0 1};
print X[format=5.0] y[format=8.2];
kn = nrow(y); 
p = ncol(X); 
k = round(trace(X*ginv(X))); * Calcula o posto de XlinhaX;
print 'rank(X) =' k;
XLX = t(X)*X;
XLy = t(X)*y;
print 'Sistema de equações normais:' XLX XLy;

Beta = ginv(XLX)*XLy;   * Usando inversa generalizada de Moore-Penrose;

igXLX = (1/7)*{0 0 0 0,
               0 1 0 0, 
               0 0 1 0, 
               0 0 0 1}; 
Beta2 = igXLX*XLy;
print 'Duas soluções:' Beta[format=12.4] Beta2[format=12.4];

L = {0, 1,-1,0}; * Lambda da função alfa1-alfa2 (estimável!);
LBeta = t(L)*Beta;
LBeta2 = t(L)*Beta2;
print 'Função estimável: alfa1-alfa2:',,LBeta[format=12.4] LBeta2[format=12.4];

L = {0,1,1,1}; * Lambda da função alfa1+alfa2+alfa3 (NÃO estimável!);
LBeta = t(L)*Beta;
LBeta2 = t(L)*Beta2;
print 'Função NÃO estimável: alfa1+alfa2+alfa3:',,LBeta[format=12.4] LBeta2[format=12.4];

12.2 Teste de hipótese

A seguir, testaremos a seguinte hipótese, seguindo o modelo reparametrizado de médias de caselas, \(y_{ij} = \mu_{i} + \epsilon_{ij}\), em que \(\mu_i = \mu + \alpha_i\):

\[ \begin{align} H_0&: \mu_1 = \mu_2 = \mu_3 \\ H_a&: \text{Pelo menos duas médias diferem entre si} \end{align} \]

Como \(\mu_i = \mu + \alpha_i\), a hipótese anterior é equivalente a seguinte:

\[ \begin{align} H_0&: \alpha_1 = \alpha_2 = \alpha_3 \\ H_a&: \text{Pelo menos dois tratamentos diferem entre si} \end{align} \]

Utilizaremos duas abordagens para a realização do teste de hipótese: modelo completo x modelo reduzido e contrastes.

12.2.1 Modelo Completo x Modelo Reduzido

In <- diag(kn)
Jn <- matrix(1, nrow = kn, ncol = kn)

A soma de quadrados total (SQTotal) é dada por:

\[ SQTotal = \mathbf{y'} \left(\mathbf{I} - \frac{1}n \mathbf{J}\right) \mathbf{y} \]

E seus graus de liberdade, pelo traço do produto de sua matriz núcleo pela inversa generalizada.

# Total
Tot <- In - (1 / kn) * Jn

SQTotal <- t(y) %*% Tot %*% y
SQTotal

         [,1]
[1,] 202.3126

gl_total <- round(sum(diag(Tot %*% ginv(Tot))))
gl_total

[1] 20

A soma de quadrados para os \(\alpha's\) ajustada para \(\mu\) (\(SQ(\alpha \mid \mu)\)) pode ser expressa como uma forma quadrática de \(\mathbf{y}\):

\[ SQ(\alpha \mid \mu) = \mathbf{y'} \left[\mathbf{X(X'X)^-X'} - \frac{1}n \mathbf{J}\right] \mathbf{y} \]

# Tratamentos
A <- X %*% ginv(t(X) %*% X) %*% t(X) - (1 / kn) * Jn

SQTrat <- t(y) %*% A %*% y
SQTrat

         [,1]
[1,] 147.3456

gl_trat <- round(sum(diag(A %*% ginv(A))))
gl_trat

[1] 2

QMTrat <- SQTrat / gl_trat
QMTrat

        [,1]
[1,] 73.6728

A soma de quadrados dos resíduos (SQRes) é dada por:

\[ SQRes = \mathbf{y'} \left[\mathbf{I} - \mathbf{X(X'X)^-X'}\right] \mathbf{y} \]

# Resíduo
B <- In - X %*% ginv(t(X) %*% X) %*% t(X)

SQRes <- t(y) %*% B %*% y
SQRes

         [,1]
[1,] 54.96697

gl_res <- round(sum(diag(B %*% ginv(B))))
gl_res

[1] 18

QMRes <- SQRes / gl_res
QMRes

         [,1]
[1,] 3.053721

Com os quadrados médios de tratamento (QMTrat) e dos resíduos (QMRes), calculamos a estatística F, bem como o p-valor.

Fcalc <- QMTrat / QMRes
Fcalc

         [,1]
[1,] 24.12559

ftab <- qf(0.95, gl_trat, gl_res)
ftab

[1] 3.554557

p_valor <- 1 - pf(Fcalc, gl_trat, gl_res)
p_valor

               [,1]
[1,] 0.000008066967

A seguir, os resultado estão descritos no quadro de ANOVA:

ANOVA para os dados de ácido ascórbico. H0: mi1 = mi2 = mi3
FV	gl	SQ	QM	Fcal	Ftab	p-valor (5%)
Tratamento	2	147.346	73.673	24.126	3.555	0
Resíduo	18	54.967	3.054
Total	20	202.313

Dado que \(F_{cal} > F_{tab}\), para \(F_{(0,05;2;18)}\), rejeita-se a hipótese \(H_0: \mu_1 = \mu_2 = \mu_3\), indicando que as médias de pelo menos dois métodos de congelamento diferem entre si.

T = I(kn)-(1/kn)*J(kn,kn,1);    * Matriz núcleo;
SQTotal = t(y)*T*y;             * Calcula SQTotal corrigida pela média;
gl_total = round(trace(T*ginv(T))); 

A = X*ginv(t(X)*X)*t(X) -(1/kn)*J(kn,kn,1);     * Matriz núcleo;
SQTrat = t(y)*A*y;                              * Calcula SQTrat;
gl_trat = round(trace(A*ginv(A)));
QMTrat = SQTrat/gl_trat;

B = I(kn) - X*ginv(t(X)*X)*t(X);    * Matriz núcleo;
SQRes = t(y)*B*y;                   * Calcula SQResiduo;
gl_res = round(trace(B*ginv(B)));
QMRes = SQRes/gl_res;

Fcalc = QMTrat/QMRes;
p_valor = 1-cdf('F',Fcalc,gl_trat,gl_Res);

print 'TABELA 12.3 - ANOVA para os dados de ácido ascórbico da Tabela 13.2',; 
print 'Método ' gl_trat  SQTrat[format=12.4] QMTrat[format=12.4] Fcalc[format=10.4] p_valor[format=12.4],,
      'Resíduo' gl_res   SQRes[format=12.4]  QMRes[format=12.4],,
      'Total  ' gl_total SQTotal[format=12.4];

12.2.2 Contrastes

Caso 1: As linhas são linearmente independentes (l.i.) e ortogonais

Considerando os constrastes ortogonais \(2\mu_1 - \mu_2 - \mu_3\) e \(\mu_2 - \mu_3\), podemos representá-los da seguinte maneira:

\[ \begin{align} H_{01} &:2\mu_1 - \mu_2 - \mu_3 = 2\alpha_1 - \alpha_2 - \alpha_3 = [0,2,-1,-1]\boldsymbol{\beta} = \mathbf{c'_1}\boldsymbol{\beta}\\ H_{02} &:\mu_2 - \mu_3 = \alpha_2 - \alpha_3 = [0,0,1,-1]\boldsymbol{\beta} = \mathbf{c'_2}\boldsymbol{\beta} \end{align} \]

As hipóteses \(H_{01}: \mathbf{c'_1}\boldsymbol{\beta} = 0\) e \(H_{02}: \mathbf{c'_2}\boldsymbol{\beta} = 0\) comparam a média do primeiro tratamento com a dos outros dois e a média do segundo tratamento com a do terceiro, respectivamente.

\[ \begin{align} H_{01} &:\mathbf{c'_1}\boldsymbol{\beta} = 0 \Leftrightarrow \mu_1 = \frac{\mu_2 + \mu_3}2 \\ H_{02} &: \mathbf{c'_2}\boldsymbol{\beta} = 0 \Leftrightarrow \mu_2 = \mu_3 \end{align} \]

C1 <- c(0, 2, -1, -1)
C1Beta <- C1 %*% Beta
C1Beta

          [,1]
[1,] -8.538571

C2 <- c(0, 0, 1, -1)
C2Beta <- C2 %*% Beta
C2Beta

          [,1]
[1,] -4.218571

# Partição da SQTrat
SQC1Beta <- t(C1Beta) %*% solve(t(C1) %*% ginv(t(X) %*% X) %*% C1) %*% (C1Beta)
SQC1Beta

        [,1]
[1,] 85.0584

SQC2Beta <- t(C2Beta) %*% solve(t(C2) %*% ginv(t(X) %*% X) %*% C2) %*% (C2Beta)
SQC2Beta

         [,1]
[1,] 62.28721

SomaSQ <- SQC1Beta + SQC2Beta

SQTrat

         [,1]
[1,] 147.3456

Soma das SQ de contrastes ortogonais = 147.3456 
 SQTrat = 147.3456

Aqui, note que a soma das somas de quadrados dos contrastes (SomaSQ) são iguais a soma de quadrados de tratamentos (SQTrat).

# QM
QMC1Beta <- SQC1Beta / 1
QMC1Beta

        [,1]
[1,] 85.0584

QMC2Beta <- SQC2Beta / 1
QMC2Beta

         [,1]
[1,] 62.28721

# F e p-valor
FC1 <- QMC1Beta / QMRes
FC1

         [,1]
[1,] 27.85402

FC2 <- QMC2Beta / QMRes
FC2

         [,1]
[1,] 20.39715

ftab <- qf(0.95, 1, 18)
ftab

[1] 4.413873

p_valorC1 <- 1 - pf(FC1, 1, gl_res)
p_valorC1

              [,1]
[1,] 0.00005110228

p_valorC2 <- 1 - pf(FC2, 1, gl_res)
p_valorC2

            [,1]
[1,] 0.000267208

ANOVA para contrastes ortogonais dos dados de ácido ascórbico
FV	gl	SQ	QM	Fcal	Ftab	p-valor (5%)
Contraste c'1	1	85.058	85.058	27.854	4.414	0.00005
Contraste c'2	1	62.287	62.287	20.397	4.414	0.00027
Resíduo	18	54.967	3.054
Total	20	202.313

Ambos os \(F_{cal}\) são superiores ao valor tabelado \(F_{(0,05;1,18)}\). Assim, ambas as hipóteses \(H_{01}\) e \(H_{02}\) são rejeitadas e se conclui que:

\[ \mu_1 \ne \frac{\mu_2 + \mu_3}2 \quad e \quad \mu_2 \ne \mu_3 \]

C1 ={0 2 -1 -1};
C2 ={0 0  1 -1};
C1Beta = C1*Beta;
C2Beta = C2*Beta;
SQC1Beta = t(C1*Beta)*inv(C1*ginv(t(X)*X)*t(C1))*C1*Beta;
SQC2Beta = t(C2*Beta)*inv(C2*ginv(t(X)*X)*t(C2))*C2*Beta;
Soma = SQC1Beta + SQC2Beta;
print 'Partição da SQTrat usando k=2 contrastes l.i. e ortogonais',,
      'SQcontraste1 = ' SQC1Beta[format=8.4],, 'SQcontraste2 = ' SQC2Beta[format=8.4],,
      'Soma SQs     = ' Soma[format=8.4] '    SQTrat = ' SQTrat[format=8.4];

Caso 2: As linhas são l.i. e não ortogonais

Agora, consideraremos os constrastes não ortogonais \(2\mu_1 - \mu_2 - \mu_3\) e \(\mu_1 - \mu_3\), podemos representá-los da seguinte maneira:

\[ \begin{align} 2\mu_1 - \mu_2 - \mu_3 &= 2\alpha_1 - \alpha_2 - \alpha_3 = [0,2,-1,-1]\boldsymbol{\beta} = \mathbf{c'_1}\boldsymbol{\beta}\\ \mu_1 - \mu_3 &= \alpha_1 - \alpha_3 = [0,1,0,-1]\boldsymbol{\beta} = \mathbf{c'_3}\boldsymbol{\beta} \end{align} \]

As hipóteses \(H_{01}: \mathbf{c'_1}\boldsymbol{\beta} = 0\) e \(H_{03}: \mathbf{c'_3}\boldsymbol{\beta} = 0\) são:

\[ \begin{align} H_{01} &:\mathbf{c'_1}\boldsymbol{\beta} = 0 \Leftrightarrow \mu_1 = \frac{\mu_2 + \mu_3}2 \\ H_{03} &: \mathbf{c'_3}\boldsymbol{\beta} = 0 \Leftrightarrow \mu_1 = \mu_3 \end{align} \]

C1 <- c(0, 2, -1, -1)
C1Beta <- C1 %*% Beta
C1Beta

          [,1]
[1,] -8.538571

C3 <- c(0, 1, 0, -1)
C3Beta <- C3 %*% Beta
C3Beta

          [,1]
[1,] -6.378571

# Partição da SQTrat
SQC1Beta <- t(C1Beta) %*% solve(t(C1) %*% ginv(t(X) %*% X) %*% C1) %*% (C1Beta)
SQC1Beta

        [,1]
[1,] 85.0584

SQC3Beta <- t(C3Beta) %*% solve(t(C3)%*% ginv(t(X) %*% X) %*% C3) %*% (C3Beta)
SQC3Beta

         [,1]
[1,] 142.4016

SomaSQ <- SQC1Beta + SQC3Beta

SQTrat

         [,1]
[1,] 147.3456

Soma das SQ de contrastes não ortogonais = 227.46 
 SQTrat = 147.3456

Ao contrário do caso dos constrastes ortogonais, a soma das somas de quadrados dos contrastes não ortogonais (SomaSQ) não resulta na soma de quadrados dos tratamentos (SQTrat), o que prejudica a interpretação da análise de variância, dada a não independência entre os contrastes.

C1 ={0  2 -1 -1};
C2 ={0  1  0 -1};
C1Beta = C1*Beta;
C2Beta = C2*Beta;
SQC1Beta = t(C1*Beta)*inv(C1*ginv(t(X)*X)*t(C1))*C1*Beta;
SQC2Beta = t(C2*Beta)*inv(C2*ginv(t(X)*X)*t(C2))*C2*Beta;
Soma = SQC1Beta + SQC2Beta;
print 'Partição da SQTrat usando contrastes l.i. mas não ortogonais',,
      'SQcontraste1 = ' SQC1Beta[format=8.4],, 'SQcontraste2 = 'SQC2Beta[format=8.4],,
      'Soma SQs     = ' Soma[format=8.4] '    SQTrat = ' SQTrat[format=8.4];
quit;