Previo

Siempre: revisar el directorio. Cargamos nuestro ambiente donde está la ehpm 2017

setwd("/Users/anaescoto/Dropbox/2019/Cursos ESA/UCA_R")
load("./datos/ehpm2017.RData")

Librerías

#install.packages(c("estout", "car", "robustbase","RCurl"), dependencies = TRUE)

#install.packages("lm.beta", dependencies = TRUE)

library(tidyverse) #porque usaremos dplyr y por cualquier cosa
## ── Attaching packages ─────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.2.0     ✔ purrr   0.3.2
## ✔ tibble  2.1.3     ✔ dplyr   0.8.1
## ✔ tidyr   0.8.3     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0
## ── Conflicts ────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ 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
library(DescTools)

Sub-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<- ehpm2017[ehpm2017$actpr2012==10, c("money", "region","r106", "r104", "aproba1", "fac00")]
tail(mydata)
## # A tibble: 6 x 6
##   money        region  r106        r104 aproba1 fac00
##   <dbl>     <dbl+lbl> <dbl>   <dbl+lbl>   <dbl> <dbl>
## 1    12  4 [Oriental]    84  2 [Mujer]        0    71
## 2    NA NA               NA NA               NA    NA
## 3     0  4 [Oriental]    50  2 [Mujer]        3    71
## 4     0  4 [Oriental]    75  1 [Hombre]       0    71
## 5    52  4 [Oriental]    24  1 [Hombre]       7    71
## 6     3  4 [Oriental]    29  2 [Mujer]        2    71
gg <- ggplot(mydata, aes(aproba1, log(money+1)))
gg +  geom_jitter()
## Warning: Removed 6286 rows containing missing values (geom_point).

mydata$log_money<-log(mydata$money+1)
cor(mydata$log_money, mydata$aproba1,  use = "pairwise")
## [1] 0.3140178

Por default está la correlación de Pearson, pero en realidad podemos obtener otros tipos

#Pearson - default
cor(mydata$log_money, mydata$aproba1, 
    use = "pairwise", method = "pearson")
## [1] 0.3140178
#Tau-Kendall
cor(mydata$log_money, mydata$aproba1, 
    use = "pairwise", method = "kendall")
## [1] 0.3229262
#Rho-Spearman
cor(mydata$log_money, mydata$aproba1, 
    use = "pairwise", method = "spearman")
## [1] 0.4395571

Una prueba de hipotésis

cor.test(mydata$log_money, mydata$aproba1, use="pairwise.complete.obs")
## 
##  Pearson's product-moment correlation
## 
## data:  mydata$log_money and mydata$aproba1
## t = 58.996, df = 31816, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3040791 0.3238881
## sample estimates:
##       cor 
## 0.3140178

Regresió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$money+1))

Una vez transformada nuestra variable, corremos el modelo

cor(mydata$log_money, mydata$aproba1, 
    use = "pairwise")
## [1] 0.3140178
modelo <-lm(log_money ~aproba1, data=mydata, 
            na.action=na.exclude)

summary(modelo) # show results
## 
## Call:
## lm(formula = log_money ~ aproba1, data = mydata, na.action = na.exclude)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.3611 -0.1050  0.7217  1.2636  5.0401 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.36309    0.02165   155.3   <2e-16 ***
## aproba1      0.14276    0.00242    59.0   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.074 on 31816 degrees of freedom
##   (6286 observations deleted due to missingness)
## Multiple R-squared:  0.09861,    Adjusted R-squared:  0.09858 
## F-statistic:  3480 on 1 and 31816 DF,  p-value: < 2.2e-16
plot(modelo)

Diagnósticos

library(car)
## 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
  1. Outliers y Normalidad
# Assessing Outliers
outlierTest(modelo) # Bonferonni p-value for most extreme obs
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##        rstudent unadjusted p-value Bonferroni p
## 13271 -3.067833           0.002158           NA
qqPlot(modelo, main="QQ Plot") #qq plot for studentized resid 

## [1]  9970 13271
  1. Homocedasticidad
# non-constant error variance test
ncvTest(modelo)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 947.6233, Df = 1, p = < 2.22e-16
# plot studentized residuals vs. fitted values 
spreadLevelPlot(modelo)

## 
## Suggested power transformation:  3.517753

Volvemos a correr nuestro modelo, hoy con esta base. Como es nuestro modelo original, le pondremos cero

modelo0<-lm(log_money ~aproba1, data=mydata, na.action=na.exclude)
summary(modelo0)
## 
## Call:
## lm(formula = log_money ~ aproba1, data = mydata, na.action = na.exclude)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.3611 -0.1050  0.7217  1.2636  5.0401 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.36309    0.02165   155.3   <2e-16 ***
## aproba1      0.14276    0.00242    59.0   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.074 on 31816 degrees of freedom
##   (6286 observations deleted due to missingness)
## Multiple R-squared:  0.09861,    Adjusted R-squared:  0.09858 
## F-statistic:  3480 on 1 and 31816 DF,  p-value: < 2.2e-16

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_money ~aproba1 + as_label(r104), data=mydata, na.action=na.exclude)
summary(modelo1)
## 
## Call:
## lm(formula = log_money ~ aproba1 + as_label(r104), data = mydata, 
##     na.action = na.exclude)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2221 -0.1820  0.7306  1.2756  5.1674 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         3.259322   0.023182  140.60   <2e-16 ***
## aproba1             0.141083   0.002418   58.35   <2e-16 ***
## as_label(r104)Mujer 0.292429   0.023734   12.32   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.069 on 31815 degrees of freedom
##   (6286 observations deleted due to missingness)
## Multiple R-squared:  0.1029, Adjusted R-squared:  0.1028 
## F-statistic:  1824 on 2 and 31815 DF,  p-value: < 2.2e-16

Este modelo tiene coeficientes que deben leerse “condicionados”. Es decir, en este caso tenemos que el coeficiente asociado a “aproba1”, 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(modelo0, modelo1)
pruebaf0
## Analysis of Variance Table
## 
## Model 1: log_money ~ aproba1
## Model 2: log_money ~ aproba1 + as_label(r104)
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1  31816 136863                                  
## 2  31815 136213  1    649.96 151.81 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Como 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_money ~aproba1 + as_label(r104) + r106, data=mydata, na.action=na.exclude)
summary(modelo2)
## 
## Call:
## lm(formula = log_money ~ aproba1 + as_label(r104) + r106, data = mydata, 
##     na.action = na.exclude)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.7412 -0.1955  0.7071  1.2749  5.1352 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         2.3740924  0.0416089   57.06   <2e-16 ***
## aproba1             0.1597568  0.0025032   63.82   <2e-16 ***
## as_label(r104)Mujer 0.2627884  0.0235241   11.17   <2e-16 ***
## r106                0.0198769  0.0007793   25.50   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.048 on 31814 degrees of freedom
##   (6286 observations deleted due to missingness)
## Multiple R-squared:  0.1209, Adjusted R-squared:  0.1208 
## F-statistic:  1458 on 3 and 31814 DF,  p-value: < 2.2e-16

Y podemos ver si introducir esta variable afectó al ajuste global del modelo

pruebaf1<-anova(modelo1, modelo2)
pruebaf1
## Analysis of Variance Table
## 
## Model 1: log_money ~ aproba1 + as_label(r104)
## Model 2: log_money ~ aproba1 + as_label(r104) + r106
##   Res.Df    RSS Df Sum of Sq     F    Pr(>F)    
## 1  31815 136213                                 
## 2  31814 133483  1    2729.3 650.5 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hoy 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”,

  1. Linealidad en los parámetros (será más díficil entre más variables tengamos)

  2. La normalidad también, porque debe ser multivariada

  3. 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)
##        aproba1 as_label(r104)           r106 
##       1.097044       1.005645       1.094596

Tabla 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(modelo0, modelo1,modelo2, type = 'text', header=FALSE)
## 
## ==========================================================================================================
##                                                      Dependent variable:                                  
##                     --------------------------------------------------------------------------------------
##                                                           log_money                                       
##                                 (1)                          (2)                          (3)             
## ----------------------------------------------------------------------------------------------------------
## aproba1                       0.143***                     0.141***                     0.160***          
##                               (0.002)                      (0.002)                      (0.003)           
##                                                                                                           
## as_label(r104)Mujer                                        0.292***                     0.263***          
##                                                            (0.024)                      (0.024)           
##                                                                                                           
## r106                                                                                    0.020***          
##                                                                                         (0.001)           
##                                                                                                           
## Constant                      3.363***                     3.259***                     2.374***          
##                               (0.022)                      (0.023)                      (0.042)           
##                                                                                                           
## ----------------------------------------------------------------------------------------------------------
## Observations                   31,818                       31,818                       31,818           
## R2                             0.099                        0.103                        0.121            
## Adjusted R2                    0.099                        0.103                        0.121            
## Residual Std. Error      2.074 (df = 31816)           2.069 (df = 31815)           2.048 (df = 31814)     
## F Statistic         3,480.487*** (df = 1; 31816) 1,824.397*** (df = 2; 31815) 1,457.927*** (df = 3; 31814)
## ==========================================================================================================
## Note:                                                                          *p<0.1; **p<0.05; ***p<0.01

Tambié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
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 r106 (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_money ~ aproba1 + as_label(r104) + r106, data = mydata, 
##     na.action = na.exclude)
## 
## Standardized Coefficients::
##         (Intercept)             aproba1 as_label(r104)Mujer 
##          0.00000000          0.35139589          0.05888904 
##                r106 
##          0.14027123

Hoy 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_money ~ aproba1 + as_label(r104) + r106, data = mydata, 
##     na.action = na.exclude)
## 
## Standardized Coefficients::
##         (Intercept)             aproba1 as_label(r104)Mujer 
##          0.00000000          0.35139589          0.05888904 
##                r106 
##          0.14027123

Para 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)
## Loading required package: bitops
## 
## Attaching package: 'RCurl'
## The following object is masked from 'package:tidyr':
## 
##     complete
url_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_money ~ aproba1 + as_label(r104) + r106, data = mydata, 
##     na.action = na.exclude)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.7412 -0.1955  0.7071  1.2749  5.1352 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)
## (Intercept)          2.37409         NA      NA       NA
## aproba1              0.15976         NA      NA       NA
## as_label(r104)Mujer  0.26279         NA      NA       NA
## r106                 0.01988         NA      NA       NA
## 
## Residual standard error: 2.048 on 31814 degrees of freedom
##   (6286 observations deleted due to missingness)
## Multiple R-squared:  0.1209, Adjusted R-squared:  0.1208 
## F-statistic:    NA on 3 and 31814 DF,  p-value: NA
summary(modelo2, robust=F)
## 
## Call:
## lm(formula = log_money ~ aproba1 + as_label(r104) + r106, data = mydata, 
##     na.action = na.exclude)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.7412 -0.1955  0.7071  1.2749  5.1352 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         2.3740924  0.0416089   57.06   <2e-16 ***
## aproba1             0.1597568  0.0025032   63.82   <2e-16 ***
## as_label(r104)Mujer 0.2627884  0.0235241   11.17   <2e-16 ***
## r106                0.0198769  0.0007793   25.50   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.048 on 31814 degrees of freedom
##   (6286 observations deleted due to missingness)
## Multiple R-squared:  0.1209, Adjusted R-squared:  0.1208 
## F-statistic:  1458 on 3 and 31814 DF,  p-value: < 2.2e-16

Matriz varianza-covarianza

library(robustbase)
modelo2rob<-lmrob(log_money ~aproba1 + r104 + r106, data=mydata, na.action=na.exclude)
summary(modelo2rob)
## 
## Call:
## lmrob(formula = log_money ~ aproba1 + r104 + r106, data = mydata, na.action = na.exclude)
##  \--> method = "MM"
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.7248 -0.9249 -0.1086  0.2957  4.4330 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.7928560  0.0178866  267.96   <2e-16 ***
## aproba1      0.0778324  0.0008162   95.36   <2e-16 ***
## r104        -0.2695597  0.0077621  -34.73   <2e-16 ***
## r106         0.0107857  0.0003060   35.25   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 0.7188 
##   (6286 observations deleted due to missingness)
## Multiple R-squared:  0.279,  Adjusted R-squared:  0.2789 
## Convergence in 9 IRWLS iterations
## 
## Robustness weights: 
##  5704 observations c(2,3,11,12,20,21,42,43,44,70,83,93,99,102,103,113,114,130,140,152,161,167,180,182,184,185,186,196,197,207,209,213,216,229,232,233,238,254,263,286,298,313,315,316,318,319,323,325,326,327,328,332,346,349,355,374,379,381,383,385,386,390,394,406,416,418,421,422,428,429,447,449,457,467,477,480,497,519,542,551,557,561,565,574,585,586,587,638,639,649,661,664,665,677,691,692,707,708,712,721,723,745,762,767,775,793,794,797,799,802,823,828,852,867,868,875,885,897,898,899,906,907,909,913,919,922,927,940,942,948,967,976,980,985,1011,1015,1031,1040,1042,1057,1102,1104,1107,1110,1112,1117,1120,1130,1131,1132,1133,1134,1135,1136,1137,1141,1143,1154,1173,1178,1179,1206,1207,1208,1209,1215,1216,1217,1238,1242,1243,1257,1262,1266,1275,1288,1290,1291,1292,1296,1302,1303,1307,1308,1310,1313,1314,1315,1316,1317,1318,1325,1328,1330,1343,1344,1346,1350,1358,1361,1363,1369,1370,1373,1374,1376,1377,1379,1385,1387,1388,1389,1393,1395,1397,1400,1405,1408,1411,1412,1413,1418,1419,1420,1421,1423,1424,1425,1426,1432,1433,1464,1466,1475,1476,1477,1498,1500,1501,1507,1525,1526,1531,1549,1556,1562,1563,1571,1580,1601,1602,1604,1607,1608,1611,1621,1651,1652,1661,1663,1665,1666,1667,1669,1670,1672,1673,1674,1675,1676,1677,1678,1679,1690,1697,1702,1708,1710,1715,1719,1724,1728,1732,1763,1764,1766,1770,1771,1772,1774,1776,1778,1782,1785,1794,1795,1797,1799,1800,1801,1803,1805,1806,1807,1809,1815,1820,1825,1834,1836,1841,1844,1846,1847,1851,1852,1853,1854,1855,1856,1857,1858,1859,1862,1864,1865,1866,1868,1872,1873,1874,1875,1876,1879,1881,1882,1892,1894,1895,1896,1897,1906,1910,1920,1921,1934,1936,1939,1959,1995,1997,2014,2025,2028,2030,2041,2043,2051,2069,2074,2077,2082,2083,2105,2110,2112,2114,2115,2117,2124,2126,2131,2136,2147,2148,2160,2162,2174,2176,2177,2178,2182,2187,2191,2196,2197,2199,2201,2202,2208,2209,2210,2211,2216,2221,2222,2223,2224,2227,2230,2231,2233,2237,2239,2240,2241,2244,2245,2246,2248,2255,2256,2260,2262,2266,2267,2274,2276,2278,2286,2318,2346,2347,2355,2367,2377,2391,2393,2394,2396,2397,2398,2400,2401,2402,2411,2412,2414,2415,2416,2419,2421,2426,2427,2429,2431,2432,2433,2439,2440,2443,2445,2446,2448,2450,2454,2456,2457,2464,2465,2466,2469,2470,2472,2473,2475,2478,2480,2481,2482,2483,2484,2485,2487,2489,2490,2493,2501,2503,2507,2508,2509,2519,2521,2522,2523,2524,2525,2526,2528,2529,2530,2531,2534,2536,2539,2540,2541,2542,2545,2551,2554,2564,2569,2571,2573,2574,2575,2576,2625,2637,2646,2647,2648,2650,2651,2652,2653,2654,2656,2657,2665,2666,2668,2670,2671,2675,2682,2713,2722,2739,2748,2750,2751,2755,2756,2792,2802,2814,2833,2845,2880,2882,2883,2894,2910,2945,2985,2986,2987,2990,2993,2996,2997,2999,3003,3005,3016,3020,3021,3024,3025,3032,3033,3034,3035,3037,3039,3041,3042,3043,3045,3047,3048,3049,3051,3052,3053,3054,3058,3060,3062,3063,3067,3082,3087,3102,3107,3108,3110,3119,3120,3126,3127,3153,3179,3180,3183,3185,3186,3188,3189,3200,3201,3222,3233,3234,3236,3264,3267,3276,3284,3286,3290,3306,3307,3308,3309,3312,3313,3314,3315,3316,3317,3320,3322,3324,3325,3329,3332,3334,3335,3342,3343,3344,3345,3364,3387,3392,3393,3442,3446,3450,3454,3455,3456,3466,3506,3516,3527,3528,3537,3541,3549,3559,3562,3563,3565,3569,3571,3572,3573,3574,3575,3576,3577,3579,3581,3589,3591,3593,3602,3603,3608,3609,3610,3611,3613,3618,3619,3620,3621,3623,3624,3626,3627,3630,3633,3634,3635,3640,3641,3642,3649,3654,3673,3674,3679,3681,3684,3685,3686,3688,3700,3701,3702,3703,3704,3706,3707,3709,3710,3711,3726,3732,3737,3739,3740,3741,3746,3748,3752,3753,3755,3756,3757,3758,3759,3762,3763,3764,3765,3770,3772,3773,3774,3775,3776,3778,3780,3788,3789,3791,3800,3801,3805,3806,3821,3822,3823,3826,3827,3828,3829,3830,3831,3832,3833,3834,3835,3836,3837,3838,3839,3840,3842,3843,3844,3845,3846,3849,3850,3851,3853,3855,3856,3857,3858,3859,3865,3866,3867,3868,3870,3872,3873,3875,3876,3877,3879,3881,3884,3885,3886,3887,3888,3890,3892,3893,3898,3899,3905,3906,3907,3908,3917,3934,3936,3937,3938,3939,3941,3946,3947,3948,3954,3964,3966,3967,3968,3971,3973,3975,3976,3977,3978,3979,3980,3981,3983,3984,3986,3987,3988,3989,3992,3995,3996,3997,4000,4001,4002,4004,4005,4006,4007,4008,4009,4011,4016,4020,4021,4026,4030,4032,4033,4034,4036,4043,4044,4045,4046,4047,4052,4053,4075,4077,4079,4081,4085,4088,4099,4100,4102,4104,4110,4112,4122,4124,4130,4135,4136,4138,4140,4143,4144,4145,4146,4148,4150,4151,4173,4210,4211,4213,4217,4227,4234,4236,4240,4287,4294,4337,4340,4342,4349,4356,4359,4424,4433,4434,4435,4436,4446,4466,4477,4483,4495,4514,4521,4523,4524,4529,4542,4547,4552,4563,4565,4576,4578,4584,4585,4586,4593,4594,4596,4604,4605,4607,4608,4609,4610,4611,4612,4615,4616,4617,4618,4619,4620,4621,4623,4624,4625,4630,4631,4632,4633,4638,4640,4641,4656,4657,4660,4666,4668,4669,4671,4676,4678,4689,4693,4696,4698,4699,4700,4702,4703,4704,4705,4708,4710,4717,4718,4719,4720,4722,4725,4726,4730,4732,4733,4734,4735,4736,4737,4742,4744,4745,4747,4748,4749,4751,4754,4756,4757,4758,4760,4764,4765,4766,4769,4770,4773,4775,4790,4797,4800,4807,4808,4809,4810,4811,4813,4814,4816,4817,4818,4819,4820,4821,4822,4823,4824,4827,4828,4829,4831,4832,4834,4838,4842,4843,4844,4845,4847,4853,4856,4859,4860,4862,4863,4864,4865,4867,4868,4869,4872,4873,4874,4875,4876,4877,4878,4880,4882,4886,4887,4889,4919,4924,4932,4952,4956,4977,4999,5000,5006,5007,5012,5074,5076,5077,5087,5093,5094,5096,5099,5100,5102,5103,5106,5107,5109,5110,5111,5112,5115,5123,5129,5158,5159,5160,5161,5164,5172,5175,5204,5205,5210,5231,5232,5233,5235,5236,5237,5238,5244,5245,5270,5275,5282,5284,5285,5286,5287,5303,5342,5356,5358,5402,5405,5412,5414,5416,5421,5426,5428,5434,5436,5450,5454,5474,5481,5483,5492,5493,5494,5498,5518,5523,5524,5525,5526,5528,5530,5534,5535,5536,5542,5565,5568,5570,5571,5572,5581,5584,5587,5593,5596,5615,5633,5636,5637,5640,5647,5657,5660,5678,5692,5695,5713,5716,5717,5719,5722,5724,5725,5726,5733,5736,5737,5738,5739,5741,5742,5743,5744,5748,5749,5750,5755,5757,5763,5772,5774,5775,5778,5787,5815,5837,5844,5845,5855,5859,5860,5864,5868,5870,5872,5874,5875,5880,5881,5882,5883,5884,5885,5886,5888,5894,5896,5903,5904,5947,5951,5952,5953,5967,5975,5979,5994,5995,5997,6005,6013,6043,6049,6050,6063,6067,6070,6081,6132,6133,6162,6165,6166,6179,6180,6181,6182,6183,6184,6190,6192,6194,6196,6197,6199,6201,6202,6203,6204,6205,6206,6207,6208,6209,6226,6227,6228,6233,6273,6285,6291,6292,6327,6331,6344,6356,6359,6362,6363,6364,6365,6366,6376,6379,6401,6407,6435,6477,6479,6481,6484,6490,6501,6514,6543,6545,6564,6570,6586,6604,6609,6615,6622,6633,6667,6670,6698,6708,6713,6743,6747,6749,6760,6785,6787,6810,6816,6828,6838,6848,6849,6850,6854,6900,6901,6923,6925,6927,6930,6931,6932,6945,6946,6947,6948,6953,6954,6955,6957,6958,6961,6962,6963,6964,6965,6966,6967,6969,6986,6994,7014,7035,7048,7120,7135,7201,7216,7217,7223,7237,7248,7258,7266,7273,7274,7277,7288,7295,7298,7305,7308,7310,7316,7317,7322,7323,7327,7344,7346,7352,7353,7357,7362,7380,7383,7393,7421,7482,7489,7501,7508,7511,7512,7517,7518,7525,7537,7538,7542,7544,7546,7547,7548,7549,7551,7553,7554,7556,7559,7560,7561,7562,7570,7573,7581,7586,7617,7620,7621,7650,7651,7654,7655,7658,7664,7665,7668,7683,7684,7685,7688,7702,7705,7731,7754,7755,7759,7765,7771,7772,7786,7787,7825,7832,7838,7839,7848,7881,7883,7884,7885,7887,7895,7905,7907,7909,7910,7913,7916,7918,7921,7922,7923,7924,7927,7937,7963,7967,7980,7981,7982,7984,7991,7993,7994,7995,7999,8000,8005,8009,8010,8012,8014,8017,8018,8022,8023,8026,8027,8029,8038,8065,8081,8089,8095,8121,8124,8125,8126,8127,8130,8134,8138,8143,8146,8151,8157,8161,8165,8168,8169,8178,8192,8193,8214,8215,8221,8222,8225,8227,8229,8230,8236,8254,8259,8260,8271,8274,8298,8299,8300,8309,8316,8317,8319,8320,8321,8341,8355,8356,8360,8361,8362,8378,8379,8380,8385,8387,8390,8392,8398,8399,8411,8414,8422,8423,8424,8425,8426,8427,8428,8430,8432,8433,8436,8437,8438,8441,8448,8449,8450,8455,8456,8458,8462,8463,8464,8465,8466,8467,8469,8470,8472,8473,8476,8477,8481,8482,8487,8490,8493,8500,8508,8510,8511,8512,8513,8514,8516,8517,8518,8519,8520,8521,8524,8525,8526,8531,8532,8533,8535,8536,8544,8545,8546,8550,8554,8556,8558,8566,8585,8587,8588,8590,8591,8594,8595,8597,8599,8602,8603,8604,8605,8606,8607,8608,8609,8611,8612,8613,8640,8651,8654,8655,8656,8657,8661,8662,8663,8664,8665,8666,8667,8669,8672,8673,8676,8678,8680,8681,8683,8686,8687,8688,8689,8691,8692,8693,8694,8695,8697,8698,8705,8710,8714,8716,8717,8718,8720,8722,8723,8725,8726,8729,8730,8731,8734,8736,8737,8740,8741,8742,8744,8746,8748,8749,8750,8752,8753,8755,8756,8757,8759,8761,8763,8764,8765,8766,8767,8768,8769,8771,8772,8773,8774,8782,8784,8785,8788,8790,8791,8793,8794,8795,8798,8806,8813,8819,8821,8824,8827,8828,8832,8833,8836,8838,8843,8844,8845,8846,8847,8848,8850,8852,8853,8854,8855,8856,8860,8868,8881,8889,8891,8893,8894,8896,8898,8906,8908,8911,8913,8917,8919,8920,8924,8925,8961,8991,8994,9086,9088,9096,9100,9101,9102,9124,9145,9146,9160,9163,9164,9215,9219,9231,9233,9252,9263,9270,9273,9296,9338,9358,9366,9371,9383,9384,9385,9395,9396,9416,9419,9422,9425,9429,9431,9443,9447,9448,9449,9464,9473,9474,9477,9491,9492,9493,9506,9512,9549,9550,9551,9565,9573,9574,9575,9587,9598,9601,9614,9682,9685,9686,9704,9709,9716,9717,9721,9723,9749,9757,9779,9806,9823,9835,9836,9876,9905,9916,9947,9960,9961,9964,9965,9970,9979,9980,9996,9999,10001,10005,10006,10031,10061,10065,10071,10076,10077,10130,10171,10180,10188,10225,10243,10247,10250,10253,10267,10288,10339,10368,10369,10380,10387,10399,10400,10401,10414,10416,10426,10432,10440,10457,10466,10473,10479,10527,10546,10554,10556,10559,10563,10573,10582,10583,10610,10617,10618,10624,10632,10634,10635,10640,10652,10654,10657,10673,10676,10681,10690,10694,10697,10713,10722,10735,10743,10778,10822,10846,10847,10860,10910,10925,10926,10927,10957,10967,10985,11017,11019,11041,11046,11047,11048,11049,11050,11054,11056,11060,11061,11062,11063,11064,11065,11066,11067,11072,11073,11076,11077,11079,11080,11093,11098,11103,11104,11113,11114,11115,11117,11129,11138,11139,11143,11156,11190,11204,11252,11253,11255,11272,11294,11328,11332,11358,11368,11391,11392,11411,11450,11475,11500,11503,11509,11518,11519,11522,11534,11540,11562,11564,11593,11594,11607,11608,11626,11656,11657,11658,11663,11664,11666,11669,11671,11689,11690,11721,11724,11725,11728,11740,11742,11769,11770,11773,11775,11777,11786,11791,11793,11804,11819,11831,11836,11838,11848,11849,11855,11858,11860,11863,11874,11880,11881,11893,11903,11906,11907,11912,11932,11936,11937,11941,11942,11943,11944,11945,11949,11950,11952,11957,11958,11970,11976,11977,11985,11999,12000,12002,12003,12006,12007,12008,12009,12010,12011,12012,12014,12017,12018,12028,12029,12034,12035,12052,12066,12068,12069,12084,12096,12109,12113,12115,12126,12130,12149,12155,12157,12158,12161,12162,12166,12180,12190,12198,12199,12200,12201,12203,12207,12213,12215,12216,12218,12220,12221,12222,12223,12226,12227,12234,12235,12237,12238,12239,12240,12241,12243,12245,12247,12248,12249,12253,12255,12266,12267,12269,12274,12275,12276,12277,12278,12279,12280,12281,12282,12283,12284,12286,12287,12288,12289,12290,12292,12293,12294,12295,12296,12297,12298,12299,12302,12307,12309,12311,12317,12318,12325,12326,12328,12331,12332,12333,12336,12344,12345,12347,12348,12350,12353,12355,12359,12360,12361,12362,12364,12365,12366,12374,12378,12383,12384,12387,12389,12390,12396,12397,12410,12411,12412,12415,12416,12417,12418,12420,12422,12424,12425,12428,12431,12432,12450,12482,12483,12492,12494,12496,12511,12520,12527,12530,12535,12537,12538,12548,12573,12576,12580,12590,12600,12602,12654,12657,12667,12671,12688,12694,12695,12700,12701,12715,12720,12733,12738,12742,12767,12772,12773,12799,12808,12809,12825,12830,12831,12840,12851,12852,12861,12889,12898,12901,12902,12907,12914,12923,12927,12940,12953,12978,12979,13010,13011,13035,13038,13042,13043,13045,13047,13072,13076,13078,13079,13108,13185,13206,13220,13226,13227,13228,13244,13295,13296,13297,13453,13516,13554,13561,13609,13627,13637,13644,13663,13714,13771,13781,13782,13784,13788,13811,13815,13820,13832,13836,13841,13844,13856,13857,13860,13861,13865,13866,13867,13868,13871,13872,13880,13881,13898,13903,13909,13910,13911,13916,13917,13929,13937,13967,13973,13978,13986,14003,14026,14029,14031,14034,14037,14041,14042,14047,14082,14112,14120,14125,14126,14141,14142,14143,14149,14168,14192,14209,14232,14251,14278,14283,14299,14321,14351,14356,14399,14432,14469,14505,14523,14527,14595,14602,14704,14717,14737,14788,14795,14824,14832,14843,14844,14848,14876,14884,14885,14920,14922,14923,14990,14992,14993,14995,15002,15003,15004,15005,15006,15008,15017,15018,15021,15024,15025,15029,15031,15035,15059,15075,15077,15078,15084,15097,15111,15112,15113,15119,15134,15136,15137,15139,15140,15141,15142,15147,15148,15151,15157,15158,15159,15160,15161,15164,15168,15170,15173,15174,15185,15186,15187,15211,15212,15213,15214,15245,15249,15259,15260,15261,15277,15292,15301,15326,15350,15362,15371,15388,15408,15426,15427,15475,15486,15491,15494,15502,15506,15514,15522,15529,15532,15534,15540,15545,15546,15549,15550,15553,15562,15571,15590,15593,15608,15609,15619,15621,15622,15623,15625,15626,15650,15652,15665,15667,15668,15675,15680,15681,15682,15685,15697,15701,15706,15712,15731,15732,15747,15765,15781,15786,15794,15795,15796,15807,15874,15887,15901,15912,15924,15934,15945,15957,15958,15959,15960,15976,15979,15981,16001,16006,16017,16018,16019,16040,16059,16064,16067,16070,16071,16072,16075,16082,16089,16090,16102,16103,16104,16106,16109,16119,16163,16175,16176,16195,16206,16208,16209,16212,16236,16250,16262,16267,16272,16282,16284,16314,16325,16331,16334,16339,16356,16360,16372,16392,16413,16424,16425,16442,16448,16453,16466,16471,16482,16483,16488,16495,16497,16502,16520,16545,16567,16596,16604,16615,16617,16619,16626,16627,16694,16702,16705,16747,16760,16764,16765,16772,16781,16793,16796,16822,16828,16840,16841,16842,16852,16862,16873,16888,16910,16913,16951,16952,16953,16971,16984,16990,16992,16997,17031,17041,17051,17071,17072,17093,17094,17096,17097,17098,17101,17103,17118,17119,17140,17153,17191,17200,17211,17213,17227,17335,17357,17359,17360,17384,17386,17399,17453,17472,17473,17476,17480,17481,17490,17493,17495,17498,17499,17505,17512,17541,17547,17557,17567,17568,17569,17570,17586,17589,17647,17657,17664,17669,17675,17678,17703,17709,17721,17749,17762,17768,17817,17818,17867,17872,17900,17967,17975,17977,17978,17981,17984,17988,17991,17993,17998,17999,18015,18027,18079,18091,18093,18094,18102,18106,18115,18118,18128,18130,18132,18134,18135,18145,18156,18198,18280,18286,18310,18336,18362,18371,18383,18386,18393,18394,18415,18416,18423,18441,18450,18452,18458,18459,18464,18465,18469,18475,18476,18484,18486,18492,18494,18514,18550,18580,18621,18623,18629,18652,18657,18675,18676,18684,18707,18711,18712,18713,18714,18727,18736,18739,18740,18743,18751,18752,18757,18758,18760,18766,18769,18782,18783,18797,18798,18799,18805,18812,18813,18814,18815,18819,18824,18829,18830,18836,18841,18845,18859,18864,18871,18877,18888,18930,18976,18977,18979,18980,18981,18982,18983,18991,18993,18994,18995,18996,18998,18999,19000,19009,19010,19015,19016,19017,19020,19021,19026,19050,19054,19064,19065,19066,19080,19100,19101,19103,19105,19106,19109,19110,19127,19128,19132,19136,19149,19154,19198,19209,19213,19215,19218,19223,19229,19231,19244,19261,19262,19310,19327,19335,19342,19346,19349,19364,19365,19366,19385,19387,19407,19418,19419,19423,19427,19434,19439,19442,19452,19454,19455,19457,19458,19459,19460,19461,19462,19463,19464,19466,19467,19468,19469,19470,19471,19472,19473,19474,19478,19479,19481,19498,19502,19503,19505,19507,19509,19510,19511,19516,19523,19524,19525,19526,19527,19529,19535,19537,19544,19545,19546,19547,19548,19549,19553,19554,19556,19560,19561,19563,19564,19567,19568,19569,19570,19571,19572,19573,19574,19575,19576,19577,19579,19580,19582,19583,19587,19588,19589,19590,19591,19592,19594,19595,19597,19599,19600,19603,19604,19605,19608,19609,19628,19638,19639,19641,19654,19674,19689,19695,19698,19756,19772,19776,19790,19792,19797,19800,19804,19828,19836,19837,19842,19874,19908,19929,19942,19951,19976,19993,20004,20011,20038,20049,20056,20058,20075,20087,20104,20107,20126,20127,20149,20174,20210,20219,20222,20239,20241,20242,20243,20257,20258,20261,20269,20283,20286,20288,20289,20290,20293,20294,20321,20322,20330,20331,20332,20347,20353,20356,20359,20371,20374,20388,20391,20393,20394,20406,20436,20438,20448,20449,20455,20464,20465,20466,20467,20469,20470,20471,20472,20474,20475,20483,20488,20498,20506,20507,20516,20519,20520,20522,20523,20524,20525,20527,20528,20529,20531,20533,20536,20544,20567,20570,20590,20607,20608,20613,20623,20624,20626,20630,20631,20632,20634,20635,20644,20645,20656,20657,20658,20659,20667,20693,20716,20718,20719,20723,20725,20741,20742,20764,20776,20781,20785,20786,20808,20820,20835,20837,20845,20854,20855,20869,20874,20906,20910,20912,20937,20938,20940,20948,20958,20959,20960,20962,20963,20964,20992,20994,20996,21002,21003,21027,21055,21072,21110,21111,21151,21157,21199,21201,21215,21218,21219,21220,21229,21230,21231,21242,21257,21258,21268,21270,21281,21290,21291,21297,21299,21327,21337,21343,21344,21345,21355,21359,21362,21366,21367,21368,21371,21379,21395,21412,21419,21423,21424,21429,21431,21438,21439,21442,21444,21445,21453,21454,21460,21463,21471,21472,21477,21487,21488,21504,21521,21532,21540,21546,21547,21550,21552,21553,21554,21555,21557,21559,21563,21576,21586,21587,21589,21596,21606,21615,21618,21620,21621,21622,21627,21638,21640,21647,21648,21649,21650,21653,21654,21655,21656,21659,21661,21676,21682,21691,21693,21697,21698,21706,21734,21736,21737,21740,21745,21749,21750,21762,21785,21792,21802,21825,21836,21841,21860,21864,21875,21876,21878,21879,21882,21884,21887,21889,21899,21909,21918,21921,21923,21926,21933,21957,21971,21973,21975,21980,21989,22000,22002,22005,22006,22026,22030,22052,22065,22068,22069,22079,22083,22088,22089,22090,22092,22094,22095,22096,22102,22103,22104,22105,22106,22110,22111,22116,22117,22118,22119,22121,22122,22123,22124,22125,22127,22129,22130,22131,22132,22138,22140,22142,22144,22146,22149,22150,22151,22152,22153,22155,22156,22158,22159,22160,22161,22162,22164,22170,22173,22174,22176,22177,22180,22181,22182,22195,22196,22203,22211,22212,22214,22227,22231,22236,22237,22238,22239,22240,22246,22251,22257,22259,22260,22269,22270,22277,22278,22279,22280,22281,22288,22295,22298,22302,22303,22306,22307,22310,22314,22316,22318,22327,22341,22342,22344,22345,22348,22349,22351,22352,22353,22355,22356,22357,22358,22359,22360,22361,22362,22363,22365,22366,22367,22369,22373,22375,22380,22386,22387,22391,22392,22393,22394,22399,22400,22401,22403,22404,22405,22406,22407,22409,22410,22411,22416,22417,22419,22420,22422,22423,22424,22427,22428,22429,22430,22437,22439,22442,22443,22445,22456,22457,22458,22459,22460,22464,22469,22475,22478,22479,22483,22484,22485,22486,22487,22492,22493,22497,22499,22500,22505,22507,22513,22514,22516,22519,22521,22523,22525,22526,22528,22529,22536,22538,22539,22540,22542,22543,22545,22546,22547,22548,22549,22550,22551,22558,22561,22569,22571,22575,22576,22579,22587,22588,22589,22593,22594,22596,22599,22601,22603,22607,22625,22626,22627,22632,22633,22634,22635,22636,22638,22639,22640,22641,22642,22644,22645,22652,22654,22655,22656,22661,22665,22666,22668,22669,22671,22672,22675,22676,22677,22678,22684,22722,22725,22730,22731,22735,22736,22781,22802,22804,22806,22808,22813,22814,22815,22865,22868,22886,22907,22909,22911,22921,22924,22928,22932,22936,22959,22963,22973,22976,22977,22978,22979,22980,22983,22984,22985,22986,22987,22989,22990,22992,23007,23013,23017,23018,23019,23024,23027,23028,23029,23030,23031,23033,23034,23035,23036,23037,23039,23040,23041,23042,23043,23044,23045,23046,23047,23049,23052,23054,23055,23056,23059,23061,23063,23064,23065,23074,23075,23076,23077,23078,23079,23080,23081,23085,23086,23087,23088,23089,23090,23092,23093,23094,23095,23096,23097,23098,23099,23101,23109,23110,23118,23128,23133,23135,23136,23139,23142,23144,23156,23157,23161,23167,23168,23170,23171,23172,23176,23180,23191,23193,23200,23201,23202,23204,23206,23207,23208,23210,23211,23212,23214,23215,23216,23218,23219,23221,23222,23223,23224,23225,23226,23227,23228,23229,23238,23249,23250,23251,23252,23253,23254,23256,23257,23258,23262,23263,23266,23267,23268,23271,23272,23273,23274,23275,23276,23277,23278,23279,23281,23282,23285,23288,23289,23290,23291,23292,23293,23295,23296,23297,23299,23301,23303,23304,23305,23308,23309,23310,23311,23312,23315,23316,23318,23323,23326,23336,23337,23339,23341,23343,23344,23345,23352,23379,23381,23382,23383,23384,23385,23386,23387,23388,23390,23396,23403,23405,23408,23414,23415,23416,23422,23431,23432,23433,23434,23438,23439,23441,23442,23443,23444,23445,23446,23447,23449,23455,23458,23464,23465,23466,23467,23468,23470,23471,23473,23474,23475,23486,23487,23489,23494,23497,23499,23504,23507,23509,23510,23511,23512,23513,23514,23524,23526,23527,23531,23532,23533,23534,23535,23536,23537,23541,23542,23543,23544,23545,23546,23547,23548,23558,23559,23560,23563,23564,23569,23571,23572,23581,23587,23593,23596,23598,23610,23611,23612,23613,23616,23639,23643,23644,23645,23646,23648,23650,23654,23655,23656,23657,23658,23660,23661,23662,23663,23664,23668,23669,23670,23671,23672,23673,23675,23676,23678,23679,23681,23682,23683,23702,23703,23704,23705,23707,23709,23712,23740,23744,23745,23748,23749,23751,23753,23755,23757,23776,23778,23780,23781,23782,23784,23785,23786,23791,23794,23796,23804,23841,23848,23863,23865,23866,23879,23882,23883,23885,23898,23899,23905,23906,23907,23908,23909,23910,23913,23914,23915,23916,23917,23920,23921,23923,23924,23927,23928,23930,23931,23951,23954,23955,23956,23957,23958,23960,23964,23965,23966,23968,23972,23973,23974,23975,23977,23978,23979,23980,23981,23982,23998,24028,24031,24032,24039,24043,24047,24050,24052,24054,24056,24061,24065,24066,24068,24075,24080,24081,24082,24083,24090,24094,24098,24099,24131,24147,24166,24193,24195,24208,24221,24246,24278,24303,24313,24356,24365,24367,24375,24377,24378,24379,24380,24381,24382,24383,24384,24385,24386,24387,24390,24401,24405,24408,24409,24412,24424,24448,24459,24460,24462,24468,24469,24472,24476,24481,24484,24485,24486,24488,24491,24492,24493,24494,24495,24496,24498,24499,24501,24502,24503,24504,24505,24506,24510,24511,24512,24516,24517,24518,24520,24521,24523,24528,24534,24540,24550,24552,24560,24561,24582,24585,24587,24588,24590,24591,24595,24597,24600,24602,24606,24608,24609,24610,24612,24615,24622,24629,24631,24632,24634,24635,24637,24639,24641,24642,24644,24646,24647,24648,24653,24655,24656,24659,24660,24669,24673,24674,24684,24699,24700,24715,24730,24734,24740,24742,24743,24748,24751,24756,24762,24764,24770,24786,24796,24801,24803,24819,24831,24832,24833,24834,24836,24843,24851,24853,24854,24857,24890,24907,24910,24915,24923,24926,24930,24932,24936,24951,24954,24955,24965,24966,24968,24979,24995,24999,25000,25023,25056,25070,25073,25075,25076,25077,25078,25079,25083,25085,25086,25087,25105,25117,25126,25135,25139,25145,25146,25150,25155,25171,25176,25177,25178,25179,25185,25204,25205,25230,25269,25280,25283,25300,25301,25310,25313,25314,25317,25318,25320,25321,25322,25325,25333,25335,25346,25355,25358,25453,25459,25460,25461,25462,25464,25471,25489,25492,25494,25498,25509,25510,25514,25522,25523,25535,25540,25541,25542,25544,25545,25547,25548,25549,25550,25552,25569,25601,25607,25617,25618,25619,25620,25621,25622,25624,25628,25629,25630,25631,25639,25648,25650,25653,25654,25655,25657,25660,25662,25666,25669,25671,25672,25675,25677,25679,25680,25681,25682,25684,25687,25692,25693,25711,25714,25718,25731,25745,25752,25756,25771,25772,25792,25796,25800,25807,25842,25843,25849,25857,25859,25872,25873,25875,25890,25909,25930,25943,25954,25960,25961,25962,25968,25972,25973,25994,26014,26025,26028,26037,26049,26050,26052,26054,26055,26056,26057,26062,26063,26064,26072,26076,26077,26082,26084,26085,26089,26092,26101,26108,26111,26113,26119,26121,26123,26124,26125,26136,26153,26156,26157,26165,26166,26167,26174,26182,26183,26184,26185,26189,26198,26206,26208,26228,26233,26235,26248,26258,26259,26263,26272,26275,26277,26280,26281,26285,26294,26309,26316,26361,26362,26365,26369,26370,26372,26373,26374,26378,26379,26381,26385,26392,26393,26396,26406,26411,26414,26415,26420,26423,26425,26426,26428,26430,26434,26439,26440,26441,26442,26443,26447,26448,26449,26450,26451,26452,26453,26454,26457,26458,26459,26477,26492,26514,26525,26549,26567,26573,26578,26579,26580,26584,26585,26594,26597,26601,26604,26606,26612,26614,26615,26616,26617,26618,26625,26629,26634,26635,26637,26638,26643,26648,26652,26658,26668,26671,26675,26676,26677,26684,26687,26688,26694,26696,26700,26702,26703,26704,26709,26755,26757,26762,26763,26768,26771,26799,26815,26818,26828,26829,26831,26832,26833,26834,26836,26849,26850,26851,26854,26855,26856,26859,26864,26869,26879,26880,26889,26891,26893,26898,26901,26904,26908,26911,26915,26916,26917,26918,26922,26929,26934,26941,26949,26971,26980,26982,26983,26985,26989,26991,26993,26994,26995,27012,27017,27019,27020,27022,27024,27028,27031,27039,27058,27060,27063,27069,27072,27078,27079,27083,27086,27101,27103,27109,27111,27117,27122,27123,27130,27131,27156,27167,27168,27172,27175,27178,27179,27180,27182,27183,27184,27185,27186,27187,27188,27189,27190,27191,27192,27193,27194,27195,27196,27197,27198,27200,27201,27204,27205,27206,27207,27212,27214,27216,27218,27221,27223,27227,27230,27237,27238,27241,27245,27251,27255,27257,27259,27260,27261,27262,27265,27266,27267,27270,27271,27272,27273,27274,27276,27278,27279,27280,27285,27286,27290,27291,27294,27296,27297,27302,27305,27306,27308,27309,27310,27312,27313,27316,27317,27319,27323,27341,27343,27345,27347,27350,27354,27355,27356,27357,27358,27359,27360,27364,27366,27368,27369,27370,27372,27374,27376,27377,27378,27382,27383,27384,27385,27387,27403,27404,27431,27448,27469,27475,27479,27483,27507,27512,27514,27515,27527,27528,27535,27563,27570,27581,27598,27633,27668,27676,27689,27703,27704,27705,27708,27709,27710,27711,27712,27715,27720,27724,27727,27728,27730,27731,27734,27735,27749,27750,27761,27762,27763,27775,27777,27778,27791,27805,27810,27813,27814,27824,27826,27827,27843,27847,27848,27849,27850,27851,27853,27854,27855,27856,27858,27860,27861,27862,27864,27865,27866,27869,27870,27871,27872,27873,27874,27876,27877,27881,27882,27885,27886,27887,27890,27892,27895,27896,27897,27898,27900,27901,27902,27906,27908,27911,27912,27913,27914,27915,27916,27917,27918,27919,27921,27936,27937,27945,27953,27968,27969,27976,27977,27985,27986,27987,27989,27994,27995,27998,27999,28009,28012,28013,28014,28016,28019,28020,28021,28023,28024,28025,28026,28029,28030,28032,28033,28036,28037,28038,28039,28043,28046,28047,28048,28062,28064,28065,28066,28071,28072,28073,28075,28076,28077,28078,28079,28082,28088,28108,28109,28120,28122,28128,28130,28131,28136,28142,28143,28149,28153,28155,28174,28178,28180,28181,28183,28187,28188,28189,28191,28197,28204,28214,28228,28231,28233,28234,28235,28236,28237,28238,28239,28241,28243,28245,28246,28249,28251,28255,28259,28260,28262,28263,28264,28265,28267,28270,28274,28281,28283,28284,28288,28308,28327,28335,28336,28337,28338,28339,28347,28360,28361,28374,28375,28376,28382,28389,28390,28401,28412,28414,28415,28416,28417,28420,28421,28422,28434,28435,28437,28441,28445,28451,28460,28462,28463,28472,28473,28477,28488,28489,28490,28496,28505,28517,28518,28519,28521,28522,28524,28525,28526,28530,28531,28536,28542,28543,28544,28545,28546,28550,28551,28552,28553,28554,28557,28558,28559,28560,28561,28563,28564,28565,28567,28568,28569,28570,28571,28572,28573,28576,28580,28582,28586,28587,28590,28591,28593,28600,28605,28609,28611,28612,28613,28614,28615,28616,28618,28619,28624,28629,28632,28633,28634,28647,28648,28650,28657,28658,28670,28673,28677,28711,28723,28724,28734,28737,28743,28745,28747,28748,28754,28757,28759,28763,28764,28770,28771,28779,28781,28782,28801,28804,28810,28819,28820,28823,28828,28832,28833,28834,28835,28836,28837,28838,28844,28850,28851,28852,28858,28859,28861,28866,28869,28877,28888,28891,28898,28899,28900,28904,28905,28918,28921,28937,28939,28940,28945,28951,28952,28958,28960,28961,28968,28985,28986,28995,28996,29005,29016,29018,29053,29083,29086,29097,29110,29111,29126,29151,29156,29157,29162,29164,29176,29184,29190,29191,29203,29204,29226,29229,29231,29232,29241,29244,29252,29256,29257,29263,29268,29276,29278,29282,29303,29313,29315,29320,29339,29344,29345,29346,29347,29348,29368,29375,29377,29402,29405,29424,29427,29428,29431,29432,29436,29437,29438,29439,29445,29453,29454,29455,29460,29463,29465,29468,29476,29477,29478,29482,29483,29488,29492,29493,29495,29496,29497,29498,29503,29505,29507,29515,29526,29528,29537,29538,29539,29541,29552,29554,29555,29561,29563,29567,29572,29574,29578,29583,29586,29587,29593,29594,29595,29597,29598,29599,29600,29602,29603,29605,29606,29610,29613,29618,29619,29622,29625,29627,29628,29630,29631,29634,29637,29642,29643,29648,29650,29651,29652,29653,29654,29655,29659,29676,29678,29680,29683,29686,29690,29702,29705,29707,29713,29714,29715,29717,29719,29720,29722,29727,29728,29729,29731,29732,29733,29734,29735,29736,29742,29743,29747,29752,29753,29755,29756,29760,29762,29763,29764,29769,29772,29773,29774,29776,29779,29780,29781,29782,29783,29786,29787,29788,29789,29790,29792,29793,29796,29799,29800,29801,29804,29805,29807,29808,29810,29811,29813,29818,29820,29821,29822,29823,29824,29830,29835,29838,29844,29849,29850,29858,29860,29861,29876,29890,29900,29908,29912,29913,29928,29940,29943,29944,29947,29950,29951,29952,29954,29955,29956,29957,29958,29959,29960,29961,29962,29963,29966,29967,29968,29970,29971,29972,29979,29986,29989,29990,29991,29999,30000,30002,30004,30010,30012,30015,30020,30033,30034,30037,30047,30053,30062,30065,30066,30070,30071,30078,30079,30083,30084,30086,30088,30089,30097,30139,30141,30142,30143,30145,30146,30150,30152,30154,30158,30159,30163,30164,30165,30167,30168,30173,30175,30176,30179,30181,30186,30191,30195,30196,30197,30211,30225,30227,30231,30235,30239,30246,30248,30250,30251,30255,30256,30259,30263,30266,30287,30296,30301,30306,30310,30319,30338,30359,30366,30367,30376,30377,30385,30387,30390,30401,30402,30403,30404,30426,30429,30433,30438,30443,30444,30445,30446,30448,30458,30464,30465,30467,30469,30470,30471,30473,30474,30475,30476,30477,30481,30482,30483,30485,30486,30487,30491,30494,30498,30501,30510,30546,30548,30551,30552,30554,30555,30558,30559,30560,30563,30571,30573,30578,30581,30583,30584,30595,30596,30597,30600,30609,30610,30613,30639,30641,30646,30647,30648,30670,30674,30704,30705,30723,30725,30726,30743,30744,30759,30760,30774,30779,30780,30787,30803,30815,30830,30831,30834,30835,30838,30840,30841,30845,30849,30866,30870,30874,30875,30877,30879,30883,30885,30887,30888,30890,30891,30895,30897,30898,30900,30901,30902,30905,30906,30907,30909,30910,30911,30912,30913,30915,30916,30921,30922,30923,30925,30929,30930,30931,30936,30955,30956,30962,30964,30965,30966,30970,30971,30986,30996,31020,31023,31028,31031,31038,31039,31041,31043,31044,31045,31046,31047,31049,31051,31052,31055,31056,31064,31066,31067,31068,31069,31072,31074,31076,31077,31078,31079,31080,31081,31082,31084,31085,31089,31090,31092,31094,31096,31097,31098,31099,31100,31103,31106,31111,31114,31115,31116,31120,31122,31123,31127,31131,31140,31144,31154,31156,31161,31171,31174,31175,31176,31177,31178,31179,31180,31182,31183,31184,31185,31186,31187,31188,31189,31191,31192,31193,31195,31196,31199,31202,31208,31209,31210,31211,31213,31214,31217,31220,31222,31223,31225,31226,31227,31228,31229,31234,31235,31236,31237,31238,31241,31245,31246,31247,31251,31252,31253,31254,31255,31256,31258,31270,31271,31272,31274,31275,31277,31281,31285,31287,31290,31292,31293,31314,31315,31316,31317,31321,31322,31323,31327,31330,31331,31332,31342,31345,31352,31353,31371,31372,31387,31388,31398,31401,31402,31403,31405,31410,31411,31412,31413,31416,31418,31419,31435,31440,31443,31445,31446,31447,31463,31465,31467,31483,31488,31491,31494,31495,31500,31515,31538,31557,31559,31573,31587,31621,31627,31629,31630,31636,31637,31638,31639,31641,31644,31646,31648,31651,31654,31657,31658,31660,31661,31663,31664,31681,31685,31689,31698,31702,31703,31711,31717,31720,31721,31727,31729,31735,31759,31763,31771,31772,31774,31785,31788,31792,31796,31797,31801,31815,31816)
##   are outliers with |weight| = 0 ( < 3.1e-06); 
##  2766 weights are ~= 1. The remaining 23348 ones are summarized as
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0001056 0.9177000 0.9706000 0.9244000 0.9904000 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         3.143e-06 
##             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.

stargazer(modelo2, modelo2rob, type = 'text', header=FALSE)
## 
## =======================================================================
##                                           Dependent variable:          
##                                  --------------------------------------
##                                                log_money               
##                                              OLS               MM-type 
##                                                                linear  
##                                              (1)                 (2)   
## -----------------------------------------------------------------------
## aproba1                                    0.160***           0.078*** 
##                                            (0.003)             (0.001) 
##                                                                        
## as_label(r104)Mujer                        0.263***                    
##                                            (0.024)                     
##                                                                        
## r104                                                          -0.270***
##                                                                (0.008) 
##                                                                        
## r106                                       0.020***           0.011*** 
##                                            (0.001)            (0.0003) 
##                                                                        
## Constant                                   2.374***           4.793*** 
##                                            (0.042)             (0.018) 
##                                                                        
## -----------------------------------------------------------------------
## Observations                                31,818             31,818  
## R2                                          0.121               0.279  
## Adjusted R2                                 0.121               0.279  
## Residual Std. Error (df = 31814)            2.048               0.719  
## F Statistic                      1,457.927*** (df = 3; 31814)          
## =======================================================================
## Note:                                       *p<0.1; **p<0.05; ***p<0.01

Normalidad y Outliers

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.

Quitar casos

En nuestro caso podríamos quitar los gastos igual a cero. No obstante, esto resuelve un problema estadístico, pero nos supone una decisión metodológica díficil de atender.

outlierTest(modelo2)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##        rstudent unadjusted p-value Bonferroni p
## 32753 -3.292338         0.00099466           NA
mydata$rownames<-rownames(mydata)
mydata1<-subset(mydata, rownames!=43905)


modelo_no_out<-lm(log_money ~aproba1 + r104+ r106, data=mydata1, na.action=na.exclude)

outlierTest(modelo_no_out)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##        rstudent unadjusted p-value Bonferroni p
## 32753 -3.292338         0.00099466           NA
qqPlot(modelo_no_out)

## [1] 10769 32753

Quitar ceros

mydata2<-subset(mydata1, money!=0)
mydata2<-subset(mydata2, rownames!=c(26296,46120,10340,17586 ))
## Warning in rownames != c(26296, 46120, 10340, 17586): longer object length
## is not a multiple of shorter object length
modelo_normal<-lm(log_money ~aproba1 + r104+ r106, data=mydata2, na.action=na.exclude)
hist(modelo_normal$residuals)

¿Cuando parar?

qqPlot(modelo_normal)

## [1] 12888 15705
outlierTest(modelo_normal)
##        rstudent unadjusted p-value Bonferroni p
## 12888  6.809466         1.0008e-11   2.6145e-07
## 15705 -6.434858         1.2573e-10   3.2845e-06
## 20272 -6.275358         3.5431e-10   9.2556e-06
## 14352 -5.872850         4.3355e-09   1.1326e-04
## 18460 -5.546337         2.9454e-08   7.6942e-04
## 23257 -5.410154         6.3524e-08   1.6594e-03
## 18657  5.399039         6.7583e-08   1.7655e-03
## 24348 -5.367321         8.0595e-08   2.1054e-03
## 15880 -5.328613         9.9785e-08   2.6067e-03
## 18420 -5.182851         2.2014e-07   5.7508e-03

¡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.

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.

Por ello es recomendable mejor utilizar otro tipo de modelos más robustos a la presencia de outliers(regresión robusta) y menos dependientes de una distribución normal(regresión mediana).