13  ANOVA balanceada com dois fatores

13.1 Introdução

Nesta seção, será apresentada a análise de variância de modelos balanceados com dois fatores de tratamento (two-way ANOVA). Seu modelo pode ser representado como:

\[ \begin{align} y_{ijk} &= \mu + \alpha_i + \beta_j + \gamma_{ij} + \epsilon_{ijk} \\ i = 1,2,\dots,a & \quad\quad j = 1,2,\dots,b \quad\quad k = 1,2,\dots,n \end{align} \]

• \(\alpha_i\) é o efeito do \(i\)-ésimo nível do fator \(A\);

• \(\beta_j\) é o efeito do \(j\)-ésimo nível do fator \(B\);

• \(\gamma_{ij}\) é o efeito da interação entre o \(i\)-ésimo nível do fator \(A\) e o \(j\)-ésimo nível do fator \(B\).

Como suposições, assumiremos:

• \(E(\epsilon_{ijk}) = 0\) para todo \(i,j,k\);

• \(var(\epsilon_{ijk}) = \sigma^2\) para todo \(i,j,k\);

• \(cov(\epsilon_{ijk}, \epsilon_{i'j'k'}) = 0\) para todo \((i,j,k) \ne (i',j',k')\);

• \(\epsilon_{ijk} \sim N(0, \sigma^2)\) para todo \(i,j,k\).

Dessa maneira, podemos reescrever o modelo balanceado na forma reparametrizada, conhecido como modelo de médias de caselas:

\[ \begin{align} y_{ijk} &= \mu_{ij} + \epsilon_{ijk} \\ i = 1,2,\dots,a \quad\quad j &= 1,2,\dots,b \quad\quad k = 1,2,\dots,n \end{align} \]

Em que \(\mu_{ij} = E(y_{ijk})\) é a média de uma observação \(k\) na casela \((ij)\).

Para ilustrar, considere o seguinte exemplo: O conteúdo da mistura de três tipos de queijo produzidos por dois métodos foi anotado por Marcuse (1949) (formato alterado). Duas peças de queijo foram medidas para cada tipo e cada método. Designando Método como o fator A e Tipo como o fator B, então a = 2, b = 3 e n = 2.

Método \ Tipo	1	2	3
1	39,02	35,74	37,02
1	38,79	35,41	36,00
2	38,96	35,58	35,70
2	39,01	35,52	36,04

• Fator A (Tratamento A): Método de produção de queijo, com 2 níveis (a = 2);

• Fator B (Tratamento B): Tipo de queijo, com 3 níveis (b = 3);

• Repetições: 2 repetições por tratamento (caso balanceado).

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

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

13.2 Estimação dos parâmetros

R

y <- c(39.02, 38.79, 35.74, 35.41, 37.02, 36.00, 38.96, 39.01, 35.58, 35.52, 35.70, 36.04)
y

 [1] 39.02 38.79 35.74 35.41 37.02 36.00 38.96 39.01 35.58 35.52 35.70 36.04

X <- matrix(c(
  1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0,
  1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0,
  1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0,
  1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0,
  1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0,
  1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0,
  1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0,
  1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0,
  1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0,
  1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0,
  1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1,
  1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1
), 
nrow = 12, byrow = TRUE)
X

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

O modelo matricial \(\mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}\) é dado por:

\[ \begin{cases} y_{111} = \mu + \alpha_1 + \beta_1 + \gamma_{11}\\ y_{112} = \mu + \alpha_1 + \beta_1 + \gamma_{11}\\ y_{121} = \mu + \alpha_1 + \beta_2 + \gamma_{12}\\ y_{122} = \mu + \alpha_1 + \beta_2 + \gamma_{12}\\ y_{131} = \mu + \alpha_1 + \beta_3 + \gamma_{21}\\ y_{132} = \mu + \alpha_1 + \beta_3 + \gamma_{21}\\ y_{211} = \mu + \alpha_2 + \beta_1 + \gamma_{22}\\ y_{212} = \mu + \alpha_2 + \beta_1 + \gamma_{22}\\ y_{221} = \mu + \alpha_2 + \beta_2 + \gamma_{31}\\ y_{222} = \mu + \alpha_2 + \beta_2 + \gamma_{31}\\ y_{231} = \mu + \alpha_2 + \beta_3 + \gamma_{32}\\ y_{232} = \mu + \alpha_2 + \beta_3 + \gamma_{32} \end{cases} \]

\[ \begin{bmatrix} y_{111} \\ y_{112} \\ y_{121} \\ y_{122} \\ y_{131} \\ y_{132} \\ y_{211} \\ y_{212} \\ y_{221} \\ y_{222} \\ y_{231} \\ y_{232} \end{bmatrix} = \begin{bmatrix} 39,02 \\ 38,79 \\ 35,74 \\ 35,41 \\ 37,02 \\ 36,00 \\ 38,96 \\ 39,01 \\ 35,58 \\ 35,52 \\ 35,70 \\ 36,04 \end{bmatrix} = \begin{bmatrix} 1 & 1 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 1 & 0 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 \\ 1 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} \mu \\ \alpha_1 \\ \alpha_2 \\ \beta_1 \\ \beta_2 \\ \beta_3 \\ \gamma_{11} \\ \gamma_{12} \\ \gamma_{21} \\ \gamma_{22} \\ \gamma_{31} \\ \gamma_{32} \end{bmatrix} + \begin{bmatrix} \epsilon_{111} \\ \epsilon_{112} \\ \epsilon_{121} \\ \epsilon_{122} \\ \epsilon_{131} \\ \epsilon_{132} \\ \epsilon_{211} \\ \epsilon_{212} \\ \epsilon_{221} \\ \epsilon_{222} \\ \epsilon_{231} \\ \epsilon_{232} \end{bmatrix} \]

Nota-se que o modelo apresenta 12 parâmetros e o posto(\(\mathbf{X}\)) = 6, ou seja, a matriz \(\mathbf{X}\) é de posto incompleto.

npar <- ncol(X)                      # número de parâmetros
npar

[1] 12

rank_X <- sum(diag(ginv(X) %*% X))   # posto de X
rank_X

[1] 6

deficit_rank <- npar - rank_X        # déficit de rank
deficit_rank

[1] 6

a <- 2    # níveis do fator A
b <- 3    # níveis do fator B
n <- 2    # número de repeticões

abn <- a * b * n    # número total de observacões
abn

[1] 12

Para visualizar as submatrizes da matriz \(\mathbf{X}\) respectivas à constante \(\mu\), ao parâmetro do fator A (\(\alpha_i\)), ao parâmetro do fator B (\(\beta_j\)) e ao de interação entre os fatores (\(\gamma_{ij}\)), procedemos da seguinte maneira:

X0 <- X[, 1]        # constante
X0

 [1] 1 1 1 1 1 1 1 1 1 1 1 1

XA <- X[, 2:3]      # fator A
XA

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

XB <- X[, 4:6]      # fator B
XB

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

XAB <- X[, 7:12]    # combinacão dos níveis dos dois fatores
XAB

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

A seguir, estimaremos o vetor de parâmetros \(\boldsymbol{\beta}\) usando duas inversas generalizadas diferentes.

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

• Inversa generalizada de Moore-Penrose:

Beta <- ginv(t(X) %*% X) %*% t(X) %*% y
Beta

           [,1]
 [1,] 18.449583
 [2,]  9.297917
 [3,]  9.151667
 [4,]  7.513750
 [5,]  5.258750
 [6,]  5.677083
 [7,]  3.643750
 [8,]  2.568750
 [9,]  3.085417
[10,]  3.870000
[11,]  2.690000
[12,]  2.591667

• Inversa generalizada simples (Searle):

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

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

iXLX <- (1 / 2) * kronecker(
  matrix(c(0,0,0,1), byrow = TRUE, ncol = 2), diag(6)
)
fractions(iXLX)

      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
 [1,]   0    0    0    0    0    0    0    0    0    0     0     0  
 [2,]   0    0    0    0    0    0    0    0    0    0     0     0  
 [3,]   0    0    0    0    0    0    0    0    0    0     0     0  
 [4,]   0    0    0    0    0    0    0    0    0    0     0     0  
 [5,]   0    0    0    0    0    0    0    0    0    0     0     0  
 [6,]   0    0    0    0    0    0    0    0    0    0     0     0  
 [7,]   0    0    0    0    0    0  1/2    0    0    0     0     0  
 [8,]   0    0    0    0    0    0    0  1/2    0    0     0     0  
 [9,]   0    0    0    0    0    0    0    0  1/2    0     0     0  
[10,]   0    0    0    0    0    0    0    0    0  1/2     0     0  
[11,]   0    0    0    0    0    0    0    0    0    0   1/2     0  
[12,]   0    0    0    0    0    0    0    0    0    0     0   1/2  

Betag <- iXLX %*% t(X) %*% y
Betag

        [,1]
 [1,]  0.000
 [2,]  0.000
 [3,]  0.000
 [4,]  0.000
 [5,]  0.000
 [6,]  0.000
 [7,] 38.905
 [8,] 35.575
 [9,] 36.510
[10,] 38.985
[11,] 35.550
[12,] 35.870

A tabela abaixo traz as duas estimativas de \(\boldsymbol{\beta}\) calculadas anteriormente.

Parâmetros	Beta MP	Beta Simples
$\hat{\mu}$	18.4496	0.000
$\hat{\alpha_1}$	9.2979	0.000
$\hat{\alpha_2}$	9.1517	0.000
$\hat{\beta_1}$	7.5137	0.000
$\hat{\beta_2}$	5.2587	0.000
$\hat{\beta_3}$	5.6771	0.000
$\hat{\gamma_{11}}$	3.6437	38.905
$\hat{\gamma_{12}}$	2.5687	35.575
$\hat{\gamma_{21}}$	3.0854	36.510
$\hat{\gamma_{22}}$	3.8700	38.985
$\hat{\gamma_{31}}$	2.6900	35.550
$\hat{\gamma_{32}}$	2.5917	35.870

SAS

proc iml;
reset fuzz;
y = {39.02,38.79,35.74,35.41,37.02,36.00,38.96,39.01,35.58,35.52,35.70,36.04};
X = {1 1 0 1 0 0 1 0 0 0 0 0,
     1 1 0 1 0 0 1 0 0 0 0 0,
     1 1 0 0 1 0 0 1 0 0 0 0,
     1 1 0 0 1 0 0 1 0 0 0 0,
     1 1 0 0 0 1 0 0 1 0 0 0,
     1 1 0 0 0 1 0 0 1 0 0 0,
     1 0 1 1 0 0 0 0 0 1 0 0,
     1 0 1 1 0 0 0 0 0 1 0 0,
     1 0 1 0 1 0 0 0 0 0 1 0,
     1 0 1 0 1 0 0 0 0 0 1 0,
     1 0 1 0 0 1 0 0 0 0 0 1,
     1 0 1 0 0 1 0 0 0 0 0 1};
print X[format=4.0] y[format=12.2];

rank_X = round(trace(ginv(X)*X));   * Calcula o posto da matriz X;
npar = ncol(X);
a=2;    * Número de níveis do fator A;
b=3;    * Número de níveis do fator B;
n=2;    * Número de repetições;
abn = a*b*n;    * Número total de observações;
X0 = X[,1];
XA = X[,2:3];
XB = X[,4:6];
XAB = X[,7:12];
print X[format=8.0],,XA[format=8.0],,XB[format=8.0],,XAB[format=8.0];
* ----------------------------------------------------------------------------;

* Estima Beta usando inversa de Moore-Penrose;
BetaMP = ginv(t(X)*X)*t(X)*y;       * Beta usando invG (Moore Penrose));
XLX = t(X)*X;
iXLX = (1/2)*{0 0, 0 1}@I(6);       
BetaG = iXLX*t(X)*y;                * Beta usando invG mais simples (Searle);

print BetaMP BetaG;

* ----------------------------------------------------------------------------;

13.3 Teste de hipótese

O primeiro teste que realizamos em ANOVA com dois fatores é o teste de interação entre os fatores A e B, cuja hipótese é:

\[ \begin{align} H_0 &: \text{Não há efeito da interação entre os fatores A e B} \\ H_1 &: \text{Há efeito da interação entre os fatores A e B} \end{align} \]

Caso a hipótese de interação entre os fatores seja não significativa (não se rejeita \(H_0\)), ou seja, os fatores são independentes, devemos analisar os efeitos principais dos fatores A e B, presentes no mesmo quadro de ANOVA construído para avaliar a interação. Assim, as hipóteses para os fatores A e B são dadas por:

\[ \begin{align} H_0 &: \text{Não há efeito do fator A} \\ H_1 &: \text{Há efeito do fator A para algum dos níveis de A} \end{align} \]

\[ \begin{align} H_0 &: \text{Não há efeito do fator B} \\ H_1 &: \text{Há efeito do fator B para algum dos níveis de B} \end{align} \]

Por outro lado, caso a hipótese de interação entre os fatores seja significativa (rejeita-se \(H_0: \text{Não há efeito da interação entre os fatores A e B}\)), devemos avaliar os efeitos simples de cada fator, ou seja, avaliar o efeito do fator A dentro de cada nível do fator B e/ou o efeito do fator B dentro de cada nível do fator A. Para isso, devemos realizar uma segunda ANOVA.

A seguir, demonstraremos um teste de hipótese em que a interação não é significativa, a partir das abordagens de modelo completo x modelo reduzido e hipótese linear geral.

13.3.1 Modelo Completo x Modelo Reduzido

R

Novamente, utilizaremos as seguintes formas quadráticas para calcular a soma de quadrados total (SQTotal) e a soma de quadrados dos resíduos (SQRes).

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

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

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

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

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

         [,1]
[1,] 26.97869

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

[1] 11

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

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

        [,1]
[1,] 0.66195

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

[1] 6

QMRes <- SQRes / gl_res
QMRes

         [,1]
[1,] 0.110325

Para calcular a soma de quadrados da interação entre os fatores A e B, utilizamos a seguinte forma quadrática:

\[ SQ_{A\times B} = \mathbf{y'} [\boldsymbol{X(X'X)^-X' - Z(X'X)^-Z'}] \mathbf{y} \]

onde \(\mathbf{Z}\) é a submatriz da matriz de delineamento \(\mathbf{X}\) respectiva à constante \(\mu\), ao parâmetro do fator A (\(\alpha_i\), \(i = 1,2\)) e ao parâmetro do fator B (\(\beta_j\), \(j = 1,2,3\)).

# Interação AxB
X1 <- cbind(X0, XA, XB)
X1

      X0          
 [1,]  1 1 0 1 0 0
 [2,]  1 1 0 1 0 0
 [3,]  1 1 0 0 1 0
 [4,]  1 1 0 0 1 0
 [5,]  1 1 0 0 0 1
 [6,]  1 1 0 0 0 1
 [7,]  1 0 1 1 0 0
 [8,]  1 0 1 1 0 0
 [9,]  1 0 1 0 1 0
[10,]  1 0 1 0 1 0
[11,]  1 0 1 0 0 1
[12,]  1 0 1 0 0 1

PAB <- X %*% ginv(t(X) %*% X) %*% t(X) - X1 %*% ginv(t(X1) %*% X1) %*% t(X1)

SQAB <- t(y) %*% PAB %*% y
SQAB

        [,1]
[1,] 0.30255

glAB <- round(sum(diag(ginv(PAB) %*% PAB)))
glAB

[1] 2

QMAB <- SQAB / glAB
QMAB

         [,1]
[1,] 0.151275

FAB <- QMAB / QMRes
FAB

         [,1]
[1,] 1.371176

ftabAB <- qf(0.95, glAB, gl_res)
ftabAB

[1] 5.143253

p_valorAB <- 1 - pf(FAB, glAB, gl_res)
p_valorAB

          [,1]
[1,] 0.3232726

A soma de quadrados do efeito principal do fator A é dado por:

\[ SQ_A = \mathbf{y'} \left[\boldsymbol{A(A'A)^-A'} - \frac{\mathbf{J}}n \right] \mathbf{y} \]

onde \(\mathbf{A}\) é a submatriz da matriz de delineamento \(\mathbf{X}\) respectiva ao parâmetro do fator A (\(\alpha_i\), \(i = 1,2\)).

# Fator A
PA <- XA %*% ginv(t(XA) %*% XA) %*% t(XA) - Jn/abn

SQA <- t(y) %*% PA %*% y
SQA

         [,1]
[1,] 0.114075

glA <- round(sum(diag(ginv(PA) %*% PA)))
glA

[1] 1

QMA <- SQA / glA
QMA

         [,1]
[1,] 0.114075

FA <- QMA / QMRes
FA

        [,1]
[1,] 1.03399

ftabA <- qf(0.95, glA, gl_res)
ftabA

[1] 5.987378

p_valorA <- 1 - pf(FA, glA, gl_res)
p_valorA

          [,1]
[1,] 0.3484596

De modo análogo, soma de quadrados do efeito principal do fator B é dado por:

\[ SQ_B = \mathbf{y'} \left[\boldsymbol{B(B'B)^-B'} - \frac{\mathbf{J}}n \right] \mathbf{y} \]

onde \(\mathbf{A}\) é a submatriz da matriz de delineamento \(\mathbf{X}\) respectiva ao parâmetro do fator A (\(\alpha_i\), \(i = 1,2\)).

# Fator B
PB <- XB %*% ginv(t(XB) %*% XB) %*% t(XB) - Jn/abn

SQB <- t(y) %*% PB %*% y
SQB

         [,1]
[1,] 25.90012

glB <- round(sum(diag(ginv(PB) %*% PB)))
glB

[1] 2

QMB <- SQB / glB
QMB

         [,1]
[1,] 12.95006

FB <- QMB / QMRes
FB

        [,1]
[1,] 117.381

ftabB <- qf(0.95, glB, gl_res)
ftabB

[1] 5.143253

p_valorB <- 1 - pf(FB, glB, gl_res)
p_valorB

              [,1]
[1,] 0.00001547711

Os resultados do teste de interação e de efeitos principais estão descritos no quadro de ANOVA a seguir.

ANOVA: Fatorial 2x3 (caso balanceado) - Interação e efeitos principais

FV	gl	SQ	QM	Fcal	Ftab	p-valor (5%)
Método (A)	1	0.114	0.114	1.034	5.987	0.3485
Tipo (B)	2	25.900	12.950	117.381	5.143	0.0000
Interação (AxB)	2	0.303	0.151	1.371	5.143	0.3233
Resíduo	6	0.662	0.110
Total	11	26.979

Dado que o \(F_{cal} < F_{tab}\) para a interação, considerando \(F_{(0,05;2;6)}\), a interação não foi significativa, logo não se rejeita \(H_0: \text{Não há efeito da interação entre os fatores A e B}\). Dessa forma, os fatores A e B podem ser considerados independentes, sem efeito de interação.

Visto que não há interação entre os fatores, podemos analisar os efeitos principais de cada um dos fatores. No caso do fator A (Método), também não se observou significância estatística para o efeito do método na fabricação de queijo. Dessa forma, não se rejeita \(H_0: \text{Não há efeito do fator A}\). Por outro lado, o fator B (Tipo) foi significativo, assim, rejeita-se \(H_0: \text{Não há efeito do fator B}\), permitindo dizer que entre, pelo menos, dois tipos de queijo há diferenças no processo de fabricação.

SAS

* Cálculo da soma de quadrados total - SQTotal;
P = I(abn) - J(abn,abn,1)/abn;
SQTotal = t(y)*(P)*y;
glTotal = round(trace(ginv(P)*P));

* Cálculo da soma de quadrados de resíduos - SQRes;
PR = I(abn) - X*ginv(t(X)*X)*t(X);
SQRes = t(y)*(PR)*y;
glRes = round(trace(ginv(PR)*PR));
QMRes = SQRes/glRes;

* Cálculo SQAxB - forma quadrática;
X1 = X0||XA||XB;
PAB = X*ginv(t(X)*X)*t(X) - X1*ginv(t(X1)*X1)*t(X1);
SQAB = t(y)*PAB*y;                  * Calcula SQ(AB);
glAB = round(trace(ginv(PAB)*PAB)); * Calcula gl da interação AxB;
QMAB = SQAB/glAB;
FAB = QMAB/QMRes;
p_valorAB = 1-cdf('F',FAB,glAB,glRes);

* Cálculo da soma de quadrados do fator A - SQ(A);
PA = XA*ginv(t(XA)*XA)*t(XA) - J(abn,abn,1)/abn;
SQA = t(y)*PA*y;
glA = round(trace(ginv(PA)*PA));
QMA = SQA/glA;
FA = QMA/QMRes;
p_valorA = 1-cdf('F',FA,glA,glRes);

* Cálculo da soma de quadrados do fator B - SQ(B);
PB = XB*ginv(t(XB)*XB)*t(XB) - J(abn,abn,1)/abn;
SQB = t(y)*PB*y;
glB = round(trace(ginv(PB)*PB));
QMB = SQB/glB;
FB = QMB/QMRes;
p_valorB = 1-cdf('F',FB,glB,glRes);

* Imprime o quadro de análise de variância – ANOVA (pág. 252);
print 'QUADRO DE ANOVA: Exemplo 14.4.2: Fatorial 2x3 (caso balanceado)';
print 'Método    ' glA[format=8.0]     SQA[format=12.4]   QMA[format=12.4]  FA[format=12.4]  p_valorA[format=12.4];
print 'Tipo      ' glB[format=8.0]     SQB[format=12.4]   QMB[format=12.4]  FB[format=12.4]  p_valorB[format=12.4];
print 'Interação ' glAB[format=8.0]    SQAB[format=12.4]  QMAB[format=12.4] FAB[format=12.4] p_valorAB[format=12.4];
print 'Resíduo   ' glRes[format=8.0]   SQRes[format=12.4] QMRes[format=12.4];
print 'Total     ' glTotal[format=8.0] SQTotal[format=12.4] ;

13.3.2 Hipótese Linear Geral

Também podemos calcular as somas de quadrados utilizando a abordagem da hipótese linear geral. A soma de quadrados obtida a partir desta abordagem é dada por:

\[ SQ_{Hip} = (\boldsymbol{C\hat{\beta}})' \left[\mathbf{C(X'X)^-C'}\right]^{-1} (\boldsymbol{C\hat{\beta}}) \]

onde \(\mathbf{C}\) é a matriz de coeficientes que reproduz, matricialmente, a hipótese de interesse.

R

CAxB <- matrix(
  c(0, 0, 0, 0, 0, 0, 1, -1, 0, -1, 1, 0,
    0, 0, 0, 0, 0, 0, 1, 0, -1, -1, 0, 1), nrow = 2, byrow = TRUE
)
CAxB

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,]    0    0    0    0    0    0    1   -1    0    -1     1     0
[2,]    0    0    0    0    0    0    1    0   -1    -1     0     1

SQ_CAxB <- t(CAxB %*% Beta) %*% solve(CAxB %*% ginv(t(X) %*% X) %*% t(CAxB)) %*% (CAxB %*% Beta)
SQ_CAxB

        [,1]
[1,] 0.30255

CA <- 
  (1/3) * matrix(
    c(0, 3, -3, 0, 0, 0, 1, 1, 1, -1, -1, -1), nrow = 1, byrow = TRUE
  )
fractions(CA)

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,]    0    1   -1    0    0    0  1/3  1/3  1/3 -1/3  -1/3  -1/3 

SQ_CA <- t(CA %*% Beta) %*% solve(CA %*% ginv(t(X) %*% X) %*% t(CA)) %*% (CA %*% Beta)
SQ_CA

         [,1]
[1,] 0.114075

CB <- 
  (1/2) * matrix(
    c(0, 0, 0, 2, -2, 0, 1, -1, 0, 1, -1, 0,
      0, 0, 0, 2, 0, -2, 1, 0, -1, 1, 0, -1), nrow = 2, byrow = TRUE
  )
fractions(CB)

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,]    0    0    0    1   -1    0  1/2 -1/2    0  1/2  -1/2     0 
[2,]    0    0    0    1    0   -1  1/2    0 -1/2  1/2     0  -1/2 

SQ_CB <- t(CB %*% Beta) %*% solve(CB %*% ginv(t(X) %*% X) %*% t(CB)) %*% (CB %*% Beta)
SQ_CB

         [,1]
[1,] 25.90012

Note que as somas de quadrados obtidas a partir da hipótese linear geral são iguais às obtidas pelo método do modelo completo x modelo reduzido.

Somas de quadrados entre as abordagens de Modelo completo x reduzido e Hipótese linear geral

Método	SQ A	SQ B	SQ AxB
Completo x Reduzido	0.114075	25.90012	0.30255
Hip.Lin.Geral	0.114075	25.90012	0.30255

SAS

* ----------------------------------------------;
* Cálculo das SQ usando Hipótese Linear Geral;
* ----------------------------------------------;

*            mi a1 a2 b1 b2 b3 g11 g12 g13 g21 g22 g23;      
CA   = (1/3)*{0  3 -3  0  0  0   1   1   1  -1  -1  -1};

CB   = (1/2)*{0  0  0  2 -2  0   1  -1   0   1  -1   0,
              0  0  0  2  0 -2   1   0  -1   1   0  -1};

CAxB =       {0  0  0  0  0  0   1  -1   0  -1   1   0,
              0  0  0  0  0  0   1   0  -1  -1   0   1}; 

SQ_CA = t(CA*Beta)*inv(CA*ginv(t(X)*X)*t(CA))*CA*Beta;
SQ_CB = t(CB*Beta)*inv(CB*ginv(t(X)*X)*t(CB))*CB*Beta;
SQ_CAxB = t(CAxB*Beta)*inv(CAxB*ginv(t(X)*X)*t(CAxB))*CAxB*Beta;

print 'Somas de quadrados usando Hipótese Linear Geral:',,
    SQ_CA[format=12.4] SQ_CB[format=12.4] SQ_CAxB[format=12.4];

13.4 Estimabilidade no modelo superparametrizado

13.4.1 Sem restrições nos parâmetros

Verificaremos se funções são estimáveis no modelo superparametrizado sem restrições nos parâmetros.

R

Neste primeiro caso, utilizaremos a abordagem de que se

\[ (\mathbf{X'X})(\mathbf{X'X})^-\boldsymbol{\lambda} = \boldsymbol{\lambda} \]

a função é estimável.

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

L1 <- matrix(c(0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0), ncol = 1)
L1

      [,1]
 [1,]    0
 [2,]    1
 [3,]   -1
 [4,]    0
 [5,]    0
 [6,]    0
 [7,]    0
 [8,]    0
 [9,]    0
[10,]    0
[11,]    0
[12,]    0

ver <- t(X) %*% X %*% ginv(t(X) %*% X)

verL1 <- ver %*% L1 |> round(2)
verL1

       [,1]
 [1,]  0.00
 [2,]  0.75
 [3,] -0.75
 [4,]  0.00
 [5,]  0.00
 [6,]  0.00
 [7,]  0.25
 [8,]  0.25
 [9,]  0.25
[10,] -0.25
[11,] -0.25
[12,] -0.25

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

        [,1]
[1,] 0.14625

Verifica-se que \(\beta = \alpha_1 - \alpha_2\) não é estimável no modelo sem restrição nos parâmetros.

• \(\beta = \alpha_1 - \alpha_2 + (\frac{1}3(\gamma_{11} + \gamma_{12} + \gamma_{13}) - (\frac{1}3)(\gamma_{21} = \gamma_{22} + \gamma_{23})\):

L2 <- (1/3) * t(matrix(c(0, 3, -3, 0, 0, 0, 1, 1, 1, -1, -1, -1), nrow = 1))
round(L2, 2)

       [,1]
 [1,]  0.00
 [2,]  1.00
 [3,] -1.00
 [4,]  0.00
 [5,]  0.00
 [6,]  0.00
 [7,]  0.33
 [8,]  0.33
 [9,]  0.33
[10,] -0.33
[11,] -0.33
[12,] -0.33

ver <- t(X) %*% X %*% ginv(t(X) %*% X)

verL2 <- ver %*% L2 |> round(2)
verL2

       [,1]
 [1,]  0.00
 [2,]  1.00
 [3,] -1.00
 [4,]  0.00
 [5,]  0.00
 [6,]  0.00
 [7,]  0.33
 [8,]  0.33
 [9,]  0.33
[10,] -0.33
[11,] -0.33
[12,] -0.33

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

      [,1]
[1,] 0.195

Verifica-se que \(\beta = \alpha_1 - \alpha_2 + (\frac{1}3(\gamma_{11} + \gamma_{12} + \gamma_{13}) - (\frac{1}3)(\gamma_{21} = \gamma_{22} + \gamma_{23})\) é estimavel no modelo superparametrizado sem restrição nos parâmetros.

SAS

L1 = t({0 1 -1 0 0 0 0 0 0 0 0 0});
ver = t(X)*X*ginv(t(X)*X);
verL1 = ver*L1;
L1BetaMP = t(L1)*BetaMP;
print 'Mostra que L1Beta = a1-a2 NÃO É estimável no modelo', 'SEM restrição nos parâmetros';
print L1 verL1 L1BetaMP[format=12.4];

print 'Mostra que L2Beta = a1-a2 + (1/3(g11+g12+g13)-(1/3)(g21+g22+g23) É estimável no modelo',
      'SEM restrição nos parâmetros';
L2 = (1/3)*t({0 3 -3 0 0 0 1 1 1 -1 -1 -1});
verL2 = ver*L2;
L2BetaMP = t(L2)*BetaMP;
print L2 verL2 L2BetaMP[format=12.4];

13.4.2 Com restrições nos parâmetros

Agora, verificaremos se funções são estimáveis no modelo superparametrizado com restrições nos parâmetros.

R

Primeiramente, assumiremos as seguintes condições marginais (restrições):

\[ \sum_i \hat{\alpha}_i = 0 \quad , \quad \sum_j \hat{\beta}_j = 0 \quad , \quad \sum_i \hat{\gamma}_{ij} = 0 \quad e \quad \sum_j \hat{\gamma}_{ij} = 0 \]

e expressas na matriz \(T\) da seguinte forma:

\[ \boldsymbol{T\hat{\beta}} = \begin{bmatrix} 0 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} \mu \\ \alpha_1 \\ \alpha_2 \\ \beta_1 \\ \beta_2 \\ \beta_3 \\ \gamma_{11} \\ \gamma_{12} \\ \gamma_{21} \\ \gamma_{22} \\ \gamma_{31} \\ \gamma_{32} \end{bmatrix} \]

t <- matrix(
  c(0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1,
    0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0,
    0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1), 
  nrow = 7, ncol = 12, byrow = TRUE)
t

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

rank_T <- sum(diag(ginv(t) %*% t))
rank_T

[1] 6

A matriz de coeficientes \(\mathbf{T}\) tem dimensão \(7 \times 12\) e posto\((\mathbf{W}) = 6\).

Juntando a matriz \(\mathbf{X}\) (de posto incompleto) com a matriz \(\mathbf{T}\), obtemos uma matriz \(\mathbf{W}\) de posto completo (posto\((\mathbf{W}) = 12\)).

W <- rbind(X, t)
W

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

rank_W <- sum(diag(ginv(W) %*% W))
rank_W

[1] 12

Ao vetor de valores observados, também colocamos restrições.

yr <- c(y, rep(0, 7))
yr

 [1] 39.02 38.79 35.74 35.41 37.02 36.00 38.96 39.01 35.58 35.52 35.70 36.04
[13]  0.00  0.00  0.00  0.00  0.00  0.00  0.00

Calculando o \(\boldsymbol{\beta}\) sujeito às condições marginais, obtemos:

Beta_R <- solve(t(W) %*% W) %*% t(W) %*% yr
Beta_R

            [,1]
 [1,] 36.8991667
 [2,]  0.0975000
 [3,] -0.0975000
 [4,]  2.0458333
 [5,] -1.3366667
 [6,] -0.7091667
 [7,] -0.1375000
 [8,] -0.0850000
 [9,]  0.2225000
[10,]  0.1375000
[11,]  0.0850000
[12,] -0.2225000

Com o \(\boldsymbol{\beta}\) com restrições nos parâmetros, verificaremos se as mesmas funções utilizadas no caso anterior (estimação de funções sem restrição) são estimáveis.

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

L1 <- t(matrix(c(0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0), nrow = 1))
L1

      [,1]
 [1,]    0
 [2,]    1
 [3,]   -1
 [4,]    0
 [5,]    0
 [6,]    0
 [7,]    0
 [8,]    0
 [9,]    0
[10,]    0
[11,]    0
[12,]    0

ver <- t(W) %*% W %*% solve(t(W) %*% W)

verL1 <- ver %*% L1
round(verL1)

      [,1]
 [1,]    0
 [2,]    1
 [3,]   -1
 [4,]    0
 [5,]    0
 [6,]    0
 [7,]    0
 [8,]    0
 [9,]    0
[10,]    0
[11,]    0
[12,]    0

L1Beta_r <- t(L1) %*% Beta_R
L1Beta_r

      [,1]
[1,] 0.195

No caso sem restrição nos parâmetros, \(\beta = \alpha_1 - \alpha_2\) é não estimável. Contudo, ao impor restrições nos parâmetros, a mesma função passa a ser estimável.

• \(\beta = \alpha_1 - \alpha_2 + (\frac{1}3(\gamma_{11} + \gamma_{12} + \gamma_{13}) - (\frac{1}3)(\gamma_{21} = \gamma_{22} + \gamma_{23})\):

L22 <- (1/3) * t(matrix(c(0, 3, -3, 0, 0, 0, 1, 1, 1, -1, -1, -1), nrow = 1))
L22 |> round(3)

        [,1]
 [1,]  0.000
 [2,]  1.000
 [3,] -1.000
 [4,]  0.000
 [5,]  0.000
 [6,]  0.000
 [7,]  0.333
 [8,]  0.333
 [9,]  0.333
[10,] -0.333
[11,] -0.333
[12,] -0.333

verL22 <- ver %*% L22 |> round(3)
verL22

        [,1]
 [1,]  0.000
 [2,]  1.000
 [3,] -1.000
 [4,]  0.000
 [5,]  0.000
 [6,]  0.000
 [7,]  0.333
 [8,]  0.333
 [9,]  0.333
[10,] -0.333
[11,] -0.333
[12,] -0.333

L2Beta_R <- t(L22) %*% Beta_R
L2Beta_R

      [,1]
[1,] 0.195

Bem como no caso sem restrição nos parâmetros, verifica-se que \(\beta = \alpha_1 - \alpha_2 + (\frac{1}3(\gamma_{11} + \gamma_{12} + \gamma_{13}) - (\frac{1}3)(\gamma_{21} = \gamma_{22} + \gamma_{23})\) continua a ser estimável no modelo com restrição nos parâmetros.

SAS

* Matriz T de condições marginais: T*Beta = 0;
T = {0 1 1 0 0 0 0 0 0 0 0 0,
     0 0 0 1 1 1 0 0 0 0 0 0,
     0 0 0 0 0 0 1 1 1 0 0 0,
     0 0 0 0 0 0 0 0 0 1 1 1,
     0 0 0 0 0 0 1 0 0 1 0 0,
     0 0 0 0 0 0 0 1 0 0 1 0,
     0 0 0 0 0 0 0 0 1 0 0 1};
rank_T = round(trace(ginv(T)*T));   * Determina o posto da matriz T;
W = X//T;                           * Junta as matrizes X e T;
rank_W = round(trace(ginv(W)*W));   * Calcula o posto da matriz W = X//T;
print npar rank_X rank_T rank_W;
yr = y//j(7,1,0);                   * Completa o vetor y com 7 zeros;
Beta_R = inv(t(W)*W)*t(W)*yr;       * Beta sujeito às condições marginais;
print BetaMP[format=12.4] Beta_R[format=12.4];

print '------------------------------------------------------------',
      ' ESTIMABILIDADE NO MODELO SUPERPARAMETRIZADO COM RESTRIÇÕES ',
      '------------------------------------------------------------',;
L1 = t({0 1 -1 0 0 0 0 0 0 0 0 0});
ver = t(W)*W*inv(t(W)*W);
verL1 = ver*L1;
L1Beta_r = t(L1)*Beta_R;
print 'Mostra que L1Beta = a1-a2 É estimável no modelo', 'COM restrição nos parâmetros';
print L1 verL1 L1Beta_R[format=12.4];

print 'Mostra que L2Beta = a1-a2 +(1/3(g11+g12+g13)-(1/3)(g21=g22+g23) É estimável no modelo',
      'COM restrição nos parâmetros';
L2 = (1/3)*t({0 3 -3 0  0  0  1  1  1 -1 -1 -1});
verL2 = ver*L2;
L2Beta_R = t(L2)*Beta_R;
print L2 verL2 L2Beta_R[format=12.4];   

13.4.3 Hipótese linear geral

Agora, utilizaremos a abordagem da hipótese linear geral para verificar a estimabilidade de funções no modelo superparametrizado sem restrição e com restrição.

R

CA <- (1 / 3) * c(0,  3, -3,  0,  0,  0,  1,  1,  1, -1, -1, -1)
fractions(CA)

 [1]    0    1   -1    0    0    0  1/3  1/3  1/3 -1/3 -1/3 -1/3

• Sem restrição:

CABeta <- CA %*% Beta
CABeta

      [,1]
[1,] 0.195

SQ_A <- t(CABeta) %*% solve(t(CA) %*% ginv(t(X) %*% X) %*% CA) %*% CABeta
SQ_A

         [,1]
[1,] 0.114075

• Com restrição:

CABeta_R <- CA %*% Beta_R
CABeta_R

      [,1]
[1,] 0.195

SQ_A_R = t(CABeta_R) %*% solve(t(CA) %*% solve(t(W) %*% W) %*% CA) %*% CABeta_R
SQ_A_R

         [,1]
[1,] 0.114075

Nota-se que, ao utilizar a matriz de coeficientes \(\mathbf{C}\) da hipótese linear geral, obtemos o mesmo resultado de soma de quadrados tanto para o caso sem restrição, como para o caso com restrição.

SAS

* Cálculo de SQA usando hipótese linear geral;

* (1) No modelo SEM restrição;
* CA = {0  1 -1  0  0  0  0  0  0  0  0  0};
CA = (1/3)*{0  3 -3  0  0  0  1  1  1 -1 -1 -1};
CABeta = CA*Beta_R;
print '-----------------------------------------------------------',
      ' ESTIMABILIDADE NO MODELO SUPERPARAMETRIZADO SEM RESTRIÇÃO ',
      '-----------------------------------------------------------',
      CA[format=6.2],,CABeta[format=12.4],,;

SQ_A = t(CABeta)*inv(CA*ginv(t(X)*X)*t(CA))*CABeta;

* (1) No modelo COM restrição;
CABeta_R = CA*Beta_R;
print 'No modelo COM restrição nos parâmetros:',,CA[format=6.2],,CABeta_R[format=12.4],,;

SQ_A_R = t(CABeta_r)*inv(CA*inv(t(W)*W)*t(CA))*CABeta_r;
print SQ_A[format=12.4] SQ_A_R[format=12.4];

quit;