26 Análise de covariância

26.1 Exemplo 1

Os dados disponíveis em neste link são um exemplo onde a análise de covariância pode ser utilizada. Neste caso temos uma variável resposta e duas variáveis explicativas sendo uma delas qualitativa (um fator) e outra quantitativa (numérica). O arquivo com conjunto de dados para sua área de trabalho de depois importando para o R com read.table() Alternativamente, se o computador estiver conectado a internet read.table() pode ler diretamente o arquivo como mostrado a seguir.

O primeiro passo é a leitura e organização dos dados, em particular assegurando que a variável qualitativa indicadora dos tratamentos seja indicada como sendo um fator evitando que seus valores sejam interpretados como quantidades numéricas. Note que neste caso temos 2 variáveis numéricas, a resposta (resp) e a covariável (cov).

  > ex12 <- read.table("http://www.leg.ufpr.br/~paulojus/dados/exemplo12.txt",
  +     header = T)
  > ex12

     maq cov resp
  1    1  20   36
  2    1  25   41
  3    1  24   39
  4    1  25   42
  5    1  32   49
  6    2  22   40
  7    2  28   48
  8    2  22   39
  9    2  30   45
  10   2  28   44
  11   3  21   35
  12   3  23   37
  13   3  26   42
  14   3  21   34
  15   3  15   32

  > ex12$maq <- as.factor(ex12$maq)
  > str(ex12)

  'data.frame': 15 obs. of  3 variables:
   $ maq : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 2 2 2 2 2 ...
   $ cov : int  20 25 24 25 32 22 28 22 30 28 ...
   $ resp: int  36 41 39 42 49 40 48 39 45 44 ...

  > summary(ex12)

   maq        cov             resp
   1:5   Min.   :15.00   Min.   :32.0
   2:5   1st Qu.:21.50   1st Qu.:36.5
   3:5   Median :24.00   Median :40.0
         Mean   :24.13   Mean   :40.2
         3rd Qu.:27.00   3rd Qu.:43.0
         Max.   :32.00   Max.   :49.0

Os dados podem ser visualizados na Figura 26.1 produzida com o comando a seguir. Note-se que os dados de cada um dos tratamentos são indicados por cores (argumento col)] e padrões de pontos (argumento pch) indexados pema variável maq que define os níveis dos tratamentos

  > with(ex12, plot(resp ~ cov, col = c(1, 2, 4)[maq], pch = (1:3)[maq]))


PIC
Figura 68: Dados para análise de covariância.


Na análise de covariância não temos ortogonalidade entre os fatores. Desta forma os testes de significância tem que ser obtidos em ajustes separados: (i) para o efeito de covariáveis, corrigido pelo efeito dos tratamentos qualitativos e (ii) para o efeito dos tratamentos qualitativos, corrigidos pelo efeito da covariável.

Primeiro vamos testar a inclinação (coeficiente β1) da reta de regressão. Na análise de variância abaixo devemos considerar apenas o teste referente à variável cov que neste caso está corrigida para o efeito de maq. Note que para isto a variável cov tem que ser a última na especificação do modelo.

  > ex12.cov <- aov(resp ~ maq + cov, data = ex12)
  > summary(ex12.cov)

              Df  Sum Sq Mean Sq F value    Pr(>F)
  maq          2 140.400  70.200  27.593 5.170e-05 ***
  cov          1 178.014 178.014  69.969 4.264e-06 ***
  Residuals   11  27.986   2.544
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

A seguir testamos o efeito do fator maq corrigindo para o efeito da covariável. Para isto basta inverter a ordem dos termos na especificação do modelo.

  > ex12.trat <- aov(resp ~ cov + maq, data = ex12)
  > summary(ex12.trat)

              Df  Sum Sq Mean Sq  F value   Pr(>F)
  cov          1 305.130 305.130 119.9330 2.96e-07 ***
  maq          2  13.284   6.642   2.6106   0.1181
  Residuals   11  27.986   2.544
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Portanto, olhamos o primeiro quadro da ANOVA para verificar o efeito da covariável, e no segundo para verificar o efeito do tratamento. Se desejável poderia-se tomar os resultados de cada um deles para compor um quandro de análise, porém com a resalva que, devido a não ortogonalidade, a soma das somas de quadrados não corresonde a soma de quadrados total. Entretanto, há uma função Anova() no pacote car do R que já monta tal quadro automaticamente conforme ilustrado a seguir.

  > require(car)
  > Anova(ex12.cov, type = "III")

  Anova Table (Type III tests)
  
  Response: resp
               Sum Sq Df F value    Pr(>F)
  (Intercept)  87.434  1 34.3664 0.0001089 ***
  maq          13.284  2  2.6106 0.1180839
  cov         178.014  1 69.9694 4.264e-06 ***
  Residuals    27.986 11
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Note que esta função irá retornar o mesmo resultado para qualquer ordem dos termos no modelo, ou seja, no exemplo acima Anova(ex12.cov, type="III") e Anova(ex12.trat, type="III") retornam os mesmos resultados.

O argumento type="III" refere-se a um jargão consagrado pelo software SAS que corresponde a soma de quadrados do tipo III. Em geral nas funções básicas do R evita-se tal jargão e procura-se usar so conceitos ligados à parametrização do modelo através da definição dos contrastes e por isto tal terminologia está apenas em um pacote contribuído.

Neste caso a função Anova faz o mesmo que mostrado nas duas análises de variâncias iniciais, obtendo para cada termo a soma de quadrados quando este é corrigido para os demais, ou seja, colocado na última posição na especificação do modelo.