Datos

Inicialmente, se carga la base ENUT2020 completa, guardantola en una tabla con el mismo nombre:

library(haven)
library("dplyr")
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
ENUT2020 <- read_dta("/Users/victorcifuentes/Downloads/ENUT2020.dta")

Pregunta de investigación

Partiendo de la base de datos, se plantea la siguiente pregunta de investigación:

¿Qué factores influyen en el valor de subsidios recibidos por las personas encuestadas?

Variables de investigación

Dependiente: Valor total de subsidios recibidos.

Independientes: Edad, ingreso, número de habitantes en el hogar, número de horas que trabaja en la semana, sexo.

Inicialmente, se construye una base llamada “df” que contiene las variables de interés en el estudio:

## Selección de columnas de interés

df <- ENUT2020[c("p1103s1a2", "p1103s2a2", "p1103s3a2","p6040", "p6500","p6541", 
                 "p6008", "p1151", "p6020")]

A continuación, se construyen las variables “subsidios recibidos” e “ingreso”: La primera, es el resultado de la suma de los tres subsidios que puede recibir una persona, mientras que la variable “ingreso” es la suma de los ingresos que reporta una persona como independiente o asalariado.

df$total_subsidios <- rowSums(df[ , c(1,2,3)], na.rm=TRUE) # Calculo de variable y

df$ingreso <- rowSums(df[ , c(5,6)], na.rm=TRUE) # cálculo de ingreso

Limpieza de datos

Dado que la pregunta de investigación contempla a las personas encuestadas que reciben algún valor por concepto de subsidio, se excluyen aquellos que no reciben subsidio alguno:

df <- subset(df, df$total_subsidios > 0)

Además, se renombran las variables para facilitar su comprensión y manejo:

# Renombrar variables

colnames(df)[4] <-"edad"
colnames(df)[7] <-"num_habitantes"
colnames(df)[8] <-"num_horas_trabajo"
colnames(df)[9] <-"sexo"

df <- df[c("total_subsidios", "edad", "ingreso", "num_habitantes", 
           "num_horas_trabajo", "sexo")]

Por último, se evalúa el tipo de dato al que corresponde cada variable y se ajusta, si es necesario:

lapply(df,class)
## $total_subsidios
## [1] "numeric"
## 
## $edad
## [1] "numeric"
## 
## $ingreso
## [1] "numeric"
## 
## $num_habitantes
## [1] "numeric"
## 
## $num_horas_trabajo
## [1] "numeric"
## 
## $sexo
## [1] "numeric"

Cabe resaltar que la variable sexo no es númerica, sino que corresponde a un factor, por ende se corrige el tipo de dato:

df=transform(df,
             sexo=factor(sexo,levels=1:2)  #definición de la variable sexo como factor
)

lapply(df,class)
## $total_subsidios
## [1] "numeric"
## 
## $edad
## [1] "numeric"
## 
## $ingreso
## [1] "numeric"
## 
## $num_habitantes
## [1] "numeric"
## 
## $num_horas_trabajo
## [1] "numeric"
## 
## $sexo
## [1] "factor"

Análisis de datos faltantes

Ya con la base organizada, se evalúa el porcentaje de datos faltantes para cada una de las variables:

## Porcentaje de datos faltantes
visdat::vis_miss(df)

Como resultado, se evidencia que la variable #num_horas_trabajo contiene un 70% de datos ausentes. Para darle manejo a esta situación y no renunciar a esta variable, se propone imputar el dato ausente por cero, esto quiere decir que las personas que no respondieron a esta pregunta, no tienen trabajo y por ende sus horas semanales dedicadas a esa labor debería ser cero.

df <- mutate_at(df, c("num_horas_trabajo"), ~replace(., is.na(.), 0))

Análisis descriptivo

summary(df)
##  total_subsidios        edad           ingreso         num_habitantes  
##  Min.   :     99   Min.   :  0.00   Min.   :       0   Min.   : 1.000  
##  1st Qu.:  60000   1st Qu.: 11.00   1st Qu.:       0   1st Qu.: 4.000  
##  Median : 115000   Median : 21.00   Median :       0   Median : 4.000  
##  Mean   : 144243   Mean   : 26.18   Mean   :  180865   Mean   : 4.714  
##  3rd Qu.: 180000   3rd Qu.: 39.00   3rd Qu.:  130000   3rd Qu.: 5.000  
##  Max.   :1800000   Max.   :101.00   Max.   :21000000   Max.   :13.000  
##  num_horas_trabajo sexo     
##  Min.   :  0.00    1: 9443  
##  1st Qu.:  0.00    2:10762  
##  Median :  0.00             
##  Mean   : 12.81             
##  3rd Qu.: 24.00             
##  Max.   :120.00
par(mfrow=c(1,3))
boxplot(df$total_subsidios, xlab = 'Valor de subsidio', main = 'Valor recibido por concepto de subsidio' )
boxplot(df$edad, xlab = 'Edad', main = 'Edad' )
boxplot(df$num_habitantes, xlab = 'Número de habitantes', main = 'Número de habitantes' )

par(mfrow=c(1,3))
boxplot(df$ingreso, xlab = 'Ingreso', main = 'Ingreso')
boxplot(df$num_horas_trabajo, xlab = 'Número de horas', main = 'Número de horas que dedica a su trabajo')
Sexo2<-table(df$sexo)
Porcentajes<- round(Sexo2/sum(Sexo2)*100)
Etiqueta<- paste(Porcentajes, "%", sep=" ")
pie(x= Sexo2, labels= Etiqueta, col = c("lightgreen","lightpink"), main= "Sexo")
legend("topright", c("Masculino","Femenino"), cex=0.7, fill= c("lightgreen","lightpink"))

El \(75\%\) de los encuestados beneficiados con subsidios, recibieron \(\$180.000\) o menos; además, la edad promedio es de 26 años y en sus hogares habita un promedio de 4 personas. El \(50\%\) de estos participantes no reciben ingresos por trabajo y de igual manera no dedican ni una sola hora a alguna labor, por lo que se puede concluir que al menos el 50% de los encuestados que reciben subsidios se encuentran desempleados. Cabe destacar que el \(53\%\) de los encuestados que reciben subsidios son mujeres.

Relación entre el valor recibido de subsidio y cada una de las variables independientes

Valor total de subsidios vs la Edad

#install.packages("ggplot2",repos = "http://cran.us.r-project.org")
library(ggplot2)

ggplot(df, aes(y=total_subsidios , x=edad))+
  geom_point()     # Edad

No se evidencia una relación lineal directa entre el valor total de subsidios vs la Edad, lo que índica que el valor del subsidio recibido no está afectado por la edad del beneficiado

Valor total de subsidios vs el ingreso

ggplot(df, aes(y=total_subsidios , x=ingreso))+
  geom_point()     # Ingreso

La mayor concentración de ingresos está alrededor del cero, sin embargo, no se evidencia una relación lineal fuerte entre el ingreso y el valor del subsidio, lo que índica que una persona que no recibe ingresos por trabajar puede recibir montos bajos o altos de subsidios.

Valor total de subsidios vs el número de habitantes en el hogar

ggplot(df, aes(y=total_subsidios , x=num_habitantes))+
  geom_point()     # Num_habitantes

El número de habitantes en el hogar de los beneficiados con subsidios, no tiene un. impacto directo en el valor del subsidio que recibe

Valor total de subsidios vs el número de horas trabajadas

ggplot(df, aes(y=total_subsidios , x=num_horas_trabajo))+
  geom_point()     # num_horas_trabajo

En los encuestados que no tienen trabajos y por ende las horas dedicadas a esta labor es cero, se evidencia una gran variabilidad en el valor del subsidio recibido. Esta variable no tiene una relación directa con el valor del subsidio.

Valor total de subsidios vs el sexo

boxplot(df$total_subsidios~df$sexo,
        main = "Total de subsidios recibidos por sexo",
        ylab="Total subsidio",
        xlab = "sexo", las=1,
        col=c("lightgreen","lightpink"))

El sexo no es un factor influyente en el valor total de subsidios que pueda recibir una persona, ya que la distribución en ambos sexos es casi identica.

Correlación lineal

## Correlación
library(corrplot)
## corrplot 0.92 loaded
co <- cor(df[,1:5])
corrplot(co, method = "number")

Por último, evidenciamos que no existe una correlación superior a 0.11 entre el total de subsidios vs las variables dependientes.

A modo de conclusión, pudimos evidenciar que ninguna de las variables independientes tienen una relación directa sobre el valor del subsidio que recibieron los participantes de la encuesta, lo que a priori índica que el modelo de regresión tendrá un \(R^2\) bajo.

Modelo de regresión

Para responder nuestra pregunta de investigación, se ajusta un modelo de regresión lineal cuya variable dependiente es total_subsidios y las variables predictoras o independientes serán la edad, el ingreso, el número de habitantes en el hogar, número de horas semanales dedicadas a trabajo y el sexo

model=lm(total_subsidios ~ edad + ingreso + num_habitantes + num_horas_trabajo + sexo, data=df)
summary(model)
## 
## Call:
## lm(formula = total_subsidios ~ edad + ingreso + num_habitantes + 
##     num_horas_trabajo + sexo, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -270563  -82422  -30843   38805 1662794 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        1.082e+05  3.489e+03  31.018  < 2e-16 ***
## edad              -1.752e+02  5.752e+01  -3.046  0.00232 ** 
## ingreso            1.365e-02  2.521e-03   5.414 6.23e-08 ***
## num_habitantes     8.708e+03  5.835e+02  14.924  < 2e-16 ***
## num_horas_trabajo -2.230e+02  5.952e+01  -3.746  0.00018 ***
## sexo2             -1.193e+02  2.027e+03  -0.059  0.95308    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 137200 on 20199 degrees of freedom
## Multiple R-squared:  0.01332,    Adjusted R-squared:  0.01307 
## F-statistic: 54.53 on 5 and 20199 DF,  p-value: < 2.2e-16

Al análizar el valor p de la prueba t, se evidencia que sólo el sexo no es significativo para el modelo, es decir, no tiene ningun efecto en el valor total de subsidios que puede recibir una persona.

La estimación del coeficiente de la edad índica que por cada año de edad que tenga el encuestado, el valor recibido por subsidios desminuirá \(\$175\), por su parte, el ingreso influye de manera positiva, es decir, por cada peso recibido de ingreso por su labora, puede recibir $ $136$ de subsidio. Respecto al número de habitantes en el hogar, la estimación del coeficiente índica que por cada habitante, aumenta el subsidio en $ $8.708$. Por último, la estimación del coeficiente de número de horas de trabajo indica que por cada hora semanal dedicada al trabajo, el valor total recibido por concepto de subsidios disminuye en 223 pesos.

No obstante, el \(R^2\) del modelo es de solo 0.01332, lo que quiere decir que el modelo, con las variables propuestas, solo explica un \(1.3\%\) de la varianza del valor de subsidios que recibe una persona. Este valor tan bajito ya se podía deducir, dado que en el análisis descriptivo no se evidenciaba relación fuerte de ninguna de las variables independientes con la variable dependiente.

Validación de supuestos

En este punto, se llevará a cabo la validación de supuestos sobre los residuales del modelo: normalidad en los errores, homogeneidad de varianza (homocedasticidad) e independencia, a través de métodos gráficos y pruebas de hipótesis.

### Validación de supuestos de manera gráfica


summary(model$residuals) 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -270563  -82422  -30843       0   38805 1662794
hist(model$residuals, freq = F, ylim = c(0,0.000005), main="Distribución de densidad de los errores")
lines(density(model$residuals), col='blue')

# normalidad en los errores
par(mfrow=c(2,2))
plot(model)

Gráficamente se puede concluir qué:

  • Normalidad: inicialmente, se presenta el gráfico de distribución de densidad de los errores, donde se puede apreciar que aparentemente no se distibuyen normal debido a que no se tiene un efecto de campana con los datos distribuidos alrededor de la media, es decir, la distribución tiene una proporción importante de errores en la cola derecha, lo que se puede llamar cola pesada. Es claro que la interpretación de la gráfica anterior puede ser un poco subjetiva, sin embargo, se presentan algunos gráficos más sobre la normalidad de los errores, en los que se tiene la distribución de los residuales, los residuales estandarizados y el Q-Q plot.
    Al comparar los residuales vs los residuales ajustados, se observan algunos puntos considerados atípicos, ya que los errores están muy por encima del promedio (que es cero) en algunas observaciones; esto quiere decir que la estimación del valor de los subsidios para estas observaciones, presenta un error bastante importante. Estos errores atípicos se observan más claramente en el gráfico ‘Residuals vs Leverage’, donde es evidente la presencia de estos al lado derecho del gráfico y alejados de la concentración de errores.
    En el Q-Q plot, se observa que los residuales no tienen una distribución normal, ya que gráficamente se observa que en la cola derecha de la distribución, se aleja de los cuantiles teoricos de la distribución normal. Por todo lo anterior, se puede concluir que no se cumple el supuesto de normalidad en los errores.

Cabe resaltar que para la validación de este supuesto, sólo se utilizará la prueba gráfica, ya que la prueba Shapiro Wilk está limitada para 5.000 errores o menos, y en este ejercicio se supera este valor.

  • Homogeneidad de la varianza: Para evaluar gráficamente el supuesto de homocedasticidad, podemos observar el gráfico ‘Scale-Location’, que presenta la distribución de la raíz de los errores estandarizados, aquí es evidente que los errores presentan una fuerte dispersión y una tendencia creciente en algunos puntos que pueden ser considerados atípicos, por lo que se puede concluir que no se cumple este supuesto.

Conclusión: De acuerdo a los gráficos presentados anteriormente, el modelo planteado no cumple ninguno de los supuestos sobre los errores. Para confirmar esta conclusión, a continuación se realizarán las puebas de hipótesis apropiadas para validar cada uno de los supuestos:

### Pruebas de Hipotésis

# Homocedasticidad (Breush-Pagan)
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
bptest(model)
## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 26.466, df = 5, p-value = 7.246e-05
# Independencia de errores - Durbin-Watson
dwtest(model)
## 
##  Durbin-Watson test
## 
## data:  model
## DW = 0.46694, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0
# Multicolinealidad

#install.packages("regclass",repos = "http://cran.us.r-project.org")
library(regclass)
## Loading required package: bestglm
## Loading required package: leaps
## Loading required package: VGAM
## Loading required package: stats4
## Loading required package: splines
## 
## Attaching package: 'VGAM'
## The following object is masked from 'package:lmtest':
## 
##     lrtest
## Loading required package: rpart
## Loading required package: randomForest
## randomForest 4.7-1.1
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:ggplot2':
## 
##     margin
## The following object is masked from 'package:dplyr':
## 
##     combine
## Important regclass change from 1.3:
## All functions that had a . in the name now have an _
## all.correlations -> all_correlations, cor.demo -> cor_demo, etc.
VIF(mod = model)
##              edad           ingreso    num_habitantes num_horas_trabajo 
##          1.213756          1.559201          1.007198          1.822924 
##              sexo 
##          1.097896

Para evaluar la homocedasticidad, se realiza la prueba Breush-Pagan, cuya hipótesis nula es que los residuos se distribuyen con varianza constante. Los resultados demuestran que el valor-p es aproximadamente cero, lo que permite concluir que los errores no cumplen con el supuesto de homocedasticidad.

A continuación, se lleva a cabo la prueba Durbin-Watson para evaluar la independencia de los errores, en los que el valor-p permite concluir que se rechaza la hipótesis nula, es decir, que la autocorrelación de los errores es distinta de cero y al ser el estadistico BP positivo, quiere decir que la autocorrelación es positiva. Finalmente se utiliza el VIF para evaluar el supuesto de multicolinealidad, concluyendo que como los valores obtenidos no superan el puntaje de 3, se puede afirma que se cumple con el supuesto de multicolinealidad.