Inferência sobre variâncias caso de 1 e 2 populações

1 Introdução:

Nas seções seguintes veremos como utilizar o software R para

• construir testes de hipóteses para a variância de uma variável aleatória com distribuição Normal.

• Construir testes de hipóteses e intervalos de confiança para a comparação de variâncias de 2 distribuições Normais.

Para compreensão desta seção é necessário que você já tenha estudado o conteúdo referente a inferência sobre variâncias.

2 Caso de uma população

Para exemplificarmos o uso do R para a construção de inferências para uma variância populacional considere o seguinte problema.

Exemplo 1: Uma máquina é utilizada para encher garrafas d’água. Se a variância do volume de água das garrafas exceder \(0.01\), uma proporção inaceitável de garrafas estará ou muito cheia, ou muito vazia. Para uma amostra aleatória de 20 garrafas, encontrou-se um valor da variância do volume igual a \(s^2 = 0.0153\). Há evidências nos dados amostrais que sugerem um mal funcionamento na máquina de encher garrafas?

Para responder esta pergunta vamos construir um teste de hipóteses. Antes de mostrar como fazer isso no R, vamos relembrar o resultado teórico utilizado para construção das inferências sobre uma variância populacional quando a amostra observada provém de uma distribuição Normal.

Quando a amostra é formada por observações independentes de uma distribuição \(N(\mu,\sigma)\), temos que a estatistica \(\chi^2 = \frac{(n-1)S^2}{\sigma^2}\) possui distribuição Qui-Quadrado com \(n-1\) graus de liberdade. Observe que essa estatistica é função da variância populacional \(\sigma^2\), do seu estimador, a variância amostral \(S^{2}\) e do tamanho da amostra \(n\).

Retomando o exemplo, formulamos as hipóteses nula e alternativa para esta teste como: \[H_{0}: \sigma^{2}=0,01\] \[H_{a}: \sigma^{2}>0.01\]

Para testar essas hipóteses vamos utilizar a estatistica de teste \(\chi_0^2 = \frac{(n-1)S^2}{\sigma_0^2}\), onde \(\sigma_{0}^2\) é o valor de \(\sigma^2\) estabelecido na hipótese nula, e que possui distribuição Qui-Quadrado com \(n - 1 = 19\) graus de liberdade quando essa hipótese é verdadeira.

Para tomar nossa decisão quanto a aceitar ou rejeitar a hipótese nula, vamos calcular a estatistica do teste

est_teste <- 19*(0.0153)/0.01
est_teste

## [1] 29.07

Considerando nível de significância \(\alpha=0,05\), nós devemos rejeitar a hipótese nula quando o valor da estatistica de teste pertencer à região crítica do teste, isto é, quando ele se encontrar na região à direita do quantil de ordem \(0,95\) da sua distribuição quando a hipótese nula é verdadeira. O valor deste quantil é obtido no R como:

valor_critico <- qchisq(0.95, df=19)
valor_critico

## [1] 30.14353

Como pode ser visto no gráfico seguinte, o valor observado da estatistica deste (destacado pela linha vertical azul) encontra-se fora da região critica do teste (destacada em vermelho). Portanto, considerando nível de significãncia de 5%, concluimos os dados não contém evidências de que a hipótese nula \(H_{0}: \sigma^{2}=0,01\) seja falsa.

curve(dchisq(x,df=19),0, qchisq(0.999,19), xlab=bquote(chi^{2}),yaxs="i")
title(bquote("Distribuição da estatistica de teste quando"~sigma^{2}==0.01))
abline(v=est_teste, col=4)
a=qchisq(0.95,df=19)
b=qchisq(0.999,df=19)
polygon(x=c(a,seq(a,b,0.01),b),y=c(0,dchisq(seq(a,b,0.01),df=19),0),col=2)
mtext(side=1,at=a, round(a,2),col=2, cex=0.8)
text(x=est_teste,y=0.02,round(est_teste,2),col=4,cex=0.8)

Entretanto, a maneira mais usual de decidir sobre a rejeição ou não da hipótese nula é através do calculo do p-valor, que é a probabilidade da estatistica de teste assumir valores tão ou mais extremos que o seu valor observado na amostra. Para o exemplo acima, ele corresponde à àrea abaixo da curva e à direita da linha vertical azul. Lembre-se que sempre que o p-valor for menor do que o nível de significância \(\alpha\), nós decidimos a favor da hipótese alternativa. Para o exemplo, o p-valor é dado pela \(P(\chi^{2} \ge 29.07)\), e é obtido no R fazendo:

p_valor=pchisq(est_teste, df=19,lower.tail=FALSE)
p_valor

## [1] 0.064892

Observe que, como esperado, o p-valor igual a 0.064892 é maior do que o nível de significância fixado a priori, \(\alpha=0.05\).

O teste realizado acima pode ser feito no R usando a função varTest do pacote ENVStats. A utilização dessa função requer a especificação do vetor os valores observados na amostra para a variável de interesse. Considere o seguinte exemplo:

Exemplo 2: Um fabricante de capacetes de segurança para trabalhadores da construção está preocupado com a média e a variação das forças que seus capacetes transmitem aos usuários quando submetidos a uma força externa. O fabricante projetou os capacetes de modo que a força média transmitida pelos capacetes aos trabalhadores seja de 800 libras (ou menos) com um desvio padrão inferior a 40 libras. Testes foram executados em uma amostra aleatória de n = 40 capacetes, obtendo-se os seguintes valores para a força transmitida aos usuários:

813,32 804,24 809,02 855,08 728,26 857,22 790,03 830,69 824,92 892,53 837,84 900,60 807,73 802,75 846,06 792,61 709,28 750,61 747,40 770,41 822,09 841,61 885,86 772,67 865,12 789,18 921,31 864,74 768,36 767,30 825,17 777,22 771.39 794.24 818.10 845.36 760.91 801.58 855.05 783.22

Considerando nível de significância de 5%, os resultados do teste fornecem evidências para concluir que o desvio da força transmitida ao usuário para a população de capacetes produzidos pela empresa é diferente 40 libras?

Para testar as hipóteses \(H_{0}: \sigma^{2}=40^2 \quad \times \quad H_{a}: \sigma^{2} \neq 40^2\), primeiro carregamos o pacote EnvStats de depois usamos a função VarTest, especificando o vetor com as observações amostrais e o valor da variância populacional estabelecido na hipótese nula usando o argumento sigma.squared, como segue:

library(EnvStats)

## Warning: package 'EnvStats' was built under R version 4.0.5

## 
## Attaching package: 'EnvStats'

## The following objects are masked from 'package:stats':
## 
##     predict, predict.lm

## The following object is masked from 'package:base':
## 
##     print.default

x=c(813.32, 804.24, 809.02, 855.08, 728.26, 857.22, 790.03, 830.69, 824.92, 892.53, 837.84, 900.60, 807.73, 802.75, 846.06, 792.61, 709.28, 750.61, 747.40, 770.41, 822.09, 841.61, 885.86, 772.67, 865.12, 789.18, 921.31, 864.74, 768.36, 767.30, 825.17, 777.22, 771.39, 794.24, 818.10, 845.36, 760.91, 801.58, 855.05, 783.22)
varTest(x,alternative="two.sided", sigma.squared=1600)

## $statistic
## Chi-Squared 
##    55.43216 
## 
## $parameters
## df 
## 39 
## 
## $p.value
## [1] 0.08499349
## 
## $estimate
## variance 
##  2274.14 
## 
## $null.value
## variance 
##     1600 
## 
## $alternative
## [1] "two.sided"
## 
## $method
## [1] "Chi-Squared Test on Variance"
## 
## $data.name
## [1] "x"
## 
## $conf.int
##      LCL      UCL 
## 1526.004 3749.482 
## attr(,"conf.level")
## [1] 0.95
## 
## attr(,"class")
## [1] "htestEnvStats"

Para um teste bilateral, o p-valor é obtido como:

• \(2 \times P(\chi_{0}^{2} \leq A)\), se \(A\) é menor do que (n-1), o valor esperado da estatistica de teste quando a hipótese nula é verdadeira,

• \(2 \times P(\chi_{0}^{2} \geq A)\), se caso contrário.

Para o exemplo, o p-valor foi obtido como resultado do cálculo da probabilidade: \(2 \times p(\chi_{0}^{2} \geq 55.543)\), cujo valor retornado pela função varTest, igual a \(0.08499\), nos leva à aceitação da hipótese nula de que o desvio padrão populacional, \(\sigma\), é igual a 40 libras (\(\sigma^2=40^2 \space libras^2\)). Observe também que essa função retorna o intervalo de \(100(1-\alpha)\%\) de confiança para a variância populacional \(\sigma^{2}\). Esse intervalo é obtido como:
\[IC_{\sigma^{2}}^{100(1 − \alpha)\%} = \left[\frac{(n-1)S^{2}}{\chi^{2}_{1-\alpha/2}};\frac{(n-1)S^{2}}{\chi^{2}_{\alpha/2}}\right]\]
onde \(\chi^{2}_{\alpha/2}\) e \(\chi^{2}_{1 - \alpha/2}\) são respectivamente os quantis de ordem \(\alpha/2\) e \(1 - \alpha/2\) da distribuição Qui-Quadrado com n - 1 graus de liberdade.

Para o exemplo, temos que, com \(95\%\) de confiança, a variância da força transmitida pelo capacete ao usuario está entre \(1526,00\) e \(3749,48 \space libras ^2\). Consequentente o desvio padrão está entre \(\sqrt{1526,00} = 39,06\) e \(\sqrt{3749,48} = 61,23 \space libras\).

Quando o teste realizado é unilateral, a função varTest retorna o intervalo unilateral de confiança. No caso do teste unilateral esquerdo, alternative = "less", o limite inferior de confiança é igual a zero e o limite superior é calculado como no caso do intervalo bilateral, porém substituindo \(\chi^{2}_{\alpha/2}\) por \(\chi^{2}_{\alpha}\). No caso do teste unilateral direito, alternative = “greater”, o limite superior de confiança é igual a \(\infty\) e o limite inferior é calculado como no caso do intervalo bilateral, porém substituindo \(\chi^{2}_{1 - \alpha/2}\) por \(\chi^{2}_{1 - \alpha}\).

A seguir apresentamos os resultados do teste realizado com os dados do exemplo 2 para as hipóteses \(H_{0}: \sigma^{2}=40^2 \quad \times \quad H_{a}: \sigma^{2} > 40^2\):

varTest(x,alternative="greater", sigma.squared=1600)

## $statistic
## Chi-Squared 
##    55.43216 
## 
## $parameters
## df 
## 39 
## 
## $p.value
## [1] 0.04249674
## 
## $estimate
## variance 
##  2274.14 
## 
## $null.value
## variance 
##     1600 
## 
## $alternative
## [1] "greater"
## 
## $method
## [1] "Chi-Squared Test on Variance"
## 
## $data.name
## [1] "x"
## 
## $conf.int
##      LCL      UCL 
## 1625.212      Inf 
## attr(,"conf.level")
## [1] 0.95
## 
## attr(,"class")
## [1] "htestEnvStats"

3 Caso de duas populações:

Nessa seção exemplificarmos o uso do R para construção de inferências para comparação de variâncias de 2 populações Normais. Mais especificamente, veremos como realizar o teste F para comparação de 2 variâncias e também como construir intervalos de confiança para a razão de 2 variâncias populacionais. Para realização do teste F, fazemos os seguintes pressupostos:

• \(X_{1},X_{2},....,X_{n}\) é uma amostra aleatória de uma população Normal com media \(\mu_{x}\) e variância \(\sigma_{x}^{2}\),

• \(Y_{1},Y_{2},....,Y_{m}\) é uma amostra aleatória de uma população Normal com media \(\mu_{y}\) e variância \(\sigma_{y}^{2}\),

• As 2 amostras são independentes.

Quando estes pressupostos são verdadeiros, a estatistica \(F =\frac{S_{x}^{2}/\sigma_{x}^{2}}{S_{y}^2/\sigma_{y}^{2}}\) possui distribuição F com \({n-1}\) graus de liberdade realtivos ao numerador e \({m-1}\) graus de liberdade relativos ao denominador. Desse resultado derivam o teste F para comparação das 2 variâncias populacionais \(\sigma_{x}^{2}\) e \(\sigma_{y}^{2}\), e também o intervalo de confiança para a razão entre elas.

Para testar a hipotese nula \(H_{0}:\sigma_{x}^{2}=\sigma_{y}^{2}\), que pode ser escrita como \(H_{0}:\frac{\sigma_{x}^{2}}{\sigma_{y}^{2}}=1\), contra uma hipótese alternativa, bilateral ou unilateral, utilizamos a estatistica \(F_{0}=\frac{S_{x}^{2}}{S_{y}^{2}}\), que possui distribuição \(F\) com \(n-1\) graus de liberdade para o numerador e \(m-1\) graus de liberdade para o denominador quando a hipótese nula \(H_{0}\) é verdadeira.

No R, o teste F para a comparação de 2 variâncias pode ser realizado usando a função var.test, do pacote stats. Para exemplificar o uso dessa função, considere o seguinte exemplo:

Exemplo 3: Uma construtora recebe tubos de aço de dois fabricantes diferentes. Para avaliar se há diferenças entre as variâncias da densidade do tubos, medida em \(kg/m\), produzidos pelos 2 fornecedores, uma amostra de 10 tubos de cada fornecedor foi selecionada e as densidades dos tubos medidas. Os resultados obtidos foram:

Fabricante 1: 7,7 7,8 7,9 8,1 7,8 8,1 8,0 7,9 7,8 8,0

Fabricante 2: 7,9 8,1 8,0 7,8 7,7 7,9 8,1 8,0 8,0 8.1

Os valores observados nas 2 amostras indicam que as variâncias das densidades dos tubos produzidos pelos 2 fabricantes são diferentes?

Para responder essa pergunta vamos realizar o teste F para as hipóteses descritas a seguir, considerando nível de significância \(\alpha = 0,05\):

\[H_{0}:\frac{\sigma_{1}^{2}}{\sigma_{2}^{2}} = 1 \quad \quad \quad H_{a}:\frac{\sigma_{1}^{2}}{\sigma_{2}^{2}} \neq 1\]

onde \(\sigma_{1}^{2}\) e \(\sigma_{2}^{2}\) são respctivamente as variâncias da densidade dos tubos frabricados, respectivamente pelos fornecedores 1 e 2.

Ao usarmos a funçao var.test é necessário especificarmos nos argumentos x e y os valores das 2 amostras. Para o exemplo, fazemos:

#definindo variáveis
fabricante1 = c(7.7,7.8,7.9,8.1,7.8,8.1,8,7.9,7.8,8)
fabricante2 = c(7.9,8.1,8,7.8,7.7,7.9,8.1,8,8,8.1)
var.test(fabricante1,fabricante2,alternative="two.sided")

## 
##  F test to compare two variances
## 
## data:  fabricante1 and fabricante2
## F = 1.0305, num df = 9, denom df = 9, p-value = 0.9651
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.2559586 4.1487379
## sample estimates:
## ratio of variances 
##           1.030488

Ao invés de especificarmos os vetores \(x\) e \(y\) contendo os valores das 2 amostras, podemos especificá-los usando a fómula \(y \sim x\), onde \(y\) é o vetor contendo os valores das 2 amostras e \(x\) é um fator que identifica as 2 amostras. Para os dados do Exemplo 3, fazemos:

densidade=c(fabricante1, fabricante2)
fabricante=rep(1:2, each=10)
var.test(densidade~factor(fabricante))

## 
##  F test to compare two variances
## 
## data:  densidade by factor(fabricante)
## F = 1.0305, num df = 9, denom df = 9, p-value = 0.9651
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.2559586 4.1487379
## sample estimates:
## ratio of variances 
##           1.030488

Observe que o R retorna o valor da estatística \(F_{0}\), \(f_{0}=1,0305\), que está muito próximo de \(1\); os graus de liberdade associados ao seu numerador e ao seu denominador, ambos iguais a \(9\); além do p-valor igual a \(0,9651\). Portanto, como o p-valor é maior do que \(\alpha = 0,05\) concluímos que não existem diferenças significativas entre os 2 fornecedores quanto à variância da densidade dos tubos por eles fabricados.

Para um teste bilateral o p-valor corresponde a \(2 \times P(F_{0} \geq f_{0})\), se a razão \(\frac{s_{x}^{2}}{s_{y}^{2}}\) é maior do que 1, e à \(2 \times P(F_{0} \leq f_{0})\), se \(\frac{s_{x}^{2}}{s_{y}^{2}}\) é menor do que 1. Para o teste unilateral esquerdo ele é dado por \(P(F_{0}\leq f_{0})\) e para o teste unilateral direito pela probabildade \(P(F_{0}\geq f_{0})\). Para o exemplo, como o valor da estatísica de teste é maior do que 1, o p-valor do teste pode ser obtido, sem fazer uso da fução var.test como:

p_valor = 2*pf(var(fabricante1)/var(fabricante2),df1=length(fabricante2)-1,df2=length(fabricante1)-1, lower.tail=FALSE)
p_valor

## [1] 0.9650511

A função var.test retorna também o intervalo de confiança para a razão entre as 2 variâncias populacionais, obtido através da expressão:

\[IC_{\frac{\sigma_{x}^{2}}{\sigma_{y}^{2}}}{100(1 − \alpha)\%} = \left[ \frac{1}{F_{{n-1 , m-1};1-\alpha/2}}\cdot\frac{s_{x}^{2}}{s_{y}^{2}} \text{ };\text{ } \frac{1}{F_{{n-1 , m-1};\alpha/2}}\cdot\frac{s_{x}^{2}}{s_{y}^{2}} \right]\]

onde \(F_{n-1 , m-1;\alpha/2}\) e \(F_{n-1 , m-1;1-\alpha/2}\) são respectivamente os quantis de ordem \(\alpha/2\) e \(1 - \alpha/2\) da distribuição F com \(n -1\) graus de liberdade para o numerador e com \(m -1\) graus de liberdade para o denominador.

Para o exemplo, temos que com \(95\%\) de confiança a razão entre as variâncias das densidades dos tubos de aço fabricados pelos forncedores 1 e 2 está entre \(0,26\) e \(4,15\). Este intervalo inclui 1, o valor esperado da razão das variâncias estabelecido na hipótese nula, concordando, como esperado, com o resultado do teste de hipóteses realizado ao nível de siginficância de \(5\%\).

No argumento conf.level da função var.test especificamos o coeficiente de confiança desejado para o intervalo de confiança, que é assumido igual a \(95\%\) (con.level=0.95), quando não for especificado.

Ao invés de usarmos a função var.test, nós podemos escrever o código R para a obtenção dos limites de confiança. Para o exemplo, o intervalo de \(95\%\) de confiança para a razão das variâncias \(\frac{\sigma_{1}^{2}}{\sigma_{1}^{2}}\) é obtido fazendo:

razao=var(fabricante1)/var(fabricante2)
# LI = Limite Inferior   LS = Limite Superior
LI=razao/qf(0.975, df1=length(fabricante1)-1,df2=length(fabricante2)-1)
LS=razao/qf(0.025, df1=length(fabricante1)-1,df2=length(fabricante2)-1)
resultado=c(razao,LI,LS)
resultado

## [1] 1.0304878 0.2559586 4.1487379

Quando os testes realizados forem unilaterais, a função var.test retorna intervalos de confiança unilaterais obtidos usando procedimento similar ao usado para obtençao de intervalos de confiança unilaterais para uma variância populacional.

O teste F para comparação de 2 variâncias não é robusto aos afastamentos do pressuposto de que as amostras provém de distribuições Normais. Outros testes para comparação de variâncias são disponíveis no R: teste de Bartlett, teste de Levene e teste de Fligner-Killeen, estes 2 últimos mais robustos aos afastamentos do pressuposto de normalidade. Os 3 últimos desses testes podem ser realizados respectivamente com as funções bartlett.test (pacote stats), leveneTest (pacote car) e fligner.test (pacote stats) quando se deseja comparar as variâncias de 2 ou mais populações.