library(stargazer)
## 
## Please cite as:
##  Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
library(sandwich)
## Warning: package 'sandwich' was built under R version 4.5.3
library(lmtest)
## Cargando paquete requerido: zoo
## 
## Adjuntando el paquete: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
getwd()
## [1] "C:/Users/laura/OneDrive/Econometria"
load("C:/Users/laura/OneDrive/Econometria/smoke.RData")
desc
##    variable                         label
## 1      educ            years of schooling
## 2   cigpric  state cig. price, cents/pack
## 3     white                   =1 if white
## 4       age                      in years
## 5    income              annual income, $
## 6      cigs          cigs. smoked per day
## 7  restaurn =1 if rest. smk. restrictions
## 8   lincome                   log(income)
## 9     agesq                         age^2
## 10 lcigpric                 log(cigprice)
str(data)
## 'data.frame':    807 obs. of  10 variables:
##  $ educ    : num  16 16 12 13.5 10 6 12 15 12 12 ...
##  $ cigpric : num  60.5 57.9 57.7 57.9 58.3 ...
##  $ white   : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ age     : int  46 40 58 30 17 86 35 48 48 31 ...
##  $ income  : int  20000 30000 30000 20000 20000 6500 20000 30000 20000 20000 ...
##  $ cigs    : int  0 0 3 0 0 0 0 0 0 0 ...
##  $ restaurn: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ lincome : num  9.9 10.3 10.3 9.9 9.9 ...
##  $ agesq   : int  2116 1600 3364 900 289 7396 1225 2304 2304 961 ...
##  $ lcigpric: num  4.1 4.06 4.05 4.06 4.07 ...
##  - attr(*, "datalabel")= chr ""
##  - attr(*, "time.stamp")= chr "25 Jun 2011 23:03"
##  - attr(*, "formats")= chr [1:10] "%9.0g" "%9.0g" "%8.0g" "%8.0g" ...
##  - attr(*, "types")= int [1:10] 254 254 251 251 252 251 251 254 252 254
##  - attr(*, "val.labels")= chr [1:10] "" "" "" "" ...
##  - attr(*, "var.labels")= chr [1:10] "years of schooling" "state cig. price, cents/pack" "=1 if white" "in years" ...
##  - attr(*, "version")= int 10
head(data)
##   educ cigpric white age income cigs restaurn   lincome agesq lcigpric
## 1 16.0  60.506     1  46  20000    0        0  9.903487  2116 4.102743
## 2 16.0  57.883     1  40  30000    0        0 10.308952  1600 4.058424
## 3 12.0  57.664     1  58  30000    3        0 10.308952  3364 4.054633
## 4 13.5  57.883     1  30  20000    0        0  9.903487   900 4.058424
## 5 10.0  58.320     1  17  20000    0        0  9.903487   289 4.065945
## 6  6.0  59.340     1  86   6500    0        0  8.779557  7396 4.083283

Estimar el modelo lineal

modelo_original <- lm(
  cigs ~ cigpric + lcigpric + income + lincome +
    age + agesq + educ + white + restaurn,
  data = data
)

summary(modelo_original)
## 
## Call:
## lm(formula = cigs ~ cigpric + lcigpric + income + lincome + age + 
##     agesq + educ + white + restaurn, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.169  -9.357  -5.915   7.851  70.744 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.408e+02  2.600e+02   1.311  0.19033    
## cigpric      2.002e+00  1.493e+00   1.341  0.18022    
## lcigpric    -1.153e+02  8.542e+01  -1.349  0.17758    
## income      -4.619e-05  1.335e-04  -0.346  0.72940    
## lincome      1.404e+00  1.708e+00   0.822  0.41134    
## age          7.784e-01  1.606e-01   4.848 1.50e-06 ***
## agesq       -9.150e-03  1.749e-03  -5.231 2.16e-07 ***
## educ        -4.948e-01  1.682e-01  -2.942  0.00336 ** 
## white       -5.311e-01  1.461e+00  -0.364  0.71629    
## restaurn    -2.644e+00  1.130e+00  -2.340  0.01953 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.41 on 797 degrees of freedom
## Multiple R-squared:  0.05515,    Adjusted R-squared:  0.04448 
## F-statistic: 5.169 on 9 and 797 DF,  p-value: 7.735e-07
stargazer(
  modelo_original,
  type = "text",
  title = "Resumen del Modelo Lineal",
  digits = 4
)
## 
## Resumen del Modelo Lineal
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                cigs            
## -----------------------------------------------
## cigpric                       2.0023           
##                              (1.4928)          
##                                                
## lcigpric                     -115.2735         
##                              (85.4243)         
##                                                
## income                       -0.00005          
##                              (0.0001)          
##                                                
## lincome                       1.4041           
##                              (1.7082)          
##                                                
## age                          0.7784***         
##                              (0.1606)          
##                                                
## agesq                       -0.0092***         
##                              (0.0017)          
##                                                
## educ                        -0.4948***         
##                              (0.1682)          
##                                                
## white                         -0.5311          
##                              (1.4607)          
##                                                
## restaurn                     -2.6442**         
##                              (1.1300)          
##                                                
## Constant                     340.8044          
##                             (260.0156)         
##                                                
## -----------------------------------------------
## Observations                    807            
## R2                            0.0552           
## Adjusted R2                   0.0445           
## Residual Std. Error     13.4128 (df = 797)     
## F Statistic           5.1693*** (df = 9; 797)  
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

Matriz Varianza-Covarianza

# Matriz Varianza-Covarianza
vcov(modelo_original)
##               (Intercept)       cigpric      lcigpric        income
## (Intercept)  6.760811e+04  3.861092e+02 -2.217382e+04  4.361954e-04
## cigpric      3.861092e+02  2.228545e+00 -1.272304e+02 -6.112577e-06
## lcigpric    -2.217382e+04 -1.272304e+02  7.297314e+03  3.973538e-04
## income       4.361954e-04 -6.112577e-06  3.973538e-04  1.781986e-08
## lincome     -8.178466e+00  7.124578e-02 -4.860425e+00 -2.062280e-04
## age          4.955739e-01  3.684842e-03 -2.196985e-01 -2.490054e-07
## agesq       -1.237179e-02 -8.003706e-05  4.619639e-03  4.523652e-09
## educ         9.421945e-02  1.549304e-03 -1.048653e-01 -2.445251e-06
## white       -2.254192e+00  1.177281e-02 -3.119456e-01 -7.733768e-06
## restaurn     4.064836e+01  2.186489e-01 -1.340589e+01 -1.150501e-05
##                   lincome           age         agesq          educ
## (Intercept) -8.1784660668  4.955739e-01 -1.237179e-02  9.421945e-02
## cigpric      0.0712457770  3.684842e-03 -8.003706e-05  1.549304e-03
## lcigpric    -4.8604250078 -2.196985e-01  4.619639e-03 -1.048653e-01
## income      -0.0002062280 -2.490054e-07  4.523652e-09 -2.445251e-06
## lincome      2.9178305410 -2.864418e-02  2.930256e-04 -2.617669e-03
## age         -0.0286441758  2.577810e-02 -2.764371e-04 -2.949036e-03
## agesq        0.0002930256 -2.764371e-04  3.060024e-06  4.144658e-05
## educ        -0.0026176694 -2.949036e-03  4.144658e-05  2.828458e-02
## white        0.1323738372 -1.438531e-02  1.721940e-04  2.035627e-03
## restaurn     0.0828693504  3.811336e-03 -3.412736e-05 -3.821198e-03
##                     white      restaurn
## (Intercept) -2.254192e+00  4.064836e+01
## cigpric      1.177281e-02  2.186489e-01
## lcigpric    -3.119456e-01 -1.340589e+01
## income      -7.733768e-06 -1.150501e-05
## lincome      1.323738e-01  8.286935e-02
## age         -1.438531e-02  3.811336e-03
## agesq        1.721940e-04 -3.412736e-05
## educ         2.035627e-03 -3.821198e-03
## white        2.133708e+00  1.602894e-01
## restaurn     1.602894e-01  1.276897e+00

Verificar heterocedasticidad

# Prueba Breusch-Pagan
bptest(modelo_original)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_original
## BP = 33.525, df = 9, p-value = 0.0001082

Estimador HAC

# Matriz HAC
vcov_HAC <- vcovHC(modelo_original, type = "HC1")

# Ver matriz HAC
vcov_HAC
##               (Intercept)       cigpric      lcigpric        income
## (Intercept)  7.857213e+04  4.500587e+02 -2.574203e+04  1.043757e-03
## cigpric      4.500587e+02  2.600967e+00 -1.479162e+02  1.565177e-06
## lcigpric    -2.574203e+04 -1.479162e+02  8.448494e+03 -2.925266e-05
## income       1.043757e-03  1.565177e-06 -2.925265e-05  1.352933e-08
## lincome     -3.166179e+01 -1.290705e-01  6.739374e+00 -1.291708e-04
## age          3.133411e+00  2.102700e-02 -1.127132e+00 -5.415920e-07
## agesq       -4.167356e-02 -2.727534e-04  1.476830e-02  9.607681e-09
## educ        -1.437301e+00 -6.290098e-03  3.913626e-01 -2.262746e-06
## white       -2.062973e+01 -9.078995e-02  6.055433e+00  9.542122e-06
## restaurn     5.600759e+01  3.145252e-01 -1.870008e+01 -2.382794e-05
##                   lincome           age         agesq          educ
## (Intercept) -3.166179e+01  3.133411e+00 -4.167356e-02 -1.437301e+00
## cigpric     -1.290705e-01  2.102700e-02 -2.727534e-04 -6.290098e-03
## lcigpric     6.739374e+00 -1.127132e+00  1.476830e-02  3.913626e-01
## income      -1.291708e-04 -5.415920e-07  9.607681e-09 -2.262746e-06
## lincome      1.529320e+00 -1.314408e-02  9.307414e-05 -7.226153e-03
## age         -1.314408e-02  1.898915e-02 -1.985813e-04 -7.123192e-04
## agesq        9.307414e-05 -1.985813e-04  2.133427e-06  1.537531e-05
## educ        -7.226153e-03 -7.123192e-04  1.537531e-05  2.689192e-02
## white       -4.226208e-02 -1.681454e-03  2.954330e-05 -1.122063e-02
## restaurn     1.781647e-01 -3.662344e-04 -3.904902e-06 -9.909995e-03
##                     white      restaurn
## (Intercept) -2.062973e+01  5.600759e+01
## cigpric     -9.078995e-02  3.145252e-01
## lcigpric     6.055433e+00 -1.870008e+01
## income       9.542122e-06 -2.382794e-05
## lincome     -4.226208e-02  1.781647e-01
## age         -1.681454e-03 -3.662344e-04
## agesq        2.954330e-05 -3.904902e-06
## educ        -1.122063e-02 -9.909995e-03
## white        1.878065e+00  1.707809e-01
## restaurn     1.707809e-01  1.091499e+00

Eerrores estándar robustos

# Errores estándar robustos HAC
errores_HAC <- sqrt(diag(vcov_HAC))

errores_HAC
##  (Intercept)      cigpric     lcigpric       income      lincome          age 
## 2.803072e+02 1.612751e+00 9.191569e+01 1.163156e-04 1.236657e+00 1.378011e-01 
##        agesq         educ        white     restaurn 
## 1.460626e-03 1.639876e-01 1.370425e+00 1.044748e+00

Tabla comparativa ORIGINAL vs CORREGIDO

stargazer(
  modelo_original,
  modelo_original,
  se = list(NULL, errores_HAC),
  column.labels = c("Original", "Corregido HAC"),
  type = "text",
  title = "Comparación Modelo Original vs Modelo Corregido",
  digits = 4,
  align = TRUE
)
## 
## Comparación Modelo Original vs Modelo Corregido
## ===========================================================
##                                    Dependent variable:     
##                                ----------------------------
##                                            cigs            
##                                   Original    Corregido HAC
##                                     (1)            (2)     
## -----------------------------------------------------------
## cigpric                            2.0023        2.0023    
##                                   (1.4928)      (1.6128)   
##                                                            
## lcigpric                         -115.2735      -115.2735  
##                                  (85.4243)      (91.9157)  
##                                                            
## income                            -0.00005      -0.00005   
##                                   (0.0001)      (0.0001)   
##                                                            
## lincome                            1.4041        1.4041    
##                                   (1.7082)      (1.2367)   
##                                                            
## age                              0.7784***      0.7784***  
##                                   (0.1606)      (0.1378)   
##                                                            
## agesq                            -0.0092***    -0.0092***  
##                                   (0.0017)      (0.0015)   
##                                                            
## educ                             -0.4948***    -0.4948***  
##                                   (0.1682)      (0.1640)   
##                                                            
## white                             -0.5311        -0.5311   
##                                   (1.4607)      (1.3704)   
##                                                            
## restaurn                         -2.6442**      -2.6442**  
##                                   (1.1300)      (1.0447)   
##                                                            
## Constant                          340.8044      340.8044   
##                                  (260.0156)    (280.3072)  
##                                                            
## -----------------------------------------------------------
## Observations                        807            807     
## R2                                 0.0552        0.0552    
## Adjusted R2                        0.0445        0.0445    
## Residual Std. Error (df = 797)    13.4128        13.4128   
## F Statistic (df = 9; 797)        5.1693***      5.1693***  
## ===========================================================
## Note:                           *p<0.1; **p<0.05; ***p<0.01

Corrección HAC

# Para corregir el problema se utilizó un estimador robusto HAC mediante:

vcovHC(modelo_original, type = "HC1")
##               (Intercept)       cigpric      lcigpric        income
## (Intercept)  7.857213e+04  4.500587e+02 -2.574203e+04  1.043757e-03
## cigpric      4.500587e+02  2.600967e+00 -1.479162e+02  1.565177e-06
## lcigpric    -2.574203e+04 -1.479162e+02  8.448494e+03 -2.925266e-05
## income       1.043757e-03  1.565177e-06 -2.925265e-05  1.352933e-08
## lincome     -3.166179e+01 -1.290705e-01  6.739374e+00 -1.291708e-04
## age          3.133411e+00  2.102700e-02 -1.127132e+00 -5.415920e-07
## agesq       -4.167356e-02 -2.727534e-04  1.476830e-02  9.607681e-09
## educ        -1.437301e+00 -6.290098e-03  3.913626e-01 -2.262746e-06
## white       -2.062973e+01 -9.078995e-02  6.055433e+00  9.542122e-06
## restaurn     5.600759e+01  3.145252e-01 -1.870008e+01 -2.382794e-05
##                   lincome           age         agesq          educ
## (Intercept) -3.166179e+01  3.133411e+00 -4.167356e-02 -1.437301e+00
## cigpric     -1.290705e-01  2.102700e-02 -2.727534e-04 -6.290098e-03
## lcigpric     6.739374e+00 -1.127132e+00  1.476830e-02  3.913626e-01
## income      -1.291708e-04 -5.415920e-07  9.607681e-09 -2.262746e-06
## lincome      1.529320e+00 -1.314408e-02  9.307414e-05 -7.226153e-03
## age         -1.314408e-02  1.898915e-02 -1.985813e-04 -7.123192e-04
## agesq        9.307414e-05 -1.985813e-04  2.133427e-06  1.537531e-05
## educ        -7.226153e-03 -7.123192e-04  1.537531e-05  2.689192e-02
## white       -4.226208e-02 -1.681454e-03  2.954330e-05 -1.122063e-02
## restaurn     1.781647e-01 -3.662344e-04 -3.904902e-06 -9.909995e-03
##                     white      restaurn
## (Intercept) -2.062973e+01  5.600759e+01
## cigpric     -9.078995e-02  3.145252e-01
## lcigpric     6.055433e+00 -1.870008e+01
## income       9.542122e-06 -2.382794e-05
## lincome     -4.226208e-02  1.781647e-01
## age         -1.681454e-03 -3.662344e-04
## agesq        2.954330e-05 -3.904902e-06
## educ        -1.122063e-02 -9.909995e-03
## white        1.878065e+00  1.707809e-01
## restaurn     1.707809e-01  1.091499e+00