Análise de Variância - ANOVA
Lia Hanna Martins Morita
2020
Pacotes
require(tidyverse)
require(dplyr)
require(DT)
require(kableExtra)
require(captioner)
require(interactions)
require(jtools)
require(car)1 Anova com 1 fator - One way Analise of Variance
dados fictícios: https://www.scribbr.com/statistics/anova-in-r/
dados=read.csv("crop.data.csv")
dados = dados %>%
mutate(fertilizer=as.factor(fertilizer),
density=as.factor(density),
ind=as.factor(seq(1,nrow(dados),1)))Descritivas
descritivas=dados %>%
group_by(fertilizer,density) %>%
summarise(media=mean(yield),sd=sd(yield),n=n())## `summarise()` regrouping output by 'fertilizer' (override with `.groups` argument)tab_nums=captioner(prefix="Tabela")
cap=tab_nums("tab1","Estatíticas descritivas")
datatable(descritivas,caption=cap)One Way Anova: Fertilizante influencia na produção?
- Descritiva considerando apenas fertilizante:
- Pela tabela 2, o fertilizante 3 possui média superior em produção;
- Pelo Gráfico de Boxplot, verificamos que o fertilizante 3 em geral possui valores superiores de produção.
y_grupos=dados %>%
group_by(fertilizer) %>%
summarise(media=mean(yield),var=var(yield),n=n())## `summarise()` ungrouping output (override with `.groups` argument)cap=tab_nums("tab2","Descritivas por tipo de fertilizante")
datatable(y_grupos,cap=cap)dados %>%
ggplot(aes(x=fertilizer,y=yield,fill=fertilizer)) +
geom_boxplot()
Figura 1: Boxplot de produção por tipos de fertilizante
- Pela tabela ANOVA verificamos que o fertilizante é um fator que influencia na produção
one.way <- aov(yield ~ fertilizer, data = dados)
summary(one.way)## Df Sum Sq Mean Sq F value Pr(>F)
## fertilizer 2 6.07 3.0340 7.863 7e-04 ***
## Residuals 93 35.89 0.3859
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Os pressupostos da Anova foram atendidos?
- Pelo teste de normalidade shapiro e análises gráficas:
- verifica-se normalidade dos resíduos;
- verifica-se independência dos resíduos;
- verifica-se que os resíduos são aproximadamente normais com média igual a zero e variância constante.
par(mfrow=c(2,2))
plot(one.way)
Figura 2: Gráficos de diagnóstico Anova
par(mfrow=c(1,1))
shapiro.test(one.way$residuals)##
## Shapiro-Wilk normality test
##
## data: one.way$residuals
## W = 0.99088, p-value = 0.7594Mas e sobre as somas de quadrados? Como são calculadas?
- Seja k o número de grupos
- Seja \(N\) o número total de indivíduos
tabela <- data.frame(
col1 = c("Fonte de Variação","Entre Grupos (Fator A)","Dentro dos Grupos (Erro)","Total"),
col2 = c("Soma de Quadrados","SQA","SQE","SQT"),
col3 = c("Graus de liberdade","k-1","N-k","N-1"),
col4 = c("Quadrados médios","QMA=$\\frac{SQA}{k-1}$","QME=$\\frac{SQE}{N-k}$",""),
col5 = c("Estatística F","$F=\\frac{QMA}{QME}$","",""),
col6 = c("Valor-p","àrea sob a curva à direita de F","",""))
cap=tab_nums("tab3","Tabela Anova")
tabela %>%
kbl(caption=cap,col.names = NULL,escape = FALSE) %>%
kable_classic(full_width = F, html_font = "Cambria") %>%
kable_styling(bootstrap_options = c("striped","hold_position"),
full_width = F)| Fonte de Variação | Soma de Quadrados | Graus de liberdade | Quadrados médios | Estatística F | Valor-p |
| Entre Grupos (Fator A) | SQA | k-1 | QMA=\(\frac{SQA}{k-1}\) | \(F=\frac{QMA}{QME}\) | àrea sob a curva à direita de F |
| Dentro dos Grupos (Erro) | SQE | N-k | QME=\(\frac{SQE}{N-k}\) | ||
| Total | SQT | N-1 |
onde \[ \begin{array}{l} SQA = \sum_{i=1}^k \sum_{j=1}^{n_i} (\bar{y}_{i.}-\bar{y}_{..})^2\\ SQE = \sum_{i=1}^k \sum_{j=1}^{n_i} (y_{ij}-\bar{y}_{i.})^2 \text{ ou também: } =\sum_{i=1}^k (n_i-1) s_i^2 \\ SQT = \sum_{i=1}^k \sum_{j=1}^{n_i} (y_{ij}-\bar{y}_{..})^2, \text{ ou também: }=SQA + SQE \end{array} \] Interpretação:
- Soma de Quadrados do Fator A (SQA) é o desvio das médias estimadas em cada tratamento (nível) em torno da média geral dos dados. Representa a variabilidade devido aos diferentes níveis do fator A;
- Soma de Quadrados do Erro (SQE) é o desvio das observações em torno da média estimada do seu nível (tratamento). Representa a variabilidade dentro de cada nível do fator A.
Aplicação: Cálculos das somas de quadrados e estatística do Teste.
k=3
N=96
y = dados$yield
y_bar = mean(dados$yield)
y_grupos=dados %>%
group_by(fertilizer) %>%
summarise(media=mean(yield),var=var(yield),n=n())## `summarise()` ungrouping output (override with `.groups` argument)SQA = sum(y_grupos$n*(y_grupos$media-y_bar)^2)
SQE = sum((y_grupos$n-1)*y_grupos$var)
SQT = SQA + SQE
QMA = SQA/(k-1)
QME = SQE/(N-k)
F = QMA/QME
tabela <- data.frame(
col1 = c("Fonte de Variação","Entre Grupos (Fator A)","Dentro dos Grupos (Erro)","Total"),
col2 = c("Soma de Quadrados",round(SQA,4),round(SQE,4),round(SQT,4)),
col3 = c("Graus de liberdade",k-1,N-k,N-1),
col4 =
c("Quadrados médios",round(SQA/(k-1),4),round(SQE/(N-k),4),""),
col5 = c("Estatística F",round(QMA/QME,4),"",""),
col6 = c("Valor-p",round(1-pf(F,k-1,N-k),4),"",""))
cap=tab_nums("tab4","Tabela Anova - aplicação")
tabela %>%
kbl(caption=cap,col.names = NULL,escape = FALSE) %>%
kable_classic(full_width = F, html_font = "Cambria") %>%
kable_styling(bootstrap_options = c("striped","hold_position"),
full_width = F)| Fonte de Variação | Soma de Quadrados | Graus de liberdade | Quadrados médios | Estatística F | Valor-p |
| Entre Grupos (Fator A) | 6.068 | 2 | 3.034 | 7.8628 | 7e-04 |
| Dentro dos Grupos (Erro) | 35.8862 | 93 | 0.3859 | ||
| Total | 41.9542 | 95 |
Post-hoc tests: Teste de Tukey para comparação de médias:
- Qual (quais) fertilizante(s) difere(m) dos demais?
- Pelo teste de Tukey, o fertilizante 3 difere dos fertilizantes 1 e 2 (valor-p<0,05) nas duas comparações. Por outro lado, os fertilizantes 1 e 2 não diferem entre si (valor-p>0,05).
tukey.one.way<-TukeyHSD(one.way)
tukey.one.way## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = yield ~ fertilizer, data = dados)
##
## $fertilizer
## diff lwr upr p adj
## 2-1 0.1761687 -0.19371896 0.5460564 0.4954705
## 3-1 0.5991256 0.22923789 0.9690133 0.0006125
## 3-2 0.4229569 0.05306916 0.7928445 0.0208735Anova & modelo de regressão linear: Linear model
- Modelo de regressão linear com variável categórica:
- consiste na seguinte expressão para a resposta do indivíduo \(i\): \[ y_j = \left\{ \begin{array}{lll} \beta_1 + \epsilon_j, \text{ se foi aplicado o fertilizante 1},\\ \beta_1 + \beta_2 + \epsilon_j, \text{ se foi aplicado o fertilizante 2},\\ \beta_1 + \beta_3 + \epsilon_j, \text{ se foi aplicado o fertilizante 3}, \end{array} \right. \] onde \(\epsilon_i \sim N(0,\sigma^2)\)
Alternativamente, podemos escrever o modelo de regressão na seguinte forma: \[ y_j = \left\{ \begin{array}{lll} \underbrace{\bar{x}_{1.}}_{\beta_1: \text{ média do grupo 1}} + \epsilon_j, \text{ se foi aplicado o fertilizante 1},\\ \underbrace{\bar{x}_{2.}}_{\beta_1 + \beta_2: \text{ média do grupo 2}} + \epsilon_j, \text{ se foi aplicado o fertilizante 2},\\ \underbrace{\bar{x}_{3.}}_{\beta_1 + \beta_3: \text{ média do grupo 3}} + \epsilon_j, \text{ se foi aplicado o fertilizante 3}, \end{array} \right. \]
model.1=lm(yield ~ fertilizer,data=dados)
summary(model.1)##
## Call:
## lm(formula = yield ~ fertilizer, data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.39620 -0.45303 0.02268 0.41076 1.70473
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 176.7570 0.1098 1609.643 < 2e-16 ***
## fertilizer2 0.1762 0.1553 1.134 0.259542
## fertilizer3 0.5991 0.1553 3.858 0.000211 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6212 on 93 degrees of freedom
## Multiple R-squared: 0.1446, Adjusted R-squared: 0.1262
## F-statistic: 7.863 on 2 and 93 DF, p-value: 0.0006999Interpretação Gráfica: Modelo linear com somente 1 variável categórica
- Os indivíduos de cada grupo de fertilizante apresentam desvios com relação à média do seu grupo => contribuição para a soma de quadrados dos erros;
- As médias dos grupos apresentam desvios com relação à média geral => contribuição para a soma de quadrados do fator fertilizante.
dados1 = dados %>%
mutate(predito=model.1$fitted.values)cat_plot(model.1, pred = fertilizer, geom = "line")## Warning: fertilizer are not included in an interaction with one another in the
## model.
dados1 %>%
ggplot(aes(fertilizer,yield))+
geom_point()+
geom_point(mapping=aes(fertilizer,predito),color="red")
Figura 3: Gráfico de pontos dos valores observados e preditos por grupo
2 Anova com 2 fatores - Two way Analise of Variance
dados fictícios: https://www.scribbr.com/statistics/anova-in-r/
Two Way Anova: Fertilizante e densidade influenciam na produção?
- Descritiva considerando apenas fertilizante:
- Pela tabela 5, o fertilizante 3 com densidade igual a 2 possui média superior em produção;
- Pelos Gráficos de Boxplot, verificamos que o fertilizante 3 em geral possui valores superiores de produção e a densidade 2 em geral possui valores superiores de produção.
y_grupos=dados %>%
group_by(fertilizer,density) %>%
summarise(media=mean(yield),var=var(yield),n=n())## `summarise()` regrouping output by 'fertilizer' (override with `.groups` argument)cap=tab_nums("tab5","Descritivas por tipo de fertilizante e densidade")
datatable(y_grupos,cap=cap)dados %>%
ggplot(aes(x=fertilizer,y=yield,fill=density)) +
geom_boxplot() +
facet_wrap(vars(density),nrow=3)
Figura 4: Boxplots de produção por tipos de fertilizante e densidade
dados %>%
ggplot(aes(x=density,y=yield,fill=fertilizer)) +
geom_boxplot() +
facet_wrap(vars(fertilizer),nrow=3)
O gráfico de interação:
- O gráfico sugere que existe interação entre fertilizante e densidade;
- os segmentos de reta não são paralelos, sugerindo que o fertilizante age de forma diferente para cada densidade;
- o próximo passo é aplicar a Análise de Variância (ANOVA), ou equivalentemente, o modelo linear com variáveis categóricas (densidade e fertilizante), para testarmos a hipótese de interação.
interaction.plot(x.factor = dados$density,
trace.factor = dados$fertilizer,
response = dados$yield,
fun = mean,
type="b",
col=c("black","red","green"), ### Colors for levels of trace var.
pch=c(19, 17, 15), ### Symbols for levels of trace var.
fixed=TRUE, ### Order by factor order in data
leg.bty = "o",
ylab="Produção",
xlab="Densidade")
- Pela tabela ANOVA verificamos a interação entre fertilizante e densidade não é significativa (ao nível de 5%).
two.way <- aov(yield ~ density*fertilizer, data = dados)
summary(two.way)## Df Sum Sq Mean Sq F value Pr(>F)
## density 1 5.122 5.122 15.195 0.000186 ***
## fertilizer 2 6.068 3.034 9.001 0.000273 ***
## density:fertilizer 2 0.428 0.214 0.635 0.532500
## Residuals 90 30.337 0.337
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Os pressupostos da Anova foram atendidos?
- Pelo teste de normalidade shapiro e análises gráficas:
- verifica-se normalidade dos resíduos;
- verifica-se independência dos resíduos;
- verifica-se que os resíduos são aproximadamente normais com média igual a zero e variância constante.
par(mfrow=c(2,2))
plot(two.way)
par(mfrow=c(1,1))
shapiro.test(two.way$residuals)##
## Shapiro-Wilk normality test
##
## data: two.way$residuals
## W = 0.98527, p-value = 0.3601Post-hoc tests: Teste de Tukey para comparação de médias:
- Qual (quais) fertilizante(s) difere(m) dos demais?
- Pelo teste de Tukey, o fertilizante 3 difere dos fertilizantes 1 e 2 (valor-p<0,05) nas duas comparações. Por outro lado, os fertilizantes 1 e 2 não diferem entre si (valor-p>0,05);
- Com relação a densidade, vemos que as densidades 1 e 2 diferem entre si.
tukey.two.way<-TukeyHSD(two.way)
tukey.two.way## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = yield ~ density * fertilizer, data = dados)
##
## $density
## diff lwr upr p adj
## 2-1 0.461956 0.2265141 0.697398 0.0001864
##
## $fertilizer
## diff lwr upr p adj
## 2-1 0.1761687 -0.16972711 0.5220646 0.4482026
## 3-1 0.5991256 0.25322974 0.9450214 0.0002393
## 3-2 0.4229569 0.07706102 0.7688527 0.0123951
##
## $`density:fertilizer`
## diff lwr upr p adj
## 2:1-1:1 0.63489351 0.03715141 1.2326356 0.0306553
## 1:2-1:1 0.33868995 -0.25905215 0.9364320 0.5680611
## 2:2-1:1 0.64854101 0.05079891 1.2462831 0.0254221
## 1:3-1:1 0.69601055 0.09826845 1.2937526 0.0128711
## 2:3-1:1 1.13713411 0.53939201 1.7348762 0.0000044
## 1:2-2:1 -0.29620356 -0.89394566 0.3015385 0.7007889
## 2:2-2:1 0.01364750 -0.58409459 0.6113896 0.9999998
## 1:3-2:1 0.06111704 -0.53662505 0.6588591 0.9996758
## 2:3-2:1 0.50224060 -0.09550149 1.0999827 0.1515870
## 2:2-1:2 0.30985106 -0.28789103 0.9075932 0.6591096
## 1:3-1:2 0.35732060 -0.24042149 0.9550627 0.5089535
## 2:3-1:2 0.79844416 0.20070206 1.3961863 0.0025700
## 1:3-2:2 0.04746954 -0.55027256 0.6452116 0.9999065
## 2:3-2:2 0.48859310 -0.10914900 1.0863352 0.1742960
## 2:3-1:3 0.44112356 -0.15661854 1.0388657 0.2721714Anova & modelo de regressão linear: Linear model
- Modelo de regressão linear com duas variáveis categóricas + efeito da interação:
- Podemos obter os coeficientes de cada fator, além dos coeficientes da interação.
- Verificamos que os coeficientes da interação fertilizante x densidade não são significativos ao nível 5%.
model.2=lm(yield ~ fertilizer + density + fertilizer*density,data=dados)
summary(model.2)##
## Call:
## lm(formula = yield ~ fertilizer + density + fertilizer * density,
## data = dados)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.1917 -0.3338 -0.0402 0.3497 1.4842
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 176.4396 0.1451 1215.607 < 2e-16 ***
## fertilizer2 0.3387 0.2053 1.650 0.10243
## fertilizer3 0.6960 0.2053 3.391 0.00104 **
## density2 0.6349 0.2053 3.093 0.00264 **
## fertilizer2:density2 -0.3250 0.2903 -1.120 0.26581
## fertilizer3:density2 -0.1938 0.2903 -0.668 0.50616
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5806 on 90 degrees of freedom
## Multiple R-squared: 0.2769, Adjusted R-squared: 0.2367
## F-statistic: 6.893 on 5 and 90 DF, p-value: 1.728e-05anova(model.2)## Analysis of Variance Table
##
## Response: yield
## Df Sum Sq Mean Sq F value Pr(>F)
## fertilizer 2 6.0680 3.0340 9.0011 0.0002732 ***
## density 1 5.1217 5.1217 15.1945 0.0001864 ***
## fertilizer:density 2 0.4278 0.2139 0.6346 0.5325001
## Residuals 90 30.3367 0.3371
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13 Exemplo 2: Anova com 2 fatores - Two way Analysis of Variance
dados fictícios: https://rcompanion.org/handbook/G_09.html
- Dados de ganho de peso baseado em dieta e país de origem
Input =("
Diet Country Weight_change
A USA 0.120
A USA 0.125
A USA 0.112
A UK 0.052
A UK 0.055
A UK 0.044
B USA 0.096
B USA 0.100
B USA 0.089
B UK 0.025
B UK 0.029
B UK 0.019
C USA 0.149
C USA 0.150
C USA 0.142
C UK 0.077
C UK 0.080
C UK 0.066
")
dados2 = read.table(textConnection(Input),header=TRUE)
datatable(dados2)- Descritiva considerando apenas fertilizante:
- Pela tabela 6, a dieta C juntamente com país de origem USA possui maior média de ganho de peso;
- Pelos Gráficos de Boxplot, verificamos que o fertilizante 3 em geral possui valores superiores de produção e a densidade 2 em geral possui valores superiores de produção.
dados2_grupos=dados2 %>%
group_by(Diet,Country) %>%
summarise(media=mean(Weight_change),var=var(Weight_change),n=n())## `summarise()` regrouping output by 'Diet' (override with `.groups` argument)cap=tab_nums("tab6","Descritivas por tipo de dieta e país de origem")
datatable(dados2_grupos,cap=cap)Gráfico de interação: Já vamos direto para este gráfico que mostra a relação entre as duas variáveis na interação.
- O gráfico sugere que não existe interação entre dieta e país de origem;
- os segmentos de reta são paralelos, sugerindo que a dieta age de forma semelhante para cada país de origem;
- o próximo passo é aplicar a Análise de Variância (ANOVA), ou equivalentemente, o modelo linear com variáveis categóricas (dieta e país de origem), para testarmos a hipótese de interação.
interaction.plot(x.factor = dados2$Country,
trace.factor = dados2$Diet,
response = dados2$Weight_change,
fun = mean,
type="b",
col=c("black","red","green"), ### Colors for levels of trace var.
pch=c(19, 17, 15), ### Symbols for levels of trace var.
fixed=TRUE, ### Order by factor order in data
leg.bty = "o",
ylab="weight change",
xlab="Diet")
- Pela tabela ANOVA verificamos a interação entre país de origem e tipo de dieta não é significativa (ao nível de 5%); - O tipo de dieta e o país de origem são significativos.
two.way.2 <- aov(Weight_change ~ Diet*Country, data = dados2)
summary(two.way.2)## Df Sum Sq Mean Sq F value Pr(>F)
## Diet 2 0.007804 0.003902 114.205 1.55e-08 ***
## Country 1 0.022472 0.022472 657.717 7.52e-12 ***
## Diet:Country 2 0.000012 0.000006 0.176 0.841
## Residuals 12 0.000410 0.000034
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Os pressupostos da Anova foram atendidos?
- Pelo teste de normalidade shapiro e análises gráficas:
- NÂO verifica-se normalidade dos resíduos;
- verifica-se independência dos resíduos e variância constante;
Diagnóstico
par(mfrow=c(2,2))
plot(two.way.2,1:2)
#par(mfrow=c(1,1))
shapiro.test(two.way.2$residuals)##
## Shapiro-Wilk normality test
##
## data: two.way.2$residuals
## W = 0.87589, p-value = 0.02228
Alternativas:
Transformação de Box Cox;
Testes não paramétricos.
Pela transformação de Box-Cox, o valor de \(\lambda=1\) - significa que a variável resposta em sua forma transformada será:
\[ Y^t = \frac{Y^1-1}{1}= Y-1, \text{ diferindo de Y apenas de uma constante, então não fará diferença da variável original} \]
boxCox(two.way.2)
- Referências: Box, G. E. P. and Cox, D. R. (1964) An analysis of transformations. Journal of the Royal Statisistical Society, Series B. 26 211-46.
- <a href=“https://www.ime.usp.br/~abe/lista/pdfm9cJKUmFZp.pdf” < Material sobre Box Cox
Curiosidades: Teste para homogeneidade dos resíduos em um modelo linear:
- Pelo teste de Breusch-Pagan, não rejeitamos a variância constante dos resíduos. Sendo assim, o problema é a normalidade!
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericbptest(two.way.2)##
## studentized Breusch-Pagan test
##
## data: two.way.2
## BP = 3.08, df = 5, p-value = 0.6877Verificar outras alterativas como transformação potência Utilizada em modelos lineares, disponível no pacote car (Companion to Applied Regression).