Una empresa sin fines de lucro desea tener un modelo de los ingresos en la sociedad mexicana, esto para diferentes regiones ya que dicha empresa considera que cada estado tiene diferentes preocupaciones presupuestales. Esto debe hacerse mediante La Encuesta Nacional de Ingreso-Gasto de hogares 2018 del INEGI la cual contiene los montos que cada hogar gasta en diferentes conceptos.

El objetivo de dicho modelo es determinar que conceptos impactan más en el gasto de los hogares de dichas regiones, así como una buena estimación del ingreso en años posteriores a los de la captura de información.

En estre trabbajo se analizará al estado de Oaxaca

Oaxaca <- read_excel("~/Downloads/Oaxaca.xlsx") #Cargo mus datos
#Oaxaca <- read_excel("~/Downloads/Nuevo_Leon.xlsx") #Cargo mus datos
names(Oaxaca)[names(Oaxaca)=="entidad"]="num_familia" #cambio el nombre de la variable entidad por num_familia
Oaxaca$num_familia = row_number(Oaxaca$num_familia) #cambio mis datos de la entidad y a cada familia le asigno un número para poder identificarlos mas facilmente
Oaxaca

Análisis descriptivo de la base de datos:

fancy_summary <- skim_with(
  numeric = sfl(
    Min = min,
    Max = max,
    Q1 = ~ quantile(., probs = .25),
    Median = ~quantile(., probs = .50),
    Q3 = ~ quantile(., probs = .75),
    Mean = mean,
    Sd = sd,
  #  hist = ~ inline_hist(., 5)
  ),
  append = FALSE,  
)
fancy_summary(Oaxaca)
Data summary
Name Oaxaca
Number of rows 2279
Number of columns 12
_______________________
Column type frequency:
numeric 12
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate Min Max Q1 Median Q3 Mean Sd
num_familia 0 1 1 2279.00 570.50 1140.00 1709.50 1140.00 658.03
edad_jefe 0 1 17 110.00 39.00 51.00 65.00 52.14 16.60
tot_integ 0 1 1 18.00 2.00 3.00 5.00 3.57 1.87
ing_cor 0 1 0 2233942.19 11329.25 19713.52 35069.65 29966.98 56036.65
alimentos 0 1 0 64285.64 3889.19 6711.31 10626.24 8047.56 6104.59
vesti_calz 0 1 0 16409.30 0.00 440.20 1164.11 915.53 1423.38
vivienda 0 1 0 33098.00 206.25 549.66 1590.55 1141.32 1677.63
limpieza 0 1 0 33707.24 423.43 759.88 1314.00 1229.96 2012.50
salud 0 1 0 41391.18 0.00 58.69 538.02 675.07 2167.84
transporte 0 1 0 87108.23 385.71 1592.61 4203.68 3292.94 5370.11
educa_espa 0 1 0 127434.00 0.00 539.99 1819.16 1926.95 5240.28
personales 0 1 0 39632.03 421.30 865.15 1606.67 1321.74 1903.98

Realizo un boxplot de todos mis datos para ver mis valores atípicos:

boxplot(Oaxaca)

  • Se observa en particular los valores super alejados en ingreso corriente, procedo a ver cuales son estos dos datos para asi ver si puedo eliminarlos de mis datos
Oaxaca %>% filter(ing_cor > 500000)
  • Tambien buscare todas las personas que su ingreso sea igual a 0
Oaxaca %>% filter(ing_cor <= 0) #como solo es una peprsona, procedere a eliminarlo

Procedo a eliminar los 3 datos obtenidos anteriormente

Oaxaca = Oaxaca %>% distinct() %>% filter(num_familia != 5 & num_familia != 380 & num_familia != 967)
Oaxaca
boxplot(Oaxaca)

Regresión Lineal Simple

Análisis de correlacción con el Coeficiente de Spearman

ggcorrplot(cor(Oaxaca,method = "spearman"), 
           #method = "circle"  #Método de visualización, "square" por default.
           hc.order = TRUE, #Orden jerárquico sobre los valores obtenidos.
           outline.col = "white", #Color del margen de los cuadrádos o círculos.
           type = "lower", #Elementos a desplegar en la gráfica
           lab = TRUE # Anotaciones de los valores obtenidos en cada celda de la matriz.
           #colors = c("green", "black", "red") #Colores 
           )

  • Checare la cantidad de 0 que hay en cada variable para checar si más adelante seria posible no contarlos en mis siguientes pruebas. (Si mi cantadad de 0 pasan del 5% no es posible eliminarlos)
var1<-colnames(Oaxaca)
casos=0
for (i in 1:length(var1)){
 casos[i]=length(which(Oaxaca[,i]==0))
}
cbind(var1,casos)
##       var1          casos
##  [1,] "num_familia" "0"  
##  [2,] "edad_jefe"   "0"  
##  [3,] "tot_integ"   "0"  
##  [4,] "ing_cor"     "0"  
##  [5,] "alimentos"   "15" 
##  [6,] "vesti_calz"  "569"
##  [7,] "vivienda"    "58" 
##  [8,] "limpieza"    "28" 
##  [9,] "salud"       "891"
## [10,] "transporte"  "297"
## [11,] "educa_espa"  "879"
## [12,] "personales"  "30"

Aquí consideramos un modelo lineal por cada variable y mostramos el resultado de la regresión entre el ingreso y mis otros datos

all_linear_models <- Oaxaca %>% 
  dplyr::select(-ing_cor,-num_familia) %>% map(~summary(lm(Oaxaca$ing_cor~.x)))
#all_linear_models

Para comparar los distintos modelos, podemos hacerlo bajo diferentes criterios;en este caso con el \(R^2\)

all_linear_models %>% map_df(~`[[`(.x, "r.squared")) %>% 
  gather("Variable predictora", "r", 1:10) %>% arrange(desc(r))

Checando los valores de mis \(R^2\) mis variables predictoras que ocupare para analizarlos es alimentos hasta educación y esparcimiento. Con ayuda de mi correlacción con el Coeficiente de Spearman aplicare transformaciones a mis variables.

Ingresos corrientes ~ Alimentos

Modelo lineal

modelo_alimentos <- lm(ing_cor~alimentos,data=Oaxaca)
summary(modelo_alimentos)
## 
## Call:
## lm(formula = ing_cor ~ alimentos, data = Oaxaca)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -69687 -11719  -4370   5595 212497 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 7.396e+03  8.266e+02   8.947   <2e-16 ***
## alimentos   2.654e+00  8.194e-02  32.396   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23830 on 2274 degrees of freedom
## Multiple R-squared:  0.3158, Adjusted R-squared:  0.3155 
## F-statistic:  1050 on 1 and 2274 DF,  p-value: < 2.2e-16

Diagrama de dispersión

Oaxaca %>% select(ing_cor,alimentos) %>% ggplot(aes(x=alimentos,y=ing_cor))  + geom_point(color = "brown", size = 1)+labs(title  =  'Ingresos corrientes ~ Alimentos', x  =  'Alimentos', y ='Ingresos corrientes' )+theme_bw() + theme(plot.title = element_text(hjust = 0.5))

Observando la gráfica podemos ver que nuestros datos necesitan alguna transformación

Modelo lineal con transformaciones sobre las variables

modelo_alimentos <-lm(log(ing_cor)~sqrt(alimentos),data=Oaxaca)
summary(modelo_alimentos)
## 
## Call:
## lm(formula = log(ing_cor) ~ sqrt(alimentos), data = Oaxaca)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.93681 -0.46327  0.00085  0.43875  2.87571 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     8.6401782  0.0405602  213.02   <2e-16 ***
## sqrt(alimentos) 0.0151630  0.0004524   33.52   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6821 on 2274 degrees of freedom
## Multiple R-squared:  0.3307, Adjusted R-squared:  0.3304 
## F-statistic:  1123 on 1 and 2274 DF,  p-value: < 2.2e-16

Gráfica con trasformaciones

Oaxaca %>% ggplot(aes(y = log(ing_cor), x = sqrt(alimentos))) + geom_point(size = 1,color = "brown") + geom_smooth(method="lm", color = "black")+labs(title  =  'Log(Ingresos corrientes) ~ Sqrt(Alimentos)', x  =  'Alimentos', y ='Ingresos corrientes' )+theme_bw() + theme(plot.title = element_text(hjust = 0.5))
## `geom_smooth()` using formula 'y ~ x'

par(mfrow = c(2, 2))
plot(modelo_alimentos, col = "darkred")

Pruebas de hipótesis y ANOVA e intervalos

broom::tidy(modelo_alimentos)

Intervalo de confianza

confint(modelo_alimentos)
##                      2.5 %     97.5 %
## (Intercept)     8.56063931 8.71971714
## sqrt(alimentos) 0.01427581 0.01605014

Supuestos

Normalidad

shapiro.test(modelo_alimentos$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_alimentos$residuals
## W = 0.99745, p-value = 0.00087

Homocedastidad

bptest(modelo_alimentos)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_alimentos
## BP = 3.0449, df = 1, p-value = 0.08099

Los resultados de mis supuestos no son lo suficientemente buenos, de echo no cumple con el supuesto de normalidad, por lo tanto procedo hacer una nueva limpieza sobre mi varriable alimentos.

Eliminación de ceros

alimentos_sin_0 = Oaxaca %>% distinct() %>% filter(alimentos != 0) #Elimino los ceros en alimentos

Modelo lineal 2

modelo_alimentos_2 <-lm(log(ing_cor)~sqrt(alimentos),data= alimentos_sin_0)
summary(modelo_alimentos_2)
## 
## Call:
## lm(formula = log(ing_cor) ~ sqrt(alimentos), data = alimentos_sin_0)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.92598 -0.46022  0.00468  0.43814  2.83682 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     8.5903742  0.0413927  207.53   <2e-16 ***
## sqrt(alimentos) 0.0156828  0.0004602   34.08   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6774 on 2259 degrees of freedom
## Multiple R-squared:  0.3396, Adjusted R-squared:  0.3393 
## F-statistic:  1162 on 1 and 2259 DF,  p-value: < 2.2e-16

Gráficas

alimentos_sin_0 %>% ggplot(aes(y = log(ing_cor), x = sqrt(alimentos))) + geom_point(size = 1,color = "brown") + geom_smooth(method="lm", color = "black")+labs(title  =  'Log(Ingresos corrientes) ~ Sqrt(Alimentos)', x  =  'Alimentos', y ='Ingresos corrientes' )+theme_bw() + theme(plot.title = element_text(hjust = 0.5))
## `geom_smooth()` using formula 'y ~ x'

par(mfrow = c(2, 2))
plot(modelo_alimentos_2, col = "darkred")

alimentos_3 = Oaxaca %>% distinct() %>% filter(num_familia != 551 & num_familia != 786 & num_familia != 442)
modelo_alimentos_3 <-lm(log(ing_cor)~sqrt(alimentos),data= alimentos_3)
summary(modelo_alimentos_3)
## 
## Call:
## lm(formula = log(ing_cor) ~ sqrt(alimentos), data = alimentos_3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.93712 -0.46333  0.00078  0.43863  2.87551 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     8.6403793  0.0405801   212.9   <2e-16 ***
## sqrt(alimentos) 0.0151644  0.0004527    33.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6824 on 2271 degrees of freedom
## Multiple R-squared:  0.3307, Adjusted R-squared:  0.3304 
## F-statistic:  1122 on 1 and 2271 DF,  p-value: < 2.2e-16

Pruebas de hipótesis y ANOVA e intervalos

broom::tidy(modelo_alimentos_2)

Intervalo de confianza

confint(modelo_alimentos_2)
##                      2.5 %     97.5 %
## (Intercept)     8.50920243 8.67154593
## sqrt(alimentos) 0.01478044 0.01658522

Supuestos

Normalidad

shapiro.test(modelo_alimentos_2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_alimentos_2$residuals
## W = 0.99797, p-value = 0.005713

Homocedastidad

bptest(modelo_alimentos_2)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_alimentos_2
## BP = 0.43808, df = 1, p-value = 0.508

Podemos observar que mis supuestos son más optimos y mejores que el anterior.

Ingresos corrientes ~ Trasporte

Modelo lineal

modelo_trasporte =lm(ing_cor~transporte,data=Oaxaca)
summary(modelo_trasporte)
## 
## Call:
## lm(formula = ing_cor ~ transporte, data = Oaxaca)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -196250  -12717   -5665    5380  236445 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.002e+04  6.159e+02   32.50   <2e-16 ***
## transporte  2.646e+00  9.775e-02   27.07   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 25050 on 2274 degrees of freedom
## Multiple R-squared:  0.2437, Adjusted R-squared:  0.2434 
## F-statistic: 732.9 on 1 and 2274 DF,  p-value: < 2.2e-16

Diagrama de dispersión

Oaxaca %>% select(ing_cor,transporte) %>% ggplot(aes(x=transporte,y=ing_cor)) + geom_point(color = "deeppink4", size = 1)+labs(title  =  'Ingresos corrientes ~ Trasporte', x  =  'Trasporte', y ='Ingresos corrientes' )+theme_bw() + theme(plot.title = element_text(hjust = 0.5))

Modelo lineal con transformaciones sobre las variables

modelo_trasporte = lm(log(ing_cor)~sqrt(transporte),data=Oaxaca)
summary(modelo_trasporte)
## 
## Call:
## lm(formula = log(ing_cor) ~ sqrt(transporte), data = Oaxaca)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.06179 -0.41973  0.03213  0.43326  2.33912 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      9.2478112  0.0221995  416.58   <2e-16 ***
## sqrt(transporte) 0.0146959  0.0003868   37.99   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6521 on 2274 degrees of freedom
## Multiple R-squared:  0.3883, Adjusted R-squared:  0.388 
## F-statistic:  1443 on 1 and 2274 DF,  p-value: < 2.2e-16

Gráficas

Oaxaca %>% ggplot(aes(y = log(ing_cor), x = sqrt(transporte))) + geom_point(color = "deeppink4", size = 1) + geom_smooth(method="lm", color = "black")+labs(title  =  'Log(Ingresos corrientes) ~ Sqrt(Transporte)', x  =  'Transporte', y ='Ingresos corrientes' )+theme_bw() + theme(plot.title = element_text(hjust = 0.5))
## `geom_smooth()` using formula 'y ~ x'

par(mfrow = c(2, 2))
plot(modelo_trasporte, col = "deeppink4")

Pruebas de hipótesis y ANOVA e intervalos

broom::tidy(modelo_trasporte)

Intervalo de confianza

confint(modelo_trasporte)
##                       2.5 %     97.5 %
## (Intercept)      9.20427780 9.29134459
## sqrt(transporte) 0.01393727 0.01545445

Supuestos

Normalidad

shapiro.test(modelo_trasporte$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_trasporte$residuals
## W = 0.99518, p-value = 9.781e-07

Homocedastidad

bptest(modelo_trasporte)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_trasporte
## BP = 51.525, df = 1, p-value = 7.07e-13

Mis supuestos no se cumple por lo que este modelo queda descartado

Ingresos corrientes ~ Limpieza

Modelo lineal

modelo_limpieza = lm(ing_cor~limpieza,data=Oaxaca)
summary(modelo_limpieza)
## 
## Call:
## lm(formula = ing_cor ~ limpieza, data = Oaxaca)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -208704  -14393   -7034    6184  239011 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.028e+04  6.265e+02   32.38   <2e-16 ***
## limpieza    6.919e+00  2.707e-01   25.56   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 25390 on 2274 degrees of freedom
## Multiple R-squared:  0.2231, Adjusted R-squared:  0.2228 
## F-statistic: 653.2 on 1 and 2274 DF,  p-value: < 2.2e-16

Diagrama de dispersión

Oaxaca %>% select(ing_cor,limpieza) %>% ggplot(aes(x=limpieza,y=ing_cor)) + geom_point(color = "green4", size = 1)+labs(title  =  'Ingresos corrientes ~ Limpieza', x  =  'Limpieza', y ='Ingresos corrientes' )+theme_bw() + theme(plot.title = element_text(hjust = 0.5))

Modelo lineal con trasformaciones sobre las variables

modelo_limpieza =lm(log(ing_cor)~sqrt(limpieza),data=Oaxaca)
summary(modelo_limpieza)
## 
## Call:
## lm(formula = log(ing_cor) ~ sqrt(limpieza), data = Oaxaca)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.83072 -0.49618 -0.00735  0.49592  2.54069 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    9.2192919  0.0319689  288.38   <2e-16 ***
## sqrt(limpieza) 0.0226814  0.0009149   24.79   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7397 on 2274 degrees of freedom
## Multiple R-squared:  0.2128, Adjusted R-squared:  0.2124 
## F-statistic: 614.6 on 1 and 2274 DF,  p-value: < 2.2e-16

Gráficas

Oaxaca %>% ggplot(aes(y = log(ing_cor), x = sqrt(limpieza)))+ geom_point(color = "green4", size = 1) + geom_smooth(method="lm", color = "black")+labs(title  =  'Log(Ingresos corrientes) ~ Sqrt(Limpieza)', x  =  'Limpieza', y ='Ingresos corrientes' )+theme_bw() + theme(plot.title = element_text(hjust = 0.5))
## `geom_smooth()` using formula 'y ~ x'

par(mfrow = c(2, 2))
plot(modelo_limpieza, col = "green4")

Pruebas de hipótesis y ANOVA e intervalos

broom::tidy(modelo_limpieza)

Intervalo de confianza

confint(modelo_limpieza)
##                     2.5 %     97.5 %
## (Intercept)    9.15660057 9.28198324
## sqrt(limpieza) 0.02088719 0.02447553

Supuestos

Normalidad

shapiro.test(modelo_limpieza$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_limpieza$residuals
## W = 0.99919, p-value = 0.4235

Homocedastidad

bptest(modelo_limpieza)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_limpieza
## BP = 0.080977, df = 1, p-value = 0.776

Por el momento este es mi mejor modelo

Ingresos corrientes ~ Vivienda

Modelo lineal

modelo_vivienda =lm(ing_cor~vivienda,data=Oaxaca)
summary(modelo_vivienda)
## 
## Call:
## lm(formula = ing_cor ~ vivienda, data = Oaxaca)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -94754 -13907  -6457   6482 262733 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.012e+04  6.555e+02   30.69   <2e-16 ***
## vivienda    7.579e+00  3.241e-01   23.39   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 25860 on 2274 degrees of freedom
## Multiple R-squared:  0.1939, Adjusted R-squared:  0.1935 
## F-statistic: 546.9 on 1 and 2274 DF,  p-value: < 2.2e-16

Diagrama de dispersión

Oaxaca %>% select(ing_cor,vivienda) %>% ggplot(aes(x=vivienda,y=ing_cor)) + geom_point(color = "yellow4", size = 1)+labs(title  =  'Ingresos corrientes ~ Vivienda', x  =  'Vivienda', y ='Ingresos corrientes' )+theme_bw() + theme(plot.title = element_text(hjust = 0.5))

Modelo lineal con trasformaciones sobre las variables

modelo_vivienda =lm(log(ing_cor)~sqrt(vivienda),data=Oaxaca)
summary(modelo_vivienda)
## 
## Call:
## lm(formula = log(ing_cor) ~ sqrt(vivienda), data = Oaxaca)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.88522 -0.49756  0.00779  0.51463  2.35943 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    9.2873364  0.0280795  330.75   <2e-16 ***
## sqrt(vivienda) 0.0220910  0.0008328   26.52   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7286 on 2274 degrees of freedom
## Multiple R-squared:  0.2363, Adjusted R-squared:  0.236 
## F-statistic: 703.6 on 1 and 2274 DF,  p-value: < 2.2e-16

Gráficas

Oaxaca %>% ggplot(aes(y = log(ing_cor), x = sqrt(vivienda))) + geom_point(color = "yellow4", size = 1) + geom_smooth(method="lm", color = "black")+labs(title  =  'Log(Ingresos corrientes) ~ Sqrt(Vivienda)', x  =  'Vivienda', y ='Ingresos corrientes' )+theme_bw() + theme(plot.title = element_text(hjust = 0.5))
## `geom_smooth()` using formula 'y ~ x'

par(mfrow = c(2, 2))
plot(modelo_vivienda, col = "yellow4")

Pruebas de hipótesis y ANOVA e intervalos

broom::tidy(modelo_vivienda)

Intervalo de confianza

confint(modelo_vivienda)
##                    2.5 %     97.5 %
## (Intercept)    9.2322722 9.34240052
## sqrt(vivienda) 0.0204578 0.02372417

Supuestos

Normalidad

shapiro.test(modelo_vivienda$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_vivienda$residuals
## W = 0.99842, p-value = 0.02835

Homocedastidad

bptest(modelo_vivienda)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_vivienda
## BP = 19.875, df = 1, p-value = 8.266e-06

#Ingresos corrientes ~ Personales

Modelo lineal

modelo_personales =lm(ing_cor~personales,data=Oaxaca)
summary(modelo_personales)
## 
## Call:
## lm(formula = ing_cor ~ personales, data = Oaxaca)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -144864  -14022   -7240    6038  237629 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.050e+04  6.700e+02   30.60   <2e-16 ***
## personales  6.239e+00  2.898e-01   21.53   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 26250 on 2274 degrees of freedom
## Multiple R-squared:  0.1694, Adjusted R-squared:  0.169 
## F-statistic: 463.6 on 1 and 2274 DF,  p-value: < 2.2e-16

Diagrama de dispersión

Oaxaca %>% select(ing_cor,personales) %>% ggplot(aes(x=personales,y=ing_cor)) + geom_point(color = "green4", size = 1)+labs(title  =  'Ingresos corrientes ~ Personales', x  =  'Personales', y ='Ingresos corrientes' )+theme_bw() + theme(plot.title = element_text(hjust = 0.5))

Pruebas de hipótesis y ANOVA e intervalos

broom::tidy(modelo_personales)

Intervalo de confianza

confint(modelo_personales)
##                    2.5 %       97.5 %
## (Intercept) 19189.388027 21817.195121
## personales      5.671083     6.807554

Modelo lineal con transformaciones sobre las variables

modelo_personales=lm(log(ing_cor)~sqrt(personales),data=Oaxaca)
summary(modelo_personales)
## 
## Call:
## lm(formula = log(ing_cor) ~ sqrt(personales), data = Oaxaca)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0619 -0.4622 -0.0067  0.4746  2.7552 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      9.1264877  0.0314231  290.44   <2e-16 ***
## sqrt(personales) 0.0246289  0.0008652   28.46   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7159 on 2274 degrees of freedom
## Multiple R-squared:  0.2627, Adjusted R-squared:  0.2624 
## F-statistic: 810.2 on 1 and 2274 DF,  p-value: < 2.2e-16

Gráficas

Oaxaca %>% ggplot(aes(y = log(ing_cor), x = sqrt(personales))) + geom_point(color = "green4", size = 1) + geom_smooth(method="lm", color = "black")+labs(title  =  'Log(Ingresos corrientes) ~ Sqrt(Personales)', x  =  'Personales', y ='Ingresos corrientes' )+theme_bw() + theme(plot.title = element_text(hjust = 0.5))
## `geom_smooth()` using formula 'y ~ x'

par(mfrow = c(2, 2))
plot(modelo_personales, col = "green4")

Supuestos

Normalidad

shapiro.test(modelo_personales$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_personales$residuals
## W = 0.99601, p-value = 9.728e-06

Homocedastidad

bptest(modelo_personales)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_personales
## BP = 53.826, df = 1, p-value = 2.19e-13

Este modelo no cumple con los supuestos por lo que no me es útil

Ingresos corrientes ~ Vestimenta y Calzado

Modelo lineal

modelo_vesti_calz =lm(ing_cor~vesti_calz,data=Oaxaca)
summary(modelo_vesti_calz)
## 
## Call:
## lm(formula = ing_cor ~ vesti_calz, data = Oaxaca)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -75624 -14518  -6800   5897 209678 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.114e+04  6.549e+02   32.27   <2e-16 ***
## vesti_calz  8.291e+00  3.868e-01   21.43   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 26280 on 2274 degrees of freedom
## Multiple R-squared:  0.1681, Adjusted R-squared:  0.1677 
## F-statistic: 459.3 on 1 and 2274 DF,  p-value: < 2.2e-16

Diagrama de dispersión

Oaxaca %>% select(ing_cor,vesti_calz) %>% ggplot(aes(x=vesti_calz,y=ing_cor)) + geom_point(color = "red", size = 1)+labs(title  =  'Ingresos corrientes ~ Vestimenta y Calzado', x  =  'Vestimenta y Calzado', y ='Ingresos corrientes' )+theme_bw() + theme(plot.title = element_text(hjust = 0.5))

Modelo lineal con trasformaciones sobre las variables

modelo_vesti_calz =lm(log(ing_cor)~sqrt(vesti_calz),data=Oaxaca)
summary(modelo_vesti_calz)
## 
## Call:
## lm(formula = log(ing_cor) ~ sqrt(vesti_calz), data = Oaxaca)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.01620 -0.52090 -0.00388  0.49200  2.82575 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      9.515974   0.024095  394.93   <2e-16 ***
## sqrt(vesti_calz) 0.017415   0.000796   21.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7578 on 2274 degrees of freedom
## Multiple R-squared:  0.1739, Adjusted R-squared:  0.1735 
## F-statistic: 478.6 on 1 and 2274 DF,  p-value: < 2.2e-16

Gráficas

Oaxaca %>% ggplot(aes(y = log(ing_cor), x = sqrt(vesti_calz))) + geom_point(color = "red", size = 1) + geom_smooth(method="lm", color = "black")+labs(title  =  'Log(Ingresos corrientes) ~ Sqrt(Vestimenta y Calzado)', x  =  'Vestimenta y Calzado', y ='Ingresos corrientes' )+theme_bw() + theme(plot.title = element_text(hjust = 0.5))
## `geom_smooth()` using formula 'y ~ x'

par(mfrow = c(2, 2))
plot(modelo_vesti_calz, col = "red")

Supuestos

Normalidad

shapiro.test(modelo_vesti_calz$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_vesti_calz$residuals
## W = 0.99764, p-value = 0.001692

Homocedastidad

bptest(modelo_vesti_calz)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_vesti_calz
## BP = 20.469, df = 1, p-value = 6.062e-06

Ingresos corrientes ~ Educación y Esparcimiento

Modelo lineal

modelo_educa_espa =lm(ing_cor~educa_espa,data=Oaxaca)
summary(modelo_educa_espa)
## 
## Call:
## lm(formula = ing_cor ~ educa_espa, data = Oaxaca)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -126453  -15166   -7466    6166  240725 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.448e+04  5.890e+02   41.56   <2e-16 ***
## educa_espa  2.210e+00  1.055e-01   20.95   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 26370 on 2274 degrees of freedom
## Multiple R-squared:  0.1618, Adjusted R-squared:  0.1615 
## F-statistic:   439 on 1 and 2274 DF,  p-value: < 2.2e-16

Diagrama de dispersión

Oaxaca %>% select(ing_cor,educa_espa) %>% ggplot(aes(x=educa_espa,y=ing_cor)) + geom_point(color = "aquamarine3", size = 1)+labs(title  =  'Ingresos corrientes ~ Educación y Esparcimiento', x  =  'Educación y Esparcimiento', y ='Ingresos corrientes' )+theme_bw() + theme(plot.title = element_text(hjust = 0.5))

Modelo lineal con transformaciones sobre las variables

modelo_educa_espa =lm(log(ing_cor)~sqrt(educa_espa),data=Oaxaca)
summary(modelo_educa_espa)
## 
## Call:
## lm(formula = log(ing_cor) ~ sqrt(educa_espa), data = Oaxaca)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.75787 -0.51505 -0.02792  0.47833  2.68202 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      9.5977926  0.0200676  478.27   <2e-16 ***
## sqrt(educa_espa) 0.0113024  0.0004573   24.72   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7402 on 2274 degrees of freedom
## Multiple R-squared:  0.2118, Adjusted R-squared:  0.2114 
## F-statistic: 610.9 on 1 and 2274 DF,  p-value: < 2.2e-16

Gráficas

Oaxaca %>% ggplot(aes(y = log(ing_cor), x = sqrt(educa_espa))) + geom_point(color = "aquamarine3", size = 1) + geom_smooth(method="lm", color = "black")+labs(title  =  'Log(Ingresos corrientes) ~ Sqrt(Educación y Esparcimiento)', x  =  'Educación y Esparcimiento', y ='Ingresos corrientes' )+theme_bw() + theme(plot.title = element_text(hjust = 0.5))
## `geom_smooth()` using formula 'y ~ x'

par(mfrow = c(2, 2))
plot(modelo_educa_espa, col = "aquamarine3")

Supuestos

Normalidad

shapiro.test(modelo_educa_espa$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_educa_espa$residuals
## W = 0.99801, p-value = 0.006223

Homocedastidad

bptest(modelo_educa_espa)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_educa_espa
## BP = 5.5408, df = 1, p-value = 0.01858

Conclusiones

Mis 2 mejores modelos es limpieza y el segundo de alimentos ya que cumple con todos los supestos.

Regresión Lineal Multiple

Método Forward

Las variables que ocupare para mi análisis serán las mismas que ocupe en la simple

  • Primero creo una tabla con solo las variables que me interesan y que sus datos ya esten transformados con logaritmo y raíz
#install.packages("GGally")
#library(GGally)
Oaxaca_2 <- Oaxaca %>% filter(alimentos != 0 & limpieza != 0)%>%  
  dplyr::select(ing_cor, alimentos,vivienda,limpieza,transporte,vesti_calz,educa_espa, personales) %>% 
  mutate(ing_cor = log(ing_cor),
         alimentos = sqrt(alimentos),
         vivienda = sqrt(vivienda),
         limpieza = sqrt(limpieza),         
         transporte = sqrt(transporte),
         vesti_calz = sqrt(vesti_calz),
         educa_espa = sqrt(educa_espa),
         personales = sqrt(personales))

Relación entre Variables

Oaxaca_2 %>% GGally::ggscatmat() 

### Generación de Modelos y elección del Mejor Modelo

# La función regsubsets() nos permite generar los mejores modelos de acuerdo al números de regresores y tomando como criterio RSS, además nos permite elegir el método que queramos: forward o backward.
mejores_modelos=regsubsets(ing_cor~., data= Oaxaca_2, method="forward",nvmax=7)
mejores=summary(mejores_modelos)

mejores
## Subset selection object
## Call: regsubsets.formula(ing_cor ~ ., data = Oaxaca_2, method = "forward", 
##     nvmax = 7)
## 7 Variables  (and intercept)
##            Forced in Forced out
## alimentos      FALSE      FALSE
## vivienda       FALSE      FALSE
## limpieza       FALSE      FALSE
## transporte     FALSE      FALSE
## vesti_calz     FALSE      FALSE
## educa_espa     FALSE      FALSE
## personales     FALSE      FALSE
## 1 subsets of each size up to 7
## Selection Algorithm: forward
##          alimentos vivienda limpieza transporte vesti_calz educa_espa
## 1  ( 1 ) " "       " "      " "      "*"        " "        " "       
## 2  ( 1 ) "*"       " "      " "      "*"        " "        " "       
## 3  ( 1 ) "*"       "*"      " "      "*"        " "        " "       
## 4  ( 1 ) "*"       "*"      " "      "*"        " "        " "       
## 5  ( 1 ) "*"       "*"      " "      "*"        " "        "*"       
## 6  ( 1 ) "*"       "*"      "*"      "*"        " "        "*"       
## 7  ( 1 ) "*"       "*"      "*"      "*"        "*"        "*"       
##          personales
## 1  ( 1 ) " "       
## 2  ( 1 ) " "       
## 3  ( 1 ) " "       
## 4  ( 1 ) "*"       
## 5  ( 1 ) "*"       
## 6  ( 1 ) "*"       
## 7  ( 1 ) "*"

El mejor modelo de acuerdo al BIC

Si tomamos como criterio el BIC entonces tenmos lo siguiente:

ggplot(data=data.frame(n_predictores=1:7, BIC= mejores$bic),
       aes(x = n_predictores, y = BIC)) +
  geom_line() +
  geom_point(size=3) +
  geom_point(aes(x = n_predictores[which.min(BIC)],
                 y = BIC[which.min(BIC)]),
             colour = "red", size = 4) +
  scale_x_continuous(breaks = c(0:19))+
  theme_bw() +
  labs(title = "Método Forward bajo el criterio BIC",
       x = "Número Predictores")

El mejor modelo de acuerdo al \(C_{p}\)

si tomamos como criterio el Cp de Mallow entonces tenemos lo siguiente:

ggplot(data = data.frame(n_predictores = 1:7, Cp= mejores$cp),
       aes(x = n_predictores, y = Cp)) +
  geom_line() +
  geom_point(size=3) +
  geom_point(aes(x = n_predictores[which.min(Cp)],
                 y = Cp[which.min(Cp)]),
             colour = "blue", size = 4) +
  scale_x_continuous(breaks = c(0:19))+
  theme_bw() +
  labs(title = "Método Forward bajo el criterio Cp",
       x = "Número Predictores")

El mejor modelo de acuerdo al Coeficiente de Determinación Ajustado

Si tomamos como criterio el R2 ajustado entonces tenemos lo siguiente:

ggplot(data = data.frame(n_predictores = 1:7, R2_ajustado= mejores$adjr2),
       aes(x = n_predictores, y = R2_ajustado)) +
  geom_line() +
  geom_point(size=3) +
  geom_point(aes(x = n_predictores[which.max(R2_ajustado)],
                 y = R2_ajustado[which.max(R2_ajustado)]),
             colour = "pink", size = 4) +
  scale_x_continuous(breaks = c(0:19))+
  theme_bw() +
  labs(title = "Método Forward bajo el criterio R2 ajustado",
       x = "Número Predictores")

Mejores modelos ocupando el método Forward

De acuerdo a lo que hemos visto en la gráficas anteriores hay un top2: el modelo de 6 regresores (alimentos,vivienda,limpieza,transporte,educación y esparcimiento,personales) y el modelo con 7 regresores(alimentos,vivienda,limpieza,transporte,vestimenta y calzado,educación y esparcimiento,personales).

Modelo con 6 regresores

modelo_6_regresores=lm(ing_cor~alimentos+vivienda+limpieza+transporte+educa_espa+personales,data=Oaxaca_2)
summary(modelo_6_regresores)
## 
## Call:
## lm(formula = ing_cor ~ alimentos + vivienda + limpieza + transporte + 
##     educa_espa + personales, data = Oaxaca_2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.75784 -0.36585  0.02498  0.36638  2.26439 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 8.5533287  0.0374331 228.497  < 2e-16 ***
## alimentos   0.0052432  0.0005348   9.804  < 2e-16 ***
## vivienda    0.0074953  0.0007704   9.729  < 2e-16 ***
## limpieza    0.0049171  0.0009011   5.457 5.39e-08 ***
## transporte  0.0074569  0.0004435  16.815  < 2e-16 ***
## educa_espa  0.0023743  0.0004304   5.516 3.86e-08 ***
## personales  0.0045742  0.0009268   4.935 8.59e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5723 on 2229 degrees of freedom
## Multiple R-squared:  0.5293, Adjusted R-squared:  0.528 
## F-statistic: 417.7 on 6 and 2229 DF,  p-value: < 2.2e-16

Pruebas de hipótesis y ANOVA e intervalos

broom::tidy(modelo_6_regresores)#Las variables tienen transformaciones de raíz cuadrado

Intervalo de confianza

confint(modelo_6_regresores)
##                   2.5 %      97.5 %
## (Intercept) 8.479921397 8.626736005
## alimentos   0.004194396 0.006292021
## vivienda    0.005984551 0.009006128
## limpieza    0.003149976 0.006684173
## transporte  0.006587256 0.008326582
## educa_espa  0.001530243 0.003218393
## personales  0.002756657 0.006391646
GGally::ggcoef_model(modelo_6_regresores) 

par(mfrow = c(2, 2))
plot(modelo_6_regresores, col = "violet")

Supuestos

Normalidad
shapiro.test(modelo_6_regresores$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_6_regresores$residuals
## W = 0.99281, p-value = 4.911e-09
Homocedasticidad
bptest(modelo_6_regresores)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_6_regresores
## BP = 24.516, df = 6, p-value = 0.0004195

Este modelo no es adecuado ya que no cumple con mis supuestos

Modelo con 7 regresores

modelo_7_regresores=lm(ing_cor~alimentos+vivienda+limpieza+transporte+vesti_calz+educa_espa+personales,data=Oaxaca_2)
summary(modelo_7_regresores)
## 
## Call:
## lm(formula = ing_cor ~ alimentos + vivienda + limpieza + transporte + 
##     vesti_calz + educa_espa + personales, data = Oaxaca_2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.76741 -0.36066  0.02583  0.37439  2.24868 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 8.5592915  0.0374920 228.297  < 2e-16 ***
## alimentos   0.0050516  0.0005410   9.337  < 2e-16 ***
## vivienda    0.0075864  0.0007708   9.843  < 2e-16 ***
## limpieza    0.0047478  0.0009034   5.256 1.62e-07 ***
## transporte  0.0073780  0.0004444  16.601  < 2e-16 ***
## vesti_calz  0.0017031  0.0007547   2.257   0.0241 *  
## educa_espa  0.0022012  0.0004368   5.039 5.05e-07 ***
## personales  0.0040207  0.0009579   4.197 2.81e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5718 on 2228 degrees of freedom
## Multiple R-squared:  0.5304, Adjusted R-squared:  0.5289 
## F-statistic: 359.4 on 7 and 2228 DF,  p-value: < 2.2e-16
par(mfrow = c(2, 2))
plot(modelo_7_regresores, col = "blueviolet")

Supuestos

Normalidad
shapiro.test(modelo_7_regresores$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_7_regresores$residuals
## W = 0.99281, p-value = 4.935e-09
Homocedasticidad
bptest(modelo_7_regresores)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_7_regresores
## BP = 25.022, df = 7, p-value = 0.000752

Este modelo no es adecuado ya que no cumple con mis supuestos

Mejores modelos con 2 regresores

Modelo 1

modelo_1=lm(log(ing_cor)~sqrt(alimentos)+sqrt(limpieza),
                              data = alimentos_sin_0) #s se eliminaron los 0 de alimentos
summary(modelo_1)
## 
## Call:
## lm(formula = log(ing_cor) ~ sqrt(alimentos) + sqrt(limpieza), 
##     data = alimentos_sin_0)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.83716 -0.43869  0.01005  0.42125  2.70353 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     8.5113214  0.0408136  208.54   <2e-16 ***
## sqrt(alimentos) 0.0126486  0.0005189   24.38   <2e-16 ***
## sqrt(limpieza)  0.0109387  0.0009485   11.53   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6584 on 2258 degrees of freedom
## Multiple R-squared:  0.3763, Adjusted R-squared:  0.3758 
## F-statistic: 681.2 on 2 and 2258 DF,  p-value: < 2.2e-16

Gráficas

alimentos_sin_0 %>% ggplot(aes(y = log(ing_cor), x = sqrt(alimentos)+sqrt(limpieza))) + geom_point(color = "SkyBlue", size = 1) + geom_smooth(method="lm", color = "black")+labs(title  =  'Log(Ingresos corrientes) ~ Sqrt(Alimentos + Limpieza)', x  =  'Alimentos + Limpieza', y ='Ingresos corrientes' )+theme_bw() + theme(plot.title = element_text(hjust = 0.5))
## `geom_smooth()` using formula 'y ~ x'

par(mfrow = c(2, 2))
plot(modelo_1, col = "SkyBlue")

Pruebas de hipótesis y ANOVA e intervalos

broom::tidy(modelo_1)#Las variables tienen transformaciones de raíz cuadrado

Intervalo de confianza

confint(modelo_1)
##                       2.5 %     97.5 %
## (Intercept)     8.431285415 8.59135748
## sqrt(alimentos) 0.011630968 0.01366619
## sqrt(limpieza)  0.009078722 0.01279867
GGally::ggcoef_model(modelo_1) 

Supuestos

Normalidad

shapiro.test(modelo_1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_1$residuals
## W = 0.99847, p-value = 0.0354

Homocedastidad

bptest(modelo_1)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_1
## BP = 0.14871, df = 2, p-value = 0.9283

Multicolinearidad

modelo_1 %>% vif()
## sqrt(alimentos)  sqrt(limpieza) 
##         1.34599         1.34599

Predicción

new_data <-  data.frame(alimentos = 100,limpieza = 100)
simple_model <-lm(ing_cor~(limpieza+alimentos), data = Oaxaca_2)
  
simple_model %>% predict(new_data)
##        1 
## 10.91919

Modelo 2

modelo_2=lm(log(ing_cor)~sqrt(alimentos)+sqrt(educa_espa),
                              data = alimentos_sin_0) #s se eliminaron los 0 de alimentos
summary(modelo_2)
## 
## Call:
## lm(formula = log(ing_cor) ~ sqrt(alimentos) + sqrt(educa_espa), 
##     data = alimentos_sin_0)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.81271 -0.44071 -0.00784  0.42677  2.45360 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      8.6787264  0.0404985  214.30   <2e-16 ***
## sqrt(alimentos)  0.0126679  0.0005005   25.31   <2e-16 ***
## sqrt(educa_espa) 0.0059445  0.0004561   13.03   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6534 on 2258 degrees of freedom
## Multiple R-squared:  0.3858, Adjusted R-squared:  0.3852 
## F-statistic: 709.1 on 2 and 2258 DF,  p-value: < 2.2e-16

Gráficas

alimentos_sin_0 %>% ggplot(aes(y = log(ing_cor), x = sqrt(alimentos)+sqrt(educa_espa))) + geom_point(color = "Peru", size = 1) + geom_smooth(method="lm", color = "black")+labs(title  =  'Log(Ingresos corrientes) ~ Sqrt(Alimentos + Educación y Esparcimiento)', x  = 'Alimentos + Educación y Esparcimiento', y ='Ingresos corrientes' )+theme_bw() + theme(plot.title = element_text(hjust = 0.5))
## `geom_smooth()` using formula 'y ~ x'

par(mfrow = c(2, 2))
plot(modelo_2, col = "Peru")

Pruebas de hipótesis y ANOVA e intervalos

broom::tidy(modelo_2)#Las variables tienen transformaciones de raíz cuadrado

Intervalo de confianza

confint(modelo_2)
##                        2.5 %      97.5 %
## (Intercept)      8.599308183 8.758144648
## sqrt(alimentos)  0.011686375 0.013649463
## sqrt(educa_espa) 0.005050117 0.006838898
GGally::ggcoef_model(modelo_2) 

Supuestos

Normalidad

shapiro.test(modelo_2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_2$residuals
## W = 0.99776, p-value = 0.002746

Homocedastidad

bptest(modelo_2)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_2
## BP = 2.5891, df = 2, p-value = 0.274

Multicolinearidad

modelo_2 %>% vif()
##  sqrt(alimentos) sqrt(educa_espa) 
##         1.271578         1.271578

Predicción

new_data <-  data.frame(alimentos = 100,educa_espa = 100)
multi_model <-lm(ing_cor~(alimentos+educa_espa), data = Oaxaca_2)
  
multi_model %>% predict(new_data)
##        1 
## 10.54172

Mejores modelos con 3 regresores

Modelo 3

modelo_3=lm(log(ing_cor)~sqrt(alimentos)+sqrt(educa_espa)+sqrt(limpieza),
                              data = alimentos_sin_0) #solo los 0 de alimentos se eliminaron
summary(modelo_3)
## 
## Call:
## lm(formula = log(ing_cor) ~ sqrt(alimentos) + sqrt(educa_espa) + 
##     sqrt(limpieza), data = alimentos_sin_0)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.72994 -0.41733  0.00023  0.40391  2.21108 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      8.6015283  0.0404228 212.789   <2e-16 ***
## sqrt(alimentos)  0.0104710  0.0005383  19.452   <2e-16 ***
## sqrt(educa_espa) 0.0052336  0.0004524  11.568   <2e-16 ***
## sqrt(limpieza)   0.0092200  0.0009336   9.875   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6399 on 2257 degrees of freedom
## Multiple R-squared:  0.4112, Adjusted R-squared:  0.4104 
## F-statistic: 525.5 on 3 and 2257 DF,  p-value: < 2.2e-16

Gráficas

alimentos_sin_0 %>% ggplot(aes(y = log(ing_cor), x = sqrt(alimentos)+sqrt(educa_espa) +sqrt(limpieza))) + geom_point(color = "Red", size = 1) + geom_smooth(method="lm", color = "black")+labs(title  =  'Log(Ingresos corrientes) ~ Sqrt(Alimentos + Limpieza + Educación)', x  =  'Alimentos + Limpieza + Educación', y ='Ingresos corrientes' )+theme_bw() + theme(plot.title = element_text(hjust = 0.5))
## `geom_smooth()` using formula 'y ~ x'

par(mfrow = c(2, 2))
plot(modelo_3, col = "Red")

Supuestos

Normalidad

shapiro.test(modelo_3$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_3$residuals
## W = 0.99786, p-value = 0.003815

Homocedastidad

bptest(modelo_3)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_3
## BP = 4.5423, df = 3, p-value = 0.2086

Multicolinearidad

modelo_2 %>% vif()
##  sqrt(alimentos) sqrt(educa_espa) 
##         1.271578         1.271578

Predicción

new_data <-  data.frame(alimentos = 100,educa_espa = 100)
multi_model <-lm(ing_cor~(alimentos+educa_espa), data = Oaxaca_2)
  
multi_model %>% predict(new_data)
##        1 
## 10.54172

Modelo 4

ing_limpieza_datos = Oaxaca %>% distinct() %>% filter(alimentos != 0 & limpieza != 0)
modelo_4 =lm(log(ing_cor)~sqrt(alimentos)+ sqrt(vivienda)+sqrt(limpieza),
                              data = ing_limpieza_datos)
summary(modelo_4)
## 
## Call:
## lm(formula = log(ing_cor) ~ sqrt(alimentos) + sqrt(vivienda) + 
##     sqrt(limpieza), data = ing_limpieza_datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.83604 -0.42779  0.01824  0.42066  2.29668 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     8.4475803  0.0396759 212.915   <2e-16 ***
## sqrt(alimentos) 0.0101908  0.0005282  19.295   <2e-16 ***
## sqrt(vivienda)  0.0112591  0.0008272  13.611   <2e-16 ***
## sqrt(limpieza)  0.0092118  0.0009470   9.728   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6306 on 2232 degrees of freedom
## Multiple R-squared:  0.4278, Adjusted R-squared:  0.427 
## F-statistic: 556.2 on 3 and 2232 DF,  p-value: < 2.2e-16

Gráficas

ing_limpieza_datos %>% ggplot(aes(y = log(ing_cor), x = sqrt(alimentos)+ sqrt(vivienda) +sqrt(limpieza))) + geom_point(color = "VioletRed", size = 1) + geom_smooth(method="lm", color = "black")+labs(title  =  'Log(Ingresos corrientes) ~ Sqrt(Alimentos + Vivienda + Limpieza)', x  = 'Alimentos + Vivienda + Limpieza', y ='Ingresos corrientes' )+theme_bw() + theme(plot.title = element_text(hjust = 0.5))
## `geom_smooth()` using formula 'y ~ x'

par(mfrow = c(2, 2))
plot(modelo_4, col = "VioletRed")

Supuestos

Normalidad

shapiro.test(modelo_4$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_4$residuals
## W = 0.99705, p-value = 0.0002685

Homocedastidad

bptest(modelo_4)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_4
## BP = 4.6921, df = 3, p-value = 0.1958

Observaciones

Estan son las 4 mejores pruebas de regresión lineal multiple que se encontraron

Predicción para mis modelos simples

new_data <-  data.frame(limpieza = 100)
simple_model <- lm(ing_cor~limpieza, data = Oaxaca_2)
simple_model %>% predict(new_data)
##        1 
## 11.53812
new_data <-  data.frame(alimentos = 100)
simple_model <- lm(ing_cor~alimentos, data = Oaxaca_2)
simple_model %>% predict(new_data)
##        1 
## 10.15622