Tutorial ANCOVA
October 23, 2019
Temos um problema
Quando se trabalha como ecologo uma das tarefas possíveis é a avaliação de impacto ambiental.Mas como fazer essa avaliação? Um caso que tem chamado atenção é o derramamento de petróleo crú nas praias do nordeste brasileiro. Muito se tem falado sobre o possível impacto desse evento na fauna e flora marinha da região. Os impactos econômicos normalmente se medem pela perda financeira evidente das famílias que subsistem da pesca e coleta de mariscos na região. O dano social vem atrelado ao dano econômico e cultural caso não haja um plano de contigência e restabelecimento das áreas afetadas. Mas a pergunta é e um ecólogo como mede o impacto ambiental?
Uma pergunta e uma hipótese
Digamos então que você tem como tarefa avaliar o impacto da contaminação por óleo cru nas praias do nordeste brasileiro na reprodução das tartarugas da espécie Caretta caretta. A primeira pergunta seria: qual o efeito do óleo na reprodução? Para formular uma hipótese biológica recorremos ao conhecimento que temos sobre o ecossistema marinho e a própria espécie analisada. Ainda recorremos a estudos prévios sobre contaminação com óleo, já que, infelizmente, não é a primeira vez que óleo é derramado no mar (embora seja o maior evento ocorrido no Brasil). Sabemos que a presença do óleo no mar diminui a produção de alga,um vez que o óleo cria uma barreira para a passagem da luz solar, diminuindo assim a taxa fotossintética. Mas e as tartarugas? As tartarugas se alimentam das algas e com menor disponibilidade de alimento acabam entrando em um “regime”, o corpo então passa a economizar energia para garantir a sobrevivência do animal e umas das áreas de corte é a de produção de ovos.
Dessa forma, nossa hipótese é que o número de ovos por ninho de tartaruga seja menor em praias contaminadas por petróleo do que em praias não contaminadas.
Já tenho uma pergunta e uma hipótese, e agora?
O próximo passo é pensar sobre qual métrica e quais dados usar para responder nossa pergunta. Bem, queremos investigar se existe diferença entre a quantidade de ovos colocadas por fêmeas que se alimentaram em regiões expostas à contaminação por petróleo. Para saber o efeito da contaminação temos que ter dados também de quando não há contaminação. Só assim podemos comparar e saber se há uma diferença entre as duas condições. Mas não basta ter apenas uma diferença, ela tem que ser estatísticamente significante, ou seja, tem que ser uma diferença maior do que a encontrada ao acaso.
Como saber se o efeito visto é real?
No mundo biológico sempre temos muitas possibilidades que explicam um determinado comportamento observado, por exemplo: Lendo sobre a biologia das tartarugas, sabemos que o tamanho corporal da fêmea também influência a quantidade de ovos colocados por ela. Logo, só pela diferença biológica uma fêmea pode colocar mais ou menos ovos e isso não ter nada haver com a presença ou ausência de óleo. Portanto, temos que pensar quais variáveis podem covariar com nossos dados e levar elas em consideração em nossas análises estatísticas.
Vamos ver como ficou nossa tabela de campo:
data <- read.csv("Datos.csv", stringsAsFactors = F)
head(data)## Sample Tratamento Ovos Tamanho
## 1 1 SEM 95 92.74622
## 2 2 SEM 80 78.80903
## 3 3 SEM 100 82.95129
## 4 4 SEM 94 109.94159
## 5 5 SEM 94 83.35462
## 6 6 SEM 69 92.43718Carregando os pacotes para data wrangling:
O que temos até agora?
Temos então a variável da presença ou ausência de óleo, o tamanho de corpo e ambas as variáveis explicando a quantidade de ovos. A presença ou ausência de óleo é o que chamamos de variável categórica, ou seja, temos dois estados possíveis. Podemos dizer assim que é uma variável categórica de dois níveis. Já tamanho de corpo temos um contínuo de possíveis tamanhos, assim como quantidade de ovos temos um contínuo de possibilidade de números de ovos, logo, ambas são variáveis contínuas.
Qual teste estatístico utilizar?
Já temos nossa pergunta, hipótese, variáveis, covariáveis e sua natureza (contínua ou categórica). A partir daqui, temos como olhar na literatura e buscar o teste estatístico que melhor se encaixa. Caso não haja, ainda é possível criar seu próprio teste estatístico, uma vez que sabemos que queremos medir o efeito do óleo apenas sobre a produção de ovos e eliminar de nossos dados o efeito da covariável tamanho do corpo. Para cada teste estatístico temos ainda os pressupostos que nossos dados precisam atender para se encaixar no teste.
Análise de Covariancia (ANCOVA)
Como queremos nos livrar de efeitos de uma covariável a análise de covariância (ANCOVA) é a mais indicada, já que é usada para comparar duas ou mais retas de regressão linear, e investigar se existe diferença entre elas. As retas são comparadas estudando a interação de uma variável categórica (presença de petróleo) com a variável dependente (quantidade de ovos), enquanto controla o efeito de uma covariável contínua (o tamanho do corpo da tartaruga). Na ANCOVA podemos também considerar a interação entre a variável categórica (presença de petróleo) e a variável independente (tamanho do corpo da tartaruga) e isso é determinado pela fórrmula utilizada no cálculo.
Pressupostos do teste
A ANCOVA combina os testes de regressão linear e ANOVA (Análise de Variância) e por isso os dados utilizados no teste de ANCOVA deve obedecer todos os pressupostos desses testes.
De maneira geral, os pressupostos são:
1. A distribuíção dos dados contínuos deve ser normal
Podemos usar um histograma para analisar visualmente a distribuíção dos valores para quantidade de ovos e tamanho do corpo. Podemos também fazer um teste de normalidade, chamado de Shapiro Teste que, de maneira geral, vai nos dar o p-valor para ver se a distribuição não é significativamente diferente da distribuição normal. A hipótese nula do Shapiro Test é que a distribuição dos dados é normal. Portanto, quando o p-value for maior que 0.05 indica que a hipótese nula foi aceita e que os dados possuem distribuíção normal.
Para a distribuíção de quantidade de ovos, considerando os níveis da variável categórica temos:
data_sem <- filter(data, data$Tratamento == "SEM")
hist(data_sem$Ovos, main = "Quantidade de ovos sem tratamento")
shapiro.test(data_sem$Ovos)##
## Shapiro-Wilk normality test
##
## data: data_sem$Ovos
## W = 0.99719, p-value = 0.555data_com <- filter(data, data$Tratamento == "COM")
hist(data_com$Ovos, main = "Quantidade de ovos com tratamento")
shapiro.test(data_com$Ovos)##
## Shapiro-Wilk normality test
##
## data: data_com$Ovos
## W = 0.99659, p-value = 0.3683Para a distribuição de tamanho de corpo temos:
hist(data_sem$Tamanho, main = "Tamanho do corpo (cm) sem tratamento")
shapiro.test(data_sem$Tamanho)##
## Shapiro-Wilk normality test
##
## data: data_sem$Tamanho
## W = 0.99456, p-value = 0.07405hist(data_com$Tamanho, main = "Tamanho do corpo (cm) com tratamento")
shapiro.test(data_com$Tamanho)##
## Shapiro-Wilk normality test
##
## data: data_com$Tamanho
## W = 0.9975, p-value = 0.66112. Homogeneidade da Variância entre as categorias da variável categórica
Podemos fazer um boxplot e ver se a variância dos valores assumidos por cada variável contínua entre os níveis da variável categórica é homogênea, ou podemos fazer um teste de homogeneidade. Nesse caso usaremos o Barlett test, que é adequado apenas para dados contínuos que seguem a distribuíção normal. Como já fizemos o Shapiro test que confirma a normalidade dos dados, podemos prosseguir com o Barlett test. Para saber se a variância é homogênea, o p-valor tem que ser maior que 0.05, para que a hipótese nula de homogeneidade seja aceita.
data$Tratamento <- as.factor(data$Tratamento)
plot(data$Ovos ~ data$Tratamento)
bartlett.test(data$Ovos ~ data$Tratamento)##
## Bartlett test of homogeneity of variances
##
## data: data$Ovos by data$Tratamento
## Bartlett's K-squared = 41.541, df = 1, p-value = 1.154e-10data$Tratamento <- as.factor(data$Tratamento)
plot(data$Tamanho ~ data$Tratamento)
bartlett.test(data$Tamanho ~ data$Tratamento)##
## Bartlett test of homogeneity of variances
##
## data: data$Tamanho by data$Tratamento
## Bartlett's K-squared = 1.039, df = 1, p-value = 0.30813. A relação entre as variáveis contínuas deve obedecer o comportamento linear
plot(data_sem$Ovos ~ data_sem$Tamanho,
main = "Qtide de ovos x tamanho (cm) sem tratamento")
sem_t <- lm(data_sem$Ovos ~ data_sem$Tamanho)
summary(sem_t)##
## Call:
## lm(formula = data_sem$Ovos ~ data_sem$Tamanho)
##
## Residuals:
## Min 1Q Median 3Q Max
## -28.2926 -6.4356 -0.0608 6.2103 26.8137
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 96.3960 3.4461 27.973 <2e-16 ***
## data_sem$Tamanho -0.0111 0.0359 -0.309 0.757
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.705 on 497 degrees of freedom
## Multiple R-squared: 0.0001924, Adjusted R-squared: -0.001819
## F-statistic: 0.09567 on 1 and 497 DF, p-value: 0.7572plot(data_com$Ovos ~ data_com$Tamanho,
main = "Qtide de ovos x tamanho (cm) com tratamento")
com_t <- lm(data_com$Ovos ~ data_com$Tamanho)
summary(com_t)##
## Call:
## lm(formula = data_com$Ovos ~ data_com$Tamanho)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.3020 -5.0825 -0.2113 4.8823 23.8941
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 43.20055 2.43608 17.734 <2e-16 ***
## data_com$Tamanho -0.01055 0.02559 -0.412 0.68
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.257 on 499 degrees of freedom
## Multiple R-squared: 0.0003402, Adjusted R-squared: -0.001663
## F-statistic: 0.1698 on 1 and 499 DF, p-value: 0.68044. Homogeneidade na distribuição de resíduos da Regressão Linear
resid_sem_t <- resid(sem_t)
plot(resid_sem_t ~ data_sem$Tamanho,
ylab="Resíduos", xlab="Quantidade de ovos",
main="Tamanho do corpo (cm) sem tratamento")
resid_com_t <- resid(com_t)
plot(resid_com_t ~ data_com$Tamanho,
ylab="Resíduos", xlab="Quantidade de ovos",
main="Tamanho do corpo (cm) com tratamento")
Para mais informações consulte os tutoriais de ANOVA e Regressão Linear.
Ancova no R
Uma vez que nossos dados passaram em todos os testes acima, podemos então fazer a Ancova. A fórmula usada para fazer a ANCOVA no R é simples e combina elementos da regressão linear em um teste de variância, como descrito abaixo:
aov(variavel_y ~ variavel_x * variavel_categorica)
No nosso caso, a variável_y será tamanho do corpo, variavel_x será número de ovos e a variavel_categórica será presença de petróleo Podemos também decidir não considerar o efeito da interação da variável categórica na variavel x, e, nesse caso, substituímos o “*" pelo “+”, como descrito abaixo:
aov(variavel_y ~ variavel_x + variavel_categorica)
Interpretação dos resultados
Com interação
aov(variavel_y ~ variavel_x * variavel_categorica)
summary(aov(data$Ovos ~ data$Tamanho * data$Tratamento))## Df Sum Sq Mean Sq F value Pr(>F)
## data$Tamanho 1 673 673 9.166 0.00253 **
## data$Tratamento 1 705124 705124 9608.531 < 2e-16 ***
## data$Tamanho:data$Tratamento 1 0 0 0.000 0.98986
## Residuals 996 73092 73
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Sem interação
aov(variavel_y ~ variavel_x + variavel_categorica)
summary(aov(data$Ovos ~ data$Tamanho + data$Tratamento))## Df Sum Sq Mean Sq F value Pr(>F)
## data$Tamanho 1 673 673 9.176 0.00252 **
## data$Tratamento 1 705124 705124 9618.176 < 2e-16 ***
## Residuals 997 73092 73
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpretando a tabela da ANCOVA
O teste de ANCOVA imprime como resultado uma tabela. A primeira coisa que a ANCOVA testa é a relação entre as duas variáveis contínuas para cada nível da variável categórica. Em outras palavras, ela testa se os coeficientes das regressões lineares entre Tamanho do Corpo da tartaruga (x) e a Quantidade de ovos (y) é igual para a situação Com petróleo e Sem petróleo (tratamento). Podemos olhar diretamente para a coluna Pr(>F) para avaliar o nível de significância entre o Tamanho do Corpo da tartaruga (y) e a Quantidade de ovos (x), e vemos que o p-valor é menor que 0.05 para a variável contínua “Tamanho”. Dessa forma, observamos que há efeito do tamanho do corpo na quantidade de ovos colocados pela tartaruga. Logo abaixo, vemos também que tem efeito do ‘Tratamento’ sobre a quatidade de ovos produzida, onde o p-valor é menor que 0.05. Podemos também acessar o efeito do tratamento sobre o tamanho do corpo e seu reflexo na quantidade de ovos colocados pela fêmea da tartaruga. Isso é avaliado pelo teste somente se usarmos um “*" e não um “+” na fórmula. No caso em que colocamos o asterisco, vemos que não existe efeito significativo da presença do petróleo no tamanho da tartaruga.
Dessa forma, podemos concluir que existe relação significativa entre o tamanho do corpo e a quantidade de ovos, e que essa relação é afetada pela presença de petróleo.
data$Ovos <- as.numeric(data$Ovos)
plot(data$Ovos ~ data$Tamanho, type='n')
points(data_sem$Tamanho, data_sem$Ovos, pch=20, col="blue")
points(data_com$Tamanho, data_com$Ovos, pch=20, col='red')
abline(sem_t, lty=2, col='blue')
abline(com_t, lty=2, col='red')
legend("bottomright", c("Sem Petróleo","Com Petróleo"), lty=c(2,2),
pch=c(20,1), col=c('blue', 'red'))