Análise de Variância
Marcelo Franco
2023-08-20
Objetivo
Aplicação de modelo linear para inferir largura da sépala de uma flor em função de sua espécie.
Metodologia
O modelo linear usado será o modelo de Análise de Variância. Derivado do modelo de regressão linear simples, é usado quando temos uma variável dependente quantitativa e uma variável independente qualitativa ou categórica com três ou mais grupos. Possue equação e interpretação semelhante ao modelo de regressão linear simples. De tal modo que,
\[Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + erro\]
Pacotes
Pacotes usados durante o script para manipução de dados, visualização de dados e testes estatísticos.
library(tidyverse)
library(gvlma)
library(car)
library(broom)Dados
Os dados usados fornecem informações sobre as medidas em centímetros de pétalas e sépalas de 50 flores de 3 espécies de iris (setosa, versicolor e virginica).
head(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 setosatail(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 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 virginicaAnálise Exploratória
O objetivo dessa análise será investigar possíveis condições presentes nos dados que possam gerar problemas com as premissas do modelo de regressão.
iris %>%
ggplot(aes(x = Species, y = Sepal.Width, fill = Species)) +
geom_boxplot()
O boxplot nos indica que há presença de outliers nas espécies setosa e verginica. Todavia, nem todo valor extremo tem poder estatístico suficiente para prejudicar um modelo de regressão. Como a distância de tais outliers não é muito grande em relação ao valor mínimo/máximo, iremos mantê - los. Após estimar o modelo, será possível verificar se os valores extremos devem ser removidos.
Outra informação relevante que percebemos ao analisarmos o boxplot é que um dos principais valores de tendência central, a mediana, difere entre as três espécies. Isso pode implicar em resíduos agrupados ou, em outras palavras, isso pode gerar depência nos resíduos de regressão. Mas como essa diferença não é, pelo menos visualmente, tão grande e o tamanho do corpo/caixa (distância interquantil) do grafíco é semelhante entre as três espécies, isso indica que a variância entre as espécies não é muito grande, consequentemente iremos prosseguir com a modelagem.
Para que o leitor possa visualizar e comparar melhor a variância dos dados entre as espécies, podemos refazer o gráfico boxplot com algumas alterações.
iris %>%
ggplot(aes(x = Species, y = Sepal.Width)) +
geom_boxplot(outlier.colour = NA) +
geom_jitter()
Note que o espalhamento dos pontos para as três espécies é semelhante. Indicativo de que não teremos problemas na modelagem.
Modelagem
Estimando o modelo linear da largura da sépala em função da espécie.
modelo <- lm(Sepal.Width ~ Species, data = iris)Antes de prosseguirmos com a inferência ou interpretação do modelo estatísco, precimos avaliar se suas premissas foram atendidas. Faremos isso através do metódo grafico e dos testes de hipótese.
Validação do Modelo
par(mfrow=c(2,2))
plot(modelo) #(1) método gráfico
Como podemos observar no primeiro gráfico Residuals vs Fitted, não há tendência nos resíduos. Para chegar a tal conclusão, devemos observar a linha vermelha praticamente reta/constante em zero e que os resíduos possuem distância vertical semelhante ao longo de toda a reta vermelha. Embora possa parecer que há agrupamento de dados (o que seria um problema) não é esse o caso. Isso é causado e esperado pela própria estrutura do modelo que tem como variável independente uma var categórica. Mais a frente iremos plotar um gráfico semelhante mostrando mais claramente a homocedasticidade.
No segundo gráfico Q-Q Residuals, observamos resíduos aproximadante normais. Caso isso não tenha ficado claro ao leitor, iremos, mais a frente, plotar o mesmo gráfico acrescido de um intervalo de confiança para a normalidade.
O gráfico Scale-Location pode ser interpretado de forma análoga ao primeiro gráfico, observando se os pontos estão distribuídos verticalmente de forma semelhante ao longo da reta vermelha. Observamos a homocedasticidade (condição necessária para inferência).
No quarto e último gráfico Constant Lavarage: Residuals vs Factor Levels, avaliamos através da distância de Cook se há observações influentes (valores extremos que podem distorcer o modelo). Tais pontos, quando existem, são indicados pela sua posição destacada em relação aos demais e delimitados por uma curva nas extremidades do gráfico. Corroborando a avalição de outliers feita durante a análise exploratória, não observamos qualquer evidência de valores extremos no gráfico.
Como mencionado anteriormente, iremos destacar algumas avaliações refazendo dois plots para que o leitor tenha uma perspectiva mais clara em relação às conclusões feitas a partir do primeiro e do segundo gráfico.
Para isso, vamos extrair os resíduos do modelo estimado juntamente com outras informações necessárias para o plot do primeiro gráfico.
resultado <- augment(modelo)head(resultado)## # A tibble: 6 × 8
## Sepal.Width Species .fitted .resid .hat .sigma .cooksd .std.resid
## <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3.5 setosa 3.43 0.0720 0.0200 0.341 0.000312 0.214
## 2 3 setosa 3.43 -0.428 0.0200 0.339 0.0110 -1.27
## 3 3.2 setosa 3.43 -0.228 0.0200 0.340 0.00313 -0.678
## 4 3.1 setosa 3.43 -0.328 0.02 0.340 0.00647 -0.975
## 5 3.6 setosa 3.43 0.172 0.02 0.341 0.00178 0.511
## 6 3.9 setosa 3.43 0.472 0.02 0.339 0.0134 1.40resultado %>%
ggplot(aes(x = Species, y = .resid)) +
geom_jitter() +
geom_hline(yintercept = 0, color = 'red', type = "deshad")
Agora, de modo mais claro, podemos verificar a homocedasticidade. Note que os pontos estão distribuidos de forma semelhante em torno da linha vermelha que é constante em zero.
Em relação à normalidade dos resíduos, vemos refazer o segundo gráfico Q-Q Residuals acrescido do intervalo de confiança para a distribuição normal.
qqPlot(modelo)
## [1] 16 42Nitidamente observamos que a grande maioria dos pontos, exceto por uma ínfima parcela, está dentro do intervalo de confiança. Mesma conclusão que chegamos anteriormente.
Para finalizarmos a avaliação do modelo e partirmos para a inferência, iremos testar as premissão do modelo atrevés do teste de hipótese.
gvlma(modelo)##
## Call:
## lm(formula = Sepal.Width ~ Species, data = iris)
##
## Coefficients:
## (Intercept) Speciesversicolor Speciesvirginica
## 3.428 -0.658 -0.454
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = modelo)
##
## Value p-value Decision
## Global Stat 3.661e+00 0.4538 Assumptions acceptable.
## Skewness 1.857e-02 0.8916 Assumptions acceptable.
## Kurtosis 1.613e+00 0.2041 Assumptions acceptable.
## Link Function 5.545e-15 1.0000 Assumptions acceptable.
## Heteroscedasticity 2.029e+00 0.1543 Assumptions acceptable.Considerando um \(\alpha\) de \(5\)% não rejeitamos a \(H_0\). Em outras palavras, a avalição dos p-valores, de forma resumida, nos indica que os resíduos têm distribuição normal com média zero e a variância de tais resíduos é constante (homocedasticidade). Mesma conclusão que chegamos através da avaliação dos gráficos.
Inferência
Com as premissas do modelo atendidas, podemos interpretar os parâmetros e fazer a inferência estatística.
Primeiro iremos verificar se existe diferença estatística entre as médias das especíes. Como queremos testar médias entre três grupos (setosa, versicolar e virginica), o teste estatístico apropriado é a ANOVA. De tal modo que,
\[H_o: \mu_1 = \mu_2 = \mu_3\]
anova(modelo)## Analysis of Variance Table
##
## Response: Sepal.Width
## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 11.345 5.6725 49.16 < 2.2e-16 ***
## Residuals 147 16.962 0.1154
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Dado o p-valor obtido rejeitamos a \(H_0\). Em outras palavras, podemos inferir que pelo menos uma das espécies tem média diferente das demais.
Seguindo para a interpretação dos coeficientes.
coef(modelo)## (Intercept) Speciesversicolor Speciesvirginica
## 3.428 -0.658 -0.454Os coficientes do modelo indicam que a setosa foi jogada na baseline. Traduzindo, o valor médio da espécie setosa é dado pelo próprio \(\beta_0\) Intercept (3.428) e a interpretação dos demais coeficientes deve ser feita sempre em relação ao próprio \(\beta_0\).
O coeficiente \(\beta_1\) Speciesversicolor nos diz que o valor médio da espécie versicolor é 0.658cm menor do que a espécie de referência que está no baseline, a setosa.
O coeficiente \(\beta_2\) Speciesvirginica nos diz que o valor médio da espécie virginica é 0.454cm menor que a espécie de referência que está no baseline, a setosa.
summary(modelo)##
## Call:
## lm(formula = Sepal.Width ~ Species, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.128 -0.228 0.026 0.226 0.972
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.42800 0.04804 71.359 < 2e-16 ***
## Speciesversicolor -0.65800 0.06794 -9.685 < 2e-16 ***
## Speciesvirginica -0.45400 0.06794 -6.683 4.54e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3397 on 147 degrees of freedom
## Multiple R-squared: 0.4008, Adjusted R-squared: 0.3926
## F-statistic: 49.16 on 2 and 147 DF, p-value: < 2.2e-16\(R^2\) de 0.40, indicando que o modelo estimado é capaz de explicar 40% das variações em \(Y\). Também observamos, com base nos p-valores, que todos os coeficientes são significativos bem como o coeficiente de determinação \(R^2\).
Breve Comentário Sobre Intervalo de Confiança
Durante a interpretação dos coeficientes, o termo valor médio foi usado algumas vezes. Porém, da onde vem esse valor? Por que ele é o valor médio? Tentaremos responder isso observando o intervalo de confiança para os valores estimados.
confint(modelo)## 2.5 % 97.5 %
## (Intercept) 3.3330635 3.5229365
## Speciesversicolor -0.7922604 -0.5237396
## Speciesvirginica -0.5882604 -0.3197396A função confint nos retorna o intervalo para os valores dos parâmetros do modelo estimado com 95% de confiança.
Observe que o valor zero não está presente em nenhum intervalo. Isso faz com que exista uma tendência de que os p-valores dos coeficientes estimados sejam próximos de zero e ainda, que a diferença entre o tamanho médio entre as espécies não é zero. Em outros termos, isso indica que é muito provável que existe, sim, diferença estatística de tamanho médio para as diferentes espécies.
Note que o valor médio é dado pela média aritimética das extremidades do intervalo. De tal modo que se o valor médio para \(\beta_0 = 3.42800\) podemos chegar nesse número com:
(3.3330635 + 3.5229365)/2 #extremidades do intervalo de confiança de beta0 dividido por 2.## [1] 3.428Da mesma forma, se o valor médio para \(\beta_1 = -0.65800\) podemos encontrar esse número com a seguinte expressão:
(-0.7922604 -0.5237396)/2 #extremidades do intervalo de confiança de beta1 dividido por 2.## [1] -0.658Por último, se \(\beta_2 = -0.454\) podemos encontrar esse valor de modo análogo.
(-0.5882604 -0.3197396)/2 #extremidades do intervalo de confiança de beta2 dividido por 2.## [1] -0.454Conclusão
Atrevés de um modelo linear, derivado da mesma equação do modelo de regressão linear simples, podemos duduzir o modelo de análise de variância. Modelo também linear usado quando temos uma variável dependente quantitativa e uma variável independente qualitativa com três ou mais níveis.
Percebemos a importância da análise exploratória para identificar possíveis problemas na hora da inferência. Entende - se que nessa etapa, mesmo que considerados outilers, os pontos mais discrepantes nem sempre prejudicam o modelo e suas interpretações.
Destaca - se a relação entre valores e aborgens na interpretação e avaliação de modelos. Por exemplo, os valores dos coeficientes estimados estão relacionados com o intervalo de confiança dos mesmos. Embora não mencionado, o valor do coeficiente de determinação \(R^2\) está relacionado com valores de covariancia. Esses são apenas alguns exemplos em meio a tantos outros. No caso da avaliação dos modelos, entende - se que a melhor abordagem é feita através da visualização com auxílio de diferentes gráficos. Todavia, a análise númerica via teste de hipótese é uma ferrementa interessante que não deve ser ignorada.