##################################################### ################## AULA INFERENCIA ############ ################## CAPÍTULO 4 WOOLDRIGDE ############ ##################################################### # 29/05/2015 # Monitor: Wander Plassa # Para deletar tudo que está aberto rm(list=ls()) # 1) Definir diretório getwd() setwd("") # 2) chamar o comando "foreign" para obter bases de dados de outros programas, por exemplo Stata install.packages("foreign") library(foreign) # 3) Carregar uma base de dados do Wooldridge # Abrindo arquivo com extensão Dta 4 edição setwd("C:/Users/Wander Plassa/Google Drive/Doutorado/PAE/Econometria1_201701/Bases_wooldridge_4ed_STATA") CEOSAL2 = read.dta("CEOSAL2.dta") TWOYEAR = read.dta("twoyear.dta") ## Caso queira salvar o banco no formato do R use o seguinte comando: setwd("C:/Users/Wander Plassa/Google Drive/Doutorado/PAE/Econometria1_201701/Aula_Cap4") save(CEOSAL2,file="CEOSAL2.Rda") save(TWOYEAR,file="TWOYEAR.Rda") ## Para carregar esse banco de dados no novo formato em uma nova ocasião basta usar: load("CEOSAL2.Rda") load("TWOYEAR.Rda") # 4) Realizando regressões, com a base de dados # "CEOSAL2", em que o salário é a variável dependente. # A primeira é uma regressão simples, com apenas uma variável independente (vendas). reg1 = lm(log(CEOSAL2$salary) ~ log(CEOSAL2$sales), data=CEOSAL2) summary(reg1) ## A segunda é a regressão adicionando a variável "mktval" ## Que é o valor de mercado da firma e "profmarg" que é o lucro como porcentagem das vendas reg2 = lm(log(CEOSAL2$salary) ~ log(CEOSAL2$sales) + log(CEOSAL2$mktval) + CEOSAL2$profmarg, data=CEOSAL2) summary(reg2) ## A terceira é a regressão completa com "ceoten" que é anos como CEO na firma ## e "comten" é o total de anos com a firma reg3 = lm(log(CEOSAL2$salary) ~ log(CEOSAL2$sales) + log(CEOSAL2$mktval) + CEOSAL2$profmarg + CEOSAL2$ceoten + CEOSAL2$comten, data=CEOSAL2) summary(reg3) ###################### ##### TESTE T ######## ###################### # Teste para UM [UNICO] parâmetro populacional ######################## # 1) Calculo teste de T ######################## # No resultado da regressão o teste de t já aparece para cada parâmetro, # contudo caso queria calcular manualmente basta fazer a seguinte conta # tcal = â - 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"] # 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 # Valores de t calculado 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. # Mas quão distante de zero ele deve estar? ############################ # 2) Nivel de significância: ############################ # Para tal, precisamos definir níveis de significância, # Isto é, a probabilidade de reijeitar a H0 quando ela é verdadeira. # Portanto, quanto menor o nível de significância mais robustos # serão resultados obtidos, uma vez que a chance de falarmos # que a covariada tem efeito sobre a variável dependente # quando na verdade ela não tem, se torna pequena. # Valores comuns do nível de significância: 1%, 5% e 10%. # Com valor t calculado em mãos e 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) # 177 Graus de Liberdade e 5% de significância para um teste unilateral t_tab_bi_5 = qt(.975, df=171) # 177 Graus de Liberdade e 5% de significância para um teste bilateral # No nosso caso, a 5%, considerando uma comparação bilateral (usado quando # o sinal de "a" não é bem determinado pela teoria - mais prudente), teremos # um t tabelado com 177 - 5 - 1 ( graus de liberdade) = 1,97. # Isto é, t calculado 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. ######################################## # 3) 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. # (â - 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. ############################################################ # 4) Computando o valor de "p" ou p-valor para o teste de t ########################################################### # outro valor que o programa R já nos apresenta no sumário da regressão # é o p-value. Ele contorna o problema de arbitrariedade # que fazemos 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) # 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. ############################ # 5) Intervalo de Confiança ############################ # fornecem uma gama de valores prováveis para o parâmetro populacional, # e não apenas uma estimativa pontual. # O calculo do intervalo se dá pela seguinte conta: # â +/- (t tabelado*sd(â)) # Onde â é o parametro estimado, t tabelado e 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), # o intervalo seria: # 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) # Há uma chance de 95% de que o parâmetro verdadeiro de "sales" # esteja nesse intervalo (não há uma garantia). # Podemos facilitar a conta, usando o comando a seguir. # Para ver todos os intervalos, basta usar o seguinte comando: confint(reg3, level = 0.90) confint(reg3, level = 0.95) confint(reg3, level = 0.99) confint(reg3, level = 0.10) ######################### # 6) Multiplas restrições ######################### # Conisdere a seguinte regressão, COM A BASE DE DADOS "TWOYEAR" reg4 = lm(TWOYEAR$lwage ~ TWOYEAR$jc + TWOYEAR$univ + TWOYEAR$exper, data = TWOYEAR) summary(reg4) # # Queremos testar se a diferença dos coeficientes relativos # a universidade e cp é estatisticamente signficativa. # O teste de t para esse tipo de teste será # H0: B1_est = B2_est ou B1_est - B2_est = 0 # (B1_est - B2_est)/EP(B1_est - B2_est) # Isto é, O numerador é fácil obter: numer = 0.0666967 - 0.0768762 # e nesse caso será igual a -0,0102. # contudo o denominador não é dado diretamente. # Podemos obter da seguinte forma o denominador # EP(B1_est - B2_est) = {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((TWOYEAR$univ) - mean(TWOYEAR$univ)*((TWOYEAR$jc) - mean(TWOYEAR$jc))) EP_b1_b2 = ((0.0068288)^2 + (0.0023087)^2 - 2*S12)^1/2 tcal = numer/EP_b1_b2 ##### 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 controlada. # HA é que H0 é falsa # A fórmula do teste de F é dado por: # F = [(SQRm - SQR)/q]/[SQR/(n-k-1)] # Onde SQR é a soma dos Quadrados residuais obtido na ANOVA # SQRm é referente a soma dos quadrados residuais da regressão # sem o conjunto de covariadas que queremos testar # q é o número de 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 summary(reg3) # Para calcular a Soma dos Quadrados Residuais ao quadrado devemos proceder # da seguinte forma: resid = (reg3$residuals) resid_quadrado = resid*resid SQR = 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 menor seria (Sem as variáveis que estou tentando testar): reg3menor = lm(log(CEOSAL2$salary) ~ CEOSAL2$profmarg + CEOSAL2$ceoten + CEOSAL2$comten, data=CEOSAL2) summary(reg3menor) # Calculando o SQRmenor, teremos residmenor = (reg3menor$residuals) residmenor_quadrado = residmenor*residmenor SQRmenor = sum(residmenor_quadrado) # Conferindo novamente: anova(reg3menor) # Portanto, nosso SQRm = 63.427 # Portanto o F calculado seria: # F calculado = [(63.427 - 41.856)/2] / [41.856/(177 - 5 - 1)] Fcalc_sal_mkt = ((SQRmenor - SQR)/2) /( SQR/(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, obtermos o valor 3.0718. Podemos também obter o # F tabelado diretamente: Ftab_sal_mkt = qf(.95, df1=2, df2=173) # F tabelado com 2 restrições e 173 graus de liberdade # 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 ################### # BASE PNAD2015 # ################### rm(list=ls()) # Apagando todas as bases antigas setwd("C:/Users/Wander Plassa/Google Drive/Doutorado/PAE/Econometria1_201701/Aula_Cap4") load("pes2015.Rda") # Carregando a base de dados da PNAD 2015 # Renomeando as variáveis names(pes2015)=c("Sexo","Idade","Cor","Renda","horas_trab","Experiencia","Anos_Estudos") # Retirando da amostra rendas maiores que 50 mil reais ou que não informaram sub_pes2015 = subset(pes2015,pes2015$Renda <= 50000 & pes2015$Cor ==2 & pes2015$Renda != "NA") # Rodando uma regressão multipla, com a variável dependente em log regsal = lm(sub_pes2015$Renda ~ sub_pes2015$Experiencia + sub_pes2015$Anos_Estudos) summary(regsal) tcalexper_0 = (7.077 - 0)/ 2.707 tcalexper = (7.077 - 7)/ 2.707 t_tab_bi_5 = pt(tcalexper_0, df=65264) pvalor = (1 - t_tab_bi_5) * 2