options(scipen = 999999)
library (stargazer)
library(foreign)
datos_regresion <- read.dta("https://stats.idre.ucla.edu/stat/data/crime.dta")
modelo_crime<-lm(crime~poverty+single,data=datos_regresion)
stargazer(modelo_crime,type = "html",title = "Modelo Crimen")
| 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 |
Prueba de Heterocedasticidad
Usando la libreria “lmtest”
library(lmtest)
white_test<-bptest(modelo_crime,~I(poverty^2)+I(single^2)+poverty*single,data = datos_regresion)
print(white_test)
studentized Breusch-Pagan test
data: modelo_crime BP = 10.73, df = 5, p-value = 0.057
Hay evidencia de heterocedasticidad ya que pvalue < 0.05
Prueba del Multiplicador de Lagrange (Breusch Godfrey)
Autocorrelación de 2º Orden:
library(lmtest)
prueba_LM<-bgtest(modelo_crime,order = 2)
print(prueba_LM)
Breusch-Godfrey test for serial correlation of order up to 2
data: modelo_crime LM test = 0.27165, df = 2, p-value = 0.873
No hay evidencia de Autocorrelación de 2º orden ya que pvalue > 0.05
Autocorrelación de 1º orden (prueba de Durbin Watson)
library(car)
durbinWatsonTest(model = modelo_crime)
lag Autocorrelation D-W Statistic p-value 1 -0.07014421 2.040007 0.904 Alternative hypothesis: rho != 0
No hay evidencia de Autocorrelación de 1º orden ya que pvalue > 0.05
Nota: No hay evidencia de autocorrelación, pero para la simulación del ejemplo se supondra que existe presencia de autocorrelación de primer orden, lo que se usa el Estimador HAC de NeweyWest.
Sin Corregir
options(scipen = 99999)
library(lmtest)
#Sin corregir:
coeftest(modelo_crime)
##
## 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
Corregido usando HAC
library(lmtest)
library(sandwich)
#Corregido:
# Para Autocorrelacion de Orden 1
estimacion_omega1<-vcovHC(modelo_crime,type = "HC1")
coeftest(modelo_crime,vcov. = estimacion_omega1)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1368.1887 284.9180 -4.8020 0.00001576624 ***
## poverty 6.7874 10.9273 0.6211 0.5374
## single 166.3727 26.2343 6.3418 0.00000007519 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
robust.se<-sqrt(diag(estimacion_omega1))
#Para Autocorrelacion de Orden 2
estimacion_omega<-NeweyWest(modelo_crime,lag = 2)
coeftest(modelo_crime,vcov. = estimacion_omega)
##
## 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
options(scipen = 999999)
library(robustbase)
library(stargazer)
modelo_crime_corregido<-lmrob(crime~poverty+single,data=datos_regresion,)
#Para ver el modelo original y el modelo corregido en dos columnas
stargazer(modelo_crime,modelo_crime,se=list(NULL, robust.se),type = "html",title = "Modelo de Regresion Lineal", column.labels = c("Modelo Original", "Modelo Corregido HAC") )
| Dependent variable: | ||
| crime | ||
| Modelo Original | Modelo Corregido HAC | |
| (1) | (2) | |
| poverty | 6.787 | 6.787 |
| (8.989) | (10.927) | |
| single | 166.373*** | 166.373*** |
| (19.423) | (26.234) | |
| Constant | -1,368.189*** | -1,368.189*** |
| (187.205) | (284.918) | |
| 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 | |