Objetivo: Modelar los ingresos de la sociedad mexicana del estado de Puebla para determinar qué conceptos impactan más en el gasto de los hogares de dicha región, así como una buena estimación del ingreso en años posteriores a los de la captura de información.
Llamamos a nuestra base de datos y hacemos un primer acercamiento con ella:
library(car)
library(readxl)
data<-data.frame(read_excel("C://Users//Nath//Documents//Semestre 2021-2//Mnpyr//Proyecto//Puebla.xlsx"))
Esta base nos muestra los ingresos totales mensuales de 1,837 hogares y la cantidad de dinero que gastan en las diferentes variables: alimentos, vestimenta y calzado, vivienda, limpieza, salud, transporte, educación y esparcimiento y gastos personales. Además nos brinda la edad de la jefa o jefe de familia así como el número de integrantes de cada hogar.
Antes de continuar, eliminaremos la variable “entidad”, pues sabemos que se trata del estado 21: Puebla.
data$entidad <- NULL
Veamos una pequeña tabla resumen de nuestra base de datos:
library(skimr)
ResEst <- 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
)
ResEst(data)
Name | data |
Number of rows | 1837 |
Number of columns | 11 |
_______________________ | |
Column type frequency: | |
numeric | 11 |
________________________ | |
Group variables | None |
Variable type: numeric
skim_variable | n_missing | complete_rate | Min | Max | Q1 | Median | Q3 | Mean | Sd | hist |
---|---|---|---|---|---|---|---|---|---|---|
edad_jefe | 0 | 1 | 16.00 | 108.00 | 37.00 | 48.00 | 60.00 | 49.31 | 16.16 | ▃▇▆▂▁ |
tot_integ | 0 | 1 | 1.00 | 16.00 | 3.00 | 4.00 | 5.00 | 3.88 | 1.94 | ▇▃▁▁▁ |
ing_cor | 0 | 1 | 3686.55 | 574207.55 | 18580.30 | 28921.53 | 44919.18 | 37536.71 | 35114.31 | ▇▁▁▁▁ |
alimentos | 0 | 1 | 0.00 | 46246.78 | 5682.77 | 8858.38 | 12661.47 | 10029.62 | 6342.12 | ▇▆▁▁▁ |
vesti_calz | 0 | 1 | 0.00 | 18308.12 | 60.65 | 704.34 | 1764.76 | 1318.06 | 1887.21 | ▇▁▁▁▁ |
vivienda | 0 | 1 | 0.00 | 29650.64 | 867.00 | 1575.00 | 2527.83 | 2198.09 | 2488.74 | ▇▁▁▁▁ |
limpieza | 0 | 1 | 0.00 | 24138.67 | 555.00 | 953.29 | 1657.06 | 1500.72 | 1948.88 | ▇▁▁▁▁ |
salud | 0 | 1 | 0.00 | 68478.26 | 0.00 | 11.73 | 489.11 | 732.52 | 2912.94 | ▇▁▁▁▁ |
transporte | 0 | 1 | 0.00 | 86399.60 | 904.28 | 2755.86 | 5725.68 | 4356.88 | 5500.09 | ▇▁▁▁▁ |
educa_espa | 0 | 1 | 0.00 | 70215.32 | 0.00 | 952.24 | 3919.34 | 3467.23 | 6664.46 | ▇▁▁▁▁ |
personales | 0 | 1 | 0.00 | 41299.42 | 804.15 | 1425.87 | 2412.52 | 2085.21 | 2750.45 | ▇▁▁▁▁ |
Como podemos observar, esta tabla nos muestra los valores mínimos y máximos que toma nuestra variable; además, nos muestra sus cuartiles, y otra información estadística importante (media, mediana y varianza de los datos). En la última columna también podemos observar cómo se ven los datos en una gráfica de histograma.
Esta tabla nos ayudará con los parámetros de la función OutlierReplace que vendrá más adelante, pues nos pide como argumentos el primer y tercer cuartil de los datos de nuestra variable a corregir. Por mientras veamos el siguiente boxplot:
boxplot(data, main = "Visualización de Outliers", col = "palevioletred1", xlab = "Variables")
Dado el gráfico anterior, podemos observar que hay datos atípicos en todas las variables, así que corregiremos esa información. Para ello crearemos una función que reemplazará los outliers por la media de los datos (sin outliers), y nos devolverá la variable con los nuevos datos.
OutliersReplace <- function(variable, Q1, Q3){ #Q1 y Q3 son el primer y tercer cuartil obtenidos de la tabla anterior.
IQR <- Q3-Q1
lowLimit <- Q1-1.5*IQR
highLimit <- Q3+1.5*IQR
variable[variable < lowLimit] <- trunc(mean(variable))
variable[variable > highLimit] <- trunc(mean(variable))
variable
}
Aplicamos la función a nuestras diferentes variables y guardamos los datos sin atípicos a nuestra misma tablita data
data$edad_jefe <- OutliersReplace(data$edad_jefe,37,60)
data$tot_integ <- OutliersReplace(data$tot_integ,3,5)
data$ing_cor <- OutliersReplace(data$ing_cor, 18580, 44919)
data$alimentos <- OutliersReplace(data$alimentos, 5683, 12661)
data$vesti_calz <- OutliersReplace(data$vesti_calz, 60.65, 1764.76)
data$vivienda <- OutliersReplace(data$vivienda, 867, 2528)
data$limpieza <- OutliersReplace(data$limpieza, 555, 1657.1)
data$salud <- OutliersReplace(data$salud, 0, 489.11)
data$transporte <- OutliersReplace(data$transporte, 904.3, 5725.7)
data$educa_espa <- OutliersReplace(data$educa_espa, 0, 3919.3)
data$personales <- OutliersReplace(data$personales, 804.1, 2412.5)
De esta forma nuestra base de datos ya no tendrán datos atípicos y nuestros nuevos datos se resumen de la siguiente manera
ResEst(data)
Name | data |
Number of rows | 1837 |
Number of columns | 11 |
_______________________ | |
Column type frequency: | |
numeric | 11 |
________________________ | |
Group variables | None |
Variable type: numeric
skim_variable | n_missing | complete_rate | Min | Max | Q1 | Median | Q3 | Mean | Sd | hist |
---|---|---|---|---|---|---|---|---|---|---|
edad_jefe | 0 | 1 | 16.00 | 94.00 | 37.00 | 48.00 | 60.00 | 49.20 | 15.97 | ▃▇▆▃▁ |
tot_integ | 0 | 1 | 1.00 | 8.00 | 3.00 | 4.00 | 5.00 | 3.71 | 1.63 | ▅▅▇▂▁ |
ing_cor | 0 | 1 | 3686.55 | 83903.29 | 18580.30 | 28921.53 | 39700.84 | 31638.85 | 16970.53 | ▇▇▆▂▁ |
alimentos | 0 | 1 | 0.00 | 22975.51 | 5682.77 | 8858.38 | 11957.12 | 9242.16 | 4768.29 | ▃▇▇▂▁ |
vesti_calz | 0 | 1 | 0.00 | 4304.31 | 60.65 | 704.34 | 1408.66 | 982.38 | 1014.19 | ▇▃▂▁▁ |
vivienda | 0 | 1 | 0.00 | 4968.00 | 867.00 | 1575.00 | 2198.00 | 1661.37 | 1050.09 | ▆▇▆▂▁ |
limpieza | 0 | 1 | 0.00 | 3309.63 | 555.00 | 953.29 | 1500.00 | 1064.78 | 664.83 | ▇▇▆▂▁ |
salud | 0 | 1 | 0.00 | 1209.10 | 0.00 | 11.73 | 489.11 | 236.06 | 319.05 | ▇▁▁▂▁ |
transporte | 0 | 1 | 0.00 | 12916.17 | 904.28 | 2755.86 | 4814.04 | 3418.54 | 2979.86 | ▇▅▂▁▁ |
educa_espa | 0 | 1 | 0.00 | 9580.63 | 0.00 | 952.24 | 3467.00 | 1942.90 | 2300.57 | ▇▃▁▁▁ |
personales | 0 | 1 | 0.00 | 4790.27 | 804.15 | 1425.87 | 2085.00 | 1581.80 | 1008.75 | ▇▇▆▂▁ |
Observemos que tienen datos un poco más centrados a la media, sin embargo más adelante haremos un contraste con nuestros datos sin ser corregidos para ver qué tanto peso aportan esos Outliers.
Por mientras, démosle paso a ver qué tan relacionadas están nuestras varibales respuesta con nuestra variable regresora (Ingresos corrientes). Para ello echaremos un primer vistazo a la siguiente gráfica que nos muestra la correlación de Spearman que existe entre nuestras variables.
Recordemos que la variable controlada en este proyecto de regresión lineal, son los Ingresos Corrientes, y queremos ver qué tanto influye esta en las variables restantes.
library(corrplot)
library(GGally)
library(Hmisc)
library(PerformanceAnalytics)
corrplot(cor(data, method = "spearman"), method = "ellipse", type = "lower", diag = F, tl.col="mediumvioletred", tl.cex=0.8, tl.srt=70)
Dicho correlograma nos brinda información de que la variable Alimentos, Limpieza, Transporte y Gastos personales se correlacionan fuertemente con los ingresos corrientes. Es decir, si el ingreso corriente aumenta, entonces los gastos de cada una de las variables mencionadas (en realidad todas, pero por ahora nos estamos centrando en las que tienen mayor correlación) también lo hacen. Pero observemos también las diferentes gráficas de dispersión de estas variables con relación a los Ingresos Corrientes:
library(dplyr)
library(ggplot2)
library(gridExtra)
gEdad <- data %>% ggplot(aes(x = data$ing_cor , y = data$edad_jefe)) + labs(x = "Ingreso corriente", y = "Edad Jefe o Jefa") + geom_point(colour = "palevioletred1")
gInt <- data %>% ggplot(aes(x = data$ing_cor , y = data$tot_integ)) + labs(x = "Ingreso corriente", y = "Total de integrantes") + geom_point(colour = "palevioletred1")
gAlim<- data %>% ggplot(aes(x = data$ing_cor , y = data$alimentos)) + labs(x = "Ingreso corriente", y = "Egresos en alimentos") + geom_point(colour = "palevioletred1")
gVesyCal <- data %>% ggplot(aes(x = data$ing_cor , y = data$vesti_calz)) + labs(x = "Ingreso corriente", y = "Egresos en Vestimenta y Calzado") + geom_point(colour = "palevioletred1")
gViv <- data %>% ggplot(aes(x = data$ing_cor , y = data$vivienda)) + labs(x = "Ingreso corriente", y = "Egresos en Vivienda") + geom_point(colour = "palevioletred1")
gLim <- data %>% ggplot(aes(x = data$ing_cor , y = data$limpieza)) + labs(x = "Ingreso corriente", y = "Egresos en Limpieza") + geom_point(colour = "palevioletred1")
gSal <- data %>% ggplot(aes(x = data$ing_cor , y = data$salud)) + labs(x = "Ingreso corriente", y = "Egresos en Salud") + geom_point(colour = "palevioletred1")
gTrans <- data %>% ggplot(aes(x = data$ing_cor , y = data$transporte)) + labs(x = "Ingreso corriente", y = "Egresos en Transporte") + geom_point(colour = "palevioletred1")
gEdu <- data %>% ggplot(aes(x = data$ing_cor , y = data$educa_espa)) + labs(x = "Ingreso corriente", y = "Egresos en Educación y Esparcimiento") + geom_point(colour = "palevioletred1")
gPer <- data %>% ggplot(aes(x = data$ing_cor , y = data$personales)) + labs(x = "Ingreso corriente", y = "Gastos Personales") + geom_point(colour = "palevioletred1")
grid.arrange(gEdad, gAlim, gVesyCal, gViv, gInt, nrow = 2)
grid.arrange(gLim, gSal, gTrans, gEdu, gPer, nrow = 2)
Dadas estas gráficas no es sencillo visualizar una regresión lineal, por eso mismo aplicaremos una transformación conveniente para reducir la dispersión y dar valores pequeños. Nos enfocaremos por ahora en la variable respuesta Alimentos.
library(purrr)
ModAlim1 <- lm(ing_cor~alimentos,data=data)
ModAlim2 <- lm(sqrt(ing_cor)~alimentos,data=data)
ModAlim3 <- lm(log(ing_cor)~alimentos,data=data)
Y antes de continuar, veamos cuántos ceros tiene cada una de nuestras variables, ya que nos causarán problemas al momento de ocupar una transformación con logaritmo natural:
zeros<-colnames(data)
casos=0
for (i in 1:length(zeros)){
casos[i]=length(which(data[,i]==0))
}
cbind(zeros,casos)
## zeros casos
## [1,] "edad_jefe" "0"
## [2,] "tot_integ" "0"
## [3,] "ing_cor" "0"
## [4,] "alimentos" "9"
## [5,] "vesti_calz" "433"
## [6,] "vivienda" "26"
## [7,] "limpieza" "24"
## [8,] "salud" "911"
## [9,] "transporte" "112"
## [10,] "educa_espa" "554"
## [11,] "personales" "12"
Como la variable Alimentos tiene nueve ceros (que son pocos a comparación del número de observaciones), y como entendemos que esta variable es una necesidad básica en la población que por ningún motivo podría ser cero, eliminaremos estas nueve observaciones:
AlimSinZeros <- data%>%filter(data$alimentos!=0)
Una vez eliminados los ceros, continuamos con la propuesta de nuestro modelo 4 y enseguida el resumen de nuestros modelos para determinar cuál es el más conventiete basándonos en el coeficiente de determinación:
ModAlim4 <- lm(ing_cor~log(alimentos),data=AlimSinZeros)
ModelsAlimento <- list(ModAlim1, ModAlim2, ModAlim3, ModAlim4)
library(stargazer)
stargazer(ModelsAlimento, type="text",title = "Modelos propuestos de RLS: Ingresos Corrientes ~ Alimentos", style = "ajs",align = FALSE, no.space=FALSE)
##
## Modelos propuestos de RLS: Ingresos Corrientes ~ Alimentos
##
## ===========================================================================================================================
## ING_COR SQRT(ING_COR) LOG(ING_COR) ING_COR
## 1 2 3 4
## ---------------------------------------------------------------------------------------------------------------------------
## alimentos 1.730*** .005*** 0.000***
## (.073) (0.000) (0.000)
##
## log(alimentos) 13,170.660***
## (572.091)
##
## Constant 15,647.430*** 124.829*** 9.626*** -86,584.500***
## (755.018) (2.093) (.026) (5,150.268)
##
## Observations 1,837 1,837 1,837 1,828
## R2 .236 .255 .254 .225
## Adjusted R2 .236 .254 .254 .225
## Residual Std. Error 14,834.050 (df = 1835) 41.119 (df = 1835) .511 (df = 1835) 14,949.420 (df = 1826)
## F Statistic 567.944*** (df = 1; 1835) 626.530*** (df = 1; 1835) 626.312*** (df = 1; 1835) 530.011*** (df = 1; 1826)
## ---------------------------------------------------------------------------------------------------------------------------
## Notes: *P < .05
## **P < .01
## ***P < .001
Estos modelos se ven de la siguiente forma:
plotModAlim1 <- ggplot(data, aes(x= data$ing_cor, y= data$alimentos)) +
geom_point(col = "palevioletred1") +
geom_smooth(method='lm', formula=y~x, se=FALSE, col='yellow') +
theme_light()+
labs(x = "Ingresos Corrientes", y = "Gastos en Alimentos", title = "Modelo 1")
plotModAlim2 <- ggplot(data, aes(x= sqrt(data$ing_cor), y= data$alimentos)) +
geom_point(col = "palevioletred1") +
geom_smooth(method='lm', formula=y~x, se=FALSE, col='yellow') +
theme_light()+
labs(x = "Ingresos Corrientes", y = "Gastos en Alimentos", title = "Modelo 2")
plotModAlim3 <- ggplot(data, aes(x= log(data$ing_cor), y= data$alimentos)) +
geom_point(col = "palevioletred1") +
geom_smooth(method='lm', formula=y~x, se=FALSE, col='yellow') +
theme_light()+
labs(x = "Ingresos Corrientes", y = "Gastos en Alimentos", title = "Modelo 3")
plotModAlim4 <- ggplot(AlimSinZeros, aes(x= ing_cor, y= log(alimentos))) +
geom_point(col = "palevioletred1") +
geom_smooth(method='lm', formula=y~x, se=FALSE, col='yellow') +
theme_light()+
labs(x = "Ingresos Corrientes", y = "Gastos en Alimentos", title = "Modelo 4")
grid.arrange(plotModAlim1, plotModAlim2, plotModAlim3, plotModAlim4, nrow = 2)
Por lo tanto, el modelo que más se ajusta es el Modelo 2 y el Modelo 3 (vemos que el Coeficiente de Determinación es el mismo para ambos modelos y mayor que el Modelo 1 y 4). Es decir, el 25% de la variabilidad de nuestra variable Alimentos es predicha por la variable Ingresos Corrientes. Se confieza que este porcentaje es muy bajo a comparación de lo que nos gustaría. Por eso mismo hagamos exactamente el mismo proceso para llegar a estos modelos pero con mi base de datos original, sin quitarle los datos atípicos.
DaTa<-data.frame(read_excel("C://Users//Nath//Documents//Semestre 2021-2//Mnpyr//Proyecto//Puebla.xlsx"))
DaTa$entidad <- NULL
MODELAlim1 <- lm(ing_cor~alimentos,data=DaTa)
MODELAlim2 <- lm(sqrt(ing_cor)~alimentos,data=DaTa)
MODELAlim3 <- lm(log(ing_cor)~alimentos,data=DaTa)
ReplaceZerosinAlim <- DaTa%>%filter(DaTa$alimentos!=0)
MODELAlim4 <- lm(ing_cor~log(alimentos),data=ReplaceZerosinAlim )
MODELSAlimento <- list(MODELAlim1, MODELAlim2, MODELAlim3 , MODELAlim4)
library(stargazer)
stargazer(ModelsAlimento, type="text",title = "Modelos propuestos de RLS con Outliers: Ingresos Corrientes ~ Alimentos", style = "ajs",align = FALSE, no.space=FALSE)
##
## Modelos propuestos de RLS con Outliers: Ingresos Corrientes ~ Alimentos
##
## ===========================================================================================================================
## ING_COR SQRT(ING_COR) LOG(ING_COR) ING_COR
## 1 2 3 4
## ---------------------------------------------------------------------------------------------------------------------------
## alimentos 1.730*** .005*** 0.000***
## (.073) (0.000) (0.000)
##
## log(alimentos) 13,170.660***
## (572.091)
##
## Constant 15,647.430*** 124.829*** 9.626*** -86,584.500***
## (755.018) (2.093) (.026) (5,150.268)
##
## Observations 1,837 1,837 1,837 1,828
## R2 .236 .255 .254 .225
## Adjusted R2 .236 .254 .254 .225
## Residual Std. Error 14,834.050 (df = 1835) 41.119 (df = 1835) .511 (df = 1835) 14,949.420 (df = 1826)
## F Statistic 567.944*** (df = 1; 1835) 626.530*** (df = 1; 1835) 626.312*** (df = 1; 1835) 530.011*** (df = 1; 1826)
## ---------------------------------------------------------------------------------------------------------------------------
## Notes: *P < .05
## **P < .01
## ***P < .001
PLOTModAlim1 <- ggplot(DaTa, aes(x= ing_cor, y= alimentos)) +
geom_point(col = "palevioletred1") +
geom_smooth(method='lm', formula=y~x, se=FALSE, col='yellow') +
theme_light()+
labs(x = "Ingresos Corrientes", y = "Gastos en Alimentos", title = "MODEL 1")
PLOTModAlim2 <- ggplot(DaTa, aes(x= sqrt(ing_cor), y= alimentos)) +
geom_point(col = "palevioletred1") +
geom_smooth(method='lm', formula=y~x, se=FALSE, col='yellow') +
theme_light()+
labs(x = "Ingresos Corrientes", y = "Gastos en Alimentos", title = "MODEL 2")
PLOTModAlim3 <- ggplot(DaTa, aes(x= log(ing_cor), y= alimentos)) +
geom_point(col = "palevioletred1") +
geom_smooth(method='lm', formula=y~x, se=FALSE, col='yellow') +
theme_light()+
labs(x = "Ingresos Corrientes", y = "Gastos en Alimentos", title = "MODEL 3")
PLOTModAlim4 <- ggplot(ReplaceZerosinAlim, aes(x= ing_cor, y= log(alimentos))) +
geom_point(col = "palevioletred1") +
geom_smooth(method='lm', formula=y~x, se=FALSE, col='yellow') +
theme_light()+
labs(x = "Ingresos Corrientes", y = "Gastos en Alimentos", title = "MODEL 4")
grid.arrange(PLOTModAlim1, PLOTModAlim2, PLOTModAlim3, PLOTModAlim4, nrow = 2)
Llgamos a la misma conclusión; el modelo que más se ajusta es el MODEL 2 Y MODEL 3, pero de igual forma no tenemos un Coeficiente de Determinación mayor.
No haremos todo este proceso con cada una de las variables restantes, mejor visualizaremos una tabla de regresión para que dados los datos obtenidos, podamos elegir una variable que me grite “Ey! elígeme a mí para ser estudiada!”
ModeloSimple1 <- lm(ing_cor~edad_jefe,data=data)
ModeloSimple2 <- lm(ing_cor~tot_integ,data=data)
ModeloSimple3 <- lm(ing_cor~alimentos,data=data)
ModeloSimple4 <- lm(ing_cor~vesti_calz,data=data)
ModeloSimple5 <- lm(ing_cor~vivienda,data=data)
ModeloSimple6 <- lm(ing_cor~limpieza,data=data)
ModeloSimple7 <- lm(ing_cor~salud,data=data)
ModeloSimple8 <- lm(ing_cor~transporte,data=data)
ModeloSimple9 <- lm(ing_cor~educa_espa,data=data)
ModeloSimple10 <- lm(ing_cor~personales,data=data)
MRLS <- list(ModeloSimple1, ModeloSimple2, ModeloSimple3, ModeloSimple4, ModeloSimple5, ModeloSimple6, ModeloSimple7, ModeloSimple8, ModeloSimple9, ModeloSimple10)
stargazer(MRLS , type="text",title = "Modelos propuestos de RLS SIN Outliers", style = "ajps",align = FALSE, no.space=FALSE)
##
## Modelos propuestos de RLS SIN Outliers
## -----------------------------------------------------------------------------------------------------------------------------------------------------------------
## ing_cor
## Model 1 Model 2 Model 3 Model 4 Model 5 Model 6 Model 7 Model 8 Model 9 Model 10
## -----------------------------------------------------------------------------------------------------------------------------------------------------------------
## edad_jefe -36.225
## (24.786)
## tot_integ 2843.122***
## (234.100)
## alimentos 1.730***
## (0.073)
## vesti_calz 5.948***
## (0.365)
## vivienda 4.538***
## (0.362)
## limpieza 10.560***
## (0.543)
## salud 10.757***
## (1.216)
## transporte 3.221***
## (0.110)
## educa_espa 2.039***
## (0.165)
## personales 7.450***
## (0.352)
## Constant 33421.030*** 21100.570*** 15647.430*** 25795.670*** 24100.110*** 20394.440*** 29099.640*** 20627.210*** 27677.560*** 19855.130***
## (1282.050) (947.689) (755.018) (515.452) (711.603) (680.955) (482.549) (497.153) (498.265) (660.567)
## N 1837 1837 1837 1837 1837 1837 1837 1837 1837 1837
## R-squared 0.001 0.074 0.236 0.126 0.079 0.171 0.041 0.320 0.076 0.196
## Adj. R-squared 0.001 0.074 0.236 0.126 0.078 0.171 0.040 0.320 0.076 0.196
## Residual Std. Error (df = 1835) 16965.280 16331.470 14834.050 15866.500 16292.280 15454.380 16624.430 13999.040 16313.890 15220.170
## F Statistic (df = 1; 1835) 2.136 147.499*** 567.944*** 265.395*** 157.047*** 378.911*** 78.242*** 863.155*** 151.774*** 447.570***
## -----------------------------------------------------------------------------------------------------------------------------------------------------------------
## ***p < .01; **p < .05; *p < .1
ModSimple1 <- lm(ing_cor~edad_jefe,data=DaTa)
ModSimple2 <- lm(ing_cor~tot_integ,data=DaTa)
ModSimple3 <- lm(ing_cor~alimentos,data=DaTa)
ModSimple4 <- lm(ing_cor~vesti_calz,data=DaTa)
ModSimple5 <- lm(ing_cor~vivienda,data=DaTa)
ModSimple6 <- lm(ing_cor~limpieza,data=DaTa)
ModSimple7 <- lm(ing_cor~salud,data=DaTa)
ModSimple8 <- lm(ing_cor~transporte,data=DaTa)
ModSimple9 <- lm(ing_cor~educa_espa,data=DaTa)
ModSimple10 <- lm(ing_cor~personales,data=DaTa)
MRLSOutliers <- list(ModSimple1, ModSimple2, ModSimple3, ModSimple4, ModSimple5, ModSimple6, ModSimple7, ModSimple8, ModSimple9, ModSimple10)
stargazer(MRLSOutliers, type="text",title = "Modelos propuestos de CON Outliers", style = "ajps",align = FALSE, no.space=FALSE)
##
## Modelos propuestos de CON Outliers
## ----------------------------------------------------------------------------------------------------------------------------------------------------------------
## ing_cor
## Model 1 Model 2 Model 3 Model 4 Model 5 Model 6 Model 7 Model 8 Model 9 Model 10
## ----------------------------------------------------------------------------------------------------------------------------------------------------------------
## edad_jefe -28.601
## (50.726)
## tot_integ 3714.811***
## (414.182)
## alimentos 2.951***
## (0.109)
## vesti_calz 7.579***
## (0.397)
## vivienda 3.816***
## (0.317)
## limpieza 8.572***
## (0.370)
## salud 1.780***
## (0.278)
## transporte 3.591***
## (0.123)
## educa_espa 1.998***
## (0.114)
## personales 5.879***
## (0.265)
## Constant 38947.060*** 23106.170*** 7939.625*** 27546.660*** 29147.950*** 24672.390*** 36232.640*** 21893.090*** 30609.090*** 25277.770***
## (2632.121) (1797.782) (1297.651) (912.972) (1052.736) (909.846) (835.760) (864.545) (854.823) (912.953)
## N 1837 1837 1837 1837 1837 1837 1837 1837 1837 1837
## R-squared 0.0002 0.042 0.284 0.166 0.073 0.226 0.022 0.316 0.144 0.212
## Adj. R-squared -0.0004 0.041 0.284 0.165 0.073 0.226 0.021 0.316 0.143 0.212
## Residual Std. Error (df = 1835) 35120.830 34378.410 29719.180 32077.650 33814.580 30894.080 34738.730 29042.660 32500.460 31178.210
## F Statistic (df = 1; 1835) 0.318 80.444*** 728.110*** 365.066*** 144.853*** 536.868*** 40.914*** 848.910*** 308.195*** 493.835***
## ----------------------------------------------------------------------------------------------------------------------------------------------------------------
## ***p < .01; **p < .05; *p < .1
Observamos que de ambas tablas las variables que tienen un mayor Coeficiente de determinación son Alimentos, Transporte, Limpieza y gastos Personales. Y tambén démonos cuenta que DaTa (base con Outliers) nos aporta “mejores” resultados en cuanto a este Coeficiente. Además, dado que el p-value en la variable Edad es mayor a nuestro nivel de significancia \(\alpha = 0.05\), nos está indicando que no es significativa con respecto a los Ingresos Corrientes.
Para mejorar nuestros modelos anteriormente propuestos, agregaremos más variables, es decir, ahora supondremos n predictores. Para ello nuevamente echemos ojo a un correlograma. Este nos ayudará a ver qué variables se relacionan monótamente a nuestra variable Ingreso corriente, y así poder elegir las variables que parezcan ser más significativas en nuestros próximos modelos.
Ocuparemos nuestra base de datos D1 (la base con los Outliers)
D1 <- data.frame(read_excel("C://Users//Nath//Documents//Semestre 2021-2//Mnpyr//Proyecto//Puebla.xlsx"))
D1$entidad <- NULL
corrplot(cor(D1, method = "spearman"), method = "square", type = "lower", diag = F, tl.col="mediumvioletred", tl.cex=0.8, tl.srt=70)
Como ya lo vimos anteriormente, las variables que tienen mayor correlación de Spearman con Ingresos Corrientes son Alimentos, Transporte, Limpieza y Personales. Seguramente estas variables conformarán nuestro modelo de regresión lineal que más se ajuste a nuestros datos, sin embargo lleguemos a ello o desmintámoslo paso a paso:
En principio, lineas atrás llegamos a que la variable Edad no es significativa con respecto a los Ingresos Corrientes, entonces podemos eliminarla, al igual que la variable Total de integrantes,
Propongamos algunos modelos para elegir el que más se ajuste a nuestra regresión.
D1$edad_jefe <- NULL
D1$tot_integ <- NULL
all_model <- lm(ing_cor~., data= D1)
stargazer(all_model , type="text",title = "Modelos propuestos de RLS", style = "ajps",align = FALSE, no.space=FALSE)
##
## Modelos propuestos de RLS
## ---------------------------------------------
## ing_cor
## ---------------------------------------------
## alimentos 1.093***
## (0.123)
## vesti_calz 0.938**
## (0.392)
## vivienda 0.918***
## (0.256)
## limpieza 2.897***
## (0.390)
## salud 0.635***
## (0.210)
## transporte 1.794***
## (0.139)
## educa_espa 0.594***
## (0.101)
## personales 0.898***
## (0.287)
## Constant 6757.575***
## (1164.253)
## N 1837
## R-squared 0.461
## Adj. R-squared 0.458
## Residual Std. Error 25840.180 (df = 1828)
## F Statistic 195.299*** (df = 8; 1828)
## ---------------------------------------------
## ***p < .01; **p < .05; *p < .1
En esta última tabla podemos ver que todas las variables tienen cierta significancia en nuestro modelo, pero veamos cúal será el modelo más óptimo:
Emplearemos el método Backward con el criterio AIC:
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
MBackwardAIC <- stepAIC(all_model, direction = "backward")
## Start: AIC=37335.66
## ing_cor ~ alimentos + vesti_calz + vivienda + limpieza + salud +
## transporte + educa_espa + personales
##
## Df Sum of Sq RSS AIC
## <none> 1.2206e+12 37336
## - vesti_calz 1 3.8276e+09 1.2244e+12 37339
## - salud 1 6.0826e+09 1.2267e+12 37343
## - personales 1 6.5488e+09 1.2271e+12 37343
## - vivienda 1 8.5729e+09 1.2292e+12 37347
## - educa_espa 1 2.2914e+10 1.2435e+12 37368
## - limpieza 1 3.6905e+10 1.2575e+12 37388
## - alimentos 1 5.2678e+10 1.2733e+12 37411
## - transporte 1 1.1162e+11 1.3322e+12 37494
summary(MBackwardAIC)
##
## Call:
## lm(formula = ing_cor ~ alimentos + vesti_calz + vivienda + limpieza +
## salud + transporte + educa_espa + personales, data = D1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -205781 -9510 -2979 5519 460267
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6757.5748 1164.2529 5.804 7.61e-09 ***
## alimentos 1.0931 0.1231 8.882 < 2e-16 ***
## vesti_calz 0.9379 0.3917 2.394 0.016755 *
## vivienda 0.9180 0.2562 3.583 0.000348 ***
## limpieza 2.8968 0.3896 7.434 1.61e-13 ***
## salud 0.6351 0.2104 3.018 0.002578 **
## transporte 1.7944 0.1388 12.930 < 2e-16 ***
## educa_espa 0.5940 0.1014 5.858 5.54e-09 ***
## personales 0.8977 0.2867 3.132 0.001765 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 25840 on 1828 degrees of freedom
## Multiple R-squared: 0.4608, Adjusted R-squared: 0.4585
## F-statistic: 195.3 on 8 and 1828 DF, p-value: < 2.2e-16
Según esta función, el mejor modelo contiene todas las variables propuestas, pero más adelante lo confirmaremos probando que cumpla los supuestos. Por mientras, veamos más modelos posibles usando otros métodos.
Con este método, utilizando la función regsubsets de R, que busca el mejor modelo que se ajuste a nuestra regresión tomando como criterio el RSS (suma de cuadrados de los residuos):
library(leaps)
bestmodel <- regsubsets(ing_cor~., data=D1, method="forward",nvmax=8)
MejoresModelosFordward <- summary(bestmodel)
MejoresModelosFordward
## Subset selection object
## Call: regsubsets.formula(ing_cor ~ ., data = D1, method = "forward",
## nvmax = 8)
## 8 Variables (and intercept)
## Forced in Forced out
## alimentos FALSE FALSE
## vesti_calz FALSE FALSE
## vivienda FALSE FALSE
## limpieza FALSE FALSE
## salud FALSE FALSE
## transporte FALSE FALSE
## educa_espa FALSE FALSE
## personales FALSE FALSE
## 1 subsets of each size up to 8
## Selection Algorithm: forward
## alimentos vesti_calz vivienda limpieza salud transporte educa_espa
## 1 ( 1 ) " " " " " " " " " " "*" " "
## 2 ( 1 ) "*" " " " " " " " " "*" " "
## 3 ( 1 ) "*" " " " " "*" " " "*" " "
## 4 ( 1 ) "*" " " " " "*" " " "*" "*"
## 5 ( 1 ) "*" " " " " "*" " " "*" "*"
## 6 ( 1 ) "*" " " "*" "*" " " "*" "*"
## 7 ( 1 ) "*" " " "*" "*" "*" "*" "*"
## 8 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## personales
## 1 ( 1 ) " "
## 2 ( 1 ) " "
## 3 ( 1 ) " "
## 4 ( 1 ) " "
## 5 ( 1 ) "*"
## 6 ( 1 ) "*"
## 7 ( 1 ) "*"
## 8 ( 1 ) "*"
Pero ahora elegiremos el mejor de acuerdo al BIC, \(C_{p}\) Y \(R^2\) ajustado:
ModelBIC <- ggplot(data=data.frame(n_predictores=1:8, BIC=MejoresModelosFordward
$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 = "palevioletred1", size = 4) +
scale_x_continuous(breaks = c(0:19))+
theme_bw() +
labs(title = "BAJO CIRTERIO BIC",
x = "Número Predictores")
ModelCp <- ggplot(data = data.frame(n_predictores = 1:8, Cp= MejoresModelosFordward$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 = "palevioletred1", size = 4) +
scale_x_continuous(breaks = c(0:19))+
theme_bw() +
labs(title = " BAJO CRITERIO Cp",
x = "Número Predictores")
ModelR2 <- ggplot(data = data.frame(n_predictores = 1:8, R2_ajustado= MejoresModelosFordward$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 = "palevioletred1", size = 4) +
scale_x_continuous(breaks = c(0:19))+
theme_bw() +
labs(title = " BAJO CRITERIO R2 AJUSTADO",
x = "Número Predictores")
grid.arrange(ModelBIC, ModelCp, ModelR2, nrow = 2)
Dados estos gráficos, bajo el criterio \(R^2\) y \(C_{p}\) utilizando el método Forward, nuestro mejor modelo es el que tiene las ocho variables predictoras que habíamos propuesto, de hecho coincide con el método Backward. Bajo el criterio BIC, parece ser que el mejor modelo es el que tiene 7 variables predictorias. Más adelante lo compararemos. Por mientras estudiemos el de ocho variables.
Primeramente observamos las diferentes gráficas de dispersión de estas variables predictorias con respecto a los Ingresos Corrientes para ver cuáles parecen tener una relación lineal y elegir a las que le aplicaremos alguna transformación:
dis1 <- (D1 %>%
dplyr::select(ing_cor, alimentos, educa_espa, limpieza, personales) %>%
GGally::ggscatmat() +
theme(axis.text = element_text(size = 8),
axis.text.x = element_text(angle = 45)))
dis2 <- (D1 %>%
dplyr::select(ing_cor, salud, transporte, vesti_calz, vivienda) %>%
GGally::ggscatmat() +
theme(axis.text = element_text(size = 8),
axis.text.x = element_text(angle = 45)))
grid.arrange(dis1, dis2, nrow =1)
Aplicaremos diferentes transformaciones a nuestras variables y filtraremos los datos
Disp1 <- (D1 %>% filter(ing_cor < 500000 & alimentos < 40000 & educa_espa < 40000 & limpieza < 20000
& personales < 40000) %>%
dplyr::select(ing_cor, alimentos, educa_espa, limpieza, personales) %>%
mutate(ing_cor = log(ing_cor),
alimentos = sqrt(alimentos),
educa_espa = sqrt(educa_espa),
limpieza = sqrt(limpieza),
personales = sqrt(personales)) %>%
GGally::ggscatmat() +
theme(axis.text = element_text(size = 8),
axis.text.x = element_text(angle = 45)))
Disp2 <- (D1 %>% filter(ing_cor < 500000 & salud < 50000 & transporte < 40000 & vesti_calz < 15000 & vivienda < 25000 ) %>%
dplyr::select(ing_cor, salud, transporte, vesti_calz, vivienda) %>%
mutate(ing_cor = log(ing_cor),
salud = sqrt(salud),
transporte = sqrt(transporte),
vesti_calz = sqrt(vesti_calz),
vivienda = sqrt(vivienda)) %>%
GGally::ggscatmat() +
theme(axis.text = element_text(size = 8),
axis.text.x = element_text(angle = 45)))
grid.arrange(Disp1, Disp2, nrow =1)
Pudimos elevar el coeficiente de correlación de Pearson con estas transformaciones! (:
Estamos listxs para proponer un primer modelo para ser postulado como el que mejor se ajusta a nuestros datos:
mutated_D1 <- D1 %>% filter(ing_cor < 500000 & alimentos < 40000 & educa_espa < 40000 & limpieza < 20000 &
personales < 40000 & salud < 50000 & transporte < 40000 & vesti_calz < 15000 &
vivienda < 25000 )
first_model <- mutated_D1 %>% lm(log(ing_cor)~sqrt(alimentos) + sqrt(educa_espa) + sqrt(limpieza) + sqrt(personales)+
sqrt(salud)+sqrt(transporte)+sqrt(vesti_calz)+ sqrt(vivienda), data = .)
first_model %>% summary()
##
## Call:
## lm(formula = log(ing_cor) ~ sqrt(alimentos) + sqrt(educa_espa) +
## sqrt(limpieza) + sqrt(personales) + sqrt(salud) + sqrt(transporte) +
## sqrt(vesti_calz) + sqrt(vivienda), data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.96454 -0.25462 -0.00299 0.27675 1.89511
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.9296632 0.0383084 233.099 < 2e-16 ***
## sqrt(alimentos) 0.0035099 0.0004990 7.034 2.85e-12 ***
## sqrt(educa_espa) 0.0013462 0.0003043 4.423 1.03e-05 ***
## sqrt(limpieza) 0.0043711 0.0008157 5.358 9.47e-08 ***
## sqrt(personales) 0.0045922 0.0007563 6.072 1.54e-09 ***
## sqrt(salud) 0.0014658 0.0005148 2.848 0.00446 **
## sqrt(transporte) 0.0076760 0.0004112 18.670 < 2e-16 ***
## sqrt(vesti_calz) 0.0016120 0.0005702 2.827 0.00475 **
## sqrt(vivienda) 0.0030140 0.0005651 5.334 1.08e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.449 on 1800 degrees of freedom
## Multiple R-squared: 0.5631, Adjusted R-squared: 0.5612
## F-statistic: 290 on 8 and 1800 DF, p-value: < 2.2e-16
Vale, pero vemos que la variable Salud y Vestimenta y Calzado no son tan significantes como las demás, veamos qué ocurre si las eliminamos de nuestro modelo:
secund_model <- mutated_D1 %>% lm(log(ing_cor)~sqrt(alimentos) + sqrt(educa_espa) + sqrt(limpieza) + sqrt(personales)+
sqrt(transporte)+ sqrt(vivienda), data = .)
secund_model %>% summary()
##
## Call:
## lm(formula = log(ing_cor) ~ sqrt(alimentos) + sqrt(educa_espa) +
## sqrt(limpieza) + sqrt(personales) + sqrt(transporte) + sqrt(vivienda),
## data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9185 -0.2590 -0.0010 0.2753 1.8221
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.9158248 0.0382904 232.847 < 2e-16 ***
## sqrt(alimentos) 0.0038242 0.0004935 7.749 1.53e-14 ***
## sqrt(educa_espa) 0.0014903 0.0002998 4.971 7.28e-07 ***
## sqrt(limpieza) 0.0048546 0.0008101 5.993 2.49e-09 ***
## sqrt(personales) 0.0051032 0.0007493 6.811 1.32e-11 ***
## sqrt(transporte) 0.0078064 0.0004117 18.962 < 2e-16 ***
## sqrt(vivienda) 0.0029965 0.0005674 5.281 1.44e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4509 on 1802 degrees of freedom
## Multiple R-squared: 0.5589, Adjusted R-squared: 0.5575
## F-statistic: 380.6 on 6 and 1802 DF, p-value: < 2.2e-16
Vemos que aunque eliminamos las variables con menor significancia, nuestra \(R^2\) ajustada disminuyó. Entonces en este caso nos quedaremos con nuestro primer modelo first_model:
\(ln(ing_cor) = \beta_0 + \beta_1 \sqrt {alimetos} + \beta_2 \sqrt {educa_espa} + \beta_3 \sqrt {limpieza} + \beta_4 \sqrt {personales} + \beta_5 \sqrt {salud} + \beta_6 \sqrt {transporte} + \beta_7 \sqrt {vesti_calz} + \beta_8 \sqrt {vivienda}+ \epsilon\)
Este modelo tiene una \(R^2=0.5612\) ajustada, por lo que el modelo captura un 56.12% de la variabilidad de los datos (un buen porcentaje). Pero sigamos comprobando o desmintiendo que es un buen modelo validando los supuestos de los modelos de Regresión Lineal Múltiple:
Intervalos de Confianza
Obtengamos los intervalos de confianza en los que se encuentran cada uno de nuestros parámetros de las distintas variables
confint(first_model)
## 2.5 % 97.5 %
## (Intercept) 8.8545295260 9.004796826
## sqrt(alimentos) 0.0025311509 0.004488550
## sqrt(educa_espa) 0.0007492702 0.001943046
## sqrt(limpieza) 0.0027712417 0.005971036
## sqrt(personales) 0.0031088291 0.006075528
## sqrt(salud) 0.0004561863 0.002475333
## sqrt(transporte) 0.0068696409 0.008482415
## sqrt(vesti_calz) 0.0004936834 0.002730338
## sqrt(vivienda) 0.0019057823 0.004122251
library(modelsummary)
library(GGally)
library(broom.helpers)
GGally::ggcoef_model(first_model)
Prueba ANOVA (Análisis de Varianzas)
anova(first_model)
## Analysis of Variance Table
##
## Response: log(ing_cor)
## Df Sum Sq Mean Sq F value Pr(>F)
## sqrt(alimentos) 1 284.51 284.510 1411.3958 < 2.2e-16 ***
## sqrt(educa_espa) 1 27.72 27.721 137.5188 < 2.2e-16 ***
## sqrt(limpieza) 1 42.50 42.499 210.8300 < 2.2e-16 ***
## sqrt(personales) 1 24.59 24.590 121.9875 < 2.2e-16 ***
## sqrt(salud) 1 3.29 3.286 16.3022 5.627e-05 ***
## sqrt(transporte) 1 77.80 77.804 385.9678 < 2.2e-16 ***
## sqrt(vesti_calz) 1 1.60 1.599 7.9344 0.004903 **
## sqrt(vivienda) 1 5.74 5.735 28.4517 1.082e-07 ***
## Residuals 1800 362.84 0.202
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Comprobamos que todas nuestras variables son significativas (?) <—– REVISAR
Regresemos a nuestro Método Forward bajo el criterio BIC, donde la función aseguraba que nuestro mejor modelo es con 7 variables predictoras;
MejoresModelosFordward
## Subset selection object
## Call: regsubsets.formula(ing_cor ~ ., data = D1, method = "forward",
## nvmax = 8)
## 8 Variables (and intercept)
## Forced in Forced out
## alimentos FALSE FALSE
## vesti_calz FALSE FALSE
## vivienda FALSE FALSE
## limpieza FALSE FALSE
## salud FALSE FALSE
## transporte FALSE FALSE
## educa_espa FALSE FALSE
## personales FALSE FALSE
## 1 subsets of each size up to 8
## Selection Algorithm: forward
## alimentos vesti_calz vivienda limpieza salud transporte educa_espa
## 1 ( 1 ) " " " " " " " " " " "*" " "
## 2 ( 1 ) "*" " " " " " " " " "*" " "
## 3 ( 1 ) "*" " " " " "*" " " "*" " "
## 4 ( 1 ) "*" " " " " "*" " " "*" "*"
## 5 ( 1 ) "*" " " " " "*" " " "*" "*"
## 6 ( 1 ) "*" " " "*" "*" " " "*" "*"
## 7 ( 1 ) "*" " " "*" "*" "*" "*" "*"
## 8 ( 1 ) "*" "*" "*" "*" "*" "*" "*"
## personales
## 1 ( 1 ) " "
## 2 ( 1 ) " "
## 3 ( 1 ) " "
## 4 ( 1 ) " "
## 5 ( 1 ) "*"
## 6 ( 1 ) "*"
## 7 ( 1 ) "*"
## 8 ( 1 ) "*"
Ahora nuestras variables predictoras son 7: Alimentos, Vivienda, Limpieza, Salud, Transporte, Educación y Esparcimiento y Personales.
Pasaremos a estudiar este nuevo modelo para luego compararlo con el primer modelo first_model propuesto.
mutated_D1 <- D1 %>% filter(ing_cor < 500000 & alimentos < 40000 & educa_espa < 40000 & limpieza < 20000 &
personales < 40000 & salud < 50000 & transporte < 40000 & vesti_calz < 15000 &
vivienda < 25000 )
second_best_model <- mutated_D1 %>% lm(log(ing_cor)~sqrt(alimentos) + sqrt(educa_espa) + sqrt(limpieza) + sqrt(personales)+
sqrt(salud)+sqrt(transporte)+ sqrt(vivienda), data = .)
second_best_model %>% summary()
##
## Call:
## lm(formula = log(ing_cor) ~ sqrt(alimentos) + sqrt(educa_espa) +
## sqrt(limpieza) + sqrt(personales) + sqrt(salud) + sqrt(transporte) +
## sqrt(vivienda), data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.93022 -0.26056 0.00263 0.27907 1.83987
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.9198929 0.0382262 233.345 < 2e-16 ***
## sqrt(alimentos) 0.0037436 0.0004931 7.592 5.01e-14 ***
## sqrt(educa_espa) 0.0015122 0.0002992 5.054 4.76e-07 ***
## sqrt(limpieza) 0.0046661 0.0008106 5.756 1.01e-08 ***
## sqrt(personales) 0.0048565 0.0007520 6.458 1.36e-10 ***
## sqrt(salud) 0.0015667 0.0005145 3.045 0.00236 **
## sqrt(transporte) 0.0077402 0.0004113 18.818 < 2e-16 ***
## sqrt(vivienda) 0.0030109 0.0005662 5.318 1.18e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4498 on 1801 degrees of freedom
## Multiple R-squared: 0.5612, Adjusted R-squared: 0.5595
## F-statistic: 329.1 on 7 and 1801 DF, p-value: < 2.2e-16
Este modelo está dado por:
\(ln(ing_cor) = \beta_0 + \beta_1 \sqrt {alimetos} + \beta_2 \sqrt {educa_espa} + \beta_3 \sqrt {limpieza} + \beta_4 \sqrt {personales} + \beta_5 \sqrt {salud} + \beta_6 \sqrt {transporte} + \beta_8 \sqrt {vivienda}+ \epsilon\)
Y comparamos ambos modelos:
library(tidyr)
library(tibble)
library(dplyr)
library(broom)
library(ggplot2)
first_model %>%
glance() %>% gather ("Estadística", "Modelo 1")%>%
inner_join(
second_best_model%>% glance() %>% gather ("Estadística", "Modelo 2"),
by = "Estadística"
)
## Warning: attributes are not identical across measure variables;
## they will be dropped
## Warning: attributes are not identical across measure variables;
## they will be dropped
## # A tibble: 12 x 3
## Estadística `Modelo 1` `Modelo 2`
## <chr> <dbl> <dbl>
## 1 r.squared 5.63e- 1 5.61e- 1
## 2 adj.r.squared 5.61e- 1 5.60e- 1
## 3 sigma 4.49e- 1 4.50e- 1
## 4 statistic 2.90e+ 2 3.29e+ 2
## 5 p.value 4.38e-317 1.25e-316
## 6 df 8 e+ 0 7 e+ 0
## 7 logLik -1.11e+ 3 -1.12e+ 3
## 8 AIC 2.25e+ 3 2.25e+ 3
## 9 BIC 2.30e+ 3 2.30e+ 3
## 10 deviance 3.63e+ 2 3.64e+ 2
## 11 df.residual 1.8 e+ 3 1.80e+ 3
## 12 nobs 1.81e+ 3 1.81e+ 3
Podemos ver que por peqeñas diferencias, el mejor modelo hasta ahora es el first_model, pues el coeficiente de determinación es mayor, la Loglik es más cercano a cero, el AIC y BIC son menores y el p-value es menor al segundo modelo.
Decimos que hasta ahora es el mejor porque falta comprobar supuestos y checar valores influyentes. Antes de pasar a eso consideremos otro modelo completamente nuevo que propusimos al terminar Regresión Lineal Simple; ahora nuestras variables predictorias serán Alimentos, Transporte, Limpieza y Personales, pues fueron las más correlacionadas con los IngresosCorrientes (De una vez aplicaremos tranformación):
third_model <- mutated_D1 %>% lm(log(ing_cor)~alimentos + sqrt(limpieza) + sqrt(personales)+sqrt(transporte),
data = .)
third_model %>% summary()
##
## Call:
## lm(formula = log(ing_cor) ~ alimentos + sqrt(limpieza) + sqrt(personales) +
## sqrt(transporte), data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.88181 -0.27033 0.00243 0.28593 1.91635
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.169e+00 2.788e-02 328.913 < 2e-16 ***
## alimentos 2.172e-05 2.389e-06 9.093 < 2e-16 ***
## sqrt(limpieza) 5.095e-03 8.225e-04 6.195 7.21e-10 ***
## sqrt(personales) 5.830e-03 7.554e-04 7.718 1.94e-14 ***
## sqrt(transporte) 8.551e-03 4.070e-04 21.010 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4583 on 1804 degrees of freedom
## Multiple R-squared: 0.5437, Adjusted R-squared: 0.5427
## F-statistic: 537.5 on 4 and 1804 DF, p-value: < 2.2e-16