Como comparar a variância de uma amostra

Vou supor que as amostras são de variáveis normais, com média 100 e desvio padrão 12. Para testar a homogeneidade das variâncias, há o teste de Levene, que se pode encontrar no pacote car.

library(car) set.seed(343) # Torna o exemplo reprodutível n <- 5 x <- rnorm(n, mean = 100, sd = 12) y <- rnorm(n, mean = 100, sd = 12) z <- rnorm(n, mean = 100, sd = 12)

Agora, juntamos as amostras numa variável v de uma data.frame, com um factor, a variável f que nos diz a qual amostra pertence cada linha da df.

dados <- data.frame(f = rep(letters[24:26], each = n), v = c(x, y, z)) leveneTest(v ~ f, data = dados, center = mean) #Levene's Test for Homogeneity of Variance (center = mean) # Df F value Pr(>F) #group 2 0.2198 0.8058 # 12

Uma forma alternativa, tendo em conta que referiu a tabela ANOVA, será usar as funções para o modelo linear, tanto lm como aov.

modelo <- aov(v ~ f, data = dados) leveneTest(modelo, center = mean)

Os resultados são os mesmos.

EDIÇÃO.
Há também o teste de Brown-Forsythe, que apesar de talvez ser menos utilizado, pode ter vantagens, nomeadamente quanto à robustez. Uma função que realiza esse teste pode ser encontrada no pacote onewaytests.

library(onewaytests) bf.test(v ~ f, data = dados) # # Brown-Forsythe Test #--------------------------------------------------------- # data : v and f # # statistic : 0.156811 # num df : 2 # denom df : 11.3735 # p.value : 0.8566824 # # Result : Difference is not statistically significant. #--------------------------------------------------------- #

Mais uma vez, não há evidência de heterogeneidade de variâncias. A hipótese nula não é rejeitada. Podemos portanto realizar um teste ANOVA para a diferença de médias, onde uma das assunções é precisamente a igualdade de variâncias.

Nos testes de hipóteses utilizamos os conhecimentos de estimação de parâmetros e probabilidade na tomada de decisão a favor de uma hipótese com base nas evidências (amostra) que lançamos mão.

Suponha a situação em que desejamos comparar dois treinamentos de atletas, o treinamento A e o treinamento B. Ambos tipos de treinamentos foram aplicados em atletas da mesma faixa etária e de similar performance, sendo como suposição, que qualquer melhora ou piora na performance dos atletas pode ser atribuída ao tipo de treinamento aplicado.

Para testar a hipótese de que o treinamento A faz com que a performance dos atletas seja inferior daqueles que realizaram o treinamento B testaremos as seguintes hipóteses:

$H_o:$ de que não existe diferença na performance média dos atletas treinados com o A ou B.
$H_1:$ de que a performance média dos atletas treinados com o A é inferior ao de B.

Cabe notar que a hipótese alternativa poderia ser de que a performance média dos atletas treinados com o A é superior ou ainda diferente…

Traduzindo isso em termos de parâmetros, temos:

$H_o: \mu_A = \mu_B$
$H_1: \mu_A < \mu_B$

Ou seja vamos testar a hipótese nula frente a hipótese alternativa. De certa forma verificaremos o quão inferior é o parâmetro do treinamento de A, e se é possível que essa diferença tenha vindo da hipótese alternativa.

O procedimento para se testar as hipóteses é coletar uma amostra de cada população treinada com cada uma das técnicas e verificar se as estimativas dos parâmetros são estatisticamente diferentes. Ou melhor, a diferença entre as estimativas é muito grande a ponto de não ser considerada ter vindo da hipótese nula, aquela que diz que ambos os treinamentos são iguais?

Quando dizemos que faremos o teste a partir de estimativas é por que não temos acesso à população, consequentemente desconhecemos os parâmetros populacionais. Caso conhecêssemos não precisaríamos realizar o teste, e sim somente verificar os parâmetros populacionais.

Somente como critério pedagógico vamos apresentar as populações subjacentes (porém assuma que não conhecemos nada a respeito dela).

performance_A = rnorm(100000,mean=12,sd=2.7) performance_B = rnorm(100000,mean=13,sd=2.8) hist(performance_A, freq=FALSE, breaks=100, col=rgb(0,0,1,1/4), xlim = c(0,20), ylim = c(0,0.3), main=expression(paste("População das performances de A(",mu,'=12, ',~ sigma^{2},'=2.7) '," B(",mu,'=13, ',~ sigma^{2},'=2.8)')), xlab = "A e B") hist(performance_B, freq=FALSE, breaks=100, col=rgb(1,0,0,1/4), xlim = c(0,20), add=TRUE)

Agora tomaremos a nossa amostra de cada uma das populações, treinamento A e B, ambas com parâmetros desconhecidos ($\mu_A$, $\mu_B$, $\sigma^2_A$ e $\sigma^2_B$). Para o presente caso tomaremos uma amostra aleatória de tamanho $n = 15$ para cada treinamento. Uma vez com as amostras calcularemos as estatísticas amostrais (estimadores de parâmetros) e testaremos a hipótese nula para a diferença das médias.

n_A = 15 n_B = 15 amostra_A = rnorm(n_A,mean=12,sd=2.7) amostra_B = rnorm(n_B,mean=13,sd=2.8) #cat("média A = ", mean(amostra_A), "\ndes. padrão A = ", sd(amostra_A)) #cat("\nmédia B = ", mean(amostra_B), "\ndes. padrão B = ", sd(amostra_B))

$\overline{x}_A =$ 11.8751483
$\overline{x}_B =$ 13.2034837
$S^2_A =$ 7.3349517
$S^2_B =$ 8.3419714
$S_A =$ 2.7083116
$S_B =$ 2.8882471

Para este caso faremos um Teste-T (t-student), pois não conhecemos as variâncias populacionais ($\sigma_A^2, \sigma_B^2$), dessa forma calculamos uma estimativa da mesma ($S_A^2,S_B^2$), sendo essas variáveis aleatórias com uma distribuição aproximada pela $\chi^2_{n-1}$. Considerando que $\sigma_A^2 = \sigma_B^2 = \sigma^2$, utilizaremos ambas as informações para calcular uma estimativa única de $\sigma^2$ que chamaremos de $S_p^2$ (essa abordagem é explicada nos exemplos que seguem na próxima seção).

\[ T = \frac{(x_A - x_B) - (\mu_A - \mu_B)}{S_p \sqrt{\frac{1}{n_A} + \frac{1}{n_A}}} \]

# medias amostrais x_A = mean(amostra_A) x_B = mean(amostra_B) # tamanos das amostras n_A = 15 n_B = 15 # desvio padrão amostral S_A = sd(amostra_A) S_B = sd(amostra_B) # variância amostral S2_A = var(amostra_A) S2_B = var(amostra_B) # variancia combinada, onde o peso da amostra influência na combinação S2_D = ((n_A - 1)/(n_A + n_B - 2))*S2_A + ((n_B - 1)/(n_A + n_B - 2))*S2_B S_D = sqrt(S2_D) # Calculando a diferença entre as duas amostras t = (x_A - x_B)/(S_D*sqrt(1/n_A + 1/n_B)) t # estatística do teste [1] -1.299341 T = t.test(amostra_A,amostra_B,alternative = "less")

Observe que a estatística diferença entre as médias $\overline{x}_A - \overline{x}_B = -1.3283$ é tranformada na estatística $t = -1.2993$ na distribuição $H_0$. Dessa forma podemos calcular a probabilidade de um valor de diferença entre as médias ser igual a esse valor ou um valor mais extremo que esse ($Valor-P$).

Note que o $Valor-P$ é 0.1022338, o que indica a probabilidade de cometermos o erro do tipo I a partir das evidências, ou seja a nossa amostra.

library(Plothtests) plot_T_test(alpha = 0.05, alternative = "less", var.equal = "equal", m1 = x_A, m2 = x_B, v1 = S2_A, v2 = S2_B, n1 = 15, n2 = 15)

$test.statistic [1] -1.299341 $df [1] 28 $p.value [1] 0.1022122

Para testar se de fato as variâncias são iguais ou não recorremos ao teste F, também como conhecido como razão de variâncias (pois nunca tomamos a diferença dessas quantidades).

Primeiramente precisariamos estimar as variâncias, a partir desse ponto precisamos saber se as variâncias (estimadas) são iguais ou diferentes. O método para determinar o tamanho da relação de uma variância com a outra é dada pela razão de variâncias.

Dessa forma nosso primeiro passo será determinar se as variâncias são iguais ou não.

Vamos então realizar na verdade um teste de hipóteses para a variância.

A variância amostral (ou seja a variável aleatória variância amostral) possui uma distribuição de densidade de probabilidade, que é aproximadamente uma qui-quadrado $\chi^2$, ou seja a variável aleatória variância amostral ($S^2$) possui uma distribuição qui-quadrado.

Na verdade o que podemos provar é que a quantidade $\big(\frac{n-1}{\sigma^2}\big)S^2 \sim \chi^2_{n-1}$, e ainda assim desde que a variável aleatória em questão tenha distribuição normal.

Um vez que pretendemos testar a razão de duas variáveis aleatórias (variâncias amostrais) que se distribuiem dessa forma aproximada a ($\chi^2_{n-1}$), qual será a distribuição dessa nova variável aleatória?

Temos que essa razão de variâncias é nada mais é do que uma função de variáveis aleatórias.

Existe uma distribuição, desenvolvida por Ronald A. Fischer, distribuição F, que é uma função de duas variáveis aleatórias, ou seja a variável aleatória F é dada por:

\[ F = \frac{Q_1/gl_1}{Q_2/gl_2}\\ \\ Q_1 \sim \chi^2_{n_1-1}\\ \\ Q_2 \sim \chi^2_{n_2-1} \] Ou seja a distribuição F é a razão de duas variáveis aleatórias independentes com distribuição de qui-quadrado divido pelos seus resectivos graus de liberdade ($gl$). A distribuição possui dois parâmetros, que são os graus de liberdade do numerador e denominador, $F_{gl_1,gl_2}$

George W. Snedecor em seu livro Statisical Methods (1937) descreveu o teste de razão de variâncias, chamando de Teste-F em homenagem a Ronald A. Fischer.

Utilando as quantidades amostrais $\big(\frac{n-1}{\sigma^2}\big)S^2$ que seguem uma distribuição de qui-quadrado com $n-1$ graus de liberdade temos.

\[ F = \frac{\big(\frac{n_A-1}{\sigma_A^2}\big)S_A^2 / n_A-1}{\big(\frac{n_B-1}{\sigma_B^2}\big)S_B^2 / n_B-1}\\ \\ F = \frac{S_A^2}{S_B^2} \frac{\sigma_B^2}{\sigma_A^2} \sim F_{gl_A,gl_B}\\ \]

Sabemos, que $F \sim F_{(gl_A,gl_B)}$ é uma distribuição F (Fischer) com graus de liberdade $gl_A$ do numerador e $gl_B$ do denominador, e essa distribuição é uma razão de duas variáveis aleatórias $\chi^2_{n-1}$.

Abaixo segue uma simulação do resultado da razão de duas distribuições $\chi^2_{n-1}$ e uma verificação da distribuição aproximada dessa variável aleatória, para fins didáticos. Ainda, o teste de Kolmogorov-Smirnov foi aplicado para verificar a adequação desses dados ($F$) à uma distribuição F. O resultado da estatística do teste foi $D = 0.00061$, indicando que não pode ser rejeitado que esses dados vieram de uma distribuição F com os respectivos graus de liberdade.

S2_A = rchisq(1000000,df=8)# variância A - v.a. chi-quadrado, gl = 8 S2_B = rchisq(1000000,df=8)# variância A - v.a. chi-quadrado, gl = 8 par(mfrow=c(2,2)) hist(S2_A, freq=FALSE, breaks = 100, main=expression(S^2[A])) hist(S2_B, freq=FALSE, breaks = 100, main=expression(S^2[B])) F = (S2_A/8)/(S2_B/8) hist(F, freq=FALSE, breaks = 100,xlim = c(0,70)) Fischer = rf(100000,df1 = 8, df=8) hist(Fischer, freq=FALSE, breaks = 100, xlim = c(0,70))

Note que a simulação não foi apresentada aqui, porém o codigo se encontra abaixo.

# code para verificar o ajuste dos dados da razao de chi-quadrados com a distribuição F library(fitdistrplus) fmle = fitdist(F, "f", start=list(8, 8)) # graus de liberdade df1=8 e df2=8 summary(fmle) plot(fmle) gofstat(fmle) # Goodness-of-fit

Para se verificar se as variâncias são iguais ou não é necessário recorrer ao teste F (inferência para duas variâncias de populações normais).

\[ H_0: \sigma_A^2 = \sigma_A^2 \\ H_1: \sigma_A^2 \ne \sigma_A^2 \\ \]

ou seja, a hipótese nula frente à uma hipótese alternativa, que pode ser, da diferença, ou seja a razão de variâncias ser diferente de 1, ou superior ou inferior à unidade,

\[ H_0: \sigma_A^2 = \sigma_A^2 \\ \]

sendo a hipótese alternativa uma das seguintes opções:

\[ H_1: \sigma_A^2 \ne \sigma_A^2 \\ \space \\ H_1: \sigma_A^2 > \sigma_A^2 \\ \space \\ H_1: \sigma_A^2 < \sigma_A^2 \\ \space \\ \] que são equivalentes às hipóteses:

\[ H_1: \frac{\sigma_A^2}{\sigma_A^2} \ne 1 \\ \space \\ H_1: \frac{\sigma_A^2}{\sigma_A^2} > 1 \\ \space \\ H_1: \frac{\sigma_A^2}{\sigma_A^2} < 1 \\ \space \\ \]

A estatística do teste é dada por:

Teste para diferença de variâncias

O teste de hipóteses para a diferença de variâncias de duas distribuições normais.

Hipótese nula: hipótese de nulidade, onde não há diferença entre as variâncias.

$H_0: \frac{\sigma_A^2}{\sigma_B^2} = 1 \\$

Estatística do teste: o quão a razão de variâncias se afastam da unidade \[ F = \frac{S_A^2}{S_B^2}\\ \]

Hipóteses alternativas disponíveis a testar:

$H_1: \frac{\sigma_A^2}{\sigma_B^2} \ne 1\\$
$H_1: \frac{\sigma_A^2}{\sigma_B^2} > 1\\$
$H_1: \frac{\sigma_A^2}{\sigma_B^2} < 1\\$

Voltando ao caso das amostras recolhidas A e B, testaremos a hipótese que as variâncias são diferentes entre si.

$S^2_A =$ 7.3349517
$S^2_B =$ 8.3419714

Hipóteses

Parâmetro de interesse a ser testado: Se ${\sigma_A^2}\ne {\sigma_B^2}$, consequentemente $\frac{\sigma_A^2}{\sigma_B^2} = 1$
Hipótese nula: $H_0: \frac{\sigma_A^2}{\sigma_B^2} = 1$
Hipótese alternativa: $H_1: \frac{\sigma_A^2}{\sigma_B^2} \ne 1$
Estatística do teste:

Note que por simplicidade no uso das probabilidades da F tabelada opta-se por obter o maior valor de F calculado, neste caso utilizando no numerador a maior variância amostral da comparação.

\[ F = \frac{S_A^2}{S_B^2}\\ \]

Decisão do teste: Rejeita-se $H_0$ se o $Valor-P$ (probabilidade de cometer o erro do Tipo I) for menor que $\alpha$.

Erro do Tipo I: é a probabilidade em rejeitar a $H_0$ sendo ela verdadeira.

Cálculos:

\[ F = \frac{S_A^2}{S_B^2} = \frac{8.3419}{7.3349} = 1.1372 \\ \]

library(Plothtests) plot_F_test(alpha = 0.05,alternative = "two.sided", df1 = 14, v1 = 8.3419, df2 = 14, v2 = 7.3349)

$test.statistic [1] 1.137289 $df.num [1] 14 $df.den [1] 14 $p.value [1] 0.8131681

Interpretação

Note que o $Valor-P$ obtido foi de aproximadamente $0.8131681$ (área hachurada), teste bilateral, devido a $H_1$ ser postulada como a diferença. O $Valor-P$ indica a probabilidade de comentermos o Erro do Tipo I com base na amostra, ou seja, rejeitarmos a $H_0$ sendo ela verdadeira. Dessa forma a probabilidade de cometermos tal erro na tomada de decisão é demasiado elevado, superior ao valor arbitrário $\alpha = 0.05$. Podemos dizer que não temos evidência nos dados que demonstrem diferença significativa entre as variâncias das amostras de preparações A e B.