Em econometria queremos explicar uma variável \(y\) em função de um conjunto de regressores, ou variáveis explicativas, \(x_{1}, x_{2}, \dots, x_{k}\). Isto é, criamos um modelo matemático que explica \(y\) em função destes regressores, \(y = f(x_{1}, x_{2}, \dots, x_{k})\). A forma mais simples para \(f(\dot)\) é a linear, isto é, \[\begin{equation} y_{t} = \beta_{0} + \beta_{1}x_{1,t} + \beta_{2}x_{2,t} + \dots + \beta_{k}x_{k,t} + u_{t} \end{equation}\] onde queremos estimar os valores de \(\beta_{0}, \beta_{1}, \dots, \beta_{k}\). O último termo da equação, \(u_{t}\), é geralmente denominado “ruído” ou “erro”. Os parâmetros \(\beta_{0}, \beta_{1}, \dots, \beta_{k}\) nos informam os “efeitos marginais” que a variável \(x_{i}\) respectiva tem sobre a variável de interesse \(y\). Em alguns casos, pode-se inclusive afirmar que o valor estimado de \(\beta_{i}\) é o efeito causal que \(x_{i}\) tem sobre \(y\).

Exemplo: Cigarros

Digamos que o consumo de cigarros seja linearmente relacionado ao preço do maço. A base CigarettesB traz informações do preço do maço de cigarros e do consumo (per capita) de cigarros para 46 estados dos EUA no ano de 1992. Para carregar esta precisamos primeiro do pactoe AER. Para carregar pacotes no R use a função library(), caso ainda não tenha instalado o pacote use install.packages("AER") e depois carregue o pacote. Para carregar a base de dados usamos a função data() com o nome da base de dados como argumento, isto é, data("CigarettesB").

library("AER")
data("CigarettesB")
packs price income
AL 4.96213 0.20487 4.64039
AZ 4.66312 0.16640 4.68389
AR 5.10709 0.23406 4.59435
CA 4.50449 0.36399 4.88147
CT 4.66983 0.32149 5.09472
DE 5.04705 0.21929 4.87087
DC 4.65637 0.28946 5.05960
FL 4.80081 0.28733 4.81155
GA 4.97974 0.12826 4.73299
ID 4.74902 0.17541 4.64307

Temos 46 observações de cada uma das três variáveis. Indexando cada uma delas pela letra \(i = 1, 2, \dots, 46\) queremos estimar o seguinte modelo: \[ \text{Packs}_{i} = \beta_{0} + \beta_{1}\text{Price}_{i} + u_{t} \] Isto é, queremos estimar \(\beta_{0}\) e \(\beta_{1}\) acima. Adicionalmente, como se trata de um modelo de regressão simples (i.e., com apenas uma variável explicativa) queremos também visualizar a relação linear graficamente. Vamos começar com um gráfico de dispersão. Para montar o gráfico usamos a função plot(). Além disso, é conveniente trabalhar com um nome mais simples do que “CigarettesB”. Recomenda-se, em geral, que se crie objetos com nomes curtos e com letras minúsculas. Neste caso, uma opção intuitiva seria “cigar”.

# "Renomeia" a tabela CigarettesB
cigar <- CigarettesB
# Gráfico de dispersão
plot(packs ~ price, data = cigar, main = "Gráfico de dispersão entre consumo de cigarros e seu preço", 
     xlab = "Consumo agregado", ylab = "Renda disponível")

Para estimar o modelo usamos a função lm(). Esta função devolve um objeto do tipo “lm”, que é, essencialmente, uma lista de valores. Usando o operador “$” podemos acessar os diferentes valores guardados dentro do objeto “lm” criado. Usando a função a summary em torno do “lm” retorna uma saída com informações gerais sobre a regressão.

# Regressão linear
fit <- lm(packs ~ price, data = cigar)
# Exibe os principais resultados da regressão
summary(fit)
## 
## Call:
## lm(formula = packs ~ price, data = cigar)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.45472 -0.09968  0.00612  0.11553  0.29346 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.0941     0.0627  81.247  < 2e-16 ***
## price        -1.1983     0.2818  -4.253 0.000108 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.163 on 44 degrees of freedom
## Multiple R-squared:  0.2913, Adjusted R-squared:  0.2752 
## F-statistic: 18.08 on 1 and 44 DF,  p-value: 0.0001085

Como se vê parece haver uma relação linear significante entre o preço do pacote de cigarro e o seu consumo per-capita. Em particular, esta relação é negativa, como seria esperado. Para visualizar graficamente o ajuste do modelo usamos a função abline().

# Regressão linear
plot(packs ~ price, data = cigar, main = "Gráfico de dispersão entre consumo de cigarros e seu preço", 
     xlab = "Consumo agregado", ylab = "Renda disponível")
abline(fit$coefficients, col = 'red')

residuo <- resid(fit)
hist(residuo)

coef(fit)
## (Intercept)       price 
##    5.094108   -1.198316
str(summary(fit))
## List of 11
##  $ call         : language lm(formula = packs ~ price, data = cigar)
##  $ terms        :Classes 'terms', 'formula'  language packs ~ price
##   .. ..- attr(*, "variables")= language list(packs, price)
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "packs" "price"
##   .. .. .. ..$ : chr "price"
##   .. ..- attr(*, "term.labels")= chr "price"
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(packs, price)
##   .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. ..- attr(*, "names")= chr [1:2] "packs" "price"
##  $ residuals    : Named num [1:46] 0.114 -0.232 0.293 -0.153 -0.039 ...
##   ..- attr(*, "names")= chr [1:46] "AL" "AZ" "AR" "CA" ...
##  $ coefficients : num [1:2, 1:4] 5.0941 -1.1983 0.0627 0.2818 81.2471 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "(Intercept)" "price"
##   .. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)"
##  $ aliased      : Named logi [1:2] FALSE FALSE
##   ..- attr(*, "names")= chr [1:2] "(Intercept)" "price"
##  $ sigma        : num 0.163
##  $ df           : int [1:3] 2 44 2
##  $ r.squared    : num 0.291
##  $ adj.r.squared: num 0.275
##  $ fstatistic   : Named num [1:3] 18.1 1 44
##   ..- attr(*, "names")= chr [1:3] "value" "numdf" "dendf"
##  $ cov.unscaled : num [1:2, 1:2] 0.148 -0.614 -0.614 2.989
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "(Intercept)" "price"
##   .. ..$ : chr [1:2] "(Intercept)" "price"
##  - attr(*, "class")= chr "summary.lm"
summary(fit)$r.squared
## [1] 0.2912836
str(fit)
## List of 12
##  $ coefficients : Named num [1:2] 5.09 -1.2
##   ..- attr(*, "names")= chr [1:2] "(Intercept)" "price"
##  $ residuals    : Named num [1:46] 0.114 -0.232 0.293 -0.153 -0.039 ...
##   ..- attr(*, "names")= chr [1:46] "AL" "AZ" "AR" "CA" ...
##  $ effects      : Named num [1:46] -32.8797 -0.6932 0.2675 -0.2312 -0.0999 ...
##   ..- attr(*, "names")= chr [1:46] "(Intercept)" "price" "" "" ...
##  $ rank         : int 2
##  $ fitted.values: Named num [1:46] 4.85 4.89 4.81 4.66 4.71 ...
##   ..- attr(*, "names")= chr [1:46] "AL" "AZ" "AR" "CA" ...
##  $ assign       : int [1:2] 0 1
##  $ qr           :List of 5
##   ..$ qr   : num [1:46, 1:2] -6.782 0.147 0.147 0.147 0.147 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:46] "AL" "AZ" "AR" "CA" ...
##   .. .. ..$ : chr [1:2] "(Intercept)" "price"
##   .. ..- attr(*, "assign")= int [1:2] 0 1
##   ..$ qraux: num [1:2] 1.15 1.07
##   ..$ pivot: int [1:2] 1 2
##   ..$ tol  : num 1e-07
##   ..$ rank : int 2
##   ..- attr(*, "class")= chr "qr"
##  $ df.residual  : int 44
##  $ xlevels      : Named list()
##  $ call         : language lm(formula = packs ~ price, data = cigar)
##  $ terms        :Classes 'terms', 'formula'  language packs ~ price
##   .. ..- attr(*, "variables")= language list(packs, price)
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "packs" "price"
##   .. .. .. ..$ : chr "price"
##   .. ..- attr(*, "term.labels")= chr "price"
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(packs, price)
##   .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. ..- attr(*, "names")= chr [1:2] "packs" "price"
##  $ model        :'data.frame':   46 obs. of  2 variables:
##   ..$ packs: num [1:46] 4.96 4.66 5.11 4.5 4.67 ...
##   ..$ price: num [1:46] 0.205 0.166 0.234 0.364 0.321 ...
##   ..- attr(*, "terms")=Classes 'terms', 'formula'  language packs ~ price
##   .. .. ..- attr(*, "variables")= language list(packs, price)
##   .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:2] "packs" "price"
##   .. .. .. .. ..$ : chr "price"
##   .. .. ..- attr(*, "term.labels")= chr "price"
##   .. .. ..- attr(*, "order")= int 1
##   .. .. ..- attr(*, "intercept")= int 1
##   .. .. ..- attr(*, "response")= int 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. ..- attr(*, "predvars")= language list(packs, price)
##   .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. .. ..- attr(*, "names")= chr [1:2] "packs" "price"
##  - attr(*, "class")= chr "lm"
y_estimado <- fitted(fit)
plot(y_estimado, cigar$packs)
abline(a = 0, b = 1, lty = 2)

plot(residuo, cigar$price)
abline(h = 0)
abline(v = 0)
grid()

plot(fit)

predict(fit, data.frame(price = 0.5), se.fit = TRUE)
## $fit
##       1 
## 4.49495 
## 
## $se.fit
## [1] 0.08639435
## 
## $df
## [1] 44
## 
## $residual.scale
## [1] 0.1630009

Regressão múltipla

Exemplo: cigarros (de novo)

# Regressão linear
fit <- lm(packs ~ price + income, data = cigar)
# Exibe os principais resultados da regressão
summary(fit)
## 
## Call:
## lm(formula = packs ~ price + income, data = cigar)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.41867 -0.10683  0.00757  0.11738  0.32868 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.2997     0.9089   4.730 2.43e-05 ***
## price        -1.3383     0.3246  -4.123 0.000168 ***
## income        0.1724     0.1968   0.876 0.385818    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1634 on 43 degrees of freedom
## Multiple R-squared:  0.3037, Adjusted R-squared:  0.2713 
## F-statistic: 9.378 on 2 and 43 DF,  p-value: 0.0004168

Exemplo: salário

library(wooldridge)
data("wage2")
wage hours IQ KWW educ exper tenure age married black south urban sibs brthord meduc feduc lwage
769 40 93 35 12 11 2 31 1 0 0 1 1 2 8 8 6.645091
808 50 119 41 18 11 16 37 1 0 0 1 1 NA 14 14 6.694562
825 40 108 46 14 11 9 33 1 0 0 1 1 2 14 14 6.715383
650 40 96 32 12 13 7 32 1 0 0 1 4 3 12 12 6.476973
562 40 74 27 11 14 5 34 1 0 0 1 10 6 6 11 6.331502
1400 40 116 43 16 14 2 35 1 1 0 1 1 2 8 NA 7.244227
names(wage2)
##  [1] "wage"    "hours"   "IQ"      "KWW"     "educ"    "exper"   "tenure" 
##  [8] "age"     "married" "black"   "south"   "urban"   "sibs"    "brthord"
## [15] "meduc"   "feduc"   "lwage"
fit <- lm(lwage ~ educ + exper + IQ + tenure + age + married, data = wage2)

Não-linearidades

A restrição de linearidade do modelo é somente sobre os parâmetros, isto é, sobre os betas. Pode-se incorporar não-linearidades de outras forma. O modelo abaixo, por exemplo, continua sendo um modelo de regressão linear:

\[ y_{t} = \beta_{0} + \beta_{1}x_{1t}^{3} + \beta_{2}\frac{1}{x_{2t}} + \beta_{3}\text{log}(x_{3t}) + \epsilon_{t} \]

Logaritmos e potências

library("wooldridge")
data("hprice2")
fit <- lm(log(price) ~ log(nox) + log(dist) + rooms + I(rooms^2) + stratio,
          data = hprice2)
summary(fit)
## 
## Call:
## lm(formula = log(price) ~ log(nox) + log(dist) + rooms + I(rooms^2) + 
##     stratio, data = hprice2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.04285 -0.12774  0.02038  0.12650  1.25272 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 13.385477   0.566473  23.630  < 2e-16 ***
## log(nox)    -0.901682   0.114687  -7.862 2.34e-14 ***
## log(dist)   -0.086781   0.043281  -2.005  0.04549 *  
## rooms       -0.545113   0.165454  -3.295  0.00106 ** 
## I(rooms^2)   0.062261   0.012805   4.862 1.56e-06 ***
## stratio     -0.047590   0.005854  -8.129 3.42e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2592 on 500 degrees of freedom
## Multiple R-squared:  0.6028, Adjusted R-squared:  0.5988 
## F-statistic: 151.8 on 5 and 500 DF,  p-value: < 2.2e-16

\[ \frac{\partial log(y_{t})}{\partial log(x_{1t})} \]

Interações entre variáveis

Exemplo: variáveis dummy

data("lawsch85")
d <- lawsch85
d$top10   <- ifelse(test = d$rank <= 10, yes = 1, no = 0)
d$r11_25  <- ifelse(test = d$rank > 10 & d$rank <= 25, yes = 1, no = 0)
d$r26_40  <- ifelse(test = d$rank > 25 & d$rank <= 40, yes = 1, no = 0)
d$r41_60  <- ifelse(test = d$rank > 40 & d$rank <= 60, yes = 1, no = 0)
d$r61_100 <- ifelse(test = d$rank > 60 & d$rank <= 100, yes = 1, no = 0)

fit <- lm(log(salary) ~ top10 + r11_25 + r26_40 + r41_60 + r61_100 + LSAT +
                        GPA + log(llibvol) + lcost,
          data = d)
fit
## 
## Call:
## lm(formula = log(salary) ~ top10 + r11_25 + r26_40 + r41_60 + 
##     r61_100 + LSAT + GPA + log(llibvol) + lcost, data = d)
## 
## Coefficients:
##  (Intercept)         top10        r11_25        r26_40        r41_60  
##    9.0089365     0.7005341     0.5937752     0.3747281     0.2623565  
##      r61_100          LSAT           GPA  log(llibvol)         lcost  
##    0.1312513     0.0056522     0.0147991     0.2110607     0.0007824

Testes

library(wooldridge)
data("wage2")
wage hours IQ KWW educ exper tenure age married black south urban sibs brthord meduc feduc lwage
769 40 93 35 12 11 2 31 1 0 0 1 1 2 8 8 6.645091
808 50 119 41 18 11 16 37 1 0 0 1 1 NA 14 14 6.694562
825 40 108 46 14 11 9 33 1 0 0 1 1 2 14 14 6.715383
650 40 96 32 12 13 7 32 1 0 0 1 4 3 12 12 6.476973
562 40 74 27 11 14 5 34 1 0 0 1 10 6 6 11 6.331502
1400 40 116 43 16 14 2 35 1 1 0 1 1 2 8 NA 7.244227

Usando o stargazer

O pacote stargazer facilita muito a publicação dos resultados de uma regressão. Depois de carregar o pacote simplesmente use a função stargazer() num objeto lm. Lembre-se que para verificar a classe de um objeto usamos a função class().

# Carregar pacote adicional de temas
class(fit)
## [1] "lm"

Note que abaixo chamo a função usando o argumento type = 'html', pois estou usando o RMarkdown para fazer este site. Outras opções incluem type = 'text', que retorna uma tabela simples estilo .txt, e type = 'latex', que retorna uma tabela no formato LaTeX.

# Carregar pacote adicional de temas
stargazer(fit,
          type = "html",
          title = "Resultados da regressão",
          covariate.labels = c("Preço", "Renda individual", "constante"),
          dep.var.labels = "Pacotes de cigarro"
          )
Resultados da regressão
Dependent variable:
Pacotes de cigarro
Preço 0.701***
(0.053)
Renda individual 0.594***
(0.039)
constante 0.375***
(0.034)
r41_60 0.262***
(0.028)
r61_100 0.131***
(0.021)
LSAT 0.006*
(0.003)
GPA 0.015
(0.074)
log(llibvol) 0.211
(0.152)
lcost 0.001
(0.025)
Constant 9.009***
(0.445)
Observations 136
R2 0.911
Adjusted R2 0.905
Residual Std. Error 0.086 (df = 126)
F Statistic 143.176*** (df = 9; 126)
Note: p<0.1; p<0.05; p<0.01

Usando o ggplot2

O pacote ggplot2 oferece formas alternativas de visualizar a mesma relação linear que encontramos acima. Em particular, a função qplot() tem uma sintaxe muito similar à da função plot().

# Carrega o pacote
library(ggplot2)
# Função qplot
p <- qplot(x = price, y = packs, data = cigar, geom = "point")
p

Para adicionar uma linha de regressão acrescentamos geom_smooth. Especificamos que queremos um modelo linear adicionando o argumento ‘lm’

p + geom_smooth(method = "lm")

Uma forma simples de gerar gráficos bonitos com o ggplot2 é tomando vantagem dos seus temas pré-definidos.

# Carregar pacote adicional de temas
library(ggthemes)
p + geom_smooth(method = "lm") + theme_gdocs()

p + geom_smooth(method = "lm") + theme_economist()