14  ANOVA desbalanceada com um fator

Nesta seção, trataremos da ANOVA com um fator em conjuntos de dados desbalanceados, em que o número de observações por tratamento não é o mesmo.

O modelo desbalanceado com um fator é dado por:

\[ \begin{align} y_{ij} = \mu + \tau_i + \epsilon_{ij} = \mu_{i} + \epsilon_{ij} \\ i = 1,\dots,k \space \text{ e } \space j = 1,\dots,n_i \end{align} \]

Para realizar inferências, assumiremos \(\epsilon_{ij} \sim N(0, \sigma^2)\) são independentes e identicamente distribuídos.

Como exemplo, temos o seguinte caso: Os pesos líquidos de latas enchidas por cinco máquinas de enchimento são apresentados:

A	B	C	D	E
11.95	12.18	12.16	12.25	12.10
12.00	12.11	12.15	12.3	12.04
12.25	-	12.08	12.1	12.02
12.10	-	-	-	12.02

O modelo matricial de médias de caselas \(\mathbf{y} = \mathbf{W}\boldsymbol{\mu} + \boldsymbol{\epsilon}\) é dado por:

\[ \begin{bmatrix} y_{11} \\ y_{12} \\ y_{13} \\ y_{14} \\ y_{21} \\ y_{22} \\ y_{31} \\ y_{32} \\ y_{33} \\ y_{41} \\ y_{42} \\ y_{43} \\ y_{51} \\ y_{52} \\ y_{53} \\ y_{54} \end{bmatrix} = \begin{bmatrix} 11,95 \\ 12,00 \\ 12,25 \\ 12,10 \\ 12,18 \\ 12,11 \\ 12,16 \\ 12,15 \\ 12,08 \\ 12,25 \\ 12,30 \\ 12,10 \\ 12,10 \\ 12,04 \\ 12,02 \\ 12,02 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} \mu_1 \\ \mu_2 \\ \mu_3 \\ \mu_4 \\ \mu_5 \end{bmatrix} + \begin{bmatrix} \epsilon_{11} \\ \epsilon_{12} \\ \epsilon_{13} \\ \epsilon_{14} \\ \epsilon_{21} \\ \epsilon_{22} \\ \epsilon_{31} \\ \epsilon_{32} \\ \epsilon_{33} \\ \epsilon_{41} \\ \epsilon_{42} \\ \epsilon_{43} \\ \epsilon_{51} \\ \epsilon_{52} \\ \epsilon_{53} \\ \epsilon_{54} \end{bmatrix} \]

R

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

y <- c(11.95,12.00,12.25,12.10,12.18,12.11,12.16,12.15,12.08,12.25,12.30,12.10,12.10,12.04,12.02,12.02)
y

 [1] 11.95 12.00 12.25 12.10 12.18 12.11 12.16 12.15 12.08 12.25 12.30 12.10
[13] 12.10 12.04 12.02 12.02

W <- matrix(
  c(rep(1, 4), rep(0, 12),
  rep(0, 4), rep(1, 2), rep(0, 10),
  rep(0, 6), rep(1, 3), rep(0, 7),
  rep(0, 9), rep(1, 3), rep(0, 4),
  rep(0, 12), rep(1, 4)),
  ncol = 5, byrow = FALSE
)
W

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

Note que a matriz \(\mathbf{W}\) é de posto completo (posto\((\mathbf{W}) = k = 5\)).

k <- ncol(W)                         # Número de tratamentos
k

[1] 5

rank_W <- sum(diag(ginv(W) %*% W))   # posto da matriz W
rank_W

[1] 5

deficit_rank <- k - rank_W           # déficit de rank
deficit_rank

[1] 0

Dessa forma, o sistema de equações normais \(\mathbf{(W'W)} \boldsymbol{\hat\mu} = \mathbf{(W'y)}\) tem solução única:

\[ \boldsymbol{\hat\mu} = \mathbf{(W'W)}^{-1} \mathbf{W'y} = \mathbf{\bar{y}} = \begin{bmatrix} \bar{y}_{1.} \\ \bar{y}_{2.} \\ \bar{y}_{3.} \\ \bar{y}_{4.} \\ \bar{y}_{5.} \end{bmatrix} \]

SAS

options nodate nocenter ps=1000;
proc iml;
reset fuzz;
* Exemplo 15.2.1 Os pesos líquidos de latas enchidas por cinco máquinas; 
* de enchimento (filling machines) são apresentados na Tabela 14.2;

y = {11.95,12.00,12.25,12.10,12.18,12.11,12.16,12.15,
     12.08,12.25,12.30,12.10,12.10,12.04,12.02,12.02};

Trat = {1,1,1,1,2,2,3,3,3,4,4,4,5,5,5,5};
W = design(Trat);   * W do modelo de médias de caselas;
k = ncol(W);        * Número de tratamentos;
N = nrow(W);        * Número total de repetições;
print W y[format=10.2];

14.1 Teste de hipótese

Para testar a hipótese \(H_0: \mu_1 = \mu_2 = \mu_3 = \mu_4 = \mu_5\), podemos utilizar a abordagem do modelo completo x modelo reduzido e contrastes.

14.1.1 Modelo completo x Modelo reduzido

R

N <- nrow(W)   # Número total de repetições
N

[1] 16

Jnn <- matrix(1, nrow = N, ncol = N)
In <- diag(1, nrow = N)

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

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

# Total
SQTotal <- t(y) %*% (In - Jnn / N) %*% y
SQTotal

          [,1]
[1,] 0.1441437

gl_total <- N - 1
gl_total

[1] 15

A soma de quadrados do modelo completo de média de caselas é dado por:

\[ SQ(\mu_1,\mu_2,\mu_3,\mu_4,\mu_5) = \boldsymbol{\hat\mu' W'y} \]

# Modelo completo
mi <- solve(t(W) %*% W) %*% t(W) %*% y
mi

         [,1]
[1,] 12.07500
[2,] 12.14500
[3,] 12.13000
[4,] 12.21667
[5,] 12.04500

SQcompleto <- t(mi) %*% t(W) %*% y
SQcompleto

         [,1]
[1,] 2347.704

O modelo reduzido pode ser escrito, matricialmente, da seguinte maneira:

\[ \mathbf{y} = \mu \mathbf{j} + \boldsymbol{\epsilon}^*, \quad \text{onde } \mathbf{j} \text{ é } N \times 1 \]

em que \(\mathbf{j}\) é um vetor coluna de 1’s com dimensão \(N \times 1\), sendo \(N\) o número total de observações (sem contar os valores ausentes).

Dessa forma, a soma de quadrados do modelo reduzido é dado por:

\[ SQ(\mu) = \hat\mu \mathbf{j}' \mathbf{y} \]

# Modelo reduzido
jn <- matrix(1, ncol = 1, nrow = N)
jn

      [,1]
 [1,]    1
 [2,]    1
 [3,]    1
 [4,]    1
 [5,]    1
 [6,]    1
 [7,]    1
 [8,]    1
 [9,]    1
[10,]    1
[11,]    1
[12,]    1
[13,]    1
[14,]    1
[15,]    1
[16,]    1

mi_reduzido <- solve(t(jn) %*% jn) %*% t(jn) %*% y
mi_reduzido

         [,1]
[1,] 12.11313

SQreduzido <- t(mi_reduzido) %*% t(jn) %*% y
SQreduzido

         [,1]
[1,] 2347.645

Com as somas de quadrados dos modelos completo e reduzido, obtemos a soma de quadrados entre grupos, calculada pela diferença entre a soma de quadrados do modelo completo (SQcompleto) e a soma de quadrados do modelo reduzido (SQreduzido). Além disso, apresenta \((k-1)\) graus de liberdade.

# Entre grupos
SQEntre <- SQcompleto - SQreduzido
SQEntre

           [,1]
[1,] 0.05942708

gl_entre <- k - 1
gl_entre

[1] 4

QMEntre <- SQEntre / gl_entre
QMEntre

           [,1]
[1,] 0.01485677

Já a soma de quadrados dos resíduos (SQRes) é dada por:

\[ SQRes = \mathbf{y'y} - \boldsymbol{\hat\mu' W'y} \]

# Resíduo
SQRes <- t(y) %*% y - t(mi) %*% t(W) %*% y
SQRes

           [,1]
[1,] 0.08471667

gl_res <- N - k
gl_res

[1] 11

QMRes <- SQRes / gl_res
QMRes

            [,1]
[1,] 0.007701515

Com os quadrados médios, calcularemos a estatística F e o p-valor.

Fcalc <- QMEntre / QMRes
Fcalc

         [,1]
[1,] 1.929071

Ftab <- qf(0.95, gl_entre, gl_res)
Ftab

[1] 3.35669

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

          [,1]
[1,] 0.1756589

A seguir, o quadro da ANOVA compila os resultados obtidos.

ANOVA desbalanceada: H0: m1=m2=m3=m4=m5

FV	gl	SQ	QM	Fcal	Ftab	p-valor (5%)
Hipótese	4	0.0594	0.0149	1.9291	3.3567	0.1757
Resíduo	11	0.0847	0.0077
Total	15	0.1441

Dado que o \(F_{cal} < F_{tab}\), considerando \(F_{(0,05;4;11)}\), não se rejeita \(H_0: \mu_1 = \mu_2 = \mu_3 = \mu_4 = \mu_5\), ou seja, não existem diferenças significativas entre as médias ponderadas dos pesos líquidos de latas enchidas pelas cinco máquinas. Ainda, podemos dizer que os pesos líquidos médios das latas enchidas pelas cinco máquinas não diferem entre si.

SAS

Jnn = J(N,N,1);
In = I(N);
SQTotal = t(y)*(In-Jnn/N)*y;
gl_total = N-1;

mi = inv(t(W)*W)*t(W)*y;
SQcompleto = t(mi)*t(W)*y;      * SQ do modelo completo: yij = mi(i) + eij;

jn = J(N,1,1);
mir = inv(t(jn)*jn)*t(jn)*y;
SQreduzido = t(mir)*t(jn)*y;    * SQ do modelo reduzido: yij = mi + eij;

SQEntre = SQCompleto - SQreduzido;
gl_entre = k-1;
QMEntre = SQEntre/gl_entre;

SQRes = t(y)*y - t(mi)*t(W)*y;
gl_res = N-K;
QMRes = SQRes/gl_res;

Fcalc = QMEntre/QMRes;
p_valor = 1 - cdf('F', Fcalc,gl_entre, gl_res);

print '------------------------------',
      'Exemplo 15.2.1 Quadro de ANOVA',
      '------------------------------';
print 'Ho: m1=m2=m3=m4=m5'  gl_entre SQEntre[format=10.5] QMEntre[format=10.5] Fcalc[format=8.4] p_valor[format=8.3],,
      'Resíduo           '  gl_res   SQRes[format=10.5]   QMRes[format=10.5],,
      'Total             '  gl_total SQTotal[format=10.5]; 

14.1.2 Contrastes

Definimos um contraste de \(k\) médias populacionais como:

\[ \delta = c_1\mu_1 + c_2\mu_2 + \dots + c_k\mu_k = \boldsymbol{c'\mu} \]

em que \(\sum^k_{i=1} c_i = 0\).

Quando trabalhamos com ANOVA balanceada, dois contrastes

\[ \hat\delta = \sum^k_{i=1} a_i \bar{y}_{i.} \quad \text{e} \quad \hat\gamma = \sum^k_{i=1} b_i \bar{y}_{i.} \]

são ditos ortogonais se \(\sum^k_{i=1} a_i b_i = 0\), sendo esta a condição de independência entre os contrastes.

Para o caso desbalanceado, devemos utilizar contrastes ortogonais ponderados, afim de garantir a independência entre os contrastes. Dessa forma, os contrastes \(\hat\delta\) e \(\hat\gamma\) são independentes em modelos desbalanceados se e somente se

\[ \sum^k_{i=1} \frac{a_i b_i}{n_i} = 0 \]

sendo \(n_i\) as respostas ao \(i\)-ésimo tratamento.

Contudo, na prática, os contrastes ortogonais ponderados são de menor interesse que os contrastes ortogonais não ponderados, não sendo necessário que as somas de quadrados sejam independentes para realizarmos os testes sobre os contrastes.

A seguir, veremos as diferenças entre os contrastes ortogonais não ponderados e os ponderados.

Contrastes ortogonais não ponderados

R

Os contrastes a seguir comparam as médias dos tratamentos da seguinte maneira:

\[ \begin{align} H_{a01}&: 3\mu_1 - 2\mu_2 - 2\mu_3 + 3\mu_4 - 2\mu_5 = 0 \\ H_{a02}&: \mu_2 - 2\mu_3 + \mu_5 = 0 \\ H_{a03}&: \mu_1 - \mu_4 = 0 \\ H_{a04}&: \mu_2 - \mu_5 = 0 \\ \end{align} \]

a1 <- c(3,-2,-2, 3,-2)
a2 <- c(0, 1,-2, 0, 1)
a3 <- c(1, 0, 0,-1, 0)
a4 <- c(0, 1, 0, 0,-1)

As somas de quadrados são calculadas a partir da seguinte forma quadrática:

\[ SQ(\hat\delta) = \boldsymbol{(c'\mu)}' [\mathbf{c'(W'W)c}]^{-1} \boldsymbol{(c'\mu)} \]

em que \(\mathbf{c}\) é o vetor de coeficientes do contraste.

SQa1 <- t(t(a1) %*% mi) %*% solve(t(a1) %*% solve(t(W) %*% W) %*% a1) %*% t(a1) %*% mi
F_a1 <- SQa1 / QMRes
p_valor_a1 <- 1 - pf(F_a1, 1, gl_res)

SQa2 <- t(t(a2) %*% mi) %*% solve(t(a2) %*% solve(t(W) %*% W) %*% a2) %*% t(a2) %*% mi
F_a2 <- SQa2 / QMRes
p_valor_a2 <- 1 - pf(F_a2, 1, gl_res)

SQa3 <- t(t(a3) %*% mi) %*% solve(t(a3) %*% solve(t(W) %*% W) %*% a3) %*% t(a3) %*% mi
F_a3 <- SQa3 / QMRes
p_valor_a3 <- 1 - pf(F_a3, 1, gl_res)

SQa4 <- t(t(a4) %*% mi) %*% solve(t(a4) %*% solve(t(W) %*% W) %*% a4) %*% t(a4) %*% mi
F_a4 <- SQa4 / QMRes
p_valor_a4 <- 1 - pf(F_a4, 1, gl_res)

ANOVA: Contrastes ortogonais não ponderados

Contrastes	SQ	Fcal	p-valor (5%)
A,D vs. B,C,E	0.0058	0.7482	0.4055
B,E vs. C	0.0024	0.3054	0.5916
A vs. D	0.0344	4.4673	0.0582
B vs. E	0.0133	1.7313	0.2150

Dado que nenhum p-valor é menor que \(\alpha = 0,05\), não rejeitamos qualquer uma das hipóteses \(H_0: \sum_i c_i\mu_i = 0\) associadas aos contrastes definidos anteriormente.

Ao somar as somas de quadrados dos contrastes ortogonais não ponderados, não obtemos o mesmo resultado da soma de quadrados entre grupos (SQEntre)

SQContrastes <- SQa1 + SQa2 + SQa3 + SQa4
SQContrastes

          [,1]
[1,] 0.0558527

SQEntre

           [,1]
[1,] 0.05942708

SAS

* Contrastes ortogonais do tipo t(ai)*mi;
a1 = {3,-2,-2, 3,-2};
a2 = {0, 1,-2, 0, 1};
a3 = {1, 0, 0,-1, 0};
a4 = {0, 1, 0, 0,-1};

SQA1 = t(t(a1)*mi)*inv(t(a1)*inv(t(W)*W)*a1)*t(a1)*mi;
F_A1 = SQA1/QMRes;
p_valor_A1 = 1 - cdf('F', F_A1,1, gl_res);

SQa2 = t(t(a2)*mi)*inv(t(a2)*inv(t(W)*W)*a2)*t(a2)*mi;
F_a2 = SQa2/QMRes;
p_valor_a2 = 1 - cdf('F', F_a2,1, gl_res);

SQa3 = t(t(a3)*mi)*inv(t(a3)*inv(t(W)*W)*a3)*t(a3)*mi;
F_a3 = SQa3/QMRes;
p_valor_a3 = 1 - cdf('F', F_a3,1, gl_res);

SQa4 = t(t(A4)*mi)*inv(t(A4)*inv(t(W)*W)*A4)*t(A4)*mi;
F_A4 = SQA4/QMRes;
p_valor_A4 = 1 - cdf('F', F_A4,1, gl_res);

print '-------------------------------------',
      'Contrastes ortogonais não ponderados:',
      '-------------------------------------',
      'A,D vs. B,C,E' SQA1[format=10.5] F_A1[format=12.4] p_valor_A1[format=12.3],,
      'B,E vs. C    ' SQA2[format=10.5] F_A2[format=12.4] p_valor_A2[format=12.3],,
      'A vs. D      ' SQA3[format=10.5] F_A3[format=12.4] p_valor_A3[format=12.3],,
      'B vs. E      ' SQA4[format=10.5] F_A4[format=12.4] p_valor_A4[format=12.3],,;

SQContrastes = SQA1 + SQA2 + SQA3 + SQA4;
print '----------------------------------------------------------------',
      'SQContrastes = SQA1 + SQA2 + SQA3 + SQA4 não é igual a SQEntre: ',,
       SQContrastes[format=10.5] SQEntre[format=10.5],,
      'porque os contrastes NÃO SÃO ORTOGONAIS!',
      '----------------------------------------------------------------',,,;

Contrastes ortogonais ponderados

R

Utilizaremos o primeiro contraste ortogonal do exemplo anterior, \(H_{a01}: 3\mu_1 - 2\mu_2 - 2\mu_3 + 3\mu_4 - 2\mu_5 = 0\), e o contraste ortogonal de \(H_0: 2\mu_2 - 6\mu_3 + 4\mu_5 = 0\) para ilustrar os contrastes ortogonais ponderados. Portanto:

\[ \begin{align} \delta &= \boldsymbol{a'\mu} = [3, -2, -2, 3, -2] \boldsymbol{\mu} \\ \gamma &= \boldsymbol{b'\mu} = [0, 2, -6, 0, 4] \boldsymbol{\mu} \end{align} \]

Neste caso, temos contrastes ortogonais ponderados:

\[ \sum^k_{i=1} \frac{a_ib_i}{n_i} = \frac{3(0)}4 - \frac{2(2)}2 - \frac{2(-6)}3 + \frac{3(0)}3 - \frac{2(4)}4 = 0 \]

a2p <- c(0, 2, -6, 0, 4)

SQa2p <- t(t(a2p) %*% mi) %*% solve(t(a2p) %*% solve(t(W) %*% W) %*% a2p) %*% t(a2p) %*% mi

F_a2p <- SQa2p / QMRes

p_valor_a2p <- 1 - pf(F_a2p, 1, gl_res)

Assim, as somas de quadrados associadas aos dois contrastes são independentes, cujos valores estão descritos no seguinte quadro de ANOVA.

ANOVA: Contrastes ortogonais ponderados

Contrastes	SQ	Fcal	p-valor (5%)
A,D vs. B,C,E	0.0058	0.7482	0.4055
2B + 4E vs. 6C	0.0053	0.6932	0.4228

SAS

A2p = {0, 2, -6, 0, 4};
SQA2p = t(t(a2p)*mi)*inv(t(a2p)*inv(t(W)*W)*a2p)*t(a2p)*mi;
F_a2p = SQa2p/QMRes;
p_valor_a2p = 1 - cdf('F', F_a2p,1, gl_res);

print '---------------------------------',
      'Contrastes ortogonais ponderados:',
      '---------------------------------',,
      'A,D vs. B,C,E' SQA1[format=8.6]  F_A1[format=12.4]  p_valor_A1[format=8.3],,
      '2B+4E vs. 6C ' SQA2p[format=8.6] F_A2p[format=12.4] p_valor_A2p[format=8.3],,,,;
quit;