Para esta primeira aula prática de Estatística Não Paramétrica, precisaremos dos seguintes pacotes:
require(ggplot2) # Abordagem gráfica
require(ggpubr) # Abordagem gráfica
require(nortest) # Testes de normalidade
require(tseries) # Testes de normalidade
Para fazer o teste qui-quadrado de aderência no
R, basta utilizar a função
chisq.test() do pacote base stats. Essa função
tem os seguintes argumentos:
A função retorna as seguintes saídas:
Fãs de corridas de cavalos frequentemente sustentam que uma corrida em torno de uma pista circular proporciona significante vantagem inicial para os cavalos colocados em certas posições dos postos. Cada posição do cavalo corresponde ao posto atribuído no começo do alinhamento. Em uma corrida de oito cavalos, a posição 1 é a mais próxima da raia no lado interno da pista; a posição 8 está no lado externo, mais distante da raia. Podemos testar o efeito da posição do posto analisando os resultados da corrida, dados de acordo com a posição do posto, durante o primeiro mês de corridas na estação em uma pista circular particular: (Utilize um nível \(\alpha=5\%\))
A proporção esperada em cada casela é igual. Assim, a proporção esperada é dada por 1/8. Com os valores observados, podemos escrever o seguinte comando no R:
# vetor com valores observados
Observado = c(29,19,18,25,17,10,15,11)
# vetor com as proporções esperadas
Proporcao = rep(1/8,8)
# Teste qui-quadrado de aderência
Teste = chisq.test(x= Observado,p=Proporcao)
Teste
##
## Chi-squared test for given probabilities
##
## data: Observado
## X-squared = 16.333, df = 7, p-value = 0.02224
Assim, a estatística de teste é \(X^2=16.33\). Sob \(H_0\), ela possui uma distribuição qui-quadrado com k-1=8-1=7 graus de liberdade. Além disso, essa estatística apresenta um valor-p de 0,02224 e rejeitamos a hipótese nula de que a distribuição dos postos seja uniforme para a vitória, com um nível de significância de \(\alpha=5\%\).
Um engenheiro de controle de qualidade coletou 50 amostras de tamanho 13 cada uma de um processo de produção. O número de itens defeituosos para estas amostras são registradas abaixo. Teste a hipótese nula ao nível \(\alpha=5\%\) de que o número de itens defeituosos segue uma distribuição Poisson.
Como visto em sala, inicialmente estimamos o parâmetro \(\lambda\) a partir da média amostral. Obtivemos o valor \(\hat{\lambda}=1.3\). Assim, estamos verificando se os dados seguem uma distribuição Poisson com parâmetro \(\lambda=1.3\). Além disso, precisamos agrupar as duas últimas linhas da base para atender às suposições do teste.
Para este caso, temos o seguinte comando no R:
# Vetor com valores observados
Observado = c(10,24,10,4,2)
# Vetor com proporções esperadas (oriundas da Poisson(1.3))
# Atente-se para a última probabilidade que se torna (>=4)
# Tal fato se deve à necessidade desse vetor somar 1
Proporcao = c(dpois(0,1.3),
dpois(1,1.3),
dpois(2,1.3),
dpois(3,1.3),
1-ppois(3,1.3))
# Teste Qui-Quadrado de aderência
Teste = chisq.test(x=Observado,p=Proporcao)
Teste
##
## Chi-squared test for given probabilities
##
## data: Observado
## X-squared = 3.6019, df = 4, p-value = 0.4625
Este teste aponta um valor p incoerente, uma vez que é baseado em uma distribuição sob \(H_0\) que possui 4 graus de liberdade. Como vimos em sala, estes graus de liberdade devem ser corrigidos para \(k-1-1=5-1-1=3\). A última subtração se deve ao fato de termos estimado um parâmetro. Vamos corrigir o valor p com os seguintes códigos:
# Corrigindo o valor-p:
1-pchisq(Teste$statistic,3)
## X-squared
## 0.3077789
Com um valor-p de 0,3077, não rejeita-se a hipótese nula de que os dados são distribuídos de acordo com uma Poisson(\(\lambda=1.3\))
O teste Kolmogorov-Smirnov pode ser aplicada para
uma ou duas amostras. Em ambos os
casos, utilizamos a função
ks.test() do pacote base
stats do R. Essa função tem os seguintes argumentos:
Se o argumento y contém uma distribuição cumulativa, na sequência devem ser fornecidos os valores dos parâmetros.
A função ks.test apresenta as seguintes
saídas:
A seguir, segue um exemplo para o caso de uma amostra.
Seja uma amostra de tamanho n = 20 das distâncias até um dos lados para as quais uma linha de 6cm se quebra quando sujeita a uma força:
É razoável supor que esta amostra venha de uma distribuição Uniforme no intervalo (0,6)? Faça o teste de Kolmogorov-Smirnov ao nível de 5%.
Inicialmente, vamos fazer o plot da função distribuição acumulada destes dados e construir um intervalo de 95% confiança:
# Valores da amostra
Amostra = c(0.6, 0.8, 1.1, 1.2, 1.4, 1.7, 1.8, 1.9, 2.2, 2.4, 2.5, 2.9, 3.1, 3.4, 3.6, 3.9, 4.4, 4.9, 5.2, 5.9)
# Obtenção da função de distribuição amostral
Cumulative = ecdf(Amostra)
Valores = Cumulative(Amostra)
# Valor crítico da tabela 2A para uma confiança de 95% e n=20
Critico = 0.294
# Limites superiores para F(x)
Superior <- Valores + Critico
# Limites inferiores para F(x)
Inferior <- Valores - Critico
# Truncando limites inferiores menores que 0 e limites superiores maiores que 1
Superior[Superior > 1] <- 1
Inferior[Inferior < 0] <- 0
# Função de distribuição teórica
Uniforme<-punif(Amostra,0,6)
# Criar um data frame com os resultados
KS_df <- data.frame(
x = sort(Amostra),
S_x = Valores,
Inferior = Inferior,
Superior = Superior,
F_x = Uniforme
)
# Gráfico ggplot
ggplot(data=KS_df, aes(x = x, y = S_x)) +
# Função distribuição amostral
geom_step(color = "blue", size = 1) +
# Intervalo de confiança
geom_ribbon(aes(ymin = Inferior, ymax = Superior), color="blue",alpha = 0.2) +
# Função distribuição teórica
geom_line(data=KS_df, aes(x = x, y = F_x))+
# Especificações para os eixos
labs(title = "Função de distribuição acumulada",
x = "Valores",
y = "Frequência Acumulada") +
theme_minimal()
Nota-se que a função de distribuição teórica encontra-se completamente dentro das faixas de confiança. Agora, para a realização do teste, tem-se os seguintes comandos:
# Teste de Kolmogorov Smirnov para uma amostra
# punif é a acumulada da distribuição uniforme(0,6)
Teste<-ks.test(Amostra,punif,min=0,max=6,alternative = "two.sided")
Teste
##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: Amostra
## D = 0.15, p-value = 0.7045
## alternative hypothesis: two-sided
De acordo com a saída, a estatística de teste é D=0,15 e o P-valor é dado por 0,7045. Assim, não rejeita-se a hipótese nula de que os dados são distribuídos de acordo com uma uniforme(0,6) a um nível de 5% de significância.
A seguir, segue um exemplo para o caso de duas amostras:
Os dados a seguir são taxas de produção de saliva (ml/min) para 6 mulheres e 18 homens adultos, arranjados por conveniência em ordem crescente. É razoável concluir que as amostras provêm de populações com a mesma distribuição?
Mulheres: 0.45 0.60 0.80 0.85 0.95 1.75
Homens: 0.40 0.50 0.55 0.65 0.70 0.75 0.90 1.05 1.15 1.25 1.30 1.35 1.45 1.50 1.85 1.90 2.30 2.55
Utilize o teste de Kolmogorov Smirnov para duas amostras com um nível de significância de 5%.
# Valores da amostra
Amostra1 = c(0.45, 0.60, 0.80, 0.85, 0.95, 1.75)
Amostra2 = c(0.40, 0.50, 0.55, 0.65, 0.70, 0.75, 0.90, 1.05, 1.15, 1.25, 1.30, 1.35, 1.45, 1.50, 1.85, 1.90, 2.30, 2.55)
# Função distribuição da amostra 1
CumulativeA1 = ecdf(Amostra1)
ValoresA1 = CumulativeA1(Amostra1)
# Função distribuição da amostra 2
CumulativeA2 = ecdf(Amostra2)
ValoresA2 = CumulativeA2(Amostra2)
# Criar um data frame com os resultados
KS_df2 <- data.frame(
Amostra1 = sort(Amostra1),
ValoresA1,
Amostra2 = sort(Amostra2),
ValoresA2)
# Gráfico
ggplot() +
geom_step(data=KS_df2,aes(x = Amostra1, y = ValoresA1),color = "blue", size = 1)+
geom_step(data=KS_df2,aes(x = Amostra2, y = ValoresA2),color = "red", size = 1)+
labs(title = "Função de distribuição acumulada",
x = "Valores",
y = "Frequência Acumulada") +
theme_minimal()
Para realizar o teste de Kolmogorov Smirnov para duas amostras, utilize o seguinte comando:
# Teste de Kolmogorov Smirnov para duas amostras
Teste<-ks.test(Amostra1,Amostra2,alternative = "two.sided")
Teste
##
## Exact two-sample Kolmogorov-Smirnov test
##
## data: Amostra1 and Amostra2
## D = 0.44444, p-value = 0.3023
## alternative hypothesis: two-sided
Assim, de acordo com o p-valor, não rejeita-se a hipótese nula a um nível de 5% de significância. Isto é, há evidências de que as amostras são provenientes de uma mesma distribuição.
Para a condução dos testes de normalidade, utilizaremos as seguintes bases de dados:
O conjunto de dados ToothGrowth contém o resultado de um experimento que estuda o efeito da vitamina C no crescimento dentário em 60 porquinhos-da-índia. Cada animal recebeu um dos três níveis de dose de vitamina C (0,5, 1 e 2 mg/dia) por um dos dois métodos de administração (suco de laranja ou ácido ascórbico (uma forma de vitamina C e codificada como VC).
No R, essa base é carregada a partir dos comandos:
data("ToothGrowth")
ToothGrowth
## len supp dose
## 1 4.2 VC 0.5
## 2 11.5 VC 0.5
## 3 7.3 VC 0.5
## 4 5.8 VC 0.5
## 5 6.4 VC 0.5
## 6 10.0 VC 0.5
## 7 11.2 VC 0.5
## 8 11.2 VC 0.5
## 9 5.2 VC 0.5
## 10 7.0 VC 0.5
## 11 16.5 VC 1.0
## 12 16.5 VC 1.0
## 13 15.2 VC 1.0
## 14 17.3 VC 1.0
## 15 22.5 VC 1.0
## 16 17.3 VC 1.0
## 17 13.6 VC 1.0
## 18 14.5 VC 1.0
## 19 18.8 VC 1.0
## 20 15.5 VC 1.0
## 21 23.6 VC 2.0
## 22 18.5 VC 2.0
## 23 33.9 VC 2.0
## 24 25.5 VC 2.0
## 25 26.4 VC 2.0
## 26 32.5 VC 2.0
## 27 26.7 VC 2.0
## 28 21.5 VC 2.0
## 29 23.3 VC 2.0
## 30 29.5 VC 2.0
## 31 15.2 OJ 0.5
## 32 21.5 OJ 0.5
## 33 17.6 OJ 0.5
## 34 9.7 OJ 0.5
## 35 14.5 OJ 0.5
## 36 10.0 OJ 0.5
## 37 8.2 OJ 0.5
## 38 9.4 OJ 0.5
## 39 16.5 OJ 0.5
## 40 9.7 OJ 0.5
## 41 19.7 OJ 1.0
## 42 23.3 OJ 1.0
## 43 23.6 OJ 1.0
## 44 26.4 OJ 1.0
## 45 20.0 OJ 1.0
## 46 25.2 OJ 1.0
## 47 25.8 OJ 1.0
## 48 21.2 OJ 1.0
## 49 14.5 OJ 1.0
## 50 27.3 OJ 1.0
## 51 25.5 OJ 2.0
## 52 26.4 OJ 2.0
## 53 22.4 OJ 2.0
## 54 24.5 OJ 2.0
## 55 24.8 OJ 2.0
## 56 30.9 OJ 2.0
## 57 26.4 OJ 2.0
## 58 27.3 OJ 2.0
## 59 29.4 OJ 2.0
## 60 23.0 OJ 2.0
O conjunto precip, que contém dados de precipitação anual (em polegadas) em 70 cidades dos Estados Unidos. Para carregar estes dados, utilizaremos os comandos:
data("precip")
precip
## Mobile Juneau Phoenix Little Rock
## 67.0 54.7 7.0 48.5
## Los Angeles Sacramento San Francisco Denver
## 14.0 17.2 20.7 13.0
## Hartford Wilmington Washington Jacksonville
## 43.4 40.2 38.9 54.5
## Miami Atlanta Honolulu Boise
## 59.8 48.3 22.9 11.5
## Chicago Peoria Indianapolis Des Moines
## 34.4 35.1 38.7 30.8
## Wichita Louisville New Orleans Portland
## 30.6 43.1 56.8 40.8
## Baltimore Boston Detroit Sault Ste. Marie
## 41.8 42.5 31.0 31.7
## Duluth Minneapolis/St Paul Jackson Kansas City
## 30.2 25.9 49.2 37.0
## St Louis Great Falls Omaha Reno
## 35.9 15.0 30.2 7.2
## Concord Atlantic City Albuquerque Albany
## 36.2 45.5 7.8 33.4
## Buffalo New York Charlotte Raleigh
## 36.1 40.2 42.7 42.5
## Bismark Cincinnati Cleveland Columbus
## 16.2 39.0 35.0 37.0
## Oklahoma City Portland Philadelphia Pittsburg
## 31.4 37.6 39.9 36.2
## Providence Columbia Sioux Falls Memphis
## 42.8 46.4 24.7 49.1
## Nashville Dallas El Paso Houston
## 46.0 35.9 7.8 48.2
## Salt Lake City Burlington Norfolk Richmond
## 15.2 32.5 44.7 42.6
## Seattle Tacoma Spokane Charleston Milwaukee
## 38.8 17.4 40.8 29.1
## Cheyenne San Juan
## 14.6 59.2
Antes de realizar os testes de normalidade, é extremamente útil visualizar os dados a partir de gráficos como densidade e q-q plot:
Os comandos a seguir fazem o gráfico de densidade das duas bases de dados:
#Base de dados 1
ggdensity(ToothGrowth$len, fill = "lightgray")
ggdensity(precip, fill = "lightgray")
Os comandos a seguir fazem o gráfico q-q plot das duas bases de dados:
#Base de dados 1
ggqqplot(ToothGrowth$len)
ggqqplot(precip)
O teste de Shapiro Wilk é feito no R com a função
\(shapiro.test\) da base do
R. Para a utilização da função, basta colocar um
vetor com os dados da amostra. O teste solta a
estatística de teste e apresenta o valor-p.
# Base de porquinhos da índia
shapiro.test(ToothGrowth$len)
##
## Shapiro-Wilk normality test
##
## data: ToothGrowth$len
## W = 0.96743, p-value = 0.1091
# Base de dados precipitação
shapiro.test(precip)
##
## Shapiro-Wilk normality test
##
## data: precip
## W = 0.96456, p-value = 0.04493
Segundo o teste de Shapiro Wilk, os dados de crescimento dos dentes dos porquinhos-da-índia são normalmente distribuídos e os dados de precipitação não possuem distribuição normal.
O teste de Lilliefors é feito com a função lillie.test
do pacote nortest. Para a utilização da função, basta
colocar um vetor com os dados da amostra. O teste
solta a estatística de teste e apresenta o valor-p.
# Base de dados porquinhos da índia
lillie.test(ToothGrowth$len)
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: ToothGrowth$len
## D = 0.097092, p-value = 0.172
# Base de dados precipitação
lillie.test(precip)
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: precip
## D = 0.10909, p-value = 0.03812
Segundo o teste de Lilliefors, os dados de crescimento dos dentes dos porquinhos-da-índia são normalmente distribuídos e os dados de precipitação não possuem distribuição normal.
O teste de Anderson-Darling é feito com a função ad.test
do pacote nortest. Para a utilização da função, basta
colocar um vetor com os dados da amostra. O teste
solta a estatística de teste e apresenta o valor-p.
# Base de dados porquinhos da índia
ad.test(ToothGrowth$len)
##
## Anderson-Darling normality test
##
## data: ToothGrowth$len
## A = 0.64705, p-value = 0.08709
# Base de dados precipitação
ad.test(precip)
##
## Anderson-Darling normality test
##
## data: precip
## A = 0.99894, p-value = 0.01163
Segundo o teste de Anderson Darling, os dados de crescimento dos dentes dos porquinhos-da-índia são normalmente distribuídos e os dados de precipitação não possuem distribuição normal.
O teste de Jarque-Bera é feito com a função
jarque.bera.test do pacote tseries. Para a
utilização da função, basta colocar um vetor com os dados da
amostra. O teste solta a estatística de teste e
apresenta o valor-p.
#Base de dados 1
jarque.bera.test(ToothGrowth$len)
##
## Jarque Bera Test
##
## data: ToothGrowth$len
## X-squared = 2.5931, df = 2, p-value = 0.2735
jarque.bera.test(precip)
##
## Jarque Bera Test
##
## data: precip
## X-squared = 1.2692, df = 2, p-value = 0.5302
Segundo o teste de Jarque-Bera, ambas as bases de dados possuem distribuição normal.
Consulte a lista com exercícios práticos (1 a 8)