##############################################################
#              Aula 23 - Teste de Independencia              #
##############################################################

rm(list = ls()); cat('\014')

##############################################################
###           ###  
### EXEMPLO 1 ###
###           ###
#################


# Passo 1: Especificar as hipoteses H0 e Ha
##############################################################

# Dados
dados = matrix(c(56,71,12,47,163,38,14,42,85),
               nrow = 3, byrow = T)

colnames(dados) = c("M Alta", "M Media", "M Baixa")
rownames(dados) = c("F Alta", "F Media", "F Baixa")

dados

## M Alta M Media M Baixa
## F Alta 56 71 12
## F Media 47 163 38
## F Baixa 14 42 85

# Tabela frequencia observada
obs = matrix(0, ncol = 4, nrow = 4)
obs[1:3,1:3] = dados
obs[, 4]     = rowSums(obs)
obs[4,]      = colSums(obs)
colnames(obs) = c("M Alta", "M Media", "M Baixa", "Total")
rownames(obs) = c("F Alta", "F Media", "F Baixa", "Total")

obs

##         M Alta  M Media  M Baixa  Total
## F Alta    56       71      12      139
## F Media   47      163      38      248
## F Baixa   14       42      85      141
## Total    117      276     135      528


# Passo 2: Especificar a estatistica do teste e sua
# distribuicao, sob H0
##############################################################

# Frequencia esperada
e = matrix(nrow = 3, ncol = 3)
colnames(e) = c("M Alta", "M Media", "M Baixa")
rownames(e) = c("F Alta", "F Media", "F Baixa")
for(i in 1:3)
{
  for(j in 1:3)
  {
    e[i, j] = round((obs[i,4] * obs[4,j])/obs[4,4], 1)
  }
}
e

##         M Alta  M Media  M Baixa
## F Alta   30.8     72.7    35.5
## F Media  55.0    129.6    63.4
## F Baixa  31.2     73.7    36.1

o  = dados
Q2 = round(sum( colSums((o - e)^2 / e) ), 1)

Q2

## [1] 145.5


# Passo 3: Fixar o nivel de significancia do teste (α)
##############################################################

# erro tipo 1 = alfa
alfa = 0.05

# numero de linhas e colunas
r = 3
s = 3

# hipoteses de uma cauda
# valor de qui-quadrado para o erro tipo 1
qc = round(qchisq(1 - alfa, (r-1)*(s-1)), 2)
qc

## [1] 9.49


# Passo 4: Calcular o p-valor (ou a regiao critica do teste)
##############################################################

valorp = 1 - pchisq(Q2, (r-1)*(s-1))
valorp

## [1] 0


# Passo 5: Decidir entre H0 e Ha, comparando o valor p com α
# (ou verificando se a estatistica do teste pertence ou nao 
# a regiao critica)
##############################################################

chisq.test(dados)

##
## Pearson’s Chi-squared test
##
## data: dados
## X-squared = 145.78, df = 4, p-value < 2.2e-16




##############################################################

rm(list = ls()); cat('\014')

##############################################################
###           ###  
### EXEMPLO 2 ###
###           ###
#################

x = c(176, 183, 185, 190, 191, 192, 201, 205, 214, 220)

n = length(x)
j = 1:n

# como n = 10, segue que
c = 3/8

w = round((j - c)/(n + 1 - 2*c), 2)
Prob = round(qnorm(w), 2)

Tab = cbind(j, x, w, Prob)
Tab

##        j  x   w    Prob
## [1,]   1 176 0.06 -1.55
## [2,]   2 183 0.16 -0.99
## [3,]   3 185 0.26 -0.64
## [4,]   4 190 0.35 -0.39
## [5,]   5 191 0.45 -0.13
## [6,]   6 192 0.55  0.13
## [7,]   7 201 0.65  0.39
## [8,]   8 205 0.74  0.64
## [9,]   9 214 0.84  0.99
## [10,] 10 220 0.94  1.55

par(mar = c(5,5,2,2))
plot(x, Prob, pch = 20, ylab = expression(phi^{-1}))
(z <- line(Prob~x))

abline(coef(z), col = "blue", lwd = 1.5)

# OU

# Caso nao tenha o pacote instalado
# install.packages("car", dep = T)

library(car)

qqPlot(x, distribution = "norm", pch = 20, cex = 1.5)

## [1] 10 1