15 ANOVA desbalanceada com dois fatores

Agora, trataremos da ANOVA com dois fatores em conjuntos de dados desbalanceados. O modelo desbalanceado com dois fatores é dado por:

\[ \begin{align} y_{ijk} = \mu + \alpha_i + \beta_j + \gamma{ij} + \epsilon_{ij} = \mu_{ij} + \epsilon_{ijk} \\ i = 1,\dots,a \space \text{ , } \space j = 1,\dots,b \space \text{ e } \space k = 1, \dots, n_{ij} \end{align} \]

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

Para o caso fatorial desbalanceado, utilizaremos o modelo de médias de caselas, que fornece uma abordagem mais simples e sem ambiguidades para testar hipóteses, quando comparado ao modelo superparametrizado.

Como exemplo, considere: Em um experimento de substituição do farelo de soja pelo farelo de girassol na ração de suínos, montou-se um experimento fatorial 2x5, com os fatores Sexo (1:Macho e 2:Fêmea) e Girassol (0, 25, 50, 75 e 100% de substituição). Foram utilizados 30 suínos (15 machos e 15 fêmeas) castrados da raça Duroc-Jersey, num delineamento inteiramente casualizado com 3 repetições. Na fase final do período experimental ocorreu a morte de três suínos. Os ganhos de peso dos animais aos 112 dias de experimento estão apresentados a seguir:

		Macho					Fêmea
0	25	50	75	100	0	25	50	75	100
-	94,5	99,5	93,0	83,0	77,9	71,5	67,5	71,5	89,5
86,0	96,0	98,0	96,0	80,0	83,2	73,5	-	70,8	91,8
84,0	95,8	-	90,5	78,5	83,5	70,5	65,0	72,5	92,9

O modelo matemático é dado por:

\[ \begin{align} y_{ijk}& = \mu_{ij} + \epsilon_{ijk} \\ i = 1,2 \space \text{ , } \space &j = 1,\dots,5 \space \text{ e } \space k = n_{ij} \end{align} \]

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

\[ \begin{bmatrix} y_{111} \\ y_{112} \\ y_{113} \\ y_{121} \\ y_{122} \\ y_{123} \\ y_{131} \\ y_{132} \\ y_{141} \\ y_{142} \\ y_{143} \\ y_{151} \\ y_{152} \\ y_{153} \\ y_{211} \\ y_{212} \\ y_{221} \\ y_{222} \\ y_{223} \\ y_{231} \\ y_{232} \\ y_{241} \\ y_{242} \\ y_{243} \\ y_{251} \\ y_{252} \\ y_{253} \end{bmatrix} = \begin{bmatrix} 77.9 \\ 83.2 \\ 83.5 \\ 71.5 \\ 73.5 \\ 70.5 \\ 67.5 \\ 65 \\ 71.5 \\ 70.8 \\ 72.5 \\ 89.5 \\ 91.8 \\ 92.9 \\ 86 \\ 84 \\ 94.5 \\ 96 \\ 95.8 \\ 99.5 \\ 98 \\ 93 \\ 96 \\ 90.5 \\ 83 \\ 80 \\ 78.5 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} \mu_{11} \\ \mu_{12} \\ \mu_{13} \\ \mu_{14} \\ \mu_{15} \\ \mu_{21} \\ \mu_{22} \\ \mu_{23} \\ \mu_{24} \\ \mu_{25} \end{bmatrix} \begin{bmatrix} \epsilon_{111} \\ \epsilon_{112} \\ \epsilon_{113} \\ \epsilon_{121} \\ \epsilon_{122} \\ \epsilon_{123} \\ \epsilon_{131} \\ \epsilon_{132} \\ \epsilon_{141} \\ \epsilon_{142} \\ \epsilon_{143} \\ \epsilon_{151} \\ \epsilon_{152} \\ \epsilon_{153} \\ \epsilon_{211} \\ \epsilon_{212} \\ \epsilon_{221} \\ \epsilon_{222} \\ \epsilon_{223} \\ \epsilon_{231} \\ \epsilon_{232} \\ \epsilon_{241} \\ \epsilon_{242} \\ \epsilon_{243} \\ \epsilon_{251} \\ \epsilon_{252} \\ \epsilon_{253} \end{bmatrix} \]

library(MASS)

y <- c(77.9,83.2,83.5,71.5,73.5,70.5,67.5,65,71.5,70.8,72.5,89.5,91.8,92.9,86,84,94.5,96,95.8,99.5,98,93,96,90.5,83,80,78.5)

n <- length(y)
n

[1] 27

W <- matrix(c(
  1, 0,  0,  0,  0,  0,  0,  0,  0,  0,
  1, 0,  0,  0,  0,  0,  0,  0,  0,  0,
  1, 0,  0,  0,  0,  0,  0,  0,  0,  0,
  0, 1,  0,  0,  0,  0,  0,  0,  0,  0,
  0, 1,  0,  0,  0,  0,  0,  0,  0,  0,
  0, 1,  0,  0,  0,  0,  0,  0,  0,  0,
  0, 0,  1,  0,  0,  0,  0,  0,  0,  0,
  0, 0,  1,  0,  0,  0,  0,  0,  0,  0,
  0, 0,  0,  1,  0,  0,  0,  0,  0,  0,
  0, 0,  0,  1,  0,  0,  0,  0,  0,  0,
  0, 0,  0,  1,  0,  0,  0,  0,  0,  0,
  0, 0,  0,  0,  1,  0,  0,  0,  0,  0,
  0, 0,  0,  0,  1,  0,  0,  0,  0,  0,
  0, 0,  0,  0,  1,  0,  0,  0,  0,  0,
  0, 0,  0,  0,  0,  1,  0,  0,  0,  0,
  0, 0,  0,  0,  0,  1,  0,  0,  0,  0,
  0, 0,  0,  0,  0,  0,  1,  0,  0,  0,
  0, 0,  0,  0,  0,  0,  1,  0,  0,  0,
  0, 0,  0,  0,  0,  0,  1,  0,  0,  0,
  0, 0,  0,  0,  0,  0,  0,  1,  0,  0,
  0, 0,  0,  0,  0,  0,  0,  1,  0,  0,
  0, 0,  0,  0,  0,  0,  0,  0,  1,  0,
  0, 0,  0,  0,  0,  0,  0,  0,  1,  0,
  0, 0,  0,  0,  0,  0,  0,  0,  1,  0,
  0, 0,  0,  0,  0,  0,  0,  0,  0,  1,
  0, 0,  0,  0,  0,  0,  0,  0,  0,  1,
  0, 0,  0,  0,  0,  0,  0,  0,  0,  1
), ncol = 10, byrow = TRUE
)

Nota-se que posto\((W) = k = 10\), sendo \(k\) o número de parâmetros do modelo.

# Número de parâmetros
k <- ncol(W)
k

[1] 10

# Posto de W
rank_W <- sum(diag(ginv(W) %*% W))
rank_W

[1] 10

# Déficit de rank
deficit_rank <- k - rank_W
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}_{11.} \\ \bar{y}_{12.} \\ \bar{y}_{13.} \\ \bar{y}_{14.} \\ \bar{y}_{15.} \\ \bar{y}_{21.} \\ \bar{y}_{22.} \\ \bar{y}_{23.} \\ \bar{y}_{24.} \\ \bar{y}_{25.} \end{bmatrix} \]

Onde o vetor \(\mathbf{\bar{y}}\) contém as médias amostrais das caselas.

Mi <- solve(t(W) %*% W) %*% t(W) %*% y
Mi

          [,1]
 [1,] 81.53333
 [2,] 71.83333
 [3,] 66.25000
 [4,] 71.60000
 [5,] 91.40000
 [6,] 85.00000
 [7,] 95.43333
 [8,] 98.75000
 [9,] 93.16667
[10,] 80.50000

data Girassol;
input Sexo Girassol Trat Rep GP;
* n11=2 n12=3 n13=2 n14=3 n15=3;
* n21=2 n22=3 n23=2 n24=3 n25=3;
cards;
  1   0   1   1  77.9
  1   0   1   2  83.2
  1   0   1   3  83.5
  1  25   2   1  71.5
  1  25   2   2  73.5
  1  25   2   3  70.5
  1  50   3   1  67.5
  1  50   3   3  65.0
  1  75   4   1  71.5
  1  75   4   2  70.8
  1  75   4   3  72.5
  1 100   5   1  89.5
  1 100   5   2  91.8
  1 100   5   3  92.9
  2   0   6   2  86.0
  2   0   6   3  84.0
  2  25   7   1  94.5
  2  25   7   2  96.0
  2  25   7   3  95.8
  2  50   8   1  99.5
  2  50   8   2  98.0
  2  75   9   1  93.0
  2  75   9   2  96.0
  2  75   9   3  90.5
  2 100  10   1  83.0
  2 100  10   2  80.0
  2 100  10   3  78.5
;

proc iml;
varNames = {"Sexo" "Girassol" "Trat" "Rep" "GP"};
use work.Girassol;
read all var varNames;
close work.Girassol;
print Sexo Girassol Trat Rep GP;

y = GP; 
n = nrow(y);

W = design(Trat);

Mi = inv(t(W)*W)*t(W)*GP;
print GP W Mi[format=8.2];

15.1 Teste de hipótese

15.1.1 Modelo de médias

A seguir, testaremos a hipótese de igualdade entre as médias dos tratamentos utilizando a abordagem do modelo completo x modelo reduzido.

Jnxn <- matrix(1, n, n)
Inxn <- diag(n)

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

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

# Total
AT <- Inxn - (1 / n) * Jnxn

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

         [,1]
[1,] 2896.396

gl_total <- sum(diag(AT %*% ginv(AT)))
gl_total

[1] 26

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

\[ SQRes = \mathbf{y'} [\mathbf{I} - \mathbf{W(W'W)^{-}W'}] \mathbf{y} \]

# Residuo
AR <- Inxn - W %*% ginv(t(W) %*% W) %*% t(W)

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

         [,1]
[1,] 65.23667

gl_res <- sum(diag(AR %*% ginv(AR)))
gl_res

[1] 17

QMRes <- SQRes / gl_res
QMRes

         [,1]
[1,] 3.837451

A soma de quadrados de tratamento é:

\[ SQTrat = \mathbf{y'} \left[\mathbf{W(W'W)^{-}W'} - \frac{1}n \mathbf{J}\right] \mathbf{y} \]

# Tratamento

AG <- W %*% ginv(t(W) %*% W) %*% t(W) - (1/n) * Jnxn

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

        [,1]
[1,] 2831.16

gl_trat <- sum(diag(AG %*% ginv(AG)))
gl_trat

[1] 9

QMTrat <- SQTrat / gl_trat
QMTrat

         [,1]
[1,] 314.5733

# F e p-valor

Fcalc <- QMTrat / QMRes
Fcalc

         [,1]
[1,] 81.97454

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

[1] 2.494291

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

                     [,1]
[1,] 0.000000000003195444

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

ANOVA fatorial desbalanceada - modelo de médias: H0:
FV	gl	SQ	QM	Fcal	Ftab	p-valor (5%)
Tratamentos	9	2831.1596	314.5733	81.9745	2.4943	0
Resíduo	17	65.2367	3.8375
Total	26	2896.3963

Dado que o \(F_{cal} > F_{tab}\), considerando \(F_{(0,05;9;17)}\), rejeita-se a hipótese nula de igualdade de médias, ou seja, pelo menos duas médias diferem entre si.

Com isso, precisamos analisar o efeito de interação entre os tratamentos.

AT = I(n)-(1/n)*J(n,n,1);
SQTotal = t(Y)*AT*y;
gl_total = n-1;

ARes = I(n)- W*inv(t(W)*W)*t(W);
SQRes = t(y)*ARes*y;
gl_res = round(trace(ARes*ginv(ARes)));
QMRes = SQRes/gl_Res;

ATrat = W*inv(t(W)*W)*t(W) - (1/n)*J(n,n,1);
SQTrat = t(y)*ATrat*y;
gl_trat = round(trace(ATrat*ginv(ATrat)));
QMTrat = SQTrat/gl_trat;
F_trat = QMTrat/QMRes;
p_trat = 1-cdf('F',F_trat,gl_trat,gl_res);

print 'ANOVA - modelo de médias',,
'Tratamentos ' gl_trat SQTrat[format=12.4] QMTrat[format=12.4] F_trat[format=12.4] p_trat[format=12.4],,
'Resíduo     ' gl_res SQRes[format=12.4] QMRes[format=12.4],,
'Total       ' gl_total SQTotal[format=12.4];

15.1.2 Efeito de interação

Para avaliar se há interação entre os fatores, realizaremos outra ANOVA, avaliando os efeitos principais de ambos os fatores (vide Seção 13.3).

\[ \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} \]

Para isso, utilizaremos a 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}}) \]

A matriz de coeficientes \(\mathbf{C}\) do fator Sexo (CS), do fator Girassol (CG) e da interação entre os fatores (CSxG) são dadas a seguir:

CS <- matrix(c(1,  1,  1,  1,  1, -1, -1, -1, -1, -1), ncol = 10, byrow = TRUE)
CS

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

CG <- matrix(c(
  -2, -1,  0,  1,  2, -2, -1,  0,  1,  2,
   2, -1, -2, -1,  2,  2, -1, -2, -1,  2,
  -1,  2,  0, -2,  1, -1,  2,  0, -2,  1,
   1, -4,  6, -4,  1,  1, -4,  6, -4,  1
), ncol = 10, byrow = TRUE
)
CG

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

CSxG <- rbind(
  sweep(CS, 2, CG[1,], `*`),
  sweep(CS, 2, CG[2,], `*`),
  sweep(CS, 2, CG[3,], `*`),
  sweep(CS, 2, CG[4,], `*`)
)
CSxG

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

A matriz CSxG é construída a partir da mutiplicação, elemento a elemento, da matriz CS com a primeira, segunda, terceira e quarta linhas da matriz CG. Para isso, utilizamos a função sweep() para realizar esta operação e, posteriormente, com a rbind(), juntamos os resultados por linha.

# Fator Sexo
SQSexo <- t(CS %*% Mi) %*% solve(CS %*% ginv(t(W) %*% W) %*% t(CS)) %*% (CS %*% Mi)
SQSexo

         [,1]
[1,] 1286.797

gl_sexo <- nrow(CS)
gl_sexo

[1] 1

QMSexo <- SQSexo / gl_sexo
QMSexo

         [,1]
[1,] 1286.797

F_sexo <- QMSexo / QMRes
F_sexo

         [,1]
[1,] 335.3259

p_sexo <- 1 - pf(F_sexo, gl_sexo, gl_res)
p_sexo

                     [,1]
[1,] 0.000000000001258105

# Fator Girassol
SQGirassol <- t(CG %*% Mi) %*% solve(CG %*% ginv(t(W) %*% W) %*% t(CG)) %*% (CG %*% Mi)
SQGirassol

         [,1]
[1,] 47.35935

gl_girassol <- nrow(CG)
gl_girassol

[1] 4

QMGirassol <- SQGirassol / gl_girassol
QMGirassol

         [,1]
[1,] 11.83984

F_girassol <- QMGirassol / QMRes
F_girassol

         [,1]
[1,] 3.085339

p_girassol <- 1 - pf(F_girassol, gl_girassol, gl_res)
p_girassol

          [,1]
[1,] 0.0442227

# Interação Sexo x Girassol
SQSxG <- t(CSxG %*% Mi) %*% solve(CSxG %*% ginv(t(W) %*% W) %*% t(CSxG)) %*% (CSxG %*% Mi)
SQSxG

        [,1]
[1,] 1624.61

gl_SxG <- nrow(CSxG)
gl_SxG

[1] 4

QMSxG <- SQSxG / gl_SxG
QMSxG

         [,1]
[1,] 406.1526

F_SxG <- QMSxG / QMRes
F_SxG

         [,1]
[1,] 105.8392

p_SxG <- 1 - pf(F_SxG, gl_SxG, gl_res)
p_SxG

                     [,1]
[1,] 0.000000000008890666

No caso fatorial desbalanceado, a soma das somas de quadrados não será igual a soma de quadrados dos tratamentos, pois não utilizamos contrastes ortogonais ponderados.

SSQ <- SQSexo + SQGirassol + SQSxG
SSQ

         [,1]
[1,] 2958.767

SQTrat

        [,1]
[1,] 2831.16

O quadro de ANOVA resume os resultados obtidos.

ANOVA: Fatorial Sexo x Girassol (caso desbalanceado) - Interação e efeitos principais
FV	gl	SQ	QM	Fcal	p-valor (5%)
Sexo (A)	1	1286.797	1286.797	335.326	0.0000
Girassol (B)	4	47.359	11.840	3.085	0.0442
Interação (AxB)	4	1624.610	406.153	105.839	0.0000
Resíduo	17	65.237	3.837
Total	26	2896.396

Dado que o p-valor da interação é menor que p-valor = 0,05, rejeita-se \(H_0: \text{Não há efeito da interação entre os fatores A e B}\), ou seja, há efeito de interação entre os fatores Sexo e Girassol (há dependência).

Visto que há interação entre os fatores, precisamos avaliar os efeitos simples dos tratamentos.

CS = { 1  1  1  1  1 -1 -1 -1 -1 -1};
*CG = {-4  1  1  1  1 -4  1  1  1  1,
       0 -3  1  1  1  0 -3  1  1  1,
       0  0 -2  1  1  0  0 -2  1  1,
       0  0  0 -1  1  0  0  0 -1  1}; 

CG = {-2 -1  0  1  2 -2 -1  0  1  2,
       2 -1 -2 -1  2  2 -1 -2 -1  2,
      -1  2  0 -2  1 -1  2  0 -2  1,
       1 -4  6 -4  1  1 -4  6 -4  1};
CSxG = CS#CG[1,]//CS#CG[2,]//CS#CG[3,]//CS#CG[4,];
print CS,, CG,, CSxG;

SQSexo = t(CS*Mi)*inv(CS*inv(t(W)*W)*t(CS))*(CS*Mi);
gl_Sexo = nrow(CS);
QMSexo = SQSexo/gl_Sexo;
F_sexo =QMSexo/QMRes;
p_sexo = 1-cdf('F',F_sexo,gl_sexo,gl_res);

SQGirassol = t(CG*Mi)*inv(CG*inv(t(W)*W)*t(CG))*(CG*Mi);
gl_Girassol = nrow(CG);
QMGirassol = SQGirassol/gl_Girassol;
F_Girassol = QMGirassol/QMRes;
p_girassol = 1-cdf('F',F_Girassol,gl_Girassol,gl_res);

SQSxG = t(CSxG*Mi)*inv(CSxG*inv(t(W)*W)*t(CSxG))*(CSxG*Mi);
gl_SxG = nrow(CSxG);
QMSxG = SQSxG/gl_SxG;
F_SxG = QMSxG/QMRes;
p_SxG = 1-cdf('F',F_SxG,gl_SxG,gl_res);

SQTrats = SQSexo + SQGirassol + SQSxG;

print 'ANOVA - FATORIAL SEXO x GIRASSOL',,
'Sexo            ' gl_sexo SQSexo[format=12.4] QMSexo[format=12.4] F_sexo[format=12.4] p_sexo[format=12.4],,
'Girassol      ' gl_girassol SQGirassol[format=12.4] QMGirassol[format=12.4] F_Girassol[format=12.4] p_Girassol[format=12.4],,
'Interação SxG   ' gl_SxG SQSxG[format=12.4] QMSxG[format=12.4] F_SxG[format=12.4] p_SxG[format=12.4],,
'Resíduo         ' gl_res SQRes[format=12.4] QMRes[format=12.4],,
'Total           ' gl_total SQTotal[format=12.4],,;

print 'SQTrats = SQS + SQG + SQSxG = ' SQTrats[format=12.4] SQTrat[format=12.4];

15.1.3 Efeitos Simples

Dado que há interação entre os fatores, devemos avaliar os efeitos simples de cada fator, ou seja, avaliar o efeito do fator A (Sexo) dentro de cada nível do fator B (Girassol) e/ou o efeito do fator B (Girassol) dentro de cada nível do fator A (Sexo).

Para isso, utilizaremos contrastes ortogonais não ponderados.

Sexo dentro de Girassol

Primeiramente, compararemos as médias dos dois sexos dentro de cada nível do fator Girassol.

\[ H_0: \mu_{1j} = \mu_{2j}, \quad \text{para } j = 1,2,3,4,5 \]

a1 <- matrix(c(1, 0, 0, 0, 0, -1,  0,  0,  0,  0), ncol = 1)
a2 <- matrix(c(0, 1, 0, 0, 0,  0, -1,  0,  0,  0), ncol = 1)
a3 <- matrix(c(0, 0, 1, 0, 0,  0,  0, -1,  0,  0), ncol = 1)
a4 <- matrix(c(0, 0, 0, 1, 0,  0,  0,  0, -1,  0), ncol = 1)
a5 <- matrix(c(0, 0, 0, 0, 1,  0,  0,  0,  0, -1), ncol = 1)

A soma de quadrados dos contrastes é calculada da seguinte forma:

\[ SQ_c = (\boldsymbol{\lambda' \mu})' \left[\boldsymbol{\lambda'} (\mathbf{W'W})^{-1} \boldsymbol{\lambda}\right]^{-1} (\boldsymbol{\lambda' \mu}) \]

Cada contraste tem um graus de liberdade.

SQa1 <- t(t(a1) %*% Mi) %*% solve(t(a1) %*% solve(t(W) %*% W) %*% a1) %*% (t(a1) %*% Mi)
SQa1

         [,1]
[1,] 14.42133

SQa2 <- t(t(a2) %*% Mi) %*% solve(t(a2) %*% solve(t(W) %*% W) %*% a2) %*% (t(a2) %*% Mi)
SQa2

       [,1]
[1,] 835.44

SQa3 <- t(t(a3) %*% Mi) %*% solve(t(a3) %*% solve(t(W) %*% W) %*% a3) %*% (t(a3) %*% Mi)
SQa3

        [,1]
[1,] 1056.25

SQa4 <- t(t(a4) %*% Mi) %*% solve(t(a4) %*% solve(t(W) %*% W) %*% a4) %*% (t(a4) %*% Mi)
SQa4

         [,1]
[1,] 697.6817

SQa5 <- t(t(a5) %*% Mi) %*% solve(t(a5) %*% solve(t(W) %*% W) %*% a5) %*% (t(a5) %*% Mi)
SQa5

        [,1]
[1,] 178.215

gl_a <- 1

A estatística F e o p-valor de cada contraste é dado por:

Fcalc1 <- SQa1 / QMRes
p_valor1 <- 1 - pf(Fcalc1, 1, gl_res)
p_valor1

           [,1]
[1,] 0.06934062

Fcalc2 <- SQa2 / QMRes
p_valor2 <- 1 - pf(Fcalc2, 1, gl_res)
p_valor2

                    [,1]
[1,] 0.00000000004020073

Fcalc3 <- SQa3 / QMRes
p_valor3 <- 1 - pf(Fcalc3, 1, gl_res)
p_valor3

                     [,1]
[1,] 0.000000000006192602

Fcalc4 <- SQa4 / QMRes
p_valor4 <- 1 - pf(Fcalc4, 1, gl_res)
p_valor4

                   [,1]
[1,] 0.0000000001658613

Fcalc5 <- SQa5 / QMRes
p_valor5 <- 1 - pf(Fcalc5, 1, gl_res)
p_valor5

               [,1]
[1,] 0.000003010948

O quadro de ANOVA resume os resultados obtidos.

ANOVA: Fatorial Sexo x Girassol (caso desbalanceado) - Sexo dentro de Girassol
Hipóteses	gl	SQ	Fcal	p-valor (5%)
Girassol = 0: M = F	1	14.421	3.758	0.0693
Girassol = 25: M = F	1	835.440	217.707	0.0000
Girassol = 50: M = F	1	1056.250	275.248	0.0000
Girassol = 75: M = F	1	697.682	181.809	0.0000
Girassol = 100: M = F	1	178.215	46.441	0.0000

Apenas para Girassol = 0% \(H_0: \mu_{1j} = \mu_{2j}\) não foi rejeitada. No demais níveis (25%, 50%, 75% e 100%), rejeita-se a hipótese nula, ou seja, há diferença entre o ganho de peso entre machos e fêmeas.

Para verificar a diferença do percentual de farelo de girassol entre os sexo, devemos analisar as médias dos efeitos simples.

Sexo	0	25	50	75	100
F	81.53	71.83	66.25	71.60	91.4
M	85.00	95.43	98.75	93.17	80.5

Uma vez que para girassol = 0% não houve significância, as médias entre os sexos não diferem, estatisticamente. Para os níveis 25%, 50%, 75% e 100%, o ganho médio de peso dos machos foi superior ao das fêmeas, mas com 100% de Girassol, as fêmeas tiveram maior ganho médio de peso.

Sexo	0	25	50	75	100
F	81.53 a	71.83 b	66.25 b	71.60 b	91.4 a
M	85.00 a	95.43 a	98.75 a	93.17 a	80.5 b

a1 = {1 0 0 0 0 -1  0  0  0  0};
a2 = {0 1 0 0 0  0 -1  0  0  0};
a3 = {0 0 1 0 0  0  0 -1  0  0};
a4 = {0 0 0 1 0  0  0  0 -1  0};
a5 = {0 0 0 0 1  0  0  0  0 -1};
SQa1 = t(a1*Mi)*inv(a1*inv(t(W)*W)*t(a1))*(a1*Mi);
SQa2 = t(a2*Mi)*inv(a2*inv(t(W)*W)*t(a2))*(a2*Mi);
SQa3 = t(a3*Mi)*inv(a3*inv(t(W)*W)*t(a3))*(a3*Mi);
SQa4 = t(a4*Mi)*inv(a4*inv(t(W)*W)*t(a4))*(a4*Mi);
SQa5 = t(a5*Mi)*inv(a5*inv(t(W)*W)*t(a5))*(a5*Mi);
gl_a = 1;
Fa1 = SQa1/QMRes; Fa2 = SQa2/QMRes; Fa3 = SQa3/QMRes; Fa4 = SQa4/QMRes; Fa5 = SQa5/QMRes;
p_a1 = 1-cdf('F',Fa1,1,gl_res); p_a2 = 1-cdf('F',Fa2,1,gl_res);
p_a3 = 1-cdf('F',Fa3,1,gl_res); p_a4 = 1-cdf('F',Fa4,1,gl_res);
p_a5 = 1-cdf('F',Fa5,1,gl_res);

print '------------------------------------------------------------------------',
      ' Desdobramento (1): Compara as médias de Sexo em cada nível de Girassol ',
      '------------------------------------------------------------------------',,
'Girassol=  0: M=F' gl_a SQa1[format=12.4] Fa1[format=12.4] p_a1[format=12.4],,
'Girassol= 25: M=F' gl_a SQa2[format=12.4] Fa2[format=12.4] p_a2[format=12.4],,
'Girassol= 50: M=F' gl_a SQa3[format=12.4] Fa3[format=12.4] p_a3[format=12.4],,
'Girassol= 75: M=F' gl_a SQa4[format=12.4] Fa4[format=12.4] p_a4[format=12.4],,
'Girassol=100: M=F' gl_a SQa5[format=12.4] Fa5[format=12.4] p_a5[format=12.4];

quit;

Girassol dentro de Sexo

Agora, compararemos as médias do fator Girassol dentro de cada nível de sexo. Como Girassol é um fator quantitativo e seus níveis são igualmente espaçados, vamos usar coeficientes de polinômios ortogonais para realizar os testes de tendência, separadamente, para cada Sexo.

\[ H_0: \mu_{i1} = \mu_{i2} = \mu_{i3} = \mu_{i4} = \mu_{i5}, \quad \text{para } i = 1,2 \]

Como são 5 níveis, utilizaremos os coeficientes até o 4º grau.

Para as fêmeas vamos usar:

F1 <- c(-2, -1,  0,  1, 2, 0, 0, 0, 0, 0)  # Grau 1
F2 <- c( 2, -1, -2, -1, 2, 0, 0, 0, 0, 0)  # Grau 2
F3 <- c(-1,  2,  0, -2, 1, 0, 0, 0, 0, 0)  # Grau 3
F4 <- c( 1, -4,  6, -4, 1, 0, 0, 0, 0, 0)  # Grau 4

e para os machos:

M1 <- c(0, 0, 0, 0, 0, -2, -1,  0,  1, 2)  # Grau 1
M2 <- c(0, 0, 0, 0, 0,  2, -1, -2, -1, 2)  # Grau 2
M3 <- c(0, 0, 0, 0, 0, -1,  2,  0, -2, 1)  # Grau 3
M4 <- c(0, 0, 0, 0, 0,  1, -4,  6, -4, 1)  # Grau 4

Com os contrastes ortogonais para cada sexo, realizaremos a ANOVA para cada um deles, a fim de verificar qual o grau de polinômio que melhor explica o ganho de peso de cada um dos sexos.

Fêmeas:

SQf1 <- t(t(F1) %*% Mi) %*% solve(t(F1) %*% solve(t(W) %*% W) %*% F1) %*% (t(F1) %*% Mi)

SQf2 <- t(t(F2) %*% Mi) %*% solve(t(F2) %*% solve(t(W) %*% W) %*% F2) %*% (t(F2) %*% Mi)

SQf3 <- t(t(F3) %*% Mi) %*% solve(t(F3) %*% solve(t(W) %*% W) %*% F3) %*% (t(F3) %*% Mi)

SQf4 <- t(t(F4) %*% Mi) %*% solve(t(F4) %*% solve(t(W) %*% W) %*% F4) %*% (t(F4) %*% Mi)


gl_f <- 1


Fcalc1 <- SQf1 / QMRes
p_valor1 <- 1 - pf(Fcalc1, 1, gl_res)

Fcalc2 <- SQf2 / QMRes
p_valor2 <- 1 - pf(Fcalc2, 1, gl_res)

Fcalc3 <- SQf3 / QMRes
p_valor3 <- 1 - pf(Fcalc3, 1, gl_res)

Fcalc4 <- SQf4 / QMRes
p_valor4 <- 1 - pf(Fcalc4, 1, gl_res)

ANOVA: Fatorial Sexo x Girassol (caso desbalanceado) - Girassol dentro de Sexo F
Grau	gl	SQ	Fcal	p-valor (5%)
1º	1	114.075	29.727	0.0000
2º	1	917.001	238.961	0.0000
3º	1	32.033	8.348	0.0102
4º	1	0.371	0.097	0.7596

Analisando, decrescentemente, a ordem dos graus, o 3º grau é o primeiro a apresentar significância estatística, considerando p-valor = 5%. Assim, o comportamento do ganho de peso em função do aumento da porcentagem de substituição do farelo de soja por farelo de girassol pode ser bem explicado por um polinômio de 3º grau para as fêmeas.

Machos:

SQm1 <- t(t(M1) %*% Mi) %*% solve(t(M1) %*% solve(t(W) %*% W) %*% M1) %*% (t(M1) %*% Mi)

SQm2 <- t(t(M2) %*% Mi) %*% solve(t(M2) %*% solve(t(W) %*% W) %*% M2) %*% (t(M2) %*% Mi)

SQm3 <- t(t(M3) %*% Mi) %*% solve(t(M3) %*% solve(t(W) %*% W) %*% M3) %*% (t(M3) %*% Mi)

SQm4 <- t(t(M4) %*% Mi) %*% solve(t(M4) %*% solve(t(W) %*% W) %*% M4) %*% (t(M4) %*% Mi)


gl_m <- 1


Fcalc1 <- SQm1 / QMRes
p_valor1 <- 1 - pf(Fcalc1, 1, gl_res)

Fcalc2 <- SQm2 / QMRes
p_valor2 <- 1 - pf(Fcalc2, 1, gl_res)

Fcalc3 <- SQm3 / QMRes
p_valor3 <- 1 - pf(Fcalc3, 1, gl_res)

Fcalc4 <- SQm4 / QMRes
p_valor4 <- 1 - pf(Fcalc4, 1, gl_res)

ANOVA: Fatorial Sexo x Girassol (caso desbalanceado) - Girassol dentro de Sexo M
Grau	gl	SQ	Fcal	p-valor (5%)
1º	1	31.734	8.270	0.0105
2º	1	506.002	131.859	0.0000
3º	1	0.000	0.000	0.9928
4º	1	0.439	0.114	0.7392

Enquanto isso, para os machos, o comportamento do ganho de peso em função do aumento da porcentagem de substituição do farelo de soja por farelo de girassol pode ser bem explicado por um polinômio de 2º grau.

A seguir, ajustaremos as curvas de regressão de ambos os modelos. Para isso, utilizaremos a função fat2.dic() do pacote ExpDes.pt, a fim de obter os coeficientes de regressão estimados.

install.packages("ExpDes.pt")
library(ExpDes.pt)

dados <- data.frame(
  Sexo = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2),
  Girassol = c(0, 0, 0, 25, 25, 25, 50, 50, 75, 75, 75, 100, 100, 100, 0, 0, 25, 25, 25, 50, 50, 75, 75, 75, 100, 100, 100),
  Trat = c(1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5, 6, 6, 7, 7, 7, 8, 8, 9, 9, 9, 10, 10, 10),
  Rep = c(1, 2, 3, 1, 2, 3, 1, 3, 1, 2, 3, 1, 2, 3, 2, 3, 1, 2, 3, 1, 2, 1, 2, 3, 1, 2, 3),
  GP = c(77.9, 83.2, 83.5, 71.5, 73.5, 70.5, 67.5, 65.0, 71.5, 70.8, 72.5, 89.5, 91.8, 92.9, 86.0, 84.0, 94.5, 96.0, 95.8, 99.5, 98.0, 93.0, 96.0, 90.5, 83.0, 80.0, 78.5)
)

with(
  dados,
  fat2.dic(Sexo, Girassol, resp = GP, quali = c(TRUE, FALSE))
  )

------------------------------------------------------------------------
Legenda:
FATOR 1:  F1 
FATOR 2:  F2 
------------------------------------------------------------------------


Quadro da analise de variancia
------------------------------------------------------------------------
        GL      SQ QM      Fc    Pr>Fc
F1       1 1158.91  3 302.001 0.000000
F2       4   47.63  2   3.103 0.043427
F1*F2    4 1624.61  5 105.839 0.000000
Residuo 17   65.24  4                 
Total   26 2896.40  1                 
------------------------------------------------------------------------
CV = 2.34 %

------------------------------------------------------------------------
Teste de normalidade dos residuos (Shapiro-Wilk)
valor-p:  0.9783552 
De acordo com o teste de Shapiro-Wilk a 5% de significancia, os residuos podem ser considerados normais.
------------------------------------------------------------------------



Interacao significativa: desdobrando a interacao
------------------------------------------------------------------------

Desdobrando  F1  dentro de cada nivel de  F2 
------------------------------------------------------------------------
------------------------------------------------------------------------
Quadro da analise de variancia
------------------------------------------------------------------------
          GL         SQ         QM       Fc  Pr.Fc
F2         4   47.63486   11.90872   3.1033 0.0434
F1:F2 0    1   14.42133   14.42133   3.7581 0.0693
F1:F2 25   1  835.44000  835.44000  217.707      0
F1:F2 50   1 1056.25000 1056.25000 275.2478      0
F1:F2 75   1  697.68167  697.68167 181.8086      0
F1:F2 100  1  178.21500  178.21500   46.441      0
Residuo   17   65.23667    3.83745                
Total     26 2896.39630  111.39986                
------------------------------------------------------------------------



 F1  dentro do nivel  0  de  F2 

De acordo com o teste F, as medias desse fator sao estatisticamente iguais.
------------------------------------------------------------------------
  Niveis   Medias
1      1 81.53333
2      2 85.00000
------------------------------------------------------------------------


 F1  dentro do nivel  25  de  F2 
------------------------------------------------------------------------
Teste de Tukey
------------------------------------------------------------------------
Grupos Tratamentos Medias
a    2   95.43333 
 b   1   71.83333 
------------------------------------------------------------------------


 F1  dentro do nivel  50  de  F2 
------------------------------------------------------------------------
Teste de Tukey
------------------------------------------------------------------------
Grupos Tratamentos Medias
a    2   98.75 
 b   1   66.25 
------------------------------------------------------------------------


 F1  dentro do nivel  75  de  F2 
------------------------------------------------------------------------
Teste de Tukey
------------------------------------------------------------------------
Grupos Tratamentos Medias
a    2   93.16667 
 b   1   71.6 
------------------------------------------------------------------------


 F1  dentro do nivel  100  de  F2 
------------------------------------------------------------------------
Teste de Tukey
------------------------------------------------------------------------
Grupos Tratamentos Medias
a    1   91.4 
 b   2   80.5 
------------------------------------------------------------------------



Desdobrando  F2  dentro de cada nivel de  F1 
------------------------------------------------------------------------
------------------------------------------------------------------------
Quadro da analise de variancia
------------------------------------------------------------------------
        GL         SQ         QM       Fc Pr.Fc
F1       1 1158.91432 1158.91432 302.0011     0
F2:F1 1  4 1081.49595  270.37399  70.4567     0
F2:F1 2  4  590.74936  147.68734  38.4858     0
Residuo 17   65.23667    3.83745               
Total   26 2896.39630  111.39986               
------------------------------------------------------------------------



 F2  dentro do nivel  1  de  F1 
------------------------------------------------------------------------
Ajuste de modelos polinomiais de regressao
------------------------------------------------------------------------

Modelo Linear
=========================================
   Estimativa Erro.padrao   tc    valor.p
-----------------------------------------
b0  73.3571     0.8864    82.7554    0   
b1   0.0780     0.0143    5.4522  0.00004
-----------------------------------------

R2 do modelo linear
--------
   1    
--------
0.105479
--------

Analise de variancia do modelo linear
=======================================================
                     GL    SQ       QM     Fc   valor.p
-------------------------------------------------------
Efeito linear        1  114.0750 114.0750 29.73 0.00004
Desvios de Regressao 3  967.4209 322.4737 84.03    0   
Residuos             17 65.2367   3.8374               
-------------------------------------------------------
------------------------------------------------------------------------

Modelo quadratico
==========================================
   Estimativa Erro.padrao    tc    valor.p
------------------------------------------
b0  82.6042     1.0662    77.4781     0   
b1  -0.7187     0.0530    -13.5586    0   
b2   0.0080     0.0005    15.6095     0   
------------------------------------------

R2 do modelo quadratico
--------
   1    
--------
0.970037
--------

Analise de variancia do modelo quadratico
========================================================
                     GL    SQ       QM      Fc   valor.p
--------------------------------------------------------
Efeito linear        1  114.0750 114.0750 29.73  0.00004
Efeito quadratico    1  935.0164 935.0164 243.66    0   
Desvios de Regressao 2  32.4046  16.2023   4.22  0.03246
Residuos             17 65.2367   3.8374                
--------------------------------------------------------
------------------------------------------------------------------------

Modelo cubico
=========================================
   Estimativa Erro.padrao   tc    valor.p
-----------------------------------------
b0  81.5708     1.1245    72.5364    0   
b1  -0.4224     0.1154    -3.6601 0.0019 
b2  -0.0003     0.0029    -0.1032 0.9190 
b3   0.0001     0.00002   2.8892  0.0102 
-----------------------------------------

R2 do modelo cubico
--------
   1    
--------
0.999657
--------

Analise de variancia do modelo cubico
========================================================
                     GL    SQ       QM      Fc   valor.p
--------------------------------------------------------
Efeito linear        1  114.0750 114.0750 29.73  0.00004
Efeito quadratico    1  935.0164 935.0164 243.66    0   
Efeito cubico        1  32.0333  32.0333   8.35  0.01019
Desvios de Regressao 1   0.3713   0.3713   0.1   0.75955
Residuos             17 65.2367   3.8374                
--------------------------------------------------------
------------------------------------------------------------------------


 F2  dentro do nivel  2  de  F1 
------------------------------------------------------------------------
Ajuste de modelos polinomiais de regressao
------------------------------------------------------------------------

Modelo Linear
=========================================
   Estimativa Erro.padrao   tc    valor.p
-----------------------------------------
b0  94.1030     0.9940    94.6686    0   
b1  -0.0693     0.0155    -4.4855 0.0003 
-----------------------------------------

R2 do modelo linear
--------
   2    
--------
0.130697
--------

Analise de variancia do modelo linear
=======================================================
                     GL    SQ       QM     Fc   valor.p
-------------------------------------------------------
Efeito linear        1  77.2089  77.2089  20.12 0.00033
Desvios de Regressao 3  513.5405 171.1802 44.61    0   
Residuos             17 65.2367   3.8374               
-------------------------------------------------------
------------------------------------------------------------------------

Modelo quadratico
==========================================
   Estimativa Erro.padrao    tc    valor.p
------------------------------------------
b0  84.9466     1.2709    66.8412     0   
b1   0.5824     0.0584     9.9650     0   
b2  -0.0063     0.0005    -11.5632    0   
------------------------------------------

R2 do modelo quadratico
--------
   2    
--------
0.999255
--------

Analise de variancia do modelo quadratico
========================================================
                     GL    SQ       QM      Fc   valor.p
--------------------------------------------------------
Efeito linear        1  77.2089  77.2089  20.12  0.00033
Efeito quadratico    1  513.1003 513.1003 133.71    0   
Desvios de Regressao 2   0.4402   0.2201   0.06  0.94445
Residuos             17 65.2367   3.8374                
--------------------------------------------------------
------------------------------------------------------------------------

Modelo cubico
=========================================
   Estimativa Erro.padrao   tc    valor.p
-----------------------------------------
b0  84.9390     1.3734    61.8462    0   
b1   0.5840     0.1250    4.6709  0.0002 
b2  -0.0063     0.0030    -2.0906 0.0519 
b3  0.000000    0.00002   0.0146  0.9885 
-----------------------------------------

R2 do modelo cubico
--------
   2    
--------
0.999256
--------

Analise de variancia do modelo cubico
========================================================
                     GL    SQ       QM      Fc   valor.p
--------------------------------------------------------
Efeito linear        1  77.2089  77.2089  20.12  0.00033
Efeito quadratico    1  513.1003 513.1003 133.71    0   
Efeito cubico        1   0.0008   0.0008    0    0.98848
Desvios de Regressao 1   0.4393   0.4393   0.11  0.73924
Residuos             17 65.2367   3.8374                
--------------------------------------------------------
------------------------------------------------------------------------

A seguir, contruiremos os gráficos dos modelos ajustados.

dados <- dados |> dplyr::mutate(Sexo = as.factor(Sexo))

## Fêmeas
dados_f <- dados |> dplyr::filter(Sexo == 1)

## Machos
dados_m <- dados |> dplyr::filter(Sexo == 2)

## Gráfico
library(ggplot2)

dados |> 
  ggplot() + 
  geom_point(aes(x = Girassol, y = GP, color = Sexo)) +
  stat_function(
    data = dados_f, 
    fun = function(Girassol) {81.577 - 0.425*Girassol - 0.0003*Girassol^2 + 0.00006 * Girassol^3},
    color = "#F77D7D", linewidth = 1, linetype = 2
  ) +
  stat_function(
    data = dados_m,
    fun = function(Girassol) {84.951 + 0.5847*Girassol - 0.0063*Girassol^2},
    color = "#52A7E8", linewidth = 1, linetype = 2
  ) +
  theme_minimal() +
  labs(x = NULL, y = "Ganho de peso", color = NULL) +
  scale_color_discrete(labels = c("Fêmea", "Macho")) +
  theme(legend.position = "bottom")

proc glm;
class trat;
model GP = trat / ss3;

contrast 'Sexo' trat 1  1  1  1  1 -1 -1 -1 -1 -1;
contrast 'Girassol'  trat -2 -1  0  1  2 -2 -1  0  1  2,
                     trat  2 -1 -2 -1  2  2 -1 -2 -1  2,
                     trat -1  2  0 -2  1 -1  2  0 -2  1,
                     trat  1 -4  6 -4  1  1 -4  6 -4  1;
contrast 'Interação' trat -4  1  1  1  1  4 -1 -1 -1 -1,
                     trat  0 -3  1  1  1  0  3 -1 -1 -1,
                     trat  0  0 -2  1  1  0  0  2 -1 -1,
                     trat  0  0  0 -1  1  0  0  0  1 -1;

contrast 'Girassol=  0: F = M' trat 1 0 0 0 0 -1  0  0  0  0;
contrast 'Girassol= 25: F = M' trat 0 1 0 0 0  0 -1  0  0  0;
contrast 'Girassol= 50: F = M' trat 0 0 1 0 0  0  0 -1  0  0;
contrast 'Girassol= 75: F = M' trat 0 0 0 1 0  0  0  0 -1  0;
contrast 'Girassol=100: F = M' trat 0 0 0 0 1  0  0  0  0 -1;
contrast 'Machos: Girassol grau 1' trat -2 -1  0  1  2  0  0  0  0  0;
contrast 'Machos: Girassol grau 2' trat  2 -1 -2 -1  2  0  0  0  0  0;
contrast 'Machos: Girassol grau 3' trat -1  2  0 -2  1  0  0  0  0  0;
contrast 'Machos: Girassol grau 4' trat  1 -4  6 -4  1  0  0  0  0  0;
contrast 'Fêmeas: Girassol grau 1' trat  0  0  0  0  0 -2 -1  0  1  2;
contrast 'Fêmeas: Girassol grau 2' trat  0  0  0  0  0  2 -1 -2 -1  2;
contrast 'Fêmeas: Girassol grau 3' trat  0  0  0  0  0 -1  2  0 -2  1;
contrast 'Fêmeas: Girassol grau 4' trat  0  0  0  0  0  1 -4  6 -4  1;
run;

proc glm;
class Sexo Girassol;
model GP = Sexo Girassol  / ss1 ss2 ss3;
run;