Estimación del modelo
options(scipen = 99999)
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.373Estimación Robusta
options(scipen = 99999)
library(lmtest)## Warning: package 'lmtest' was built under R version 4.0.5## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.0.5##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(sandwich)## Warning: package 'sandwich' was built under R version 4.0.5#Sin la correción
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#Usando stargazer
library(stargazer)##
## Please cite as:## Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.2. https://CRAN.R-project.org/package=stargazerstargazer(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.01Con la correción y usando un estimador HAC
options(scipen = 99999)
#Con la correción:
#HC0 corrige Heterocedasticidad
#HC1 corrige tambien Autocorrelacion orden 1
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)
#Error estandar
error <- sqrt(diag(estimacion.omega))
stargazer(modelo.estimado, modelo.estimado,
se=list(NULL, error), # 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.01Estimador HAC de NeweyWest
#Al detectar correlación de segundo orden, este se corrige de la siguiente 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#Usando stargazer
library(stargazer)
#Obtener el error estándar
error2 <- sqrt(diag(estimacion.omega2))
#Comparando los modelos
stargazer(modelo.estimado, modelo.estimado,
se=list(NULL, error2), # 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.01Estimación Robusta
options(scipen = 999999)
library(robustbase)## Warning: package 'robustbase' was built under R version 4.0.5library(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