#######################
#Aula R################
#Monitor: Maria Isabel#
#isabeltheodoro@usp.br#
#Data: 21/05/2018######
#######################

########
#AULA 3#
########

rm(list=ls()) # removendo todos os objetos ativos

#Definir diretório
getwd()
#ou
setwd("")

#Nesta aula vamos usar duas bases de dados do Wooldridge: ceosal2 e twoyear

#Primeiro vamos abrir a base ceosal2

load("ceosal2.RData")
attach(data)

############
#REGRESSÕES#
############

#reg1: salário (salary) sendo explicado pelas vendas(sales)

reg1 = lm(log(salary) ~ log(sales), data=data)
summary(reg1)

#Interpretação:
#Coeficiente log(sale) é a elasticidade estimada de salary  em relação a sales.
#Um aumento de 1% nas vendas das empresas aumentam o salário do CEOs em  cerca de 0,224% 

#reg2: salário sendo explicado por vendas mktval (valor de mercado da firma) e profmarg (lucro como porcentagem das vendas)

reg2 = lm(log(salary) ~ log(sales) + log(mktval) + profmarg, data=data)
summary(reg2)

#reg3: acrescenta ceoten (anos como CEO da firma) e comten (total de anos na firma)

reg3 = lm(log(salary) ~ log(sales) + log(mktval) + profmarg + ceoten + comten, data=data)
summary(reg3)

######################
##### TESTE T ########
######################

# Teste para UM [ÚNICO] parâmetro populacional

###############
#ESTATÍSTICA t#
###############

#No resultado da regressão a estatística t já aparece calculada para cada parâmetro.

# Cálculo manual:

#  t = â - a / sd(â)

# onde:
#"â" é o parametro estimado;
#"a" o valor a ser testado e
#sd(â) é o erro padrão estimado do parâmetro estimado. 

# Por exemplo, variáveis na reg3
# se o meu objetivo é testar a hipótese nula de que H0: a = 0, podemos fazer
# a seguinte conta:

# Podemos colocar os coeficientes direto ou considerar os coeficientes do sumário
tsales = (0.187787 - 0)/0.040003

tsales = (coef(summary(reg3))[2, "Estimate"] - 0)/coef(summary(reg3))[2, "Std. Error"]
tsales
# Para as outras variáveis

tmktval = (0.099872 - 0)/0.049214
tprofmarg = (-0.002211 - 0)/0.002105
tceoten = (0.017104 - 0)/0.005540
tcomten = (-0.009238 - 0)/0.003337

#Interpretação: Uma estatítica t suficientemente longe de zero 
#resultará na rejeição de H0.
# Isto é, rejeitará a hipótese que a variável independente em questão 
# não tem efeito sobre a variável depedente.

#Questão: Mas quão distante de zero ele deve estar?

########################
#NÍVEL DE SIGNIFICÂNCIA#
########################

# Nível de significância: probabilidade de reijeitar a H0 quando ela é verdadeira.
# Quanto menor o nível de significância, mais robustos é o resultado obtido.
# Ou seja, a chance de falarmos que a covariada tem efeito 
# sobre a variável dependente
# quando na verdade ela não tem, é pequena.

# Valores comuns do nível de significância: 1%, 5% e 10%.
# Com a estatística t calculada e definido o nível de significância,
# podemos encontrar o t tabelado (t com distribuição n-k-1).
# Onde n é o tamanho da amostra e k o número de parâmetros na regressão.
# e -1 refere-se ao intercepto

t_tab_uni_5 = qt(.95, df=171)   # 171 Graus de Liberdade e 5% de significância para um teste unilateral
t_tab_bi_5 = qt(.975, df=171)   # 171 Graus de Liberdade e 5% de significância para um teste bilateral

# No nosso caso, a 5%, considerando uma comparação bilateral 
# uma estatística t acima de |1,97|, nesse exemplo, nos indica que,
# com um nível de significância a 5%, as covariadas afetam o log(salário)

# Por exemplo, o t calculado da covariada "profmarg" é de -1.050. Em módulo,
# ela não ultrapassa o valor de 1,96. Portanto não podemos rejeitar a 
# hipótese nula a 5% para essa variável.
summary(reg3)
#####################################
#Testando outras Hipóteses sobre "a"#
#####################################

# Podemos testar outros valores específicos para "a" além de zero.
# Por exemplo, fazer uma H0 : "a" = 1. Para isso, basta substituir
# 1 no lugar de "a" e proceder da mesma forma, escolhendo os 
# níves de significância.

# t = (â - 1)/sd(â)

# Por exemplo, queremos testar se a variável "ceoten" é estatisticamente
# diferente de 1 a 5%.

tceoten_1 = (-0.009238 - 1)/0.003337

# O valor em módulo é 302.4387. Bem acima de 1.96 (valor de t tabelado a 5%).
# Podemos rejeitar a H0 que "comten" é igual a 1.

######################################
#Cálculos dos p-valores dos testes t #
######################################

# Outro valor que o programa R já nos apresenta no sumário da regressão é o p-value.
# O p-valor contorna o problema de arbitrariedade ao escolher o valor de significância.

# Ele nos responde a seguinte pergunta: Dado o valor observado da estatística t, 
# qual o menor nível de significância para o qual a hipótese nula seria rejeitado?

summary(reg3)

pvalue<- 2*(pt(-1.050,171))
pvalue

# No nosso exemplo, para a variável "profmarg", para a hipótese nula (H0: a = 0) 
# ser rejeitada, deveríamos aceitar um nível de significância maior do que 29%.

# Isso significa que os valores menores de p são evidência contra H0; 
# p-valores elevados oferecem pouca evidência contra H0.

# Ele é um teste bilateral que testa a hipótese do coeficiente ser diferente de 0.

##########################
# Intervalo de Confiança #
##########################

# Fornecem uma extensão dos valores prováveis para o parâmetro populacional, e não somente uma estimativa pontual.

# O IC é calculado da seguinte forma:

# IC de "a" a 95% = â +/- (t tabelado*sd(â))

# Onde:
# â é o parametro estimado
# t tabelado é o valor t (n-k-1) obtido na tabela de t para o nível de signifcância desejado e bilateral.
# Para a variável "log(sales)" do nosso exemplo, com 95% de confiança (ou 5% de sig), 


# Podemos calcular plotando os valores diretamente de duas formas novamente.
0.187787 + t_tab_bi_5*0.040003 #Intervalo Superior

0.187787 - t_tab_bi_5*0.040003 #Intervalo Inferior

coef(summary(reg3))[2, "Estimate"] + (t_tab_bi_5*coef(summary(reg3))[2, "Std. Error"] ) #Intervalo Superior

coef(summary(reg3))[2, "Estimate"] - (t_tab_bi_5*coef(summary(reg3))[2, "Std. Error"] ) #Intervalo Inferior


# O intervalor seria ->  0.1085811 e 0.2669929 (zero está fora do intervalo)

# Se as amostras aleatórias fossem obtidas repetidas vezes, com o IC calculado a cada vez,
#então o valor populacional desconhecido estaria dentro dos IC calculados em 95% das vezes.
#Como temos uma única amostra para calcular o IC, não sabemos se o valor verdadeiro de beta
#está realmente contido no intervalo. Esperamos ter obtido uma amostra que seja uma das 95%
#de todas as amostras em que a estimativa  do IC contém o beta verdadeiro, mas não temos 
#essa garantia.

# Podemos facilitar a conta, usando o comando a seguir.
# Para ver todos os intervalos, basta usar o seguinte comando:
confint(reg3, level = 0.95)

#Outros níveis de confiança
confint(reg3, level = 0.90)
confint(reg3, level = 0.99)
confint(reg3, level = 0.10)


########################
# Múltiplas restrições #
########################

rm(list=ls()) # removendo todos os objetos ativos

#Agora vamos usar a base 

load("twoyear.RData")
attach(data)

# Conisdere a seguinte regressão:

#reg4: salário sendo explicado jc (créditos freq. em um curso prof. de 2 anos), 
#uni (créditos freq. em um curso superior de 4 anos) e 
#exper (meses na força de trabalho)

reg4 = lm(lwage ~ jc + univ + exper, data = data)
summary(reg4)

# Queremos testar se a diferença dos coeficientes relativos 
# a univ e jc é estatisticamente signficativa. 

# H0: B1_est = B2_est ou  B1_est - B2_est = 0

# t = (B1_est - B2_est)/ep(B1_est - B2_est)

# Isto é, O numerador é fácil obter:

numerador = 0.0666967 - 0.0768762
numerador

# e nesse caso será igual  a -0,0102. 
# contudo o denominador não é dado diretamente.

# Podemos obter da seguinte forma o denominador
# denominador = {ep(B1_est)^2 + ep(B2_est)^2 - 2s12}^1/2
# Onde S12 é a covariância entre as duas variáves
# s12 = (1/n-1)*sum[(x1 - média(x1))*(x2 - média(x2))]

S12 = (1/(177-1))*sum((univ) - mean(univ)*((jc) - mean(jc)))

denominador = ((0.0068288)^2 + (0.0023087)^2 - 2*S12)^1/2                       

t = numerador/denominador
t

t_tab_uni_5 = qt(.95, df=6759)   # 171 Graus de Liberdade e 5% de significância para um teste unilateral
t_tab_bi_5 = qt(.975, df=6759)   # 171 Graus de Liberdade e 5% de significância para um teste bilateral

######################
##### TESTE F ########
######################

# testar hipóteses que envolvem mais do que um parâmetro populacional

# Testaremos várias restrições

# H0 é que um conjunto de variáveis não tem nenhum efeito sobre a var. dependente, 
# Por exemplo H0: a1 = 0, a2 = 0.
# uma vez que um outro conjunto de variáveis foi controlado.

# HA é que H0 é falsa

# A fórmula do teste de F é dado por:
# F = [(SQRr - SQRur)/q]/[SQRur/(n-k-1)]

# Onde:
# SQRr é a Soma dos Quadrados dos Resíduos do modelos restrito
# SQRur é a Soma dos Quadrados dos Resíduos do modelos irrestrito
# q é o número de restrições - variáveis retiradas do modelo

# Exemplo: suponha que queremos testar que "log(sales)" e "log(mktval)",
# tenham efeito conjunto sobre a nossa variável dependente "log(salary)".

# Para isso a H0: log(CEOSAL$sales) = 0, log(CEOSAL$mktval) = 0

# Lembrando dos valores da regressão
rm(list=ls()) # removendo todos os objetos ativos
load("ceosal2.RData")
attach(data)
reg3 = lm(log(salary) ~ log(sales) + log(mktval) + profmarg + ceoten + comten, data=data)

summary(reg3)

# Para calcular a Soma dos Quadrados Residuais ao quadrado devemos proceder
# da seguinte forma:

resid = (reg3$residuals)
resid_quadrado = resid*resid
SQRur = sum(resid_quadrado)


# So para conferir podemos solicitar a ANOVA. O ultimo valor é a SQR.
anova(reg3)

# Portanto, nosso SQR = 41.856

# A regressão restrita (sem as variáveis que estou tentando testar) seria :

reg3r = lm(log(salary) ~  profmarg + ceoten + comten, data=data)
summary(reg3r)

# Calculando o SQRr, teremos

resid_r = (reg3r$residuals)

resid_r_quadrado = resid_r*resid_r

SQRr = sum(resid_r_quadrado)

# Conferindo novamente: 

anova(reg3r)

# Portanto, nosso SQRr = 63.427

#Então:

# F  = [(63.427 - 41.856)/2] / [41.856/(177 - 5 - 1)]

Fcalc_sal_mkt = ((SQRr - SQRur)/2) /( SQRur/(177 - 5 - 1))

# Portanto, o F calculado = 44.06347 ao considerar as variáveis sales e mkt.

# Ao analisar a tabela do teste de F, considerando os
# graus de liberdade, acima de 120, com 5% de significância
# e um q=2

Ftab_sal_mkt = qf(.95, df1=2, df2=173) 

# Obtermos o valor 3.0718.
# Interpretação:Como F calculado > F tabelado, podemos rejeitar a H0
# de que as variáveis retiradas do modelo não possuem
# efeito sobre a variável dependente

###########################################################
# ESTATÍSTICA F  PARA SIGNIFICÂNCIA GERAL DE UMA REGRESÃO #
###########################################################

# F = (R2/k)/(1-R2)/(n-k-1)

# A ideia dessa estatística F é testar se nenhuma das 
# variáveis explicativas tem efeito sobre y.

F<- (0.3525/5)/((1-0.3525)/(177 - 5 - 1))

#Condferindo:
summary(reg3)