11 ANOVA balanceada com um fator - Estimação

11.1 Introdução

O modelo de análise de variância (ANOVA) consiste em comparar as médias de tratamentos para alguma variável resposta. O modelo matemático pode ser representado da seguinte maneira:

\[ y_{ij} = \mu + \tau_i + \epsilon_{ij} \space \text{ , para } i = 1,\dots,a \space \text{ e } \space j = 1,\dots,n \]

em que:

\(y_{ij}\) é o valor observado na \(j\)−ésima repetição do \(i\)−ésimo nível de tratamento;
\(\mu\) é uma constante comum a todas as observações, geralmente a média geral;
\(\tau_i\) é o efeito do \(i\)−ésimo nível de tratamento;
\(\epsilon_{ij}\) é o erro experimental observado na \(j\)−ésima repetição do \(i\)−ésimo nível de tratamento.

Quando o modelo é representado na forma matricial \(\mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}\), a matriz de delineamento \(\mathbf{X}\) não é de posto completo (posto incompleto), sendo assim um modelo dito superparametrizado. Dessa forma, é necessário adequar a metodologia aplicada nos modelos de regressão.

Neste primeiro momento, trataremos dos modelos ANOVA balanceados, ou seja, todos os tratamentos apresentam número igual de observações, com um fator de tratamento (one-way model).

Como exemplo, temos o seguinte caso: Um pesquisador deseja comparar dois aditivos usados para melhorar o desempenho da gasolina. Suponhamos que:

Sem aditivos: carro percorre \(\mu\) quilômetros por litro;
Com aditivo 1: acréscimo de \(tau_1\) quilômetros por litro;
Com aditivo 2: acréscimo de \(tau_2\) quilômetros por litro.

Considerando seis carros idênticos e escolhendo, aleatoriamente, três carros para adicionar o aditivo 1 e três carros para adicionar o aditivo 2, ou seja, três repetições para cada tipo de aditivo, o modelo matemático fica:

\[ y_{ij} = \mu + \tau_i + \epsilon_{ij}, \space i = 1,2 \space j = 1,2,3 \]

em que \(y_{ij}\) é a quilometragem por litro observada no \(j\)-ésimo carro que recebeu o \(i\)-ésimo aditivo.

Matricialmente, \(\mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}\), temos:

\[ \begin{cases} \mu + \tau_1 + \epsilon_{11} \\ \mu + \tau_1 + \epsilon_{12} \\ \mu + \tau_1 + \epsilon_{13} \\ \mu + \tau_2 + \epsilon_{21} \\ \mu + \tau_2 + \epsilon_{22} \\ \mu + \tau_2 + \epsilon_{23} \\ \end{cases} = \begin{bmatrix} y_{11} \\ y_{12} \\ y_{13} \\ y_{21} \\ y_{22} \\ y_{23} \end{bmatrix} = \begin{bmatrix} 1 & 1 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 0 \\ 1 & 0 & 1 \\ 1 & 0 & 1 \\ 1 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} \mu \\ \tau_1 \\ \tau_2 \end{bmatrix} + \begin{bmatrix} \epsilon_{11} \\ \epsilon_{12} \\ \epsilon_{13} \\ \epsilon_{21} \\ \epsilon_{22} \\ \epsilon_{23} \end{bmatrix} \]

Note que a matriz \(\mathbf{X}\), de dimensão \(6 \times 3\), é de posto incompleto (posto(\(\mathbf{X}\)) = 2), assim, não é possível estimar diretamente os parâmetros do vetor \(\boldsymbol{\beta} = [\mu,\tau_1,\tau_2]'\), pois apresentam infinitas soluções (não apresenta solução única).

A seguir, serão apresentados métodos para estimar os parâmetros do vetor \(\boldsymbol{\beta}\), considerando os seguintes valores de \(y_{ij}\):

Aditivo 1 (y1)	Aditivo 2 (y2)
14	18
16	19
15	17

11.2 Estimação dos parâmetros

Os dados descritos anteriormente são os seguintes:

y <- as.vector(c(14,16,15,18,19,17))
y

[1] 14 16 15 18 19 17

X <- matrix(c(
  1, 1, 0, 
  1, 1, 0,
  1, 1, 0, 
  1, 0, 1, 
  1, 0, 1, 
  1, 0, 1
),
ncol = 3, byrow = TRUE
)
X

     [,1] [,2] [,3]
[1,]    1    1    0
[2,]    1    1    0
[3,]    1    1    0
[4,]    1    0    1
[5,]    1    0    1
[6,]    1    0    1

A matriz \(\mathbf{X'X}\) (XLX) e o vetor \(\mathbf{X'y}\) (XLy) são calculados a seguir:

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

     [,1] [,2] [,3]
[1,]    6    3    3
[2,]    3    3    0
[3,]    3    0    3

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

     [,1]
[1,]   99
[2,]   45
[3,]   54

Com a matriz \(\mathbf{X'X}\), obtemos o número de observações totais (6) e o número de observações dos tratamentos (3 para cada aditivo). O vetor \(\mathbf{X'y}\) nos retorna o total dos valores observados por tratamento (45 km/L para o aditivo 1 e 54 km/L para o aditivo 2), além do total global (99 km/L).

Para calcular o posto das matrizes \(\mathbf{X}\) e \(\mathbf{X'X}\), utilizaremos a função ginv() do pacote MASS, que calcula a inversa generalizada de Moore-Penrose de uma matriz.

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

O posto de uma matriz idempotente é o traço do produto de uma matriz pela sua inversa generalizada, ou seja, é a soma dos elementos da diagonal.

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

[1] 2

rank_XLX <- sum(diag(XLX %*% ginv(XLX)))
rank_XLX

[1] 2

Fazendo a diferença entre o número de parâmetros do modelo (p) com o posto da matriz de delineamento \(\mathbf{X}\) (k), confirmamos que se trata de uma matriz de posto incompleto (defRank).

# Número de parâmetros
p <- ncol(X)
p

[1] 3

# Posto(X)
k <- rank_X
k

[1] 2

# Deficiência de rank - mostra que é de posto incompleto
defRank <- p - k
defRank

[1] 1

Para estimar o vetor de parâmetros \(\boldsymbol{\beta} = [\mu,\tau_1,\tau_2]'\), utilizamos o Método dos Mínimos Quadrados Ordinários (MQO) que busca o \(\boldsymbol{\hat\beta}\) que minimize a soma de quadrados dos desvios. O sistema de equações normais é dado por:

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

Aqui, considera-se que o modelo matricial \(\mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}\) tem \(E(\mathbf{y}) = \mathbf{X} \boldsymbol{\beta}\), \(cov(\mathbf{y}) = \sigma^2 \mathbf{I}\) e que \(\mathbf{X}\) é \(n \times p\) de posto \(k < p \le n\).

Para obter uma possível solução de \(\boldsymbol{\hat\beta}\), podemos utilizar qualquer inversa generalizada de \(\mathbf{(X'X)^-}\), a partir da seguinte operação:

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

Beta0 <- ginv(XLX) %*% t(X) %*% y
Beta0

     [,1]
[1,]   11
[2,]    4
[3,]    7

y_hat0 <- X %*% Beta0
y_hat0

     [,1]
[1,]   15
[2,]   15
[3,]   15
[4,]   18
[5,]   18
[6,]   18

Outro método de estimação de \(\boldsymbol{\hat\beta}\) é impondo retrições nos parâmetros (condições marginais). Abordaremos três tipos:

Solução com restrição \(\tau_1 = 0\): é a retrição imposta pelo software R.

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

[1] 0 1 0

W <- rbind(X, TR)
W

   [,1] [,2] [,3]
      1    1    0
      1    1    0
      1    1    0
      1    0    1
      1    0    1
      1    0    1
TR    0    1    0

z <- matrix(c(y, 0), ncol = 1)
z

     [,1]
[1,]   14
[2,]   16
[3,]   15
[4,]   18
[5,]   19
[6,]   17
[7,]    0

Beta1 <- round(solve(t(W) %*% W) %*% t(W) %*% z, 2)
Beta1

     [,1]
[1,]   15
[2,]    0
[3,]    3

y_hat1 <- X %*% Beta1
y_hat1

     [,1]
[1,]   15
[2,]   15
[3,]   15
[4,]   18
[5,]   18
[6,]   18

Solução com restrição \(\tau_2 = 0\): é a retrição imposta pelo software SAS.

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

[1] 0 0 1

W <- rbind(X, TS)
W

   [,1] [,2] [,3]
      1    1    0
      1    1    0
      1    1    0
      1    0    1
      1    0    1
      1    0    1
TS    0    0    1

z <- matrix(c(y, 0), ncol = 1)
z

     [,1]
[1,]   14
[2,]   16
[3,]   15
[4,]   18
[5,]   19
[6,]   17
[7,]    0

Beta2 <- round(solve(t(W) %*% W) %*% t(W) %*% z, 2)
Beta2

     [,1]
[1,]   18
[2,]   -3
[3,]    0

y_hat2 <- X %*% Beta2
y_hat2

     [,1]
[1,]   15
[2,]   15
[3,]   15
[4,]   18
[5,]   18
[6,]   18

Solução com restrição \(\tau_1 + \tau_2 = 0\): é a retrição comumente utilizada na estatística experimental.

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

[1] 0 1 1

W <- rbind(X, TE)
W

   [,1] [,2] [,3]
      1    1    0
      1    1    0
      1    1    0
      1    0    1
      1    0    1
      1    0    1
TE    0    1    1

z <- matrix(c(y, 0), ncol = 1)
z

     [,1]
[1,]   14
[2,]   16
[3,]   15
[4,]   18
[5,]   19
[6,]   17
[7,]    0

Beta3 <- round(solve(t(W) %*% W) %*% t(W) %*% z, 2)
Beta3

     [,1]
[1,] 16.5
[2,] -1.5
[3,]  1.5

y_hat3 <- X %*% Beta3
y_hat3

     [,1]
[1,]   15
[2,]   15
[3,]   15
[4,]   18
[5,]   18
[6,]   18

O resumo dos resultados estão descritos a seguir:

Estimação de beta de acordo com a restrição
Inv. Generaliz.	tau1 = 0	tau2 = 0	tau1 + tau2 = 0
11	15	18	16.5
4	0	-3	-1.5
7	3	0	1.5

Estimação de y (y_hat) de acordo com a restrição
y	Inv. Generaliz.	tau1 = 0	tau2 = 0	tau1 + tau2 = 0
14	15	15	15	15
16	15	15	15	15
15	15	15	15	15
18	18	18	18	18
19	18	18	18	18
17	18	18	18	18

Note que, independentemente do \(\boldsymbol{\beta}\) estimado, temos a mesma estimativa de \(\mathbf{y}\) (y_hat ou \(\hat{y}\)).

proc iml;
reset fuzz;
y = {14,16,15,18,19,17}; * vetor coluna y;
X = {1 1 0, 
     1 1 0,
     1 1 0, 
     1 0 1, 
     1 0 1, 
     1 0 1};
XLX = t(X)*X;
XLy = t(X)*y;
print XLX XLy;

rank_X = trace(X*ginv(X));
rank_XLX = trace(XLX*ginv(XLX));
print rank_X rank_XLX;

p = ncol(X);
k = round(trace(Ginv(XLX)*XLX));  * Cálculo do rank da matriz X;
defRank = p-k;
print 'Número de parâmetros =' p,,
      'Posto(X) =' k,,    
      'Deficiência de rank =' defRank;

Beta0 = Ginv(XLX)*t(X)*y; * Solução usando inversa-G;
y_hat0 = X*Beta0;
print 'Solução com inversa-G:       ' Beta0 y y_hat0;

TR = {0 1 0}; * Restrição do R: t1=0;
z = y//{0};
W = X//TR;
Beta1 = inv(t(W)*W)*t(W)*z;
y_hat1 = X*Beta1;
print 'Solução com restrição: t1=0 (R)        ' Beta1 y_hat1;

TS = {0 0 1}; * Restrição do SAS: t2=0;
z = y//{0};
W = X//TS;
Beta2 = inv(t(W)*W)*t(W)*z;
y_hat2 = X*Beta2;
print 'Solução com restrição: t2=0 (SAS)      ' Beta2 y_hat2;

T = {0 1 1}; * Restrição Estatística Experimental: t1+t2=0;
z = y//{0};
W = X//T;
Beta3 = inv(t(W)*W)*t(W)*z;
y_hat3 = X*Beta3;
print 'Solução com restrição: t1+t2=0 (EstExp)' Beta3 y_hat3;

11.3 Funções estimáveis

Uma função \(\boldsymbol{\lambda'\beta}\) é estimável no modelo \(\mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}\) se existe um vetor \(\mathbf{a}\) tal que \(E(\mathbf{a'y}) = \boldsymbol{\lambda'\beta}\).

A seguir, verificaremos se algumas combinações lineares são estimáveis de forma única.

\(\beta = \mu + \tau_1\):

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

[1] 1 1 0

L1Beta0 <- L1 %*% Beta0
L1Beta0

     [,1]
[1,]   15

L1Beta1 <- L1 %*% Beta1
L1Beta1

     [,1]
[1,]   15

L1Beta2 <- L1 %*% Beta2
L1Beta2

     [,1]
[1,]   15

L1Beta3 <- L1 %*% Beta3
L1Beta3

     [,1]
[1,]   15

Dado que o produto da combinação linear pelos vetores estimados de \(\boldsymbol{\beta}\) são iguais (invariantes), dizemos que \(\beta = \mu + \tau_1\) é estimável. Portanto, independentemente da estimativa do vetor de parâmetros \(\boldsymbol{\beta}\), as estimativas são iguais.

\(\beta = \tau_1 + \tau_2\):

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

[1] 0 1 1

L2Beta0 <- L2 %*% Beta0
L2Beta0

     [,1]
[1,]   11

L2Beta1 <- L2 %*% Beta1
L2Beta1

     [,1]
[1,]    3

L2Beta2 <- L2 %*% Beta2
L2Beta2

     [,1]
[1,]   -3

L2Beta3 <- L2 %*% Beta3
L2Beta3

     [,1]
[1,]    0

Dado que o produto da combinação linear pelos vetores estimados de \(\boldsymbol{\beta}\) são diferentes, dizemos que \(\beta = \tau_1 + \tau_2\) não é estimável. Portanto, funções não estimáveis não são invariantes e resultam em estimativas diferentes.

L1 = t({1,1,0}); * mi+t1 é estimável!;
L1Beta0 = L1*Beta0;
L1Beta1 = L1*Beta1;
L1Beta2 = L1*Beta2;
L1Beta3 = L1*Beta3;
print '-------------------------------',
      'Função estimável: LBeta = mi+t1',
      '-------------------------------',;
print '  Estimativa R:     ' L1Beta1,, 
      '  Estimativa SAS:   ' L1Beta2,,
      '  Estimativa EstExp:' L1Beta3;

L2 = {0 1 1}; * t1+t2 NÃO é estimável!;
L2Beta0 = L2*Beta0;
L2Beta1 = L2*Beta1;
L2Beta2 = L2*Beta2;
L2Beta3 = L2*Beta3;
print '-----------------------------------------',
      'Função NÃO estimável: L2Beta = T1 + T2   ',
      '-----------------------------------------',;
print '  Estimativa R:     ' L2Beta1,, 
      '  Estimativa SAS:   ' L2Beta2,, 
      '  Estimativa EstExp:' L2Beta3;
quit;