Teste T & ANOVA: Aplicações, gráficos e tabelas
Renata Brandt e Fábio Barros
July 16, 2015
Uma vez que aprendemos a construir/manipular os objetos, abrir arquivos em R (do tipo .xls, .txt, entre outros) e instalar pacotes, podemos pensar em analisar os dados biológicos que foram obtidos durante os experimentos e/ou observações em campo.
A idéia do curso é, gradativamente (em termos da complexidade do modelo), mostrar as análises mais adequadas para determinado tipo de pergunta e como mostrar os resultados graficamente de modo satisfatório para publicação.
1. Teste T
Para iniciarmos, suponhamos que a pergunta do meu projeto esteja relacionado ao efeito da idade sobre a produção de frutos em laranjais.
Neste exemplo, minha hipótese seria a de que indivíduos mais novos produzem proporcionalmente MAIS frutos do que indivíduos mais velhos (medido em kg de frutos/biomassa total de cada árvore).
Durante a tomada dos dados biológicos, eu dividi os indivíduos de determinada área em dois grupos: Grupo “NEW” = Laranjeiras de até 5 anos; Grupo “OLD” = Laranjeiras de mais de 5 anos; foram medido 50 indivíduos de cada grupo.
O teste mais adequado para testar minha hipótese seria um Teste T simples pois eu possuo dois grupos categorizados (i.e. Grupo NEW e Grupo OLD). A aplicação do teste “T” de Student é uma função que NÃO requer instalação de novos pacotes, já que está contido no pacote “stats”, dentro de uma série de pacotes já pré-instalados em R.
No sentido de tratarmos somente da análise, eu simulei previamente um conjunto de dados com valores para o Grupo NEW e Grupo OLD e salvei em .csv.
Simulando Dados
rnorm é a função que simula variáveis contínuas normais. A forma de usar é: x <- rnorm(n, mean, sd) em que n é o tamanho da amostra (x será um vetor caso n > 1), e mean e sd são parâmetros (opcionais) dando a média e o desvio-padrão da normal. Se mean ou sd forem omitidos, serão usados os valores, respectivamente, de 0 e 1.
Por exemplo: x <- rnorm(sd = 2, n = 10) gera um vetor com 10 valores independentes e identicamente distribuídos, com média zero e desvio-padrão 2.
No caso dos dados biológicos simulados referentes ao exemplo 1:
NEW <- rnorm(50, 200, 45); hist(NEW)
OLD <- rnorm(50, 125, 30); hist(OLD)
orange <- data.frame(cbind(NEW, OLD))
write.csv (orange, “C:/Users/F?bio Barros/Desktop/CursoR/orange.csv”, row.names=F)
Validando o Teste T
Passos:
- abrir a planilha com os dados:
# definindo o diretório de análise
# cada usuário terá o seu
# setwd("C:/Users/F?bio Barros/Desktop/CursoR")
# para abrir e nomear o conjunto de dados
# laranja<-read.csv('orange.csv')
# Para facilitar, subimos o arquivo na internet e você pode abrir direto do link:
laranja<-read.csv('http://renatabrandt.github.io/EBC2015/data/orange.csv')- executar a função para Teste T de Student:
t.laranja<-t.test(laranja)
t.laranja##
## One Sample t-test
##
## data: laranja
## t = 29.02, df = 99, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 151.2869 173.4932
## sample estimates:
## mean of x
## 162.3901- visualizar graficamente os resultados:
boxplot(laranja)
Como diria Milton Leite: Que beleza!!! Já temos as análises para nosso primeiro problema biológico… Realmente, as laranjeiras mais novas (com até 05 anos), produzem proporcionalmente mais frutos do que árvores mais velhas!
Mas esse gráfico de caixa está meio caidinho hein?! Não está auto-explicativo como uma figura deve ser! Vamos dar um upgrade no gráfico? Para isto usaremos a função ggplot, do pacote ggplot2.
Criando gráficos
A vantagem de usar o pacote ggplot2 para produzir gráficos é que a sintaxe é praticamente a mesma em qualquer gráfico que se queria produzir.
Basicamente, para produzir um gráfico em ggplot2, devemos definir qual o conjunto de dados que será utilizado na representação e ainda o tipo de gráfico (se barras, caixa, pontos, entre outros). Depois, pode-se definir os parâmetros secundários.
Lembrando do primeiro bloco, instalaremos o pacote ggplot2, e em seguida produziremos o gráfico que representa a diferença entre os laranjeiras em termos de produção de frutos:
install.packages(“ggplot2”)
require(ggplot2)## Loading required package: ggplot2Para produzir um gráfico no ggplot, devemos reorganizar nosso conjunto de dados de forma que a função consiga “ler” a tabela. Faremos isso usando a função stack():
# função stack reorganiza os dados no data.frame
laranja <- stack(laranja)
# renomeando as colunas
names(laranja) <- c("produto","idade")
# reorganizando as colunas
laranja <- laranja[c("idade", "produto")]E agora vamos de fato produzir o gráfico:
orange.plot <- ggplot(laranja, aes(x=idade, y=produto)) +
geom_boxplot()
orange.plot
# Podemos melhorar o default gráfico preenchendo cada grupo de
# laranjeiras por cores e mudando a largura do gráfico
orange.plot <- ggplot(laranja, aes(x=idade, y=produto,
fill=idade)) +
geom_boxplot(outlier.size=3,outlier.shape=21, width=0.5)
orange.plot
Pronto!!! Temos um gráfico produzino com ggplot2. Para customizá-lo existem muitas opções (muitas mesmo!!!). Como guia básico, nós aconselhamos que os usuários utilizem o R Graphics Cookbook
Deixamos esta opção aqui, inclusive com a utilização de fontes alternativas!
install.packages(“extrafont”)
# Para acessar tipos de letras alternativas
# loadfonts(device="win")
require(extrafont)## Loading required package: extrafont
## Registering fonts with R# Aqui definimos um objeto que representa a área de um gráfico em BRANCO.
# Eu incorporei este objeto no gráfico que vem a seguir
# para retirar as linhas de grade
no.grid <- element_blank()
orange.plot.full <- orange.plot +
scale_fill_manual(values=c("lightgreen",
"navajowhite3")) +
# otimizando cor das variaveis
theme_bw(base_size = 12, base_family="Times") +
# definindo tamanho e características gerais
# da fonte
xlab("Idade das laranjeiras") +
theme(axis.title.x=element_text(angle=0, size=14)) +
theme(axis.title.x=element_text(vjust=-0.50)) +
# ajustando Título no eixo X, e distância do
# texto em relação ao gráfico
ylab("Kg de fruto por biomassa") +
theme(axis.title.y=element_text(angle=90, size=14)) +
theme(axis.title.y=element_text(vjust=1)) +
# o mesmo para o eixo y
theme(panel.grid.minor=no.grid,
panel.grid.major=no.grid) +
# retirando linhas de grade
guides(fill=FALSE)
# retirando as legendas
orange.plot.full
Salvando imagens
O melhor modo de salvar os gráficos para publicação é no formato .pdf ou Postscript (.eps). Estes formatos não alteram a qualidade da figura, mesmo quando aumentadas muitas vezes.
Como fazer? Para gráficos de ggplot2 existe a possibilidade de usar a função ggsave(), que permite customização de formato, tamanho… Para quem usa RStudio, na janela inferior direita clique em “Export -> Save as PDF… -> SAVE”.
2. ANOVA
A ANOVA (Analysis of Variance) permite testar se determinados conjuntos de dados possuem médias iguais ou não, levando-se em conta a variação dos números em torno da média. Este método, diferentemente do Teste T de Student permite que vários grupos sejam comparados ao mesmo tempo. Em outras palavras, permite testes com um fator que tenha mais de 2 níveis.
Existem muitos pacotes para realizar testes de ANOVA (i.e. ez, car, entre outros), mas a base do R também já contém essa função com uma sintaxe bem simples…
aov(RESPOSTA ~ FATOR), onde assumimos que a variável RESPOSTA tem a condição de dependencia em relação ao FATOR.
Num exemplo hipotético utilizando o conjunto de dados IRIS disponível no R, poderíamos fazer desta maneira:
# para carregar o data.frame
iris## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
## 7 4.6 3.4 1.4 0.3 setosa
## 8 5.0 3.4 1.5 0.2 setosa
## 9 4.4 2.9 1.4 0.2 setosa
## 10 4.9 3.1 1.5 0.1 setosa
## 11 5.4 3.7 1.5 0.2 setosa
## 12 4.8 3.4 1.6 0.2 setosa
## 13 4.8 3.0 1.4 0.1 setosa
## 14 4.3 3.0 1.1 0.1 setosa
## 15 5.8 4.0 1.2 0.2 setosa
## 16 5.7 4.4 1.5 0.4 setosa
## 17 5.4 3.9 1.3 0.4 setosa
## 18 5.1 3.5 1.4 0.3 setosa
## 19 5.7 3.8 1.7 0.3 setosa
## 20 5.1 3.8 1.5 0.3 setosa
## 21 5.4 3.4 1.7 0.2 setosa
## 22 5.1 3.7 1.5 0.4 setosa
## 23 4.6 3.6 1.0 0.2 setosa
## 24 5.1 3.3 1.7 0.5 setosa
## 25 4.8 3.4 1.9 0.2 setosa
## 26 5.0 3.0 1.6 0.2 setosa
## 27 5.0 3.4 1.6 0.4 setosa
## 28 5.2 3.5 1.5 0.2 setosa
## 29 5.2 3.4 1.4 0.2 setosa
## 30 4.7 3.2 1.6 0.2 setosa
## 31 4.8 3.1 1.6 0.2 setosa
## 32 5.4 3.4 1.5 0.4 setosa
## 33 5.2 4.1 1.5 0.1 setosa
## 34 5.5 4.2 1.4 0.2 setosa
## 35 4.9 3.1 1.5 0.2 setosa
## 36 5.0 3.2 1.2 0.2 setosa
## 37 5.5 3.5 1.3 0.2 setosa
## 38 4.9 3.6 1.4 0.1 setosa
## 39 4.4 3.0 1.3 0.2 setosa
## 40 5.1 3.4 1.5 0.2 setosa
## 41 5.0 3.5 1.3 0.3 setosa
## 42 4.5 2.3 1.3 0.3 setosa
## 43 4.4 3.2 1.3 0.2 setosa
## 44 5.0 3.5 1.6 0.6 setosa
## 45 5.1 3.8 1.9 0.4 setosa
## 46 4.8 3.0 1.4 0.3 setosa
## 47 5.1 3.8 1.6 0.2 setosa
## 48 4.6 3.2 1.4 0.2 setosa
## 49 5.3 3.7 1.5 0.2 setosa
## 50 5.0 3.3 1.4 0.2 setosa
## 51 7.0 3.2 4.7 1.4 versicolor
## 52 6.4 3.2 4.5 1.5 versicolor
## 53 6.9 3.1 4.9 1.5 versicolor
## 54 5.5 2.3 4.0 1.3 versicolor
## 55 6.5 2.8 4.6 1.5 versicolor
## 56 5.7 2.8 4.5 1.3 versicolor
## 57 6.3 3.3 4.7 1.6 versicolor
## 58 4.9 2.4 3.3 1.0 versicolor
## 59 6.6 2.9 4.6 1.3 versicolor
## 60 5.2 2.7 3.9 1.4 versicolor
## 61 5.0 2.0 3.5 1.0 versicolor
## 62 5.9 3.0 4.2 1.5 versicolor
## 63 6.0 2.2 4.0 1.0 versicolor
## 64 6.1 2.9 4.7 1.4 versicolor
## 65 5.6 2.9 3.6 1.3 versicolor
## 66 6.7 3.1 4.4 1.4 versicolor
## 67 5.6 3.0 4.5 1.5 versicolor
## 68 5.8 2.7 4.1 1.0 versicolor
## 69 6.2 2.2 4.5 1.5 versicolor
## 70 5.6 2.5 3.9 1.1 versicolor
## 71 5.9 3.2 4.8 1.8 versicolor
## 72 6.1 2.8 4.0 1.3 versicolor
## 73 6.3 2.5 4.9 1.5 versicolor
## 74 6.1 2.8 4.7 1.2 versicolor
## 75 6.4 2.9 4.3 1.3 versicolor
## 76 6.6 3.0 4.4 1.4 versicolor
## 77 6.8 2.8 4.8 1.4 versicolor
## 78 6.7 3.0 5.0 1.7 versicolor
## 79 6.0 2.9 4.5 1.5 versicolor
## 80 5.7 2.6 3.5 1.0 versicolor
## 81 5.5 2.4 3.8 1.1 versicolor
## 82 5.5 2.4 3.7 1.0 versicolor
## 83 5.8 2.7 3.9 1.2 versicolor
## 84 6.0 2.7 5.1 1.6 versicolor
## 85 5.4 3.0 4.5 1.5 versicolor
## 86 6.0 3.4 4.5 1.6 versicolor
## 87 6.7 3.1 4.7 1.5 versicolor
## 88 6.3 2.3 4.4 1.3 versicolor
## 89 5.6 3.0 4.1 1.3 versicolor
## 90 5.5 2.5 4.0 1.3 versicolor
## 91 5.5 2.6 4.4 1.2 versicolor
## 92 6.1 3.0 4.6 1.4 versicolor
## 93 5.8 2.6 4.0 1.2 versicolor
## 94 5.0 2.3 3.3 1.0 versicolor
## 95 5.6 2.7 4.2 1.3 versicolor
## 96 5.7 3.0 4.2 1.2 versicolor
## 97 5.7 2.9 4.2 1.3 versicolor
## 98 6.2 2.9 4.3 1.3 versicolor
## 99 5.1 2.5 3.0 1.1 versicolor
## 100 5.7 2.8 4.1 1.3 versicolor
## 101 6.3 3.3 6.0 2.5 virginica
## 102 5.8 2.7 5.1 1.9 virginica
## 103 7.1 3.0 5.9 2.1 virginica
## 104 6.3 2.9 5.6 1.8 virginica
## 105 6.5 3.0 5.8 2.2 virginica
## 106 7.6 3.0 6.6 2.1 virginica
## 107 4.9 2.5 4.5 1.7 virginica
## 108 7.3 2.9 6.3 1.8 virginica
## 109 6.7 2.5 5.8 1.8 virginica
## 110 7.2 3.6 6.1 2.5 virginica
## 111 6.5 3.2 5.1 2.0 virginica
## 112 6.4 2.7 5.3 1.9 virginica
## 113 6.8 3.0 5.5 2.1 virginica
## 114 5.7 2.5 5.0 2.0 virginica
## 115 5.8 2.8 5.1 2.4 virginica
## 116 6.4 3.2 5.3 2.3 virginica
## 117 6.5 3.0 5.5 1.8 virginica
## 118 7.7 3.8 6.7 2.2 virginica
## 119 7.7 2.6 6.9 2.3 virginica
## 120 6.0 2.2 5.0 1.5 virginica
## 121 6.9 3.2 5.7 2.3 virginica
## 122 5.6 2.8 4.9 2.0 virginica
## 123 7.7 2.8 6.7 2.0 virginica
## 124 6.3 2.7 4.9 1.8 virginica
## 125 6.7 3.3 5.7 2.1 virginica
## 126 7.2 3.2 6.0 1.8 virginica
## 127 6.2 2.8 4.8 1.8 virginica
## 128 6.1 3.0 4.9 1.8 virginica
## 129 6.4 2.8 5.6 2.1 virginica
## 130 7.2 3.0 5.8 1.6 virginica
## 131 7.4 2.8 6.1 1.9 virginica
## 132 7.9 3.8 6.4 2.0 virginica
## 133 6.4 2.8 5.6 2.2 virginica
## 134 6.3 2.8 5.1 1.5 virginica
## 135 6.1 2.6 5.6 1.4 virginica
## 136 7.7 3.0 6.1 2.3 virginica
## 137 6.3 3.4 5.6 2.4 virginica
## 138 6.4 3.1 5.5 1.8 virginica
## 139 6.0 3.0 4.8 1.8 virginica
## 140 6.9 3.1 5.4 2.1 virginica
## 141 6.7 3.1 5.6 2.4 virginica
## 142 6.9 3.1 5.1 2.3 virginica
## 143 5.8 2.7 5.1 1.9 virginica
## 144 6.8 3.2 5.9 2.3 virginica
## 145 6.7 3.3 5.7 2.5 virginica
## 146 6.7 3.0 5.2 2.3 virginica
## 147 6.3 2.5 5.0 1.9 virginica
## 148 6.5 3.0 5.2 2.0 virginica
## 149 6.2 3.4 5.4 2.3 virginica
## 150 5.9 3.0 5.1 1.8 virginica# para permitir o acesso direto às colunas do conjunto de dados
attach(iris)
aov.iris <- aov(Petal.Width ~ Species)
summary (aov.iris)## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 80.41 40.21 960 <2e-16 ***
## Residuals 147 6.16 0.04
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Desafio
Uma vez que já demonstramos como realizar diversas opera??es no R, propomos aqui nosso primeiro DESAFIO do curso:
Vamos pensar novamente na produtividade dos laranjais, que abordamos no item anterior. Digamos que durante observações prévias, percebemos que indivíduos de até 2 anos talvez produzam ainda mais frutos que árvores entre 2 a 5 anos. Neste caso, decidimos por comparar 03 grupos de ?rvores distintas por meio de uma Análise de Variância.
Simule, por meio da função rnorm, 3 vetores referentes aos grupos de laranjeiras teoricamente amostradas pelo pesquisador. O Grupo 1 (de até 2 anos) deve conter 100 amostras, com média 373, e variância de 30; Grupo 2 (de 2 a 5 anos) - 75 amostras, média de 157, e variância de 50; Grupo 3 - 90 amostras, média de 145, e variância de 65.
Crie um conjunto de dados (data.frame) usando a função cbind e stack, e altere o nome das colunas usando a função names.
Aplique o teste de ANOVA em seus dados. Em termos dos laranjais, qual é a variável RESPOSTA? E qual é o FATOR (no caso, a variável independente)? Qual foi o resultado obtido?
ANOVA de 2-fatores (Two-way) e Interações
A ANOVA realizada anteriormente levou em consideração apenas 1 fator como influente para a produç˜eo de frutos em laranjeiras, no caso, a idade dos indivíduos.
Dependendo da pergunta que consideramos responder, podemos observar mais elementos potencialmente influentes na RESPOSTA. Por exemplo:
aov(RESPOSTA ~ FatorA + FatorB)
Exemplo prático
Vamos abrir o arquivo orange_iii no qual adicionamos mais um fator ao conjunto de dados. Lembre-se que as variáveis independentes numa ANOVA são categóricas e somente a variável RESPOSTA deve ser contínua.
orange_iii <- read.csv("http://renatabrandt.github.io/EBC2015/data/orange_iii.csv")
str (orange_iii)## 'data.frame': 300 obs. of 3 variables:
## $ age : Factor w/ 3 levels "interm","old",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ altitude : Factor w/ 2 levels "high","low": 2 2 2 2 2 2 2 2 2 2 ...
## $ production: num 446 435 432 430 429 ...# reparem que o fator "AGE" possui 3 níveis
# enquanto o fator "ALTITUDE" possui 2 níveis
attach (orange_iii)Agora vamos testar o modelo tanto para Age quanto para Altitude
aov.orange_iii <- aov(production ~ age+altitude)
anova (aov.orange_iii)## Analysis of Variance Table
##
## Response: production
## Df Sum Sq Mean Sq F value Pr(>F)
## age 2 3411521 1705761 885.168 < 2.2e-16 ***
## altitude 1 74599 74599 38.711 1.674e-09 ***
## Residuals 296 570406 1927
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ainda podemos considerar que existe interação entre os fatores, que refletem alterações nos índices da variável dependente, neste caso:
aov.interaction <- aov(production ~ age + altitude + age:altitude)
anova (aov.interaction)## Analysis of Variance Table
##
## Response: production
## Df Sum Sq Mean Sq F value Pr(>F)
## age 2 3411521 1705761 959.819 < 2.2e-16 ***
## altitude 1 74599 74599 41.976 3.882e-10 ***
## age:altitude 2 47918 23959 13.482 2.501e-06 ***
## Residuals 294 522488 1777
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como diria Serginho Malandro… Yeah! Yeah! Existe um efeito tanto da altitude quanto da idade em relação à produção de frutos por laranjeira!
Mas novamente, o teste falha em demonstrar com mais acurácia como os fatores influenciam a produção de frutos e, principalmente, o efeito das interações entre os fatores! Portanto, nós sugerimos que vocês utilizem um pacote mais adequado, por exemplo o PHIA (Pos-hoc Interaction Analysis)
PHIA (Pos-hoc Interaction Analysis)
O pacote phia é direcionado especificamente para teste de interaões e visualização gráfica. O default do pacote é de fácil compreensão e por isso sugerimos que vocês se aprofundem, se tiverem interesse.
Ele funciona com modelos de base lm, com sintaxe semelhante aos testes de ANOVA que demonstramos anteriormente.
O teste de ANOVA num modelo de base lm requer a instalação do pacote car. Portanto devemos instalá-lo também.
install.packages(“car”)
install.packages(“phia”)
require(car)## Loading required package: carrequire(phia)## Loading required package: phia# O "phia" nos dá opções de fixar um dos fatores
# e somente olhar o efeito no outro fator de interesse
# Ainda podemos plotar as interações, o que é muito interessante!
# Deixamos aqui um exemplo para o conjunto de dados "orange_iii"
phia.orange <- lm(production ~ age*altitude)
Anova(phia.orange, type=3)## Anova Table (Type III tests)
##
## Response: production
## Sum Sq Df F value Pr(>F)
## (Intercept) 1188927 1 669.001 < 2.2e-16 ***
## age 2386697 2 671.488 < 2.2e-16 ***
## altitude 95366 1 53.661 2.304e-12 ***
## age:altitude 47918 2 13.482 2.501e-06 ***
## Residuals 522488 294
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1testInteractions(phia.orange, fixed="age",
across="altitude", adjustment="bonferroni")## F Test:
## P-value adjustment method: bonferroni
## Value Df Sum of Sq F Pr(>F)
## interm -64.745 1 95366 53.6615 6.912e-12 ***
## old -1.720 1 74 0.0416 1.0000000
## young -38.529 1 27077 15.2363 0.0003531 ***
## Residuals 294 522488
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1means.orange<-interactionMeans(phia.orange)
plot(means.orange, abbrev.levels=FALSE)
E aí? O que acharam? Na nossa opinião as relações ficam muito mais claras de observar utilizando phia com poucas linhas de comando.
No exemplo, onde está realmente a diferença? O que ocorre com as laranjeiras mais velhas?