Nesta sessão iremos usar o R para analisar um experimento em delineamento inteiramente casualizado com apenas um fator. Tal procedimento é também chamado em alguns textos de "análise da variância de simples entrada"(one-way anova). A seguir são apresentados os comandos exemplificando alguns procedimentos usuais para a análise dos dados de um experimento deste tipo que, neste exemplo, envolve um fator com nove níveis (tratamentos). O primeiro passo é ler os dados.
> ex01 <- read.table("http://www.leg.ufpr.br/~paulojus/dados/exemplo01.txt", head=T, stringsAsFactors=TRUE)
Caso não consiga executar o comando acima diretamente com o endereço http utilize um navegador para ir até esta página e copie o arquivo exemplo1.txt para o seu diretório de trabalho. Caso o arquivo esteja em outro diretório deve-se colocar o caminho completo deste diretório no argumento de read.table() acima. A seguir vamos inspecionar o objeto que armazena os dados e seus componentes. Em particular é importante certificar-se que a variável resposta é do tipo numeric e, se os níveis de tratamentos forem qualitativos, a variável indicadora dos tratamentos é do tipo factor . Caso isto não ocorra é necessário transformar as variáveis para estes tipos antes de prosseguir com as análises.
> head(ex01)
trat resp 1 t1 385 2 t1 323 3 t1 417 4 t1 370 5 t1 437 6 t1 340
> is.numeric(ex01$resp)
[1] TRUE
> is.factor(ex01$trat)
[1] TRUE
Portanto o objeto ex01 é um data-frame com duas variáveis, sendo uma delas um fator (a variável trat) e a outra uma variável numérica (resp). Vamos iniciar obtendo um rápido resumo dos dados que mostra que este é um experimento "balanceado"com mesmo número de repetições (seis) para cada tratamento. Calculamos também as médias, variâncias e erros padrão das médias para cada tratamento separadamente.
> summary(ex01)
trat resp t1 : 6 Min. :115.0 t2 : 6 1st Qu.:307.5 t3 : 6 Median :377.5 t4 : 6 Mean :353.5 t5 : 6 3rd Qu.:417.0 t6 : 6 Max. :474.0 (Other):18
> ex01.nrep <- with(ex01, tapply(resp, trat, length)) > ex01.mds <- with(ex01, tapply(resp, trat, mean)) > ex01.var <- with(ex01, tapply(resp, trat, var)) > ex01.se <- with(ex01, tapply(resp, trat, function(x) sqrt(var(x)/length(x)))) > data.frame(Repet= ex01.nrep, Medias=ex01.mds, Variancias = ex01.var, EP=ex01.se, row.names=paste("trat", 1:9, sep="-"))
Repet Medias Variancias EP trat-1 6 378.6667 1916.267 17.87114 trat-2 6 431.5000 987.500 12.82900 trat-3 6 346.3333 3117.867 22.79571 trat-4 6 293.6667 3494.667 24.13389 trat-5 6 341.8333 1513.767 15.88378 trat-6 6 406.0000 1903.600 17.81198 trat-7 6 164.1667 2173.367 19.03228 trat-8 6 403.8333 1242.167 14.38846 trat-9 6 415.6667 1091.067 13.48497
Vamos prosseguir com a análise exploratória com gráficos gerados pelos comandos a seguir e mostrados na Figura 57. O gráfico de esquerda utiliza a função boxcox() do pacote MASS para verificar a necessidade de transformação dos dados o que neste caso não é necessária visto que o valor um está contido no intervalo definido pelas linhas tracejadas. A transformação Box-Cox é discutida em mais detalhes em uma outra Seção deste material. O gráfico do meio mostra um boxplot para os dados de cada tratamento, o que deve ser analisado com cautela lembrando que cada boxplot é produzido com apenas seis observações. Optamos aqui por indicar também neste gráfico a média de cada tratamento. O gráfico da direita produzido com stripchart() é uma alternativa ao boxplot para amostras de tamanho pequeno. Na chamada desta função optamos por alterar valores default de alguns argumentos como por exemplo para method="jitter" que provoca um pequeno deslocamento horizontal aleatório dos pontos evitando assim sobreposição de pontos com valores coincidentes ou muito próximos. Ainda neste gráfico acrescentamos as médias e barras que somam e subtraem os erros padrões da média para cada tratamento. Na função arrows() os quatro argumentos iniciais informam coordenadas para as barras, code=3 informa que as "setas"devem ser colocadas em ambas extremidades e angle=90 faz com que a "seta"se torne uma pequena barra horizontal com o tamanho controlado por length.
> require(MASS); boxcox(resp ~trat, lambda=seq(0,3,l=101), data=ex01) > plot(ex01) > points(ex01.mds, pch="x", col=2, cex=1.5) > with(ex01, stripchart(resp~trat, met="jitter", vert=T, pch=19)) > points(ex01.mds, pch=4, cex=1.5) > arrows(1:9,ex01.mds+ex01.se,1:9, ex01.mds-ex01.se, angle=90, code=3, length=0.1)
É importante notar que as barras simplesmente refletem a variância dos dados dentro de cada tratamento e não são adequadas para detectar diferenças entre tratamentos, o que será discutido mais adiante nesta sessão. Além dos gráficos acima podemos também verificar o pressuposto de homogeneidade de variâncias com o testes de igualdade de variâncias, como por exemplo, o teste de Bartlett. Neste caso o teste indica variâncias homogêneas. Caso isto não ocorresse uma possível alternativa seria usar o procedimento descrito na Sessão 24.3.
> bartlett.test(resp ~ trat, data=ex01)
Bartlett test of homogeneity of variances
data: resp by trat
Bartlett's K-squared = 3.6738, df = 8, p-value = 0.8853
Uma vez concluída a análise exploratória e verificada a adequação de alguns pressupostos o passo seguinte é ajustar o modelo usando aov() ou lm(). Neste exemplo, por se tratar da análise de um experimento, tipicamente avaliada pelo quadro de análise de variância, optamos por usar aov(). Embora aov() use lm() internamente, os resultados são organizados internamente de forma conveniente para a efetuar a análise de variância.
> ex01.mod <- aov(resp ~ trat, data = ex01) > ex01.mod
Call: aov(formula = resp ~ trat, data = ex01) Terms: trat Residuals Sum of Squares 332918.1 87201.3 Deg. of Freedom 8 45 Residual standard error: 44.02053 Estimated effects may be unbalanced
> anova(ex01.mod)
Analysis of Variance Table Response: resp Df Sum Sq Mean Sq F value Pr(>F) trat 8 332918 41615 21.475 5.445e-13 Residuals 45 87201 1938
Portanto o objeto ex01.mod é uma lista que guarda os resultados da análise para o modelo ajustado. Vamos inspecionar este objeto e seus elementos mais detalhadamente ilustrando como usá-lo para obter a análise dos resultados e extrair elementos para a análise de resíduos. A função names() mostra os elementos da lista e adicionalmente existem funções que extraem elementos do objeto. Duas tipicamente utilizadas são coef() para extrair os coeficientes, residuals() para extrair resíduos e fitted() para valores ajustados, mas há ainda várias outras como effects(), AIC() logLik(), model.tables(), entre outras.
> names(ex01.mod)
[1] "coefficients" "residuals" "effects" "rank" "fitted.values" [6] "assign" "qr" "df.residual" "contrasts" "xlevels" [11] "call" "terms" "model"
> coef(ex01.mod)
(Intercept) tratt2 tratt3 tratt4 tratt5 tratt6 tratt7 378.66667 52.83333 -32.33333 -85.00000 -36.83333 27.33333 -214.50000 tratt8 tratt9 25.16667 37.00000
> model.tables(ex01.mod)
Tables of effects trat trat t1 t2 t3 t4 t5 t6 t7 t8 t9 25.15 77.98 -7.19 -59.85 -11.69 52.48 -189.35 50.31 62.15
> model.tables(ex01.mod, type="means")
Tables of means Grand mean 353.5185 trat trat t1 t2 t3 t4 t5 t6 t7 t8 t9 378.7 431.5 346.3 293.7 341.8 406.0 164.2 403.8 415.7
O resultado de coef() vai depender da parametrização adotada e definida pelos contrastes. Os valores default e/ou correntes são dados por options()$contrasts. Para fatores qualitativos como no caso deste exemplo a parametrização default corresponde a "contr.treatment" que assinala o valor da média do primeiro tratamento (primeiro nível do fator) ao primeiro coeficiente. Os demais representam a diferença das médias de cada um dos tratamentos a este tratamento de referência.
Uma outra forma de especificar o modelo para este exemplo é mostrada a seguir com o uso -1 que, para níveis quantitativos corresponde a ajustar um modelo com intercepto igual a zero. No caso de níveis qualitativos como neste exemplo, monta uma matriz do modelo de forma a que cada coeficiente corresponda à média de cada um dos tratamentos. Note que apenas a interpretação dos coeficientes muda e a análise de variância permanece a mesma.
> ex01.mod1 <- aov(resp ~ trat - 1, data = ex01) > coef(ex01.mod1)
tratt1 tratt2 tratt3 tratt4 tratt5 tratt6 tratt7 tratt8 tratt9 378.6667 431.5000 346.3333 293.6667 341.8333 406.0000 164.1667 403.8333 415.6667
> anova(ex01.mod1)
Analysis of Variance Table Response: resp Df Sum Sq Mean Sq F value Pr(>F) trat 9 7081587 786843 406.05 < 2.2e-16 Residuals 45 87201 1938
A parametrização para os coeficientes é determinada pela matriz do modelo e é definida pelo argumento contrasts de options() ou pela função contrasts() que mostra ou atribui a matriz de contrastes a ser utilizada. Fatores são definidos como sendo unordered (por exemplo níveis qualitativos como no caso da análise vista aqui) ou ordered, o que é usado, por exemplo, no caso de níveis quantitativos.
> options()$contrasts
unordered ordered "contr.treatment" "contr.poly"
> contrasts(ex01$trat)
t2 t3 t4 t5 t6 t7 t8 t9 t1 0 0 0 0 0 0 0 0 t2 1 0 0 0 0 0 0 0 t3 0 1 0 0 0 0 0 0 t4 0 0 1 0 0 0 0 0 t5 0 0 0 1 0 0 0 0 t6 0 0 0 0 1 0 0 0 t7 0 0 0 0 0 1 0 0 t8 0 0 0 0 0 0 1 0 t9 0 0 0 0 0 0 0 1
Para definir a parametrização a ser utilizada e definida pelos contrastes, pode-se usar outras opções de contrastes já disponibilizadas pelo R tipicamente usando options(). Nos comandos a seguir alteramos a opção para fatores unordered para "contr.sum". Os coeficientes obtidos são diferentes dos obtidos anteriormente sendo o primeiro a média geral e os demais uma comparação da média de cada tratamento contra as médias dos demais. Os resultados da análise de variância permanecem inalterado.
> options(contrasts=c("contr.sum", "contr.poly")) > contrasts(ex01$trat)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] t1 1 0 0 0 0 0 0 0 t2 0 1 0 0 0 0 0 0 t3 0 0 1 0 0 0 0 0 t4 0 0 0 1 0 0 0 0 t5 0 0 0 0 1 0 0 0 t6 0 0 0 0 0 1 0 0 t7 0 0 0 0 0 0 1 0 t8 0 0 0 0 0 0 0 1 t9 -1 -1 -1 -1 -1 -1 -1 -1
> coef(lm(resp~trat, data=ex01))
(Intercept) trat1 trat2 trat3 trat4 trat5 trat6 353.518519 25.148148 77.981481 -7.185185 -59.851852 -11.685185 52.481481 trat7 trat8 -189.351852 50.314815
Os contrastes já definidos no R são listados e descritos a seguir. Além destes outros pacotes podem ter outras definições de contrastes, como em de "contr.sdif" do pacotes MASS. Estes contrastes terão efeito se o termo -1 não for incluído no modelo pois neste caso os coeficientes são sempre as médias de cada um dos tratamentos, independente da opção de contraste adotada.
Além dos contrastes pré definidos, outros contrastes definidos pelo usuário e atribuídos ao fator em estudo usando a função contrasts(). Retornamos a este tópico com um exemplo na Sessão 24.1.3.
Retornando à análise do exemplo, vamos ver agora alguns gráficos e ferramentas para avaliar o modelo ajustado. Um método associado a plot() produz automaticamente gráficos de resíduos para objetos das classes lm e aov conforme ilustrado na Figura 58 produzida com o comando plot(ex01.mod).
Além dos gráficos "pré-preparados"pelo R, o usuário pode obter outros que desejar extraindo a informação necessária do objeto que contém o ajuste do modelo. Na Figura 59 mostramos quatro gráficos: resíduos padronizados versus valores preditos, boxplot, histograma dos resíduos padronizados e qqplot() dos resíduos. Para isto obtemos os resíduos padronizados dividindo os resíduos do modelo pela raiz quadrada da variância do termo de erro.
> ex01.res <- resid(ex01.mod) > respad <- ex01.res/sqrt(sum(ex01.res^2)/ex01.mod$df.res) > plot(fitted(ex01.mod), respad, xlab="valores ajustados", ylab="resíduos") > title("Resíduos Padronizados vs \n Valores Preditos") > boxplot(respad); title("Resíduos Padronizados") > hist(respad, main="Histograma dos resíduos padronizados") > qqnorm(ex01.res,ylab="Residuos", main=NULL); qqline(ex01.res) > title("Gráfico Normal de \n Probabilidade dos Resíduos")
Um teste de normalidade dos resíduos pode ser efetuado como indicado a seguir.
> shapiro.test(ex01.res)
Shapiro-Wilk normality test
data: ex01.res
W = 0.97159, p-value = 0.2263
Uma das formas possíveis formas de interpretar os resultados no caso de efeito de tratamentos significativos é utilizar algum procedimento de comparações de tratamentos após verificar o resultado da anova, o que justifica o termo às vezes utilizado que descreve tais procedimentos como comparações post-hoc.
A questão do uso de comparações de tratamentos é polêmica no meio estatístico e não vamos aqui entrar em tal discussão. Nesta sessão vamos ilustrar três procedimentos deixando a cargo do leitor o julgamento de qual procedimento é mais adequado para o problema em questão. Os procedimentos discutidos a seguir correspondem a três possíveis abordagens ao problema de comparação de mais de duas médias, sendo eles: (i) teste-\(t\) para comparações duas a duas, (ii) teste de Tukey e (iii) contrastes e contrastes ortogonais. O primeiro caso se desdobra em mais opções uma vez que permite que os valores \(p\) sejam ou não ajustados, e caso sejam, por diferentes métodos.
Os procedimentos mostrados aqui são implementados em pacotes básicos do R. O pacote multcomp disponibiliza uma extensa lista de procedimentos adicionais de comparações múltiplas e alguns procedimentos específicos podem ainda ser encontrados em outros pacotes do R.
A função pairwise.t.test() calcula todas as possíveis comparações entre dois grupos, podendo ser vista como uma extensão ao teste-\(t\) para duas amostras, retornando o valor-\(p\) para cada comparação. A principal diferença é que o nível de significância deve ser corrigido para garantir o nível de significância conjunto para todas as comparações. O argumento p.adjust.method da função permite o usuário escolher entre diferentes métodos propostos para ajustar o nível de significância sendo o default o procedimento proposto por Holm, que é uma modificação ao ajuste de Bonferroni, que também está disponível através do argumento p.adj="bonferroni". Mais detalhes podem ser encontrados na documentação da função.
> with(ex01,pairwise.t.test(resp, trat))
Pairwise comparisons using t tests with pooled SD
data: resp and trat
t1 t2 t3 t4 t5 t6 t7 t8
t2 0.65049 - - - - - - -
t3 1.00000 0.03768 - - - - - -
t4 0.03768 6.4e-05 0.65049 - - - - -
t5 1.00000 0.02345 1.00000 0.83853 - - - -
t6 1.00000 1.00000 0.39692 0.00160 0.28829 - - -
t7 2.5e-09 3.8e-12 1.8e-07 0.00019 3.2e-07 8.2e-11 - -
t8 1.00000 1.00000 0.45676 0.00203 0.33686 1.00000 1.0e-10 -
t9 1.00000 1.00000 0.18109 0.00048 0.11918 1.00000 2.5e-11 1.00000
P value adjustment method: holm
O teste Tukey de comparações múltiplas é implementado na função TukeyHSD(). A saída em formato texto do teste de Tukey é mostrada a seguir e plot(ex01.HSD) produz o gráfico mostrado na Figura 60. As saídas da função mostram intervalos de confiança para as diferenças entre pares de médias.
> ex01.HSD <- TukeyHSD(ex01.mod, ordered=TRUE) > ex01.HSD
Tukey multiple comparisons of means 95% family-wise confidence level factor levels have been ordered Fit: aov(formula = resp ~ trat, data = ex01) $trat diff lwr upr p adj t4-t7 129.500000 46.719034 212.28097 0.0002153 t5-t7 177.666667 94.885701 260.44763 0.0000004 t3-t7 182.166667 99.385701 264.94763 0.0000002 t1-t7 214.500000 131.719034 297.28097 0.0000000 t8-t7 239.666667 156.885701 322.44763 0.0000000 t6-t7 241.833333 159.052367 324.61430 0.0000000 t9-t7 251.500000 168.719034 334.28097 0.0000000 t2-t7 267.333333 184.552367 350.11430 0.0000000 t5-t4 48.166667 -34.614299 130.94763 0.6203900 t3-t4 52.666667 -30.114299 135.44763 0.5040619 t1-t4 85.000000 2.219034 167.78097 0.0401018 t8-t4 110.166667 27.385701 192.94763 0.0024139 t6-t4 112.333333 29.552367 195.11430 0.0018566 t9-t4 122.000000 39.219034 204.78097 0.0005599 t2-t4 137.833333 55.052367 220.61430 0.0000730 t3-t5 4.500000 -78.280966 87.28097 1.0000000 t1-t5 36.833333 -45.947633 119.61430 0.8721075 t8-t5 62.000000 -20.780966 144.78097 0.2886707 t6-t5 64.166667 -18.614299 146.94763 0.2479215 t9-t5 73.833333 -8.947633 156.61430 0.1146645 t2-t5 89.666667 6.885701 172.44763 0.0247945 t1-t3 32.333333 -50.447633 115.11430 0.9342210 t8-t3 57.500000 -25.280966 140.28097 0.3855262 t6-t3 59.666667 -23.114299 142.44763 0.3369467 t9-t3 69.333333 -13.447633 152.11430 0.1671352 t2-t3 85.166667 2.385701 167.94763 0.0394343 t8-t1 25.166667 -57.614299 107.94763 0.9849417 t6-t1 27.333333 -55.447633 110.11430 0.9749062 t9-t1 37.000000 -45.780966 119.78097 0.8693183 t2-t1 52.833333 -29.947633 135.61430 0.4998060 t6-t8 2.166667 -80.614299 84.94763 1.0000000 t9-t8 11.833333 -70.947633 94.61430 0.9999286 t2-t8 27.666667 -55.114299 110.44763 0.9730043 t9-t6 9.666667 -73.114299 92.44763 0.9999849 t2-t6 25.500000 -57.280966 108.28097 0.9836416 t2-t9 15.833333 -66.947633 98.61430 0.9993743
Visualizações mais convenientes dos resultados podem ser obtidas com operações sobre o objeto resultante, tal como a usualmente adotada de listar as médias em ordem decrescente e indicar com letras as diferenças significativas ou não entre estas médias. Vamos ilustrar aqui uma possível forma de obter tal visualização. Inicialmente vamos obter a DMS (diferença mínima significativa). No caso deste experimento balanceado, isto é, o mesmo número de repetições em cada tratamento, o intervalo de confiança para cada diferença é o mesmo e a DMS é portanto comum e dada por metade da amplitude do intervalo.
> dms <- unname(0.5 * diff(ex01.HSD[[1]][1,2:3])) > dms
[1] 82.78097
O passo seguinte é ordenar as médias de forma decrescente e verificar as diferenças significativas. O código abaixo é uma (mas certamente não a única) maneira de indicar as diferenças significativas usando código de letras usual na literatura.
> ex01.mds.ord <- sort(ex01.mds, decreasing = TRUE) > i <- pos <- letra <- 1 > letras <- character(nlevels(ex01$trat)) > while(i <= nlevels(ex01$trat)){ + print(letters[letra]) + ind <- (ex01.mds.ord[i] - (ex01.mds.ord[-(1:i)])) < dms + pos.i <- i+sum(ind) + if(pos.i > pos){ + letras.vec <- rep(" ",length(letras)) + letras.vec[i:pos.i] <- letters[letra] + letras <- paste(letras, letras.vec, sep="") + pos <- pos.i + letra <- letra + 1 + } + i <- i + 1 + }
[1] "a" [1] "b" [1] "c" [1] "c" [1] "c" [1] "c" [1] "d" [1] "d" [1] "d"
> data.frame(medias = ex01.mds.ord, diferencas = letras)
medias diferencas t2 431.5000 a t9 415.6667 ab t6 406.0000 ab t8 403.8333 ab t1 378.6667 ab t3 346.3333 bc t5 341.8333 bc t4 293.6667 c t7 164.1667 d
Neste caso o procedimento é simples pois para um experimento balanceado e pelo teste de Tukey tem-se apenas um único valor de DMS. O algoritmo deve ser modificado e generalizado para outras situações ou pode-se usar funções de pacotes como multcompLetters do pacote multcompView.
Na análise de experimentos pode-se ter interesse em estudar determinadas comparações entre as médias que podem ser especificadas pelo usuário na forma de contrastes, que são um caso particular das funções estimáveis para o modelo. Vamos iniciar revendo definições.
Seja o modelo linear escrito na forma matricial \(Y = X \beta + \epsilon \) onde \(Y\) é a variável resposta, \(X\) a matriz do modelo, \(\beta \) o vetor de \(p\) parâmetros (coeficientes) e \(\epsilon \) o vetor de erros. Uma combinação linear dos coeficientes da forma \(\sum _p \lambda _p \beta _p\) onde \(\lambda =[\lambda _1, \ldots , \lambda _p]\) é um vetor de constantes é dita uma função estimável para o dado modelo se \(\lambda \) pode ser escrita como uma combinação linear das linhas da \(X\). Um contraste é um caso especial de função estimável em que a soma das constantes é nula, isto é, pode ser escrito como \(\sum _p c_p \beta _p\) onde \(\sum _p c_p = 0\).
No que se segue vamos ver como obter estimativas de contrastes de interesse no R, onde fórmulas lineares são usadas para definir as matrizes do modelo usadas no ajuste de modelos lineares e lineares generalizados. No caso de fatores (qualitativos) a matriz \(X\) do modelo não é definida unicamente para um mesmo experimento, podendo ser escrita de diversas formas alternativas que irão produzir a ajustes equivalentes. Tais formas são definidas pela escolha de contrastes ou funções estimáveis que definirão a interpretação dos coeficientes \(\beta \) do modelo. Portanto, se o interesse é apenas na análise de variância a particular forma adotada é irrelevante. Por outro lado, a escolha deve ser adequada se os coeficientes devem ser interpretados.
Ao ajustar um modelo as estimativas de contrastes podem ser obtidas de duas formas:
Vamos discutir aqui algumas ideias iniciais sobre como implementar a segunda forma. Como na análise de contrastes os coeficientes passam a ser diretamente interpretados, passamos a usar lm() no ajuste do modelo.
Uma classe especial de contrastes é a de contrastes ortogonais. Um conjunto de contrastes ortogonais tem a propriedade de que a soma dos produtos dos coeficientes de qualquer par de contrastes deste conjunto é nula. Contrastes ortogonais são particularmente interessantes pois permitem desdobrar (particionar) a soma de quadrados de tratamentos em parcelas referentes a cada um dos contrastes. Isto permite que cada contraste seja testado diretamente por um teste \(t\) (ou o equivalente teste \(F\)).
Com nove tratamentos é possível definir oito contrastes ortogonais com cada um deles sendo associado a um dos graus de liberdade dos tratamentos. A definição destes contrastes não é única e deve refletir comparações relevantes para o problema em questão, assegurando-se que a ortogonalidade seja mantida o que garante que a soma das somas de quadrados dos contrastes seja equivalente à soma de quadrados total dos tratamentos. Para obter o desdobramento abordamos a modelagem como um problema de regressão múltipla onde os contrastes definem variáveis quantitativas a serem incluídas no modelo que é ajustado com lm(). Neste exemplo vamos considerar o seguinte conjunto de contrastes entre as médias dos tratamentos que são especificados nas linhas de uma matriz como se segue.
> c1 <- rbind(c(2,2,2,-1,-1,-1,-1,-1,-1), + c(2,-1,-1, rep(0,6)), + c(0, 1,-1, rep(0,6)), + c(rep(0,3), c(2,2,-1,-1,-1,-1)), + c(rep(0,3), c(1,-1), rep(0,4)), + c(rep(0,5), c(1,1,-1,-1)), + c(rep(0,5), c(1,-1, 0,0)), + c(rep(0,5), c(0,0, 1,-1))) > c1
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [1,] 2 2 2 -1 -1 -1 -1 -1 -1 [2,] 2 -1 -1 0 0 0 0 0 0 [3,] 0 1 -1 0 0 0 0 0 0 [4,] 0 0 0 2 2 -1 -1 -1 -1 [5,] 0 0 0 1 -1 0 0 0 0 [6,] 0 0 0 0 0 1 1 -1 -1 [7,] 0 0 0 0 0 1 -1 0 0 [8,] 0 0 0 0 0 0 0 1 -1
O próximo passo é fazer com que a matriz do modelo seja montada pelo R de forma que os coeficientes reflitam os contrastes desejados. Para isto associamos ao fator que representa os tratamentos (trat no exemplo) o atributo contrast contendo a inversa generalizada obtida por ginv() do pacote MASS. A análise de variância deste modelo é a mesma obtida anteriormente. entretanto os coeficientes são agora dados pela média geral seguida pelas estimativas de cada um dos oito contrastes definidos que podem ser testadas diretamente pelo teste-\(t\) usando o comando summary().
> c1.ginv <- ginv(c1) > colnames(c1.ginv) <- paste("contr", 1:8, sep="") > contrasts(ex01$trat) <- c1.ginv > mod1 <- lm(resp~trat, data=ex01) > anova(mod1)
Analysis of Variance Table Response: resp Df Sum Sq Mean Sq F value Pr(>F) trat 8 332918 41615 21.475 5.445e-13 Residuals 45 87201 1938
> summary(mod1)
Call: lm(formula = resp ~ trat, data = ex01) Residuals: Min 1Q Median 3Q Max -85.67 -33.29 4.75 33.17 85.67 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 353.52 5.99 59.014 < 2e-16 tratcontr1 287.83 76.25 3.775 0.000466 tratcontr2 -20.50 44.02 -0.466 0.643682 tratcontr3 85.17 25.41 3.351 0.001638 tratcontr4 -118.67 62.25 -1.906 0.063029 tratcontr5 -48.17 25.41 -1.895 0.064503 tratcontr6 -249.33 35.94 -6.937 1.26e-08 tratcontr7 241.83 25.41 9.515 2.41e-12 tratcontr8 -11.83 25.41 -0.466 0.643748 Residual standard error: 44.02 on 45 degrees of freedom Multiple R-squared: 0.7924, Adjusted R-squared: 0.7555 F-statistic: 21.48 on 8 and 45 DF, p-value: 5.445e-13
Nos comandos a seguir visualizamos os mesmos resultados de uma forma alternativa, usando model.matrix() para montar a matriz de covariáveis da forma desejada, onde excluímos o intercepto (primeira coluna) e, para visualização adequada dos resultados, trocamos os nomes das colunas. A este data-frame adicionamos os dados e ajustamos o modelo de regressão com lm(). A função anova() sobre o modelo ajustado exibe a soma de quadrados decomposta entre os contrastes agora testados pelo teste \(F\) que é equivalente ao teste-\(t\) mostrado acima pois cada contraste possui um grau de liberdade. Note que a soma delas corresponde a soma de quadrados de tratamentos mostrada no ajuste inicial do modelo e os coeficientes são os mesmos.
> ex01co <- data.frame(model.matrix(resp~trat, ex01)[,-1]) > names(ex01co) <- paste("Contraste", 1:8) > ex01co$resp <- ex01$resp > mod2 <- lm(resp ~ ., data=ex01co) > av2 <- anova(mod2); av2
Analysis of Variance Table Response: resp Df Sum Sq Mean Sq F value Pr(>F) `Contraste 1` 1 27616 27616 14.2512 0.000466 `Contraste 2` 1 420 420 0.2169 0.643682 `Contraste 3` 1 21760 21760 11.2292 0.001638 `Contraste 4` 1 7041 7041 3.6334 0.063029 `Contraste 5` 1 6960 6960 3.5917 0.064503 `Contraste 6` 1 93251 93251 48.1217 1.264e-08 `Contraste 7` 1 175450 175450 90.5405 2.409e-12 `Contraste 8` 1 420 420 0.2168 0.643748 Residuals 45 87201 1938
> sum(av2$Sum[1:8])
[1] 332918.1
> coef(mod2)
(Intercept) `Contraste 1` `Contraste 2` `Contraste 3` `Contraste 4` `Contraste 5` 353.51852 287.83333 -20.50000 85.16667 -118.66667 -48.16667 `Contraste 6` `Contraste 7` `Contraste 8` -249.33333 241.83333 -11.83333
Os coeficientes retornados equivalem a aplicar os contrastes desejados sobre as médias dos tratamentos. Pode-se ainda visualizar os contrastes assinalados ao fator trat através da inversa generalizada.
> drop(c1 %*% ex01.mds)
[1] 287.83333 -20.50000 85.16667 -118.66667 -48.16667 -249.33333 241.83333 [8] -11.83333
> fractions(contrasts(ex01$trat))
contr1 contr2 contr3 contr4 contr5 contr6 contr7 contr8 t1 1/9 1/3 0 0 0 0 0 0 t2 1/9 -1/6 1/2 0 0 0 0 0 t3 1/9 -1/6 -1/2 0 0 0 0 0 t4 -1/18 0 0 1/6 1/2 0 0 0 t5 -1/18 0 0 1/6 -1/2 0 0 0 t6 -1/18 0 0 -1/12 0 1/4 1/2 0 t7 -1/18 0 0 -1/12 0 1/4 -1/2 0 t8 -1/18 0 0 -1/12 0 -1/4 0 1/2 t9 -1/18 0 0 -1/12 0 -1/4 0 -1/2
Uma forma alternativa talvez mais direta e conveniente de obter a decomposição da soma de quadrados de tratamentos entre os contrastes é com o uso de aov() e summary() utilizando o argumento split
> mod3 <- aov(resp ~ trat, data = ex01) > summary(mod3, split=list(trat=1:8))
Df Sum Sq Mean Sq F value Pr(>F) trat 8 332918 41615 21.475 5.44e-13 trat: C1 1 27616 27616 14.251 0.000466 trat: C2 1 420 420 0.217 0.643682 trat: C3 1 21760 21760 11.229 0.001638 trat: C4 1 7041 7041 3.633 0.063029 trat: C5 1 6960 6960 3.592 0.064503 trat: C6 1 93251 93251 48.122 1.26e-08 trat: C7 1 175450 175450 90.541 2.41e-12 trat: C8 1 420 420 0.217 0.643748 Residuals 45 87201 1938
Nota: A atribuição do atributo contrast ao fator não terá efeito sobre a construção da matriz do modelo caso o termo de intercepto esteja retirado na definição do modelo, por exemplo, se o modelo acima fosse definido por resp trat - 1.
Para cancelar a atribuição dos contrastes a um fator e retornar à definida por option() basta fazer atribuir a valor NULL.
> contrasts(ex01$trat) <- NULL
Interpretando os contrastes: a interpretação dos contrastes vai depender dos valores dos coeficientes definidos na matriz dos contrastes (matriz c1 no exemplo). Por exemplo suponha que desejamos obter contrastes de forma não só a obter testes de significância mas também obter estimativas dos contrastes com interpretação desejada. No exemplo, para comparar médias \(\mu _i\) (ou efeitos \(t_i\)) dos tratamentos com efeitos estimados dos contrastes na mesma escala das médias dos tratamentos poderíamos definir c1 da forma a seguir. Todos os resultados e testes de significância (estatísticas de testes e p-valores) permanecem inalterados e apenas os valores dos efeitos estimados são diferentes
> c1 <- rbind(c(1/3,1/3,1/3,-1/6,-1/6,-1/6,-1/6,-1/6,-1/6), + c(1,-0.5,-0.5, rep(0,6)), + c(0, 1,-1, rep(0,6)), + c(rep(0,3), c(0.5,0.5,-0.25,-0.25,-0.25,-0.25)), + c(rep(0,3), c(1,-1), rep(0,4)), + c(rep(0,5), c(0.5,0.5,-0.5,-0.5)), + c(rep(0,5), c(1,-1, 0,0)), + c(rep(0,5), c(0,0, 1,-1))) > fractions(c1)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [1,] 1/3 1/3 1/3 -1/6 -1/6 -1/6 -1/6 -1/6 -1/6 [2,] 1 -1/2 -1/2 0 0 0 0 0 0 [3,] 0 1 -1 0 0 0 0 0 0 [4,] 0 0 0 1/2 1/2 -1/4 -1/4 -1/4 -1/4 [5,] 0 0 0 1 -1 0 0 0 0 [6,] 0 0 0 0 0 1/2 1/2 -1/2 -1/2 [7,] 0 0 0 0 0 1 -1 0 0 [8,] 0 0 0 0 0 0 0 1 -1
> c1.ginv <- ginv(c1) > colnames(c1.ginv) <- rownames(c1) <- paste("contr", 1:8, sep="") > contrasts(ex01$trat) <- c1.ginv > mod4 <- lm(resp~trat, data=ex01) > summary(mod4)
Call: lm(formula = resp ~ trat, data = ex01) Residuals: Min 1Q Median 3Q Max -85.67 -33.29 4.75 33.17 85.67 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 353.52 5.99 59.014 < 2e-16 tratcontr1 47.97 12.71 3.775 0.000466 tratcontr2 -10.25 22.01 -0.466 0.643682 tratcontr3 85.17 25.41 3.351 0.001638 tratcontr4 -29.67 15.56 -1.906 0.063029 tratcontr5 -48.17 25.41 -1.895 0.064503 tratcontr6 -124.67 17.97 -6.937 1.26e-08 tratcontr7 241.83 25.41 9.515 2.41e-12 tratcontr8 -11.83 25.41 -0.466 0.643748 Residual standard error: 44.02 on 45 degrees of freedom Multiple R-squared: 0.7924, Adjusted R-squared: 0.7555 F-statistic: 21.48 on 8 and 45 DF, p-value: 5.445e-13
> coef(mod4)
(Intercept) tratcontr1 tratcontr2 tratcontr3 tratcontr4 tratcontr5 tratcontr6 353.51852 47.97222 -10.25000 85.16667 -29.66667 -48.16667 -124.66667 tratcontr7 tratcontr8 241.83333 -11.83333
Definindo contrastes com funções do pacote gmodels Este pacote oferece várias funções auxiliares que facilitam o uso e a apresentação dos resultados em formatos por vezes convenientes. No exemplo a seguir vamos repetir a obtenção dos resultados para os contrastes definidos na matriz c1 conforme definida anteriormente.
> require(gmodels)
Carregando pacotes exigidos: gmodels
> mod5 <- aov(resp~trat, data=ex01, contrast=list(trat=make.contrasts(c1))) > summary(mod5, split=list(trat=1:8))
Df Sum Sq Mean Sq F value Pr(>F) trat 8 332918 41615 21.475 5.44e-13 trat: C1 1 27616 27616 14.251 0.000466 trat: C2 1 420 420 0.217 0.643682 trat: C3 1 21760 21760 11.229 0.001638 trat: C4 1 7041 7041 3.633 0.063029 trat: C5 1 6960 6960 3.592 0.064503 trat: C6 1 93251 93251 48.122 1.26e-08 trat: C7 1 175450 175450 90.541 2.41e-12 trat: C8 1 420 420 0.217 0.643748 Residuals 45 87201 1938
> fit.contrast(mod5, "trat", c1)
Estimate Std. Error t value Pr(>|t|) tratcontr1 47.97222 12.70763 3.7750713 4.660336e-04 tratcontr2 -10.25000 22.01027 -0.4656918 6.436823e-01 tratcontr3 85.16667 25.41527 3.3510042 1.638352e-03 tratcontr4 -29.66667 15.56361 -1.9061560 6.302888e-02 tratcontr5 -48.16667 25.41527 -1.8951863 6.450255e-02 tratcontr6 -124.66667 17.97131 -6.9369836 1.263698e-08 tratcontr7 241.83333 25.41527 9.5152781 2.408674e-12 tratcontr8 -11.83333 25.41527 -0.4655994 6.437479e-01
Experimentos desbalanceados: finalmente vale ressaltar que o exemplo acima tratou de um experimento balanceado, isto é, com o mesmo número de repetições para cada tratamento e no caso de desbalanceamento ajustes são necessários na definição dos contrastes.
Na sessão anterior discutimos a comparação post-hoc de tratamentos utilizando funções como pairwise.t.test() e TukeyHSD implementadas no conjunto de pacotes básicos do R.
Outros procedimentos são implementados em pacotes contribuídos do R. Entre estes encontra-se os pacotes multcomp e multcompView que implementam diversos outros procedimentos e gráficos para visualizações dos resultados. Vale notar que estes pacotes devem ser instalados com a opção dependencies=TRUE para garantir plena funcionalidade pois suas funções dependem de diversos outros pacotes.
> install.packages("multcompView", dep=TRUE) > require(multcomp) > require(multcompView)
Carregando pacotes exigidos: multcomp
Carregando pacotes exigidos: mvtnorm
Carregando pacotes exigidos: survival
Carregando pacotes exigidos: TH.data
Anexando pacote: 'TH.data'
O seguinte objeto é mascarado por 'package:MASS': geyser
Carregando pacotes exigidos: multcompView
Para ilustrar o uso desta pacote vamos efetuar novamente o teste de Tukey visto acima porém agora utilizando cálculos e gráficos gerados por funções destes pacotes, cujos resultados, embora iguais, são apresentados em forma diferente do visto anteriormente. A indicação de letras para diferenças entre pares de tratamentos mostrada a seguir requer que TukeyHSD seja invocada sem a ordenação dos tratamentos e uma representação visual é dada na Figura 61.
> multcompLetters(TukeyHSD(ex01.mod)$trat[,4])
t2 t3 t4 t5 t6 t7 t8 t9 t1 "a" "bc" "b" "bc" "ac" "d" "ac" "ac" "ac"
> multcompBoxplot(resp~trat,data=ex01,compFn="TukeyHSD",decreasing=FALSE)
No caso de variâncias não homogêneas em experimentos inteiramente casualizados a função oneway.test() pode ser utilizada nas análises. Uma outra alternativa é a análise não paramétrica da Kruskal-Wallis implementada por kruskal.test().