Category

Dicas de R

Testes de Hipóteses sobre um único parâmetro: o teste t

By | Dicas de R

Um importante tipo de hipótese que estamos interessados é o da forma em que:

(1)   \begin{align*} H_0 : \beta_j = a_j, \end{align*}

onde a_j é um número dado [em geral, a_j = 0]. Para a maioria dos testes bicaudais, a hipótese alternativa implica em:

(2)   \begin{align*} H_1 : \beta_j \neq a_j, \end{align*}

e para testes unicaudais, ou temos:

(3)   \begin{align*} H_1 : \beta_j < a_j \quad ou \quad \beta_j > a_j. \end{align*}

Essas hipóteses podem ser testadas usando um teste t que é baseado na seguinte estatística:

(4)   \begin{align*} t = \frac{\hat{\beta_j} - a_j}{ep(\hat{\beta_j})}. \end{align*}

Se H_0 é verdadeira, essa estatística possui uma distribuição t com n-k-1 graus de liberdade. Para ilustrar, estimamos uma função para o log do salário-hora. Assim, temos os parâmetros dos retornos percentuais de cada entrada no modelo. Podemos avaliar se, por exemplo, depois de controlar por educação e titularidade, experiência ainda tem um efeito estatisticamente significante no salário-hora.


library(wooldridge)

data("wage1") # puxamos os dados

summary(lm(log(wage) ~ educ + exper + tenure, data=wage1))

E abaixo, a saída da regressão.

E de fato, a 5\% de significância existe um efeito para experiência. Mais especificamente, um ano a mais de experiência na média se traduz em 0,41\% de aumento salarial. Observe ainda que a estatística t pode ser calculada como sendo o parâmetro \beta_j estimado sobre o erro-padrão. Para o caso da experiência, temos 0.004121/0.001723, que é igual a 2,39. Em outras palavras, podemos rejeitar a hipótese nula que \beta_j = 0.

__________________________

(*) Isso e muito mais você aprende em nosso Curso de Introdução à Econometria usando o R.

Análise da violência no Rio com o R

By | Dicas de R

O Instituto de Segurança Pública (ISP) disponibiliza uma série de dados relacionados à violência no Rio de Janeiro. Nesse post, mostro como coletar e visualizar alguns desses dados com o R. A seguir, baixamos os dados mensais referentes à várias métricas de violência no Estado.


#####################################################
##### Segurança Pública no Rio de Janeiro ###########
#####################################################

library(readr)
library(lubridate)
library(magrittr)
library(dplyr)
library(ggplot2)
library(scales)
library(BMR)

url = 'http://www.ispdados.rj.gov.br/Arquivos/DOMensalEstadoDesde1991.csv'
download.file(url, destfile='basededados.csv', mode='wb')
data = read_csv2('basededados.csv') %>%
mutate(date = make_datetime(vano, mes))

Uma vez que baixamos e lemos os dados do ISP, nós também criamos um vetor de datas a partir do próprio dataset. A seguir, visualizamos alguns dados referentes a homicídios e outros crimes relacionados.

A seguir, visualizamos crimes associados a roubos.

Por fim, visualizamos os homicídios associados a intervenções policiais.

Como mostra o bloxplot abaixo, os dados de homicídios por intervenção policial na ponta estão bem acima da mediana histórica. É realmente um número preocupante.

____________________________

Aprenda a fazer isso e muito mais com o R com o nosso Curso de Introdução ao R para Análise de Dados.

[Dicas de R] Regressão Múltipla

By | Dicas de R

Em post anterior das Dicas de R, vimos o modelo de regressão simples, onde y pode ser explicado por uma única variável x. O problema básico desse tipo de análise é que ela faz uma suposição bastante forte, qual seja, que x não está correlacionado com o erro, dificultando a aplicação da condição ceteris paribus. A análise de regressão múltipla, por outro lado, é mais receptiva a esse tipo de condição, uma vez que ela permite que controlemos outros fatores que afetam y, adicionando os mesmos na equação. Assim, por suposto, se queremos explicar y, podemos utilizar k variáveis, como abaixo:

(1)   \begin{align*} y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + ... + \beta_k x_k + u,  \end{align*}

onde \beta_0 é o intercepto, \beta_k é o parâmetro associado a x_k. De modo a obter uma estimativa para 1, devemos observar que

(2)   \begin{align*} E(u|x_1, x_2, ..., x_k) = 0. \end{align*}

Isto é, que todos os fatores no termo de erro não observado u sejam não correlacionados com as variáveis explicativas. De modo a obter estimativas para os \beta_k parâmetros, é possível recorrer ao método de mínimos quadrados ordinários. Isto é, dado

(3)   \begin{align*} \hat{y} = \hat{\beta_0} + \hat{\beta_1} x_1 + \hat{\beta_2} x_2 + ... + \hat{\beta_k} x_k, \end{align*}

onde \hat{\beta_k} é a estimativa de \beta_k, o método de MQO escolhe as estimativas \hat{\beta_k} que minimizam a soma dos quadrados dos resíduos:

(4)   \begin{align*} \sum_{i=1}^{n} (y_i - \hat{\beta_0} - \hat{\beta_1} x_{i1} - ... - \hat{\beta_k} x_{ik})^2. \end{align*}

O problema acima pode ser resolvido por meio de cálculo multivariado, de onde obtemos as condições de primeira ordem

(5)   \begin{align*} \sum_{i=1}^{n} (y_i - \hat{\beta_0} - \hat{\beta_1} x_{i1} - ... - \hat{\beta_k} x_{ik}) = 0 \nonumber \\ \sum_{i=1}^{n} x_{i1} (y_i - \hat{\beta_0} - \hat{\beta_1} x_{i1} - ... - \hat{\beta_k} x_{ik}) = 0 \nonumber \\ \sum_{i=1}^{n} x_{i2}(y_i - \hat{\beta_0} - \hat{\beta_1} x_{i1} - ... - \hat{\beta_k} x_{ik}) = 0 \nonumber \\ \sum_{i=1}^{n} x_{ik}(y_i - \hat{\beta_0} - \hat{\beta_1} x_{i1} - ... - \hat{\beta_k} x_{ik}) = 0, \nonumber \end{align*}

ou simplesmente, E(u) = 0 e E(x_j u) = 0.

# Interpretação da equação de regressão de MQO

Suponha que tenhamos

(6)   \begin{align*} \hat{y} = \hat{\beta_0} + \hat{\beta_1} x_1 + \hat{\beta_2} x_2. \end{align*}

O intercepto \beta_0 será então o valor previsto de y quando x_1 = x_2 = 0. Já as estimativas \hat{\beta_1} e \hat{\beta_2} devem ser interpretadas como efeito parcial ou simplesmente ceteris paribus. Isto é,

(7)   \begin{align*} \Delta \hat{y} = \hat{\beta_1} \Delta x_1 + \hat{\beta_2} \Delta x_2, \nonumber \end{align*}

de modo que obtemos a variação prevista em y dadas as variações em x_1 e x_2. Em particular, quando x_2 é mantido fixo, de modo que \Delta x_2 = 0, teremos

(8)   \begin{align*} \Delta \hat{y} = \hat{\beta_1} \Delta x_1. \nonumber \end{align*}

Ou, simplesmente,

(9)   \begin{align*} \frac{\partial \hat{y}}{\partial \hat{x_1}} = \hat{\beta_1}, \nonumber \end{align*}

onde \hat{\beta_1} irá medir o efeito da variação de x_1 em y, mantido x_2 constante.

# Exemplo: equação do salário-hora

De modo a ilustrar, vamos considerar o exemplo 3.2 de Wooldridge (2003), em que o mesmo utiliza o conjunto de dados wage1, disponível no pacote wooldridge. Ele pode ser acessado como abaixo.


library(wooldridge)
data(wage1)

modelo = lm(log(wage) ~ educ+exper+tenure, data=wage1)

E abaixo, o nosso modelo.

Dependent variable:
log(wage)
educ 0.092***
(0.007)
exper 0.004**
(0.002)
tenure 0.022***
(0.003)
Constant 0.284***
(0.104)
Observations 526
R2 0.316
Adjusted R2 0.312
Residual Std. Error 0.441 (df = 522)
F Statistic 80.391*** (df = 3; 522)
Note: *p<0.1; **p<0.05; ***p<0.01

De modo a obter a seguinte reta de regressão para o log do salário-hora

(10)   \begin{align*} \hat{log(wage)} = 0.284 + 0.092 educ + 0.0041 exper + 0.022 tenure. \end{align*}

De onde se conclui, por exemplo, que o aumento de um ano na educação formal equivale a um aumento de 9.2% no salário-hora, mantidos exper e tenure fixos.

Quer aprender mais sobre econometria? Conheça nosso Curso de Introdução à Econometria usando o R.

_______________________

Wooldridge, J. M. 2013. Introductory Econometrics: A Modern Approach. Editora Cengage.

Implementando regressões simples no R

By | Dicas de R

Regredir um variável x contra uma variável y é um poderoso recurso estatístico. De modo a explicar o método, suponha que estamos interessados em estimar os parâmetros populacionais \beta_0 e \beta_1 de um modelo de regressão simples

(1)   \begin{align*} y = \beta_0 + \beta_1 x + u  \end{align*}

a partir de uma amostra aleatória de y e x. Os estimadores de Mínimos Quadrados Ordinários (MQO) serão

(2)   \begin{align*} \hat{\beta}_0 &= \hat{y} - \hat{\beta_1} \bar{x} \\ \hat{\beta_1} &= \frac{Cov(x,y)}{Var{x}}. \end{align*}

Baseado nos parâmetros estimados, a reta de regressão será

(3)   \begin{align*} \hat{y} = \hat{\beta}_0 + \hat{\beta}_1 x. \end{align*}

Para uma dada amostra, nós precisaremos calcular as quatro estatísticas \bar{y}, \bar{x}, Cov(x,y) e Var(x) e colocá-las nessas equações. Para ilustrar, vamos considerar o exemplo 2.3 de Wooldridge (2013) sobre Salários de CEOs e Retornos sobre o patrimônio. Para isso, considere o seguinte modelo

(4)   \begin{align*} salary = \beta_0 + \beta_1 roe + u \end{align*}

onde salary é o salário anual de CEO em milhares de dólares e roe é o retorno médio sobre o patrimônio em percentual. O parâmetro \beta_1 irá medir a variação no salário anual quando o retorno médio sobre o patrimônio aumentar em um ponto percentual. Para estimar esse modelo, podemos utilizar o conjunto de dados ceosal1. Podemos dar uma olhada nas variáveis do conjunto de dados cesal1 a partir do pacote wooldridge como abaixo.


data(ceosal1, package='wooldridge')
attach(ceosal1)

Uma vez que tenhamos carregado o conjunto de dados, podemos calcular manualmente os parâmetros \beta_0 e \beta_1, como abaixo.


# Cálculo manual dos parâmetros
b1hat = cov(roe,salary)/var(roe)
b0hat = mean(salary) - b1hat*mean(roe)

Isto é, a reta de regressão será dada por

(5)   \begin{align*} \hat{salary} = 963,191 + 18,501 * roe \end{align*}

o que pode ser facilmente obtido com o código abaixo:


lm(salary ~ roe)

Implicando que para um roe = 0, teremos um salário previsto de 963,19 ou US$ 963.191, que é o intercepto. Ademais, se \Delta roe = 1, então \Delta salary = 18,5 ou US$ 18.501. Podemos, por fim, desenhar a reta de regressão com o código abaixo.


CEOregress = lm(salary ~ roe)
plot(roe, salary, ylim=c(0,4000))
abline(CEOregress, col='red')

Vamos continuar nossa revisão de modelos de regressão simples com o conjunto de dados wage1. Estamos interessados agora em estudar a relação entre educação e salários, de modo que o nosso modelo de regressão será

(6)   \begin{align*} wage = \beta_0 + \beta_1 education + u. \end{align*}

O que pode ser obtido com o código abaixo.


modelo = lm(wage ~ educ, data=wage1)
modelo

Isto é, teremos a seguinte reta de regressão

(7)   \begin{align*} \hat{wage} = -0,9 + 0,54 * education \end{align*}

de modo que um ano adicional de estudo implica em mais 54 centavos à hora de trabalho. O objeto obtido com a função lm contém todas as informações relevantes de uma regressão. Abaixo, acessamos os elementos do objeto CEOregress:


names(CEOregress)
CEOregress$coefficients

Podemos obter os valores ajustados:


plot(CEOregress$fitted.values)

E os resíduos:


plot(CEOregress$residuals)

Por fim, podemos ainda obter um sumário de todas as estatísticas relevantes da regressão com a função abaixo.


summary(CEOregress)

O que podemos gerar como tabela com o pacote stargazer como abaixo.

Dependent variable:
salary
roe 18.501*
(11.123)
Constant 963.191***
(213.240)
Observations 209
R2 0.013
Adjusted R2 0.008
Residual Std. Error 1,366.555 (df = 207)
F Statistic 2.767* (df = 1; 207)
Note: *p<0.1; **p<0.05; ***p<0.01

Gostou? Isso e muito mais você aprende em nosso Curso de Introdução à Econometria usando o R.

__________________

(*) Wooldridge, J. M. Introductory Econometrics: A Modern Approach. Editora Cengage, 2013.

Regressões espúrias: PIB per capita do Chile vs. preço do cobre

By | Dicas de R

As manifestações no Chile não param de produzir textos pelas redes sociais. A última que vi foi a relação entre o preço do cobre e o pib per capita do país vizinho. Resolvi dar uma olhada nas séries e alertar para alguns problemas nessa relação. Os dados utilizados no exercício são do FMI e estão disponíveis no final do post. Eu os importo com o código abaixo para o R.


library(readxl)
library(forecast)
library(xts)
library(ggplot2)
library(stargazer)

cobre = read_excel('cobre.xlsx', col_types = c('date', 'numeric'))
pib = window(ts(read_excel('pib.xlsx')[,-1], start=1960, freq=1),
start=1980)

A série de preço do cobre que consegui é mensal, então de modo a comparar com o pib per capita do Chile, é preciso anualizá-la. Faço isso abaixo, já pegando uma janela que eu possa comparar com o pib do Chile.


cobre = xts(cobre$cobre, order.by = cobre$date)
cobre = apply.yearly(cobre, FUN=mean)
cobre = window(cobre, start='1980-12-01', end='2018-12-01')

Com os dados devidamente comparáveis, faço um gráfico de correlação como abaixo.

A correlação entre as séries é altamente positiva, de 0,91. Os gráficos das séries são colocados em seguida.

Como podemos ver, as séries não são estacionárias. Isso, entretanto, não impede que possamos trabalhar com elas em nível, desde que exista o que chamamos de relação de longo prazo entre as mesmas. No jargão econométrico, dizemos que existe uma relação de longo prazo entre duas séries não estacionárias se as mesmas cointegram. Há algumas formas de verificar isso. Talvez a mais simples seja rodar uma regressão entre elas e verificar se os resíduos são estacionários. Se for o caso, é possível dizer que as séries cointegram. Abaixo, coloco o output de uma regressão entre o pib per capita do Chile e o preço do cobre.

Dependent variable:
pib
cobre 1.916***
(0.141)
Constant -77.583
(606.886)
Observations 39
R2 0.834
Adjusted R2 0.829
Residual Std. Error 2,049.611 (df = 37)
F Statistic 185.387*** (df = 1; 37)
Note: *p<0.1; **p<0.05; ***p<0.01

E agora, coloco um gráfico dos resíduos da regressão...

Como se vê, os resíduos da regressão estão longe de serem estacionários. O leitor mais interessado pode rodar um teste de raiz unitária sobre os mesmos. Assim, não podemos acreditar no output da regressão simplesmente porque os resultados são espúrios...

Por fim, outro alerta: não estou querendo dizer nada com esse exercício simples sobre a dependência ou não da economia chilena em relação ao cobre. Apenas que não podemos confiar em uma regressão entre o pib per capita do Chile contra o preço do Cobre... 😉

___________________________________________

(*) Quer aprender mais sobre séries temporais? Veja nosso Curso de Séries Temporais usando o R.

(**) As séries que utilizei podem ser baixadas aqui.

Cadastre-se na newsletter
e receba nossas novidades em primeira mão!