A maior rede de estudos do Brasil

Grátis
224 pág.
Estatistica usado o R

Pré-visualização | Página 35 de 50

teste t. O banco
energy contém o gasto energético de mulheres obesas (obese) e magras (lean), em Megajoules
(MJ).
data(energy)
attach(energy)
Verifique a disposição dos dados nesse data frame. Repare que nós temos uma variável para
os valores e outra para a classificação dos grupos (obese e lean). Vamos então usar a nossa função
t.test() nesse banco:
> t.test(expend~stature)
 Welch Two Sample t-test
data: expend by stature 
6
t = -3.8555, df = 15.919, p-value = 0.001411
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -3.459167 -1.004081 
sample estimates:
 mean in group lean mean in group obese 
 8.066154 10.297778 
Repare que nós usamos o til para indicar que a variável expend está sendo classificada pela
variável stature, da mesma maneira que nós usamos no boxplot, lembra?
Repare também que o default do R é usar mesmo o procedimento de Welch, para variâncias
heterogêneas. As únicas diferenças dessa saída para a que nós já vimos anteriormente são:
data: expend by stature 
Indicando que é uma variável classificada pela outra.
t = -3.8555, df = 15.919, p-value = 0.001411
Os graus de liberdade que devem ser calculados segundo aquela equação lá em cima.
alternative hypothesis: true difference in means is not equal to 0 
O teste que nesse caso é a diferença das médias ser ou não igual a zero
sample estimates:
 mean in group lean mean in group obese 
 8.066154 10.297778 
E duas estimativas de médias, uma para cada grupo a ser comparado.
Vamos agora ver como o R calcula esse teste de Welch. Bem, primeiramente vamos calcular
os graus de liberdade, segundo a equação acima:
s1 <- var(expend[stature=="obese"])
s2 <- var(expend[stature=="lean"])
n1 <- length(expend[stature=="obese"])
n2 <- length(expend[stature=="lean"])
gl <- ((s1/n1)+(s2/n2))^2 / (((1/(n1-1))*(s1/n1)^2) + ((1/(n2-
1))*(s2/n2)^2))
gl
Veja se não bate com o resultado lá de cima... Se quiser pode fazer isso tudo na mão
também, mas eu recomendaria que você tentasse identificar a equação acima nesse código. Você
seria capaz de explicar o que foi feito antes de calcularmos os graus de liberdade propriamente
ditos?
Bem, agora falta ainda calcular a estatística T. Mas isso é fácil, também baseado na fórmula
mais geral lá de cima:
media1 <- mean(expend[stature=="obese"])
media2 <- mean(expend[stature=="lean"])
T <- (media1-media2)/sqrt((s1/n1)+(s2/n2))
T
Algo de errado com o valor de T? Está com o sinal trocado? Isso faz alguma diferença?
Como você obteria o resultado com o mesmo sinal?
7
Bom, já que o resultado bate com o esperado, agora ficou fácil obter o p-valor no R, não? É
a mesma coisa que para uma amostra:
> 2*pt(T, df=gl, lower.tail=F)
[1] 0.001410692
Ora, não é que bateu direitinho... Ah, se fosse fazer com o resultado de T negativo,
precisaríamos omitir o último argumento:
> 2*pt(-T, df=gl)
[1] 0.001410692
Claro que podemos também testar as nossas hipóteses baseado nesse mesmo teste. Para isso,
vamos primeiro estabelecer:
H0: 1−2=0
H1: 1−2≠0
E agora vamos usar o mesmo raciocínio que usamos na última aula: vamos calcular um
valor crítico para essa nossa t sob a hipótese nula. Lembra como é? Vamos lá:
> qt(0.025, df=gl)
[1] -2.120785
> qt(0.025, df=gl, lower.tail=F)
[1] 2.120785
Lembrou? A diferença é que tivemos que usar os graus de liberdade corrigidos. Como a
nossa T é maior que 2.12 (ou -T é menor que -2.12, não importa), vamos rejeitar a hipótese nula,
para um alfa de 5%, corroborando o resultado que obtivemos através do p-valor.
Ah, já ia me esquecendo. Tem ainda o IC 95% para a diferença dessas médias. Vamos usar
o nosso erro-padrão da diferença para a nossa notação (o denominador da equação lá em cima). É
assim:
IC 95% =  x1− x2− t ' ,0.975 sdiff , x1− x2 t ' ,0.975 sdiff
Isso para o nosso procedimento de Welch, é claro – por isso estamos usando sdif e  ' .
No R:
> media1-media2-qt(0.975, df=gl)*sqrt((s1/n1)+(s2/n2))
[1] 1.004081
> media1-media2+qt(0.975, df=gl)*sqrt((s1/n1)+(s2/n2))
[1] 3.459167
Bateu? Não? Ihhhh, esse sinal...
Mas é claro que o R permite que você use também o método clássico para variâncias
supostamente iguais. Basta um argumento para isso:
> t.test(expend~stature, var.equal=T)
 Two Sample t-test
data: expend by stature 
t = -3.9456, df = 20, p-value = 0.000799
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -3.411451 -1.051796 
sample estimates:
 mean in group lean mean in group obese 
8
 8.066154 10.297778 
Repare que agora houve uma mudança em uma das saídas:
t = -3.9456, df = 20, p-value = 0.000799
Todos esses valores são diferentes, pois T foi calculado a partir da variância conjunta, os
graus de liberdade da t são um número inteiro (usual) e iguais ao tradicional n1n2−2 e o p-
valor é evidentemente calculado a partir desta distribuição e deste valor de T e portanto também é
bastante diferente. Nossa conclusão sobre o teste de hipóteses não mudou, porém.
Já sei: você está achando que eu vou repetir toda aquela contaria que eu fiz para o
procedimento de Welch agora para o método clássico, não é? Se enganou, vai ficar para um longo
exercício...
Comparação de variâncias
Bem, de qualquer maneira, precisamos ter um meio de testar se as variâncias são
homogêneas ou não, para podermos usar o teste clássico. O R oferece como opção o teste F para
razão de variâncias. O procedimento no R é bem simples:
> var.test(expend~stature)
 F test to compare two variances
data: expend by stature 
F = 0.7844, num df = 12, denom df = 8, p-value = 0.6797
alternative hypothesis: true ratio of variances is not equal to 1 
95 percent confidence interval:
 0.1867876 2.7547991 
sample estimates:
ratio of variances 
 0.784446 
Esse teste é bem intuitivo e ele testa se a razão das variâncias das amostras é ou não igual a
1 (caso em que elas seriam iguais, obviamente.) Ele é baseado no fato de uma razão de variâncias
seguir uma distribuição F de Snedecor. Isso não é difícil de entender se você souber que na verdade
a distribuição F é uma razão de duas distribuições gl2 divididas pelos seus respectivos graus de
liberdade. Além disso, a distribuição F tem como parâmetros justamente os graus de liberdade do
numerador e do denominador, que correspondem aos graus de liberdade das gl2 do numerador e
do denominador dessa razão, respectivamente:
Fgl 1 , gl 2=
gl1
2
gl1
gl2
2
gl2
Vamos ver como isso funciona então. Como você deve se lembrar, a variância da amostra
segue uma distribuição n −12 , se for feito um pequeno ajuste. Lembra?
n−1× s2
 2
~n −12
Então, para comparar duas variâncias, podemos fazer o seguinte: vamos comparar as
variâncias s1
2 e s2
2 de duas amostras, inicialmente fazendo uma divisão de duas expressões
como a acima:
9
n1−1× s12
12
n2−1× s22
 22
Bem, por enquanto, nada podemos dizer desta expressão acima, apenas que representa a
divisão de duas distribuições Qui-quadradas. Mas se nós dividirmos o numerador por n1-1 e o
denominador por n2-1, que são os graus de liberdade dessas distribuições, teremos o que estávamos
procurando, ou seja, uma Fn 1−1, n 2−1 :
s1
2 /12
s2
2 / 22
~F n 1−1, n 2−1
Mas há um outro detalhe: o tipo de teste de hipóteses que vamos fazer neste caso é seguinte:
H0:  12 / 21=1
H1:  12 / 21≠1
Ora, isso quer dizer que sob a hipótese nula, 12= 22 e podemos então cancelá-los na
expressão acima. Sendo assim, a razão entre s1
2 e s2
2 segue uma distribuição Fn 1−1 , n 2−1 .
No nosso caso:
> s2/s1
[1] 0.784446
O mesmo resultado da saída lá em cima. Mas e o nosso p-valor? Como ele foi calculado?
Antes disso, vamos ver o jeitão dessa distribuição F
curve(df(x, df1=n2-1, df2=n1-1), from=0, to=5)
abline(v=1)
Repare que nós estamos testando uma igualdade com a unidade e não com zero. Bem, se a
estatística obtida for menor que 1, estaremos procurando um p-valor que vai desde zero (limite
inferior