######################################################
#####################################################
# AulaR #############################################
# Monitor: Wander Plassa ############################
# wanderplassa@usp.br ###############################
# Introdução ########################################

# Organização da Aula
 # 1 - Comandos básicos
 # 2 - Extração de dados da PNAD
 # 3 - Ler Dicionário da PNAD
 # 4 - Regressão Simples
 # 5 - SQE, SQR, SQT e R²





##################################################
############### Visão Geral ######################
###############     do R    ######################
##################################################



##################################################
############### Tutorial  ########################
############### Proprio R ########################
##################################################


# instalando tutorial do R ("install.packages()" serve para instalar funções 
# que não estão de inicio no R)
# A função para tutorial chama-se "swirl":

install.packages("swirl") 

# Toda vez que for usar a função se inicia chamando ela pelo comando "library"
library("swirl") 

#A partir desse ponto indica a função
swirl() 



###################################################
############### Introdução #########################
###################################################

#1) Encontrando  e escolhendo o diretório ~> get working directory (getwd)

# Se você deseja salvar ou carregar um arquivo, 
# você precisa saber o que o diretório de trabalho atual.

getwd() 

# Definir o diretório de trabalho poderia ser 
# feito através da  linha de comando: "setwd()"

setwd("C:/Users/Wander Plassa/Desktop/Pos_Graduacao/Doutorado/PAE/Econometria1_201701/Aula_01") 

# Ou manualmente em Files/More/Set as Working Directory/Escolhe a pasta
# Obs: as barras são invertidas


#2) Para começar a trabalhar verifica-se se existe objetvos abertos

# Se você deseja excluir todos os objetos atualmente disponível, você pode usar ls ():

rm(list = ls()) ## Dentro do ls você também pode especificar quais objetos serão deletados


#3) Dúvidas de funções:

# use a função "help()", entre o () o nome da função que deseja conhecer
# ou simplesmente use "?nomedafunção"
# Exemplos

help("swirl")

?lm


##################################################
############### ESTATISTICAS   ###################
##################################################

#1) Primero iremos criar um vetor, a título de exemplo:
# Podemos escrever assim:

v3 = c(1,2,3,4,5,6)

# com igual, ou

v3 <- c(1,2,3,4,5,6)

# com <-
# Obs: 
# -Ele diferencia maiscula e minuscula
# -Cuidado com nomes de bases e vetores com caracteres especiais
# O R não aceita caracteres especiais nos nomes e nem que ele c
# comece com números


# Por exemplo
3v = c(1,2,3,4,5,6)
_v = c(1,2,3,4,5,6)
v% = c(1,2,3,4,5,6)
vô = = c(1,2,3,4,5,6)
calça = = c(1,2,3,4,5,6)

# Não são aceitos. Mas

vv = c(1,2,3,4,5,6)
VV = c(1,2,3,4,5,6)
V3 = c(1,2,3,4,5,6)
V_ = c(1,2,3,4,5,6)
v_3 = c(1,2,3,4,5,6)

# São aceitos

#2) Estatísticas descritivas

v3
summary(v3) #Descrição total
mean(v3) #Media
var(v3) #Variância
quantile(v3) #quantils
sum(v3) #soma dos elementos
median(v3) #mediana
sd(v3) #erro padrão
max(v3) #Valor máximo
min(v3) #valo mínimo

# Se quiser guardar algum desses valores basta dar um nome
Media=mean(v3)

Media

##################################################
############### DADOS   ##########################
##################################################

#1) Intro

# R vem com conjuntos de dados integrados que podem ser recuperados pelo nome, 
# usando a função de "data()". Nesta classe, vamos utilizar "mtcars".


?data # conhecendo a função
data() #listanto todas as bases disponiveis
data("mtcars") 

# visualizando os dados 

View(mtcars) 

#Conhecendo a base usando a função help
? mtcars 


# Como obtemos estatísticas de um vetor, 
# podemos aplicar a função "summary()" para dados:

summary(mtcars)

# Summary de uma variável em específico 

summary(mtcars$mpg)


#1) Visualizando variáveis

mpg


# Para indicar uma variável de uma base de dados
# Ou uso a função

attach(mtcars) #Fixa nas variáveis dessa base

mpg

detach(mtcars) #Desvincula a base de dados

# Ou uso "$", como no exemplo a seguir: 

mtcars$mpg #ao inserir $ você indica que quer analisar uma variável em específico

mtcars$mpg

#2) Criando sub-bases a partir dos dados

efficients.cars = mtcars[mtcars$mpg > 20,] 

# obs: lembrando que o nome não pode começar com números ou ter caracteres
# especiais

# Você pode combinar várias condições usando & (e) ou | (Ou), 
# Por exemplo, duas condições: que o mpg seja superior a 20, 
# e câmbio automático:

efficients.autos = mtcars[mtcars$mpg> 20 & mtcars$am ==0,] 

# caracteres: & (e); | (ou); > (maior); < (menor); <= (menor ou igual);
# >= (maior ou igual); == (par exato); != (diferente)
efficients.autos


##################################################
########## Regressão #############################
##################################################


#1) Regresão Linear

# Usado para transformar essa relação em uma previsão: 
# por exemplo

# qual seria a eficiência de um carro 4500 libras???

# Podemos fazer isso através de um modelo linear, 
# ou regressão linear, o que é feito em R com a função lm. 
# Vamos salvar modelo linear em uma variável que chamamos de fit.

fit = lm(mpg ~ wt, mtcars) #o til (~) significa "explicado por" e lm signfica Linear Model

summary(fit)

# milhas por galão explicado pelo peso usando os dados mtcars. 

# Em primeiro lugar, temos a nossa estimativa do 
# intercepto: é milhas por galão hipotéticas de um carro que pesava 0 em nosso modelo linear.
# Em segundo lugar: o efeito da variável peso em milhas por galão. 
#  também chamado o coeficiente de inclinação do peso. 
# Isso mostra que há uma relação negativa, onde o aumento do peso diminui as milhas 
# por galão. Em particular, mostra que o aumento do peso por 1000 libras diminui 
# a eficiência de 5,3 milhas por galão.

# Este modelo prevê a eficiência para cada um dos carros 
# existentes, utilizando a função predict.

predict(fit) 

# Alguma ideia de como podemos responder a Pergunta Inicial?????


37.2851 + (-5.3445) * 4.5

# Se prevê uma eficiência de combustível de 13,2 milhas por galão. 
# Isto é o que um modelo linear na verdade significa:  uma 
# combinação linear do intercepto e inclinação.

#2) Regresão Multipla

# Por exemplo, digamos que você quer medir o efeito de não apenas o peso, 
# mas também o número de cilindros e volume do cilindro dos carros? 

# Que combinação de preditores vai prever melhor eficiência de combustível? 
# Quais preditores aumentam nossa precisão  estatisticamente significativa?

mfit = lm(mpg ~ wt + disp + cyl, data=mtcars)
summary(mfit)


# Observe que a tabela de coeficiente tem agora quatro linhas: 
# uma para o intercepto e três para os preditores. 
# Cada um destes contém ainda a estimativa do coeficiente, ou inclinação. 
# Ele também contém um P-valor para cada um dos preditores 
# de forma independente. Observe que os valores de p para o peso e o número de cilindros 
# são ambas significativas. Podemos ver a classificação de importância 
# com base no número de estrelas, onde ** significa que é entre 0,001 e 0,01. 
# Mas a p-valor para o volume não é significativa, isto indica que essa 
# variável é redundante com um ou ambos os outros preditores, 
# e não fornece qualquer informação adicional significativa.

predict(mfit)

##############################################
############### Base da PNAD #################
##############################################

# Baixar dados PNAD 

# endereço: "http://www.ibge.gov.br/home/estatistica/populacao/trabalhoerendimento/pnad2013/microdados.shtm"
# baixar : Dados e Dicionários
# Colocar a base (arquivo txt) no diretorio de trabalho


# BAIXAR R 
# Arquivo Instalação do R


#1) Lendo dados, abrindo e salvando base de dados


# 1-4, 5-6, 16-17,18, 27-29, 33, 41, 325-336,667-668 
# V0101 (Ano), UF (Estado), V0301 (Número de ordem), 
# V0302 (Sexo),V8005 (Idade),V0404 (Cor), V4011 (Estado civil), V9532 (Renda), V4803 (educação)

# Para ler arquivos TXT uso read.fwf

pes2015= read.fwf(file = 'PES2015.txt', widths = c(4,2,-9,2,1,-8,3,-3,1,-7,1,-283,12,-330,2))

summary(pes2015)

# Como calcular no caso da PNAD
# Considere n1-n2, n3-n4, ...
# n3-n2-1 para achar o valor negativo


# Acrescentado nomes as variáveis
names(pes2015)=c("Ano","UF","N_ordem","Sexo","Idade","Cor","Estado_Civil","Renda","Anos_Estudo")

summary(pes2015)

# Analisando as primeiras 20 linhas
head(pes2015,20)

# Retirando valores não informado de renda
pes2015_sub=subset(pes2015,pes2015$Renda <= 100000)

summary(pes2015_sub)

# Salvando a base de dados em formato R
save(pes2015_sub, file="pes2015_sub.Rda") 


rm(list = ls())

# Carregando a base de dados salva
load("pes2015_sub.Rda")


#2) Regressão

attach(pes2015_sub)

# Regressão simples de escolaridade explicando a renda

reg_renda = lm(Renda ~ Anos_Estudo, pes2015_sub) #Recordamos que o til (~) significa "explicado por" e lm signfica Linear Model

summary(reg_renda)

# Podemos fazer essa regressão simples com o log da renda também

logrenda= log(Renda)

reg_logrenda = lm(log(Renda) ~ Anos_Estudo, pes2015_sub)

summary(reg_logrenda)

# 3) Calculo SQE, SQT, SQR e R² (PAGINA 77)

anova(reg_logrenda)

# SQE

SQE = sum((predict(reg_logrenda) - mean(logrenda))^2)

# SQR

SQR = sum(reg_logrenda$residuals^2)




# SQT

SQT = SQE + SQR


# R²

R_2 = SQE/SQT

sigmahat=sqrt(SQR/152289)


tssx=sum((Anos_Estudo-mean(Anos_Estudo))^2)

epb1=sqrt(sigmahat^2/tssx)