Estimacion del modelo
options(scipen = 999999)
library(foreign)
datos_regresion <- read.dta("https://stats.idre.ucla.edu/stat/data/crime.dta")
modelo_estimado<-lm(crime~poverty+single,data=datos_regresion)
print(modelo_estimado)
##
## Call:
## lm(formula = crime ~ poverty + single, data = datos_regresion)
##
## Coefficients:
## (Intercept) poverty single
## -1368.189 6.787 166.373
Estimacion Robusta
options(scipen = 999999)
library(lmtest)
library(sandwich)
#Sin corregir:
coeftest(modelo_estimado)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1368.1887 187.2052 -7.3085 0.00000000247861 ***
## poverty 6.7874 8.9885 0.7551 0.4539
## single 166.3727 19.4229 8.5658 0.00000000003117 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
con el uso de la biblioteca Stargazer
options(scipen = 999999)
library(stargazer)
stargazer(modelo_estimado, type = "text", title = "MODELO CRIME")
##
## MODELO CRIME
## ===============================================
## Dependent variable:
## ---------------------------
## crime
## -----------------------------------------------
## poverty 6.787
## (8.989)
##
## single 166.373***
## (19.423)
##
## Constant -1,368.189***
## (187.205)
##
## -----------------------------------------------
## Observations 51
## R2 0.707
## Adjusted R2 0.695
## Residual Std. Error 243.610 (df = 48)
## F Statistic 57.964*** (df = 2; 48)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01
usando un estimador HAC ya corregido
options(scipen = 99999)
estimacion_omega<-vcovHC(modelo_estimado,type = "HC0")
coeftest(modelo_estimado,vcov. = estimacion_omega)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1368.1887 276.4111 -4.9498 0.00000956181 ***
## poverty 6.7874 10.6010 0.6403 0.5251
## single 166.3727 25.4510 6.5370 0.00000003774 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Usando Stargazer
library(stargazer)
# Errores estandar
errores <- sqrt(diag(estimacion_omega))
stargazer(modelo_estimado, modelo_estimado,
se=list(NULL, errores), # Agregando la correción HAC
type = "text",
title = "Comparación de modelos: corrección HAC.")
##
## Comparación de modelos: corrección HAC.
## ==========================================================
## Dependent variable:
## ----------------------------
## crime
## (1) (2)
## ----------------------------------------------------------
## poverty 6.787 6.787
## (8.989) (10.601)
##
## single 166.373*** 166.373***
## (19.423) (25.451)
##
## Constant -1,368.189*** -1,368.189***
## (187.205) (276.411)
##
## ----------------------------------------------------------
## Observations 51 51
## R2 0.707 0.707
## Adjusted R2 0.695 0.695
## Residual Std. Error (df = 48) 243.610 243.610
## F Statistic (df = 2; 48) 57.964*** 57.964***
## ==========================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
Estimador HAC de NeweyWest
#Nota, si hay correlacion en el orden 2 se corrge de esta manera
library(lmtest)
library(sandwich)
estimacion_omega2<- NeweyWest(modelo_estimado, lag = 2)
coeftest(modelo_estimado, vcov. = estimacion_omega2)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1368.1887 303.8466 -4.5029 0.00004279768 ***
## poverty 6.7874 10.5943 0.6407 0.5248
## single 166.3727 25.9154 6.4198 0.00000005708 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(stargazer)
#Errores estándar
errores2 <- sqrt(diag(estimacion_omega2))
# Comparacion de Modelos
stargazer(modelo_estimado, modelo_estimado,
se=list(NULL, errores2), # Agregando la correción HAC
type = "text",
title = "Comparación de modelos: corrección HAC de NeweyWest.")
##
## Comparación de modelos: corrección HAC de NeweyWest.
## ==========================================================
## Dependent variable:
## ----------------------------
## crime
## (1) (2)
## ----------------------------------------------------------
## poverty 6.787 6.787
## (8.989) (10.594)
##
## single 166.373*** 166.373***
## (19.423) (25.915)
##
## Constant -1,368.189*** -1,368.189***
## (187.205) (303.847)
##
## ----------------------------------------------------------
## Observations 51 51
## R2 0.707 0.707
## Adjusted R2 0.695 0.695
## Residual Std. Error (df = 48) 243.610 243.610
## F Statistic (df = 2; 48) 57.964*** 57.964***
## ==========================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
Estimacion Robusta
options(scipen = 999999)
library(robustbase)
library(stargazer)
modelo_crime_robust <- lmrob(crime~poverty+single, data = datos_regresion)
stargazer(modelo_estimado, modelo_crime_robust, type = "text", title = "comparativa")
##
## comparativa
## ==================================================================
## Dependent variable:
## ------------------------------------
## crime
## OLS MM-type
## linear
## (1) (2)
## ------------------------------------------------------------------
## poverty 6.787 11.466
## (8.989) (9.263)
##
## single 166.373*** 176.569***
## (19.423) (23.223)
##
## Constant -1,368.189*** -1,539.640***
## (187.205) (235.765)
##
## ------------------------------------------------------------------
## Observations 51 51
## R2 0.707 0.795
## Adjusted R2 0.695 0.787
## Residual Std. Error (df = 48) 243.610 191.864
## F Statistic 57.964*** (df = 2; 48)
## ==================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01