Práctica 7. Regresión lineal múltiple
AE
20/01/2020
Previo
¡Recuerda poner el directorio!
setwd("C:/Users/AULA 237/Desktop/RInvierno2020")Volvemos a importar la base completa del cuestionario sociodemográfico de la Encuesta Nacional de Ocupación y Empleo, trimestre III de 2019. O si guardaste tu ambiente, vuévelo a cargar.
library(haven)## Warning: package 'haven' was built under R version 3.6.2SDEMT319 <- read_dta("SDEMT319.dta")Paquetes y librerÃas
Vamos a instalar un par de nuevas librerÃas. Por ello haremos un vector que las liste para instalarlas más fácilmente. Este es un ejemplo de cómo podemos instalar paquetes. RCurl pedirá autorización para instalar su forma binaria. Escribe “Yes”.
nuevo<-c("estout", "car", "robustbase","RCurl", "lm.beta")
#install.packages(nuevo, dependencies = TRUE)
# Aségurate de descomentar este comando para instalar los paquetesLuego cargamos algunas librerÃas. Iremos cargando más a lo largo de la práctica
library(tidyverse) #porque usaremos dplyr y por cualquier cosa## Warning: package 'tidyverse' was built under R version 3.6.2## -- Attaching packages -------------------- tidyverse 1.3.0 --## v ggplot2 3.2.1 v purrr 0.3.3
## v tibble 2.1.3 v dplyr 0.8.3
## v tidyr 1.0.0 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.4.0## Warning: package 'purrr' was built under R version 3.6.2## -- Conflicts ----------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(sjlabelled)##
## Attaching package: 'sjlabelled'## The following object is masked from 'package:forcats':
##
## as_factor## The following object is masked from 'package:dplyr':
##
## as_label## The following objects are masked from 'package:haven':
##
## as_factor, read_sas, read_spss, read_stata, write_sas,
## zap_labelslibrary(DescTools)## Warning: package 'DescTools' was built under R version 3.6.2Sub-setting para comparar modelos
Vamos a hacer una sub-base de nuestras posibles variables explicativas. Esto es importante porque sólo podemos comparar modelos con la misma cantidad de observaciones.
mydata<- SDEMT319[SDEMT319$clase2==1 & SDEMT319$anios_esc<99 & SDEMT319$ing_x_hrs>0 & SDEMT319$eda<98,
c("ing_x_hrs", "t_loc","eda", "sex", "anios_esc", "fac")]
tail(mydata)## # A tibble: 6 x 6
## ing_x_hrs t_loc eda sex anios_esc fac
## <dbl> <dbl+lbl> <dbl> <dbl+lbl> <dbl> <dbl>
## 1 21.4 4 [Localidades menores de 2 50~ 15 1 [Hombr~ 9 261
## 2 17.6 4 [Localidades menores de 2 50~ 17 2 [Mujer] 9 261
## 3 37.5 4 [Localidades menores de 2 50~ 33 1 [Hombr~ 7 228
## 4 7.41 4 [Localidades menores de 2 50~ 27 1 [Hombr~ 9 308
## 5 25.7 4 [Localidades menores de 2 50~ 57 1 [Hombr~ 12 308
## 6 26.7 4 [Localidades menores de 2 50~ 25 2 [Mujer] 10 306rm(SDEMT319) # podemos quedarnos con la base que estamos usandogg <- ggplot(mydata, aes(anios_esc, log(ing_x_hrs)))
gg + geom_jitter()
mydata$log_ing_x_hrs<-log(mydata$ing_x_hrs)
cor(mydata$log_ing_x_hrs, mydata$anios_esc, use = "pairwise")## [1] 0.3632596Una prueba de hipotésis sobe la correlación (Compara, hoy que quitamos el cero, la diferencia de las correlaciones, en relación a la práctica 4)
cor.test(mydata$log_ing_x_hrs, mydata$anios_esc, use="pairwise.complete.obs")##
## Pearson's product-moment correlation
##
## data: mydata$log_ing_x_hrs and mydata$anios_esc
## t = 135.98, df = 121636, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3583715 0.3681278
## sample estimates:
## cor
## 0.3632596Regresión lineal
Regresión lineal simple
La regresión lineal nos ayuda a describir esta relación a través de una lÃnea recta.
hist(log(mydata$ing_x_hrs))
Una vez transformada nuestra variable, corremos el modelo
cor(mydata$log_ing_x_hrs, mydata$anios_esc,
use = "pairwise")## [1] 0.3632596modelo <-lm(log_ing_x_hrs ~anios_esc, data=mydata,
na.action=na.exclude)
summary(modelo) # show results##
## Call:
## lm(formula = log_ing_x_hrs ~ anios_esc, data = mydata, na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.5290 -0.3710 -0.0075 0.3653 5.0065
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.8117217 0.0050309 558.9 <2e-16 ***
## anios_esc 0.0628322 0.0004621 136.0 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6704 on 121636 degrees of freedom
## Multiple R-squared: 0.132, Adjusted R-squared: 0.132
## F-statistic: 1.849e+04 on 1 and 121636 DF, p-value: < 2.2e-16plot(modelo)
Diagnósticos
library(car)## Warning: package 'car' was built under R version 3.6.2## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:DescTools':
##
## Recode## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## some- Outliers y Normalidad
# Assessing Outliers
outlierTest(modelo) # Bonferonni p-value for most extreme obs## rstudent unadjusted p-value Bonferroni p
## 74500 -8.249260 1.6095e-16 1.9577e-11
## 70355 7.469288 8.1166e-14 9.8728e-09
## 113456 -7.465671 8.3426e-14 1.0148e-08
## 105839 -7.413911 1.2342e-13 1.5013e-08
## 83232 -6.984532 2.8726e-12 3.4942e-07
## 14860 -6.858080 7.0124e-12 8.5297e-07
## 102404 -6.546744 5.9038e-11 7.1812e-06
## 67199 -6.486551 8.8159e-11 1.0723e-05
## 57249 6.474902 9.5233e-11 1.1584e-05
## 74506 -6.346082 2.2164e-10 2.6960e-05out<-outlierTest(modelo) # guardamos en objeto
qqPlot(modelo, main="QQ Plot") #qq plot for studentized resid
## [1] 70355 74500qqPlot<-qqPlot(modelo) #guardamos en objeto
- Homocedasticidad
# non-constant error variance test
ncvTest(modelo)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 210.9901, Df = 1, p = < 2.22e-16# plot studentized residuals vs. fitted values
spreadLevelPlot(modelo)
##
## Suggested power transformation: 0.7286022Volvemos a correr nuestro modelo, hoy con una base que nos quite estas observaciones. Como es nuestro modelo original, le pondremos cero
names(out$bonf.p)## [1] "74500" "70355" "113456" "105839" "83232" "14860" "102404"
## [8] "67199" "57249" "74506"as.integer(names(out$bonf.p))## [1] 74500 70355 113456 105839 83232 14860 102404 67199 57249 74506outliers<-rbind(as.integer(names(out$bonf.p)), qqPlot) # lista los casos Vamos a eliminar estos casos que son extremos (¡Ojo! esto tiene implicaciones de interpretación y debe ser justificado metodológicamente y ser señalado como una limitante)
Tenemos el nombre de las filas que nos dan problemas
mydata$rownames<-rownames(mydata)
#View(mydata) # verificamos que no hayamos movido el orden
mydata2<-mydata[-outliers,]Corremos un nuevo modelo
modelo0<-lm(log_ing_x_hrs ~anios_esc, data=mydata2, na.action=na.exclude)
summary(modelo0)##
## Call:
## lm(formula = log_ing_x_hrs ~ anios_esc, data = mydata2, na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.2014 -0.3712 -0.0078 0.3653 4.0350
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.8122524 0.0050210 560.1 <2e-16 ***
## anios_esc 0.0628026 0.0004611 136.2 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6691 on 121626 degrees of freedom
## Multiple R-squared: 0.1323, Adjusted R-squared: 0.1323
## F-statistic: 1.855e+04 on 1 and 121626 DF, p-value: < 2.2e-16¿Cuando parar?
qqPlot(modelo0)
## [1] 20734 44787outlierTest(modelo0)## rstudent unadjusted p-value Bonferroni p
## 44787 -6.280405 3.3883e-10 4.1211e-05
## 20734 -6.247530 4.1837e-10 5.0885e-05
## 32109 -6.212505 5.2315e-10 6.3629e-05
## 11231 -6.135966 8.4903e-10 1.0327e-04
## 53043 -6.131917 8.7092e-10 1.0593e-04
## 60401 -6.131917 8.7092e-10 1.0593e-04
## 104103 -6.129204 8.8590e-10 1.0775e-04
## 87900 6.031595 1.6282e-09 1.9803e-04
## 97680 5.988775 2.1202e-09 2.5787e-04
## 51092 5.984496 2.1766e-09 2.6474e-04¡Este puede ser un proceso infinito! Si quitamos lo anormal, esto mueve nuestros rangos y al quitar un outlier, otra variable que antes no era outlier en el ajuste se puede convertir en outlier.
Regresión Lineal múltiple
Agregando una variable categórica
Cuando nosotros tenemos una variable categórica para la condición de sexo. [nota: seguimos haciendo el ejercicio, a pesar de que ya observamos en nuestro diagnóstico el modelo no cumple con los supuestos, pero lo haremos para fines ilustrativos]
modelo1<-lm(log_ing_x_hrs ~anios_esc + as_label(sex), data=mydata, na.action=na.exclude)
summary(modelo1)##
## Call:
## lm(formula = log_ing_x_hrs ~ anios_esc + as_label(sex), data = mydata,
## na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.5608 -0.3709 -0.0087 0.3656 4.9746
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.838717 0.005191 546.82 <2e-16 ***
## anios_esc 0.063367 0.000462 137.16 <2e-16 ***
## as_label(sex)Mujer -0.080538 0.003920 -20.55 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6693 on 121635 degrees of freedom
## Multiple R-squared: 0.135, Adjusted R-squared: 0.1349
## F-statistic: 9488 on 2 and 121635 DF, p-value: < 2.2e-16Este modelo tiene coeficientes que deben leerse “condicionados”. Es decir, en este caso tenemos que el coeficiente asociado a “anios_esc”, mantiene constante el valor de sexo y viceversa.
¿Cómo saber is ha mejorado nuestro modelo? Podemos comparar el ajuste con la anova, es decir, una prueba F
pruebaf0<-anova(modelo, modelo1)
pruebaf0## Analysis of Variance Table
##
## Model 1: log_ing_x_hrs ~ anios_esc
## Model 2: log_ing_x_hrs ~ anios_esc + as_label(sex)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 121636 54672
## 2 121635 54483 1 189.07 422.1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como puedes ver, el resultado muestra un Df de 1 (lo que indica que el modelo más complejo tiene un parámetro adicional) y un valor p muy pequeño (<.001). Esto significa que agregar el sexo al modelo lleva a un ajuste significativamente mejor sobre el modelo original.
Podemos seguir añadiendo variables sólo “sumando” en la función
modelo2<-lm(log_ing_x_hrs ~anios_esc + as_label(sex) + eda, data=mydata, na.action=na.exclude)
summary(modelo2)##
## Call:
## lm(formula = log_ing_x_hrs ~ anios_esc + as_label(sex) + eda,
## data = mydata, na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.5374 -0.3519 -0.0024 0.3595 5.0358
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.5361300 0.0084352 300.66 <2e-16 ***
## anios_esc 0.0691894 0.0004759 145.40 <2e-16 ***
## as_label(sex)Mujer -0.0854705 0.0038890 -21.98 <2e-16 ***
## eda 0.0063014 0.0001392 45.28 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6637 on 121634 degrees of freedom
## Multiple R-squared: 0.1493, Adjusted R-squared: 0.1493
## F-statistic: 7116 on 3 and 121634 DF, p-value: < 2.2e-16Y podemos ver si introducir esta variable afectó al ajuste global del modelo
pruebaf1<-anova(modelo1, modelo2)
pruebaf1## Analysis of Variance Table
##
## Model 1: log_ing_x_hrs ~ anios_esc + as_label(sex)
## Model 2: log_ing_x_hrs ~ anios_esc + as_label(sex) + eda
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 121635 54483
## 2 121634 53579 1 903.29 2050.6 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Hoy que tenemos más variables podemos hablar de revisar dos supuestos más.
Otros supuestos
Además de los supuestos de la regresión simple, podemos revisar estos otros. De nuevo, usaremos la librerÃa “car”,
Linealidad en los parámetros (será más dÃficil entre más variables tengamos)
La normalidad también, porque debe ser multivariada
Multicolinealidad La prueba más común es la de Factor Influyente de la Varianza (VIF) por sus siglas en inglés. La lógica es que la multicolinealidad tendrá efectos en nuestro R2, inflándolo. De ahà que observamos de qué variable(s) proviene este problema relacionado con la multicolinealidad.
Si el valor es mayor a 5, tenemos un problema muy grave.
library(car)
vif(modelo2)## anios_esc as_label(sex) eda
## 1.082193 1.003972 1.078908Tabla de modelos estimados
Para los muy avanzados, con el paquete “stargazer” se pueden pasar a LaTeX fácilmente.
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=stargazer#stargazer(modelo0, modelo1,modelo2, type = 'latex', header=FALSE)stargazer(modelo, modelo1,modelo2, type = 'text', header=FALSE)##
## ==============================================================================================================
## Dependent variable:
## ------------------------------------------------------------------------------------------
## log_ing_x_hrs
## (1) (2) (3)
## --------------------------------------------------------------------------------------------------------------
## anios_esc 0.063*** 0.063*** 0.069***
## (0.0005) (0.0005) (0.0005)
##
## as_label(sex)Mujer -0.081*** -0.085***
## (0.004) (0.004)
##
## eda 0.006***
## (0.0001)
##
## Constant 2.812*** 2.839*** 2.536***
## (0.005) (0.005) (0.008)
##
## --------------------------------------------------------------------------------------------------------------
## Observations 121,638 121,638 121,638
## R2 0.132 0.135 0.149
## Adjusted R2 0.132 0.135 0.149
## Residual Std. Error 0.670 (df = 121636) 0.669 (df = 121635) 0.664 (df = 121634)
## F Statistic 18,490.790*** (df = 1; 121636) 9,488.450*** (df = 2; 121635) 7,115.761*** (df = 3; 121634)
## ==============================================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01También la librerÃa “sjPlot” tiene el comando “plot_model()” (instala el comando si no lo tienes)
library(sjPlot)## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod car## Learn more about sjPlot with 'browseVignettes("sjPlot")'.plot_model(modelo1)
Estandarizando
Comparar los resultados de los coeficientes es dÃficil, porque el efecto está medido en las unidades que fueron medidas. Por lo que no serÃa tan comparable el efecto que tenemos de nuestro Ãndice sumativo (proporción de lugares con inseguridad declarada) con respecto a la eda (que se mide en años). Por lo que a veces es mejor usar las medida estandarizadas (es decir, nuestra puntajes z).
Podemos hacerlo transormando nuestras variables de origen e introducirlas al modelo. O bien, podemos usar un paquete que lo hace directamente. Los coeficientes calculados se les concoe como “beta”
library("lm.beta")Simplemente aplicamos el comando a nuestros modelos ya calculados
lm.beta(modelo2)##
## Call:
## lm(formula = log_ing_x_hrs ~ anios_esc + as_label(sex) + eda,
## data = mydata, na.action = na.exclude)
##
## Standardized Coefficients::
## (Intercept) anios_esc as_label(sex)Mujer
## 0.00000000 0.40001328 -0.05823751
## eda
## 0.12439260Hoy la comparación será mucho más clara y podemos ver qué variable tiene mayor efecto en nuestra dependiente.
modelo_beta<-lm.beta(modelo2)
modelo_beta##
## Call:
## lm(formula = log_ing_x_hrs ~ anios_esc + as_label(sex) + eda,
## data = mydata, na.action = na.exclude)
##
## Standardized Coefficients::
## (Intercept) anios_esc as_label(sex)Mujer
## 0.00000000 0.40001328 -0.05823751
## eda
## 0.12439260Para graficarlos, podemos usar de nuevo el comando plot_model(), con una opción
plot_model(modelo2, type="std")
¿Qué podemos concluir de estos resultados?
No cumplo los supuestos
Heterocedasticidad
El problema de la heterocedasticidad es que los errores estándar de subestiman, por lo que si estos están en el cociente de nuestro estadÃstico de prueba t, esto implicarÃa que nuestras pruebas podrÃan estar arrojando valores significativos cuando no lo son.
Una forma muy sencilla es pedir los errores robustos, importando una función https://economictheoryblog.com/2016/08/08/robust-standard-errors-in-r/
library(RCurl)## Warning: package 'RCurl' was built under R version 3.6.2## Loading required package: bitops##
## Attaching package: 'RCurl'## The following object is masked from 'package:tidyr':
##
## completeurl_robust <- "https://raw.githubusercontent.com/IsidoreBeautrelet/economictheoryblog/master/robust_summary.R"
eval(parse(text = getURL(url_robust, ssl.verifypeer = FALSE)),
envir=.GlobalEnv)
summary(modelo2, robust=T)##
## Call:
## lm(formula = log_ing_x_hrs ~ anios_esc + as_label(sex) + eda,
## data = mydata, na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.5374 -0.3519 -0.0024 0.3595 5.0358
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.5361300 0.0086913 291.80 <2e-16 ***
## anios_esc 0.0691894 0.0005284 130.94 <2e-16 ***
## as_label(sex)Mujer -0.0854705 0.0039218 -21.79 <2e-16 ***
## eda 0.0063014 0.0001523 41.38 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6637 on 121634 degrees of freedom
## Multiple R-squared: 0.1493, Adjusted R-squared: 0.1493
## F-statistic: 5843 on 3 and 121634 DF, p-value: < 2.2e-16summary(modelo2, robust=F)##
## Call:
## lm(formula = log_ing_x_hrs ~ anios_esc + as_label(sex) + eda,
## data = mydata, na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.5374 -0.3519 -0.0024 0.3595 5.0358
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.5361300 0.0084352 300.66 <2e-16 ***
## anios_esc 0.0691894 0.0004759 145.40 <2e-16 ***
## as_label(sex)Mujer -0.0854705 0.0038890 -21.98 <2e-16 ***
## eda 0.0063014 0.0001392 45.28 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6637 on 121634 degrees of freedom
## Multiple R-squared: 0.1493, Adjusted R-squared: 0.1493
## F-statistic: 7116 on 3 and 121634 DF, p-value: < 2.2e-16Matriz varianza-covarianza
library(robustbase)## Warning: package 'robustbase' was built under R version 3.6.2modelo2rob<-lmrob(log_ing_x_hrs ~anios_esc + sex + eda, data=mydata, na.action=na.exclude)
summary(modelo2rob)##
## Call:
## lmrob(formula = log_ing_x_hrs ~ anios_esc + sex + eda, data = mydata, na.action = na.exclude)
## \--> method = "MM"
## Residuals:
## Min 1Q Median 3Q Max
## -5.548313 -0.355825 -0.008202 0.351402 5.032049
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.6393745 0.0086881 303.79 <2e-16 ***
## anios_esc 0.0643320 0.0004836 133.02 <2e-16 ***
## sex -0.0914825 0.0034217 -26.74 <2e-16 ***
## eda 0.0074911 0.0001329 56.39 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 0.5299
## Multiple R-squared: 0.179, Adjusted R-squared: 0.179
## Convergence in 11 IRWLS iterations
##
## Robustness weights:
## 755 observations c(133,419,1382,2102,2834,3121,3338,3521,3544,3563,4048,4119,4497,5156,5193,6275,6367,6407,7219,7287,7401,7730,7843,8119,8852,9081,9171,9218,9413,9780,9835,10698,10801,11214,11223,11231,11242,11245,11762,11768,12085,12118,12384,12696,12880,13164,13169,14369,14438,14582,14757,14860,15088,15094,15546,15567,15748,15764,16371,16442,17211,17409,17502,17536,17882,18179,18194,18309,18576,18618,18807,18840,18850,19131,19190,19460,19869,19871,20143,20146,20177,20408,20435,20453,20455,20471,20667,20708,20735,20866,20870,20873,20879,20883,20895,20906,20922,20924,20926,20946,21545,21564,21573,21575,21598,21608,21609,21640,21644,21649,21658,21688,21697,21809,21927,21995,22004,22029,22107,22279,22480,22541,22579,22592,22600,22650,22666,22694,22823,22943,23003,23267,23407,23698,23968,24254,24409,24694,24810,25212,26212,26231,26583,26683,27464,27540,27798,27869,27897,28075,28141,28361,28375,28615,28670,28756,28918,28991,29331,29422,29436,29462,29528,29531,29541,29547,29678,29685,29701,29827,29877,29885,29890,29892,29902,29909,29916,30156,30228,30297,30345,30725,30983,31137,31292,31482,31881,32036,32110,32326,33046,33179,33407,33568,33824,33850,33976,34450,34796,35445,35683,36223,36569,36711,36712,36934,36947,36964,37007,37025,37059,37137,37180,37194,37208,37220,37237,37406,37432,37446,37453,37458,38880,38971,39196,39323,39668,40304,40458,41093,41116,41342,41549,41740,42115,42127,42542,42674,43007,43108,43283,43423,43945,44274,44328,44341,44444,44524,44577,44589,44671,44719,44724,44788,44804,44806,44838,44839,44965,45041,45046,45048,45139,45153,45273,45295,45363,45364,45370,46360,46564,46659,46671,46678,46830,46864,48245,48344,48586,48588,49195,49286,49322,49364,50011,50147,50159,50306,50311,50389,50477,50816,50858,51093,51614,51780,51801,51814,51941,51951,51994,52083,52297,52320,52397,52406,52410,52412,52544,52560,52590,52669,52759,52761,52769,52797,52801,52899,52925,53044,53049,53054,53245,53369,53652,54252,54583,54893,55497,56108,56151,56294,56546,56982,57169,57209,57249,57432,57539,57743,58161,58434,58635,59049,59527,59552,59766,59888,59889,59891,59902,59988,60152,60167,60168,60171,60172,60177,60205,60352,60386,60403,60408,60409,60421,60782,60804,60867,62155,62161,62163,62644,62939,63260,63607,63798,64432,64801,65045,65102,65130,65142,65443,66298,66482,66761,66867,66974,67199,67255,67277,67280,67450,67587,67600,67672,67714,67720,67758,67762,67780,67967,67987,67993,67998,68001,68007,68070,68074,68115,68188,68299,68332,68343,68795,69691,69804,69827,69831,70280,70337,70355,70830,71106,71189,71327,71343,71473,71949,72088,72294,72562,72844,72928,73363,73800,73868,74148,74400,74500,74506,74528,74704,74957,75146,75160,75173,75206,75225,75306,75313,75322,75333,75334,75335,75336,75340,75363,75370,75418,75575,75577,75583,75597,75602,75893,75896,75977,76430,76469,77011,77408,78150,78156,78605,78867,79303,79607,79985,80303,80387,80726,80943,81103,81232,81449,81588,82026,82158,82309,82498,82522,82523,82547,82615,82716,82812,82873,82881,82933,82936,82942,82945,82952,82971,83004,83045,83200,83210,83219,83232,83316,83350,83568,83653,84017,84121,84136,84637,84836,85183,85714,85745,86090,86124,86828,86936,86964,86990,87000,87544,87563,87763,87907,88008,88197,88292,88337,88970,89622,89778,89928,89997,90002,90214,90329,90347,90424,90431,90432,90475,90506,90560,90579,90588,90631,90766,90851,90858,91061,91096,91159,91260,92187,92352,92517,92611,92627,92740,93613,94140,94719,95275,95455,95678,96005,96168,96194,96404,96783,97095,97133,97181,97275,97455,97524,97687,97700,97802,98009,98141,98161,98212,98398,98410,98411,98421,98428,98774,98807,98812,98966,99339,99622,99756,99872,100308,100559,100587,101146,101653,101749,101774,102251,102394,102404,102630,103568,103681,103813,104111,104143,104159,104557,104978,105069,105070,105134,105402,105405,105406,105430,105481,105508,105515,105518,105640,105650,105808,105839,105874,105884,105889,105903,105907,105909,105919,105940,106034,106142,106164,106276,106448,106491,106529,106853,106963,107713,107814,107998,108578,108629,108893,109332,109813,110051,110151,110726,111264,111279,111600,112614,112742,112748,112848,112937,112941,112942,113188,113222,113233,113259,113383,113406,113446,113456,113458,113459,113460,113463,113468,113473,113476,113477,113488,113490,113560,113653,113841,113851,113865,114203,114343,114371,114373,114414,114421,114433,114748,114978,115154,115176,115298,116296,116620,117146,117343,117706,117896,117913,117984,118081,118246,118425,118600,118761,118849,118852,118948,118990,119038,119207,119419,120042,120481,120699,120719,120733,120744,120862,120889,120902,120904,120952,120976,121189,121202,121210,121362,121453,121502,121531,121543,121567,121587)
## are outliers with |weight| <= 6.4e-07 ( < 8.2e-07);
## 10664 weights are ~= 1. The remaining 110219 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000088 0.8578000 0.9511000 0.8814000 0.9862000 0.9990000
## Algorithmic parameters:
## tuning.chi bb tuning.psi refine.tol
## 1.548e+00 5.000e-01 4.685e+00 1.000e-07
## rel.tol scale.tol solve.tol eps.outlier
## 1.000e-07 1.000e-10 1.000e-07 8.221e-07
## eps.x warn.limit.reject warn.limit.meanrw
## 1.764e-10 5.000e-01 5.000e-01
## nResample max.it best.r.s k.fast.s k.max
## 500 50 2 1 200
## maxit.scale trace.lev mts compute.rd fast.s.large.n
## 200 0 1000 0 2000
## psi subsampling cov
## "bisquare" "nonsingular" ".vcov.avar1"
## compute.outlier.stats
## "SM"
## seed : int(0)Regresión robusta
No es lo mismo la regresión robusta que los errores robustos. La regresión robusta es más robusta a los outliers. No confundir.
La regresión robusta, es esto, robusta a los outliers, porque pesa el valor de las observaciones de tal manera que los outliers tenga menor influencia.
stargazer(modelo2, modelo2rob, type = 'text', header=FALSE)##
## =========================================================================
## Dependent variable:
## ---------------------------------------
## log_ing_x_hrs
## OLS MM-type
## linear
## (1) (2)
## -------------------------------------------------------------------------
## anios_esc 0.069*** 0.064***
## (0.0005) (0.0005)
##
## as_label(sex)Mujer -0.085***
## (0.004)
##
## sex -0.091***
## (0.003)
##
## eda 0.006*** 0.007***
## (0.0001) (0.0001)
##
## Constant 2.536*** 2.639***
## (0.008) (0.009)
##
## -------------------------------------------------------------------------
## Observations 121,638 121,638
## R2 0.149 0.179
## Adjusted R2 0.149 0.179
## Residual Std. Error (df = 121634) 0.664 0.530
## F Statistic 7,115.761*** (df = 3; 121634)
## =========================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Un poquito de reflexión
Se pueden usar métodos no paramétricos, como la regresión mediana. O como ya vimos podemos transformar la variable a logaritmo, seleccionar casos.
En nuestro caso quitamos los ingresos iguales a cero. No obstante, esto resuelve un problema estadÃstico, pero nos supone una decisión metodológica dÃficil de atender.
Por ello es recomendable mejor utilizar otro tipo de modelos más robustos a la presencia de outliers (i.e. regresión robusta) y menos dependientes de una distribución normal (i.e. regresión mediana).