Estatística Básica - Bussab e Morettin (2017)
Comandos R para análises estatísticas
Capítulo 14: Análise de Aderência e Associação
Exemplo 14.1
Para realizar o teste de Qui-quadrado no R, podemos calcular o valor da estatística \(\chi_{obs}^2\) e compará-la com o valor do quantil da distribuição de Qui-quadrado:
f_obs = c(43,49,56,45,66,41)
f_esp = c(50,50,50,50,50,50)
q_obs = sum((f_obs-f_esp)^2)/50
q_c=qchisq(0.95,5)
if(q_obs > q_c){
print("Rejeito H0")
}else print("Aceito H0")## [1] "Aceito H0"Uma outra forma de realizar este teste é utilizando a função chisq.test(), que calcula o teste para tabelas de contingência:
chisq.test(f_obs)##
## Chi-squared test for given probabilities
##
## data: f_obs
## X-squared = 8.96, df = 5, p-value = 0.11Como não há um segundo fator na tabela, basta declarar na função o vetor de frequências observadas. Se houvessem dois fatores, deveríamos primeiro construir tabela de frequências cruzada e depois calcular o teste, como veremos no exemplo a seguir.
Exemplo 14.2
c_hum<-c(15,20,30,20,15)
c_bio<-c( 8,23,18,34,17)
tab14_2<-rbind(c_hum,c_bio)
test_tab14_2=chisq.test(tab14_2)
test_tab14_2##
## Pearson's Chi-squared test
##
## data: tab14_2
## X-squared = 9.09, df = 4, p-value = 0.059A Tabela 14.8 apresenta as frequências absolutas sob a suposição de homogeneidade dos dois cursos, i.e. sob a hipótese:
\[ H_0: P_1=P_2 \]
Esta tabela é gerada automaticamente ao testarmos a homogeneidade através da função chisq.test como já fizemos acima. Desta forma, temos:
tab14_8 = rbind(test_tab14_2$expected,
total=apply(test_tab14_2$expected,2,sum))| A | B | C | D | E | |
|---|---|---|---|---|---|
| c_hum | 11.5 | 21.5 | 24 | 27 | 16 |
| c_bio | 11.5 | 21.5 | 24 | 27 | 16 |
| total | 23.0 | 43.0 | 48 | 54 | 32 |
Exemplo 14.3
M<-data.frame(uso_hospital=c("usaram_hospital", "nao_usaram_hospital"),homens=c(100,900),mulheres=c(150,850), row.names = 1)
M## homens mulheres
## usaram_hospital 100 150
## nao_usaram_hospital 900 850Exemplo 14.4
A Tabela 14.5 é baseada nas frequências observada e esperada. Para calcular a frequência esperada, novamente utilizaremos a função que calcula o teste de qui-quadrado, chisq.test:
Tabela 14.5
acidentes<-c(32,40,20,25,33)
test_tab14_5<-chisq.test(acidentes)
q_i<-(acidentes-test_tab14_5$expected)^2/test_tab14_5$expected
tab14_5<- rbind(acidentes, e_i=test_tab14_5$expected, q_i)
tab14_5<-cbind(tab14_5,c(sum(tab14_5[1,]),sum(tab14_5[2,]),sum(tab14_5[3,])))| Seg | Ter | Qua | Qui | Sex | Total | |
|---|---|---|---|---|---|---|
| acidentes | 32.00000 | 40.0000 | 20.0000 | 25.00000 | 33.0 | 150.0000 |
| e_i | 30.00000 | 30.0000 | 30.0000 | 30.00000 | 30.0 | 150.0000 |
| q_i | 0.13333 | 3.3333 | 3.3333 | 0.83333 | 0.3 | 7.9333 |
O resultado do teste que já fora calculado é:
test_tab14_5##
## Chi-squared test for given probabilities
##
## data: acidentes
## X-squared = 7.93, df = 4, p-value = 0.094Exemplo 14.5
Retornando ao Exemplo 6. 17, temos a Tabela 6.13:
| x | n_k | np_k |
|---|---|---|
| 0 | 57 | 54.399 |
| 1 | 203 | 210.523 |
| 2 | 383 | 407.361 |
| 3 | 525 | 525.496 |
| 4 | 532 | 508.418 |
| 5 | 408 | 393.515 |
| 6 | 273 | 253.817 |
| 7 | 139 | 140.325 |
| 8 | 45 | 67.882 |
| 9 | 27 | 29.189 |
| 10 | 16 | 17.075 |
Escreveremos uma função para calcular o teste com dois vetores (esperado e observado):
teste_quiquad<-function(f_obs,f_esp,alpha,gl){
q_obs=sum((f_obs-f_esp)^2/f_esp)
q_c=qchisq(1-alpha,df=gl)
p_val=pchisq(q_obs,df=gl,lower.tail = FALSE)
writeLines(" *** Teste de Qui-quadrado *** ")
writeLines(paste(" Estatística observada: ",round(q_obs,2)," ",sep=""))
writeLines(paste(" Graus de liberdade: ", gl," ",sep=""))
writeLines(paste(" Valor crítico: ",round(q_c,2)," ",sep=""))
writeLines(paste(" Valor-p: ",round(p_val,4)," ",sep=""))
writeLines(paste("Nível de significância: ",alpha," ",sep=""))
if(q_obs>q_c){
writeLines(" **Rejeito H0**")
} else writeLines(" **Aceito H0**")
}
teste_quiquad(n_k,np_k,0.05,10)## *** Teste de Qui-quadrado ***
## Estatística observada: 12.88
## Graus de liberdade: 10
## Valor crítico: 18.31
## Valor-p: 0.2302
## Nível de significância: 0.05
## **Aceito H0**Exemplo 14.6
dados.ex14_6<-c( 1.04, 1.73, 3.93, 4.44, 6.37, 6.51,
7.61, 7.64, 8.18, 8.48, 8.57, 8.65,
9.71, 9.87, 9.95, 10.01, 10.52, 10.69,
11.72, 12.17, 12.61, 12.98, 13.03, 13.16,
14.11, 14.60, 14.64, 14.75, 16.68, 22.14)Para categorizar os dados do exemplo, faremos da seguinte forma:
classe<-numeric()
classe[1]<-sum(dados.ex14_6<6.63)
classe[2]<-sum(dados.ex14_6 > 6.63 & dados.ex14_6 < 10)
classe[3]<-sum(dados.ex14_6 > 10 & dados.ex14_6 < 13.37)
classe[4]<-sum(dados.ex14_6 > 13.37)
test_tab14_6<- chisq.test(classe)Tabela 14.6
tab14_6<-cbind(rbind(classe,e_i=test_tab14_6$expected),c(sum(classe),sum(test_tab14_6$expected)))
print(kable(tab14_6,caption="**Tabela 14.6**: Valores observados e esperados para os dados, sob a suposição de normalidade.", col.names = c("(-Inf;6.63]","(6.63;10]","(10;13.37]","(13.37;+Inf)","Total")))| (-Inf;6.63] | (6.63;10] | (10;13.37] | (13.37;+Inf) | Total | |
|---|---|---|---|---|---|
| classe | 6.0 | 9.0 | 9.0 | 6.0 | 30 |
| e_i | 7.5 | 7.5 | 7.5 | 7.5 | 30 |
Novamente, utilizando a função teste_quiquad, temos:
teste_quiquad(classe,test_tab14_6$expected,0.05,3)## *** Teste de Qui-quadrado ***
## Estatística observada: 1.2
## Graus de liberdade: 3
## Valor crítico: 7.81
## Valor-p: 0.753
## Nível de significância: 0.05
## **Aceito H0**chisq.test(classe)##
## Chi-squared test for given probabilities
##
## data: classe
## X-squared = 1.2, df = 3, p-value = 0.75Que é equivalente ao resultado do teste utilizando a função chisq.test.
Figura 14.1
par(mfrow=c(1,2))
hist(dados.ex14_6,col="lightblue3",border = "white",xlab="Dados Normal")
boxplot(dados.ex14_6, xlab="Dados Normal",col="lightblue", pch=16, border="darkgray")
Figura 14.1: Histograma e boxplot para os dados do Exemplo 14.6. R.
Exemplo 14.7
tab14_10<-rbind(P_1T=c(29,60,9,2),P_2C=c(37,44,13,6))| (0.4;1.0] | (1.0;1.6] | (1.6;2.2] | (2.2;2.8] | |
|---|---|---|---|---|
| P_1T | 29 | 60 | 9 | 2 |
| P_2C | 37 | 44 | 13 | 6 |
| total | 66 | 104 | 22 | 8 |
chisq.test(tab14_10)## Warning in chisq.test(tab14_10): Chi-squared approximation may be incorrect##
## Pearson's Chi-squared test
##
## data: tab14_10
## X-squared = 6.16, df = 3, p-value = 0.1Note que o aviso “Chi-squared approximation may be incorrect” aparece por conta de termos uma das caselas menor que 5.
Exemplo 14.8
# Dados das tabelas 4.8 e 4.9
library(gmodels)
dados.tab4_8<-data.frame(rbind(
matrix(rep(c("1.Consumidor","1.São Paulo"),times=214),ncol=2,byrow=T),
matrix(rep(c("1.Consumidor","2.Paraná"),times=51),ncol=2,byrow=T),
matrix(rep(c("1.Consumidor","3.Rio G. do Sul"),times=111),ncol=2,byrow=T),
matrix(rep(c("2.Produtor","1.São Paulo"),times=237),ncol=2,byrow=T),
matrix(rep(c("2.Produtor","2.Paraná"),times=102),ncol=2,byrow=T),
matrix(rep(c("2.Produtor","3.Rio G. do Sul"),times=304),ncol=2,byrow=T),
matrix(rep(c("3.Escola","1.São Paulo"),times=78),ncol=2,byrow=T),
matrix(rep(c("3.Escola","2.Paraná"),times=126),ncol=2,byrow=T),
matrix(rep(c("3.Escola","3.Rio G. do Sul"),times=139),ncol=2,byrow=T),
matrix(rep(c("4.Outras","1.São Paulo"),times=119),ncol=2,byrow=T),
matrix(rep(c("4.Outras","2.Paraná"),times=22),ncol=2,byrow=T),
matrix(rep(c("4.Outras","3.Rio G. do Sul"),times=48),ncol=2,byrow=T)))
colnames(dados.tab4_8)<-c("tipo_de_cooperativa","estado")
attach(dados.tab4_8)Tabela 14.11
Tabela 14.11: Valores observados e esperados para o Exemplo 14.8.
CrossTable(estado,tipo_de_cooperativa,
prop.r=FALSE, prop.c=FALSE, prop.t=FALSE, prop.chisq=FALSE, expected=TRUE,
digits=0)##
##
## Cell Contents
## |-------------------------|
## | N |
## | Expected N |
## |-------------------------|
##
##
## Total Observations in Table: 1551
##
##
## | tipo_de_cooperativa
## estado | 1.Consumidor | 2.Produtor | 3.Escola | 4.Outras | Row Total |
## ----------------|--------------|--------------|--------------|--------------|--------------|
## 1.São Paulo | 214 | 237 | 78 | 119 | 648 |
## | 157 | 269 | 143 | 79 | |
## ----------------|--------------|--------------|--------------|--------------|--------------|
## 2.Paraná | 51 | 102 | 126 | 22 | 301 |
## | 73 | 125 | 67 | 37 | |
## ----------------|--------------|--------------|--------------|--------------|--------------|
## 3.Rio G. do Sul | 111 | 304 | 139 | 48 | 602 |
## | 146 | 250 | 133 | 73 | |
## ----------------|--------------|--------------|--------------|--------------|--------------|
## Column Total | 376 | 643 | 343 | 189 | 1551 |
## ----------------|--------------|--------------|--------------|--------------|--------------|
##
##
## Statistics for All Table Factors
##
##
## Pearson's Chi-squared test
## ------------------------------------------------------------
## Chi^2 = 173.38 d.f. = 6 p = 8.6339e-35
##
##
## Exemplo 14.10
anos_exp<-c(2,4,5,6,8)
n_clientes<-c(48,56,64,60,72)Queremos testar:
- \(H_0:\rho=0\)
- \(H_1: \rho\ne 0\)
No R, o teste para o coeficiente de correlação é realizado com a função cor.test:
cor.test( x, y,
alternative = c("two.sided", "less", "greater"),
method = c("pearson", "kendall", "spearman"),
conf.level = 0.95 )
x, y : Vetores aos quais se deseja avaliar a correlação
alternative : "two.sided": Testa H_1: rho != 0
"greater" : testa se a associação entre x e y é positiva
"less". : testa se a associação entre x e y é negativa
method:"pearson", "kendall", ou "spearman"
conf.level: Nivel de confiança desejado.Desta forma, para testar as hipóteses desejadas, teremos:
cor.test(anos_exp,n_clientes,alternative="two.sided", method="pearson",conf.level=0.95)##
## Pearson's product-moment correlation
##
## data: anos_exp and n_clientes
## t = 5.27, df = 3, p-value = 0.013
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.4185 0.9968
## sample estimates:
## cor
## 0.9514.6 Outro teste de Aderência
Exemplo 14.6 (Continuação)
Queremos testar:
- \(H_0: F(x)=F_0(x),\forall x\).
- \(H_0: F(x)\ne F_0(x)\), para algum \(x\).
com \(F_0(x)\) a f.d.a. da v.a. \(X\sim N(10,25)\).
O teste de Kolmogorov-Smirnov (KS) e calculado, no R, através da função ks.test:
ks.test(x, y, ...,
alternative = c("two.sided", "less", "greater"))
x: Dados da amostra a qual se deseja testar a função empirica
y: Vetor de dados (valores esperados) ou o nome (entre aspas) da distribuição esperada de x, e.g., "pnorm"
...: Parâmetros da distribuição esperada
alternative = c("two.sided", "less", "greater"): Como no teste anteriorAssim, para testar se os valores do Exemplo 14.6 são aderentes à distribuição normal, faremos:
ks.test(dados.ex14_6,y="pnorm",mean=10,sd=5)##
## One-sample Kolmogorov-Smirnov test
##
## data: dados.ex14_6
## D = 0.116, p-value = 0.77
## alternative hypothesis: two-sided# Gráfico qxq da distribuição emprírica e da distribuição Normal(10,25)
qqnorm(dados.ex14_6,cex=2,pch=20,col="darkblue", xlab="Quantis da Normal Padrão", ylab="Quantis dos dados", main="")
qqline(dados.ex14_6)
Figura 14.4: Quantis da normal padrão contra quantis dos dados.