#============================================================
# Estatística Descritiva: análise bidimensional do
# comportamento/associação dos dados             (19/04/2011)
#
#                                 Professor Walmes M. Zeviani
#                                      www.leg.ufpr.br/ce003b
#============================================================

#------------------------------------------------------------

Var <- function(x) (length(x)-1)*var(x)/length(x)

#------------------------------------------------------------
# importação dos dados

da <- read.table("renda-alfab.txt",
                 header=TRUE, sep="\t", quote="")
str(da)

#------------------------------------------------------------
# criando categorias de renda e pegando os estados do sul

sul <- subset(da, estado%in%c("PR","SC","RS"))
sul <- sul[complete.cases(sul),]
sul$c.renda <- cut(sul$renda, c(0,1000,5000,10000,50000,Inf))
sul$estado <- as.factor(as.character(sul$estado))
str(sul)

#============================================================
# qualitativa vs qualitativa (classe de renda vs estado)

#------------------------------------------------------------
# tabela de frequências absolutas

tfa <- with(sul, table(estado, c.renda))
tfa

sum(tfa)

#------------------------------------------------------------
# tabela com frequências relativas ao total

tfrt <- with(sul, prop.table(tfa))
tfrt

sum(tfrt)

#------------------------------------------------------------
# tabela com frequências relativas ao total das linhas

tfrl <- with(sul, prop.table(tfa, margin=1))
tfrl

apply(tfrl, 1, sum)

#------------------------------------------------------------
# tabela com frequências relativas ao total das colunas

tfrc <- with(sul, prop.table(tfa, margin=2))
tfrc

apply(tfrc, 2, sum)

#------------------------------------------------------------
# sob a hipótese de independência entre os níveis das
# categorias qual seria a frequência esperada em cada cédula?

with(sul, table(c.renda))
with(sul, table(estado))

fel <- apply(tfrt, 1, sum) # frequencia esperada nas linhas
fel

fec <- apply(tfrt, 2, sum) # frequencia esperada nas colunas
fec

FE <- outer(fel, fec, "*") # frequencias esperadas sob ind
FO <- tfrt                 # frequencias observadas

#------------------------------------------------------------
# como medir a discrepância entre observado e esperado?

FE
FO

FO-FE

#------------------------------------------------------------
# usa-se a medida qui-quadrado de pearson

qui2 <- nrow(sul)*sum((FO-FE)^2/FE)
qui2

#------------------------------------------------------------
# adota-se a medida T, que varia de 0 a 1

T <- sqrt(
          (qui2/nrow(sul))/
          ((nlevels(sul$estado)-1)*
           (nlevels(sul$c.renda)-1))
          )
T

# quanto maior T, mas associadas as variáveis

#------------------------------------------------------------
# gráficos que representam a relação

require(lattice)

barchart(tfrl, auto.key=list(columns=5), horizontal=FALSE)
barchart(t(tfrc), auto.key=list(columns=3), horizontal=FALSE)

#============================================================
# quantitativas vs quantitativas (renda vs prop de alfabet)

#------------------------------------------------------------
# diagrama de dispersão

xyplot(alfab~log(renda), sul, pch=19)

#------------------------------------------------------------
# diagrama de dispersão separando os estados

xyplot(alfab~log(renda), groups=estado, sul, pch=19,
       auto.key=list(columns=3))

#------------------------------------------------------------
# diagrama de dispersão separando os estados com tendência

xyplot(alfab~log(renda), groups=estado, sul, pch=19,
       auto.key=list(columns=3), type=c("p","smooth"))

#------------------------------------------------------------
# tipos de padrões (independente, associação + e -)

xyplot(rnorm(1000)~c(rnorm(1000)),
       type=c("p","smooth"), pch=19)

require(MASS)

X <- mvrnorm(1000, c(0,0), matrix(c(1,0.8,0.8,1),2,2))
xyplot(X[,1]~X[,2], type=c("p","smooth"), pch=19)

Z <- mvrnorm(1000, c(0,0), matrix(c(1,-0.8,-0.8,1),2,2))
xyplot(Z[,1]~Z[,2], type=c("p","smooth"), pch=19)

#------------------------------------------------------------
# como medir a associação entre as variáveis?
# covariância/correlação

#------------------------------------------------------------
# calculando a correlação

l.renda <- log(sul$renda)
alfab <- sul$alfab

m.renda <- mean(l.renda) # média do log da renda
m.renda

m.alfab <- mean(alfab)   # média do núm de alfabetizados
m.alfab

dp.renda <- sqrt(Var(l.renda)) # desvio padrão log renda
dp.renda

dp.alfab <- sqrt(Var(alfab))   # desvio padrão num alfab

dp.alfab

sum(
    (l.renda-m.renda)*(alfab-m.alfab)/ # produto dos desvios
    (dp.renda*dp.alfab)                # pondera pelos dp
    )/nrow(sul)                        # divide pelo n

m0 <- lm(alfab~l.renda)
summary(m0)
sqrt(summary(m0)$r.squared)
plot(alfab~l.renda)
abline(m0, col="red")

#============================================================
# quantitativa vs qualitativa (alfabetizados vs estado)

boxplot(alfab~estado, sul)

#------------------------------------------------------------
# será que a distribuição dos alfab depende do estado?

xyplot(alfab~estado, sul, type=c("p","a"), jitter.x=TRUE)

#------------------------------------------------------------
# qual a variância de alfab?

v0 <- Var(alfab)
v0

#------------------------------------------------------------
# qual a variância de alfab em cada estado?

v3 <- with(sul, tapply(alfab, estado, Var))
v3

#------------------------------------------------------------
# qual a média (ponderada) das variâncias?

n <- with(sul, tapply(alfab, estado, length))
n

vp <- sum(n*v3)/sum(n)
vp

#------------------------------------------------------------
# usa-se a medida R², ganho relativo na variância, entre 0, 1

1-vp/v0

m0 <- aov(alfab~estado, sul)
summary.lm(m0)
anova(m0)
TukeyHSD(m0)

#------------------------------------------------------------