# --- Aula Prática SME0122 Introdução à Inferência Estatística ---  #
# --- 25/10/11 - Prof. Cibele Russo --- #
# --- Alguns sites úteis --- #
# --- R Project: http://www.r-project.org/  --- #
# --- Busca no R: http://www.rseek.org/ --- #
# --- Using R for Introductory Statistics: http://cran.r-project.org/doc/contrib/Verzani-SimpleR.pdf --- #

# --- Analisando dados de Acessibilidade na Web (André P. Freire) --- #

dados<-read.csv("http://www.icmc.usp.br/~cibele/Dados/acessibilidade.csv", dec=",", sep=";")

attach(dados)

names(dados)

n <- nrow(dados)
n

# --- Análise descritiva --- #

# --- Tabelas de Frequências --- #
table(genero)
table(genero)/n
round(table(genero)/n, 3)

table(idade)
round(table(idade)/n, 3)

plot(idade, main="Gráfico de Barras de idade")


rbind(table(estado), round(table(estado)/n, 3))
rbind(table(area_atuacao), round(table(area_atuacao)/n, 3))
rbind(table(formacao), round(table(formacao)/n, 3))

plot(table(formacao))

rbind(table(q1), round(table(q1)/n, 3))
rbind(table(q2), round(table(q2)/n, 3))
rbind(table(q3), round(table(q3)/n, 3))
rbind(table(q4), round(table(q4)/n, 3))
rbind(table(q5), round(table(q5)/n, 3))
rbind(table(q6), round(table(q6)/n, 3))
rbind(table(q7), round(table(q7)/n, 3))

plot(idade)
plot(genero)
plot(estado)


# --- Análise descritiva: Análise gráfica --- #

idade.valores<-summary(idade)
estado.valores<-summary(estado) 
genero.valores<-summary(genero)  

 pie(idade.valores, col = gray(seq(0.3,1,length=nlevels(idade))), main="Gráfico de setores de idade")
 pie(estado.valores, col = gray(seq(0.3,1,length=nlevels(idade))), main="Gráfico de setores de estado")
 pie(genero.valores, col = gray(seq(0.3,1,length=nlevels(idade))), main="Gráfico de setores de gênero")

 par(mfrow=c(1,3))
 pie(idade.valores, col = gray(seq(0.3,1,length=nlevels(idade))), main="Gráfico de setores de idade")
 pie(estado.valores, col = gray(seq(0.3,1,length=nlevels(idade))), main="Gráfico de setores de estado")
 pie(genero.valores, col = gray(seq(0.3,1,length=nlevels(idade))), main="Gráfico de setores de gênero")

 par(mfrow=c(1,3))
 pie(idade.valores, col = rainbow(nlevels(idade)), main="Gráfico de setores de idade")
 pie(estado.valores, col = rainbow(nlevels(estado)), main="Gráfico de setores de estado")
 pie(genero.valores, col = rainbow(nlevels(genero)), main="Gráfico de setores de gênero")



# --- Definindo regiões --- #
# --- Fonte: IBGE --- #

regiao<-matrix(0, length(estado))

 CO<-c("GO", "MT", "MS", "DF")
 NE<-c("MA", "PI", "CE", "RN", "PB", "PE", "AL", "SE", "BA")
 N<-c("AC", "AM", "RR", "RO", "PA", "AP", "TO")
 SE<-c("SP", "RJ", "ES", "MG")
 S<-c("RS", "SC", "PR")

 regiao[which(pmatch(estado, CO, dup=T) != "NA")]<-"CO"
 regiao[which(pmatch(estado, NE, dup=T) != "NA")]<-"NE"
 regiao[which(pmatch(estado, N, dup=T) != "NA")]<-"N"
 regiao[which(pmatch(estado, SE, dup=T) != "NA")]<-"SE"
 regiao[which(pmatch(estado, S, dup=T) != "NA")]<-"S"

 regiao<-factor(regiao)

 plot(regiao)
 
 regiao.valores<-summary(regiao)
 pie(regiao.valores, col = gray(seq(0.3,1,length=nlevels(idade))), main="Gráfico de setores de região")
                    
 regiao.valores<-summary(regiao)
 pie(regiao.valores, col = rainbow(length(regiao.valores)), main="Gráfico de setores de região")


 boxplot(wab)


# --- Análises gráficas de associação entre variáveis

boxplot(wab ~ q1, main="Boxplot de WAB", xlab="Q1")
boxplot(wab ~ q2, main="Boxplot de WAB", xlab="Q2")
boxplot(wab ~ q3, main="Boxplot de WAB", xlab="Q3")
boxplot(wab ~ q4, main="Boxplot de WAB", xlab="Q4")
boxplot(wab ~ q6, main="Boxplot de WAB", xlab="Q6")
boxplot(wab ~ q7, main="Boxplot de WAB", xlab="Q7")
boxplot(wab ~ idade, main="Boxplot de WAB", xlab="Idade")
boxplot(wab ~ regiao, main="Boxplot de WAB", xlab="Região")
boxplot(wab ~ formacao, main="Boxplot de WAB", xlab="Formação")
boxplot(wab ~ area_atuacao, main="Boxplot de WAB", xlab="Área de atuação")


 par(mfrow=c(1,2))
 hist(wab, main="Histograma de WAB")
 boxplot(wab, main="Boxplot de WAB", col=2)
 shapiro.test(wab) # Teste para normalidade rejeita normalidade de wab



# --- Inferência a respeito da média de WAB = mu --- #

 Xbarra <- mean(wab)
 S <- sd(wab)



# --- Estimação pontual de mu --- #

 paste("Uma estimativa não viesada da média de WAB é", round(Xbarra,3))

 paste("Uma estimativa não viesada da variância de WAB é", round(S,3))



# Intervalo de confiança para mu = média(wab): supondo sigma=1 conhecido e usando o TLC

 sigma <- 1
 z <- abs(qnorm(0.025))

 LI = round( Xbarra - z * sqrt(sigma^2/n), 3)
 LS = round( Xbarra + z * sqrt(sigma^2/n), 3)

 paste("IC(0.95,mu) = (", LI, ",", LS, ")")


# Intervalo de confiança para mu = média(wab): supondo sigma desconhecido, usando T

 t <- abs(qt(0.025, n-1))

 LI = round( Xbarra - t * sqrt(S^2/n), 3)
 LS = round( Xbarra + t * sqrt(S^2/n), 3)

 paste("IC(0.95,mu) = (", LI, ",", LS, ")") 


# Teste H0: mu=2 contra H1: mu != 2 (supondo que sigma = 1 fosse conhecido) e usando o TLC

 sigma = 1
 Z <- (Xbarra-2)/(sigma/sqrt(n))
 2*( 1- pnorm(Z) )


# Teste H0: mu=2 contra H1: mu != 2 (variância desconhecida)
 
 T <- (Xbarra-2)/(S/sqrt(n))
 2*(1 - pt(T, n-1))
 


# Modelo de regressão linear para wab 

 fit.model<-lm(wab ~ genero + idade + regiao + formacao + q1 + q2 + q3 + q4 + q6 + q7)
 summary(fit.model)


# Seleção de variáveis stepAIC
 library(MASS)  
 fit.model.AIC <- stepAIC(lm(wab ~ genero + idade + regiao + formacao + q1 + q2 + q3 + q4 + q6 + q7))
 summary(fit.model.AIC)



# Outras funções úteis

qnorm(0.05)
rnorm(100) # Gera uma amostra de N(0,1)