trabajo final data
2022-11-23
Pregunta de investigación
¿Qué factores condicionan el nivel de salario de una persona con discapacidad? ¿Qué factores determinan el salario de las personas con discapacidad?
activamos las librerías
library(rio)
library(haven)
library(dslabs)
library(MASS)
library(car)
library(dplyr)
library(tidyverse)
library(DescTools)#para el pseudo R2
library(ggfortify)
library(see)
library(patchwork)
library(performance)
library(car)
library(lmtest)#para usar breusch pagan
library(nortest)#para usar Kolmogorov-Smirnovabrimos y limpiamos la data
bd=import("Base de Datos II Estudio Nacional de la Discapacidad.sav")
bd_1=bd%>%filter(edad > 18)
subdata<- bd_1[,c("des_puntaje_adulto","educc","sexo","h4","h9","zona","c41_1","disc_grado_adulto")]
data<-na.omit(subdata)Recodificación de variables
1: pertenencia a una comunidad indígena
table(data$h9)##
## 1 2 3 4 5 6 7 9 10 96
## 161 49 17 804 18 12 4 35 10952 9class(data$h9)## [1] "numeric"data$h9=as.numeric(data$h9)
data$indigena <- car::recode(data$h9, "1=1; 2=1; 3=1; 4=1; 5=1; 6=1; 7=1; 8=1; 9=1; 10=2; else = NA")
table(data$indigena)##
## 1 2
## 1100 109522: pertenencia a zona rural o urbana
table(data$zona)##
## 1 2
## 10102 1959#urbano=1; rural=23: nivel de educación niveles: -sin educación formal -educación básica incompleta -educación básica completa -Media incompleto -Media completo -superior incompleta -Superior completa 5.115
data$educc=as.numeric(data$educc)
data$educacion <- car::recode(data$educc, "0=0; 1=1; 2=2; 3=3; 4=4; 5=5; 6=6; else = NA")
table(data$educacion)##
## 0 1 2 3 4 5 6
## 335 1903 1406 1653 3372 1249 2138class(data$educacion)## [1] "numeric"4: si la persona es casada o no
data$h4=as.numeric(data$h4)
data$casado <- car::recode(data$h4, "1=1; 2=2; 3=2; 4=2; 5=2; 6=2; 7=2; else = NA")
table(data$casado)##
## 1 2
## 4321 7740class(data$casado)## [1] "numeric"variables ya limpias
table(data$sexo)##
## 1 2
## 5203 6858#1=hombre; 2=mujer
table(data$disc_grado_adulto)##
## 0 1 2
## 9453 1518 1090#0=persona sin discapacidad; 1= persona con discapacidad leve a moderada; 2=persona con discapacidad severa
table(data$fa33)## < table of extent 0 >recodificación de las variable dependientes del primer y segundo modelo
#primer modelo: ordinal
summary(data$des_puntaje_adulto)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 25.05 38.02 35.13 46.68 100.00data$desempeño=
cut(data$des_puntaje_adulto,
breaks=c(0,25,50,75,100),
labels=c("bajo","medio-bajo","medio-alto","alto"))
table(data$desempeño)##
## bajo medio-bajo medio-alto alto
## 1992 6999 2044 14str(data$desempeño)## Factor w/ 4 levels "bajo","medio-bajo",..: 2 2 1 2 1 2 NA 1 NA 2 ...Primer modelo
data$desempeño=as.factor(data$desempeño)
modelo1<- polr(desempeño~ educacion + sexo + zona +indigena + disc_grado_adulto + c41_1, data = data, Hess=TRUE)
summary(modelo1)## Call:
## polr(formula = desempeño ~ educacion + sexo + zona + indigena +
## disc_grado_adulto + c41_1, data = data, Hess = TRUE)
##
## Coefficients:
## Value Std. Error t value
## educacion -0.12451 0.01357 -9.173
## sexo 0.42914 0.04463 9.616
## zona -0.08566 0.06127 -1.398
## indigena -0.14829 0.07694 -1.927
## disc_grado_adulto 3.10367 0.06206 50.012
## c41_1 -0.46369 0.19720 -2.351
##
## Intercepts:
## Value Std. Error t value
## bajo|medio-bajo -2.3105 0.4368 -5.2900
## medio-bajo|medio-alto 2.2611 0.4381 5.1608
## medio-alto|alto 9.8209 0.5152 19.0639
##
## Residual Deviance: 14209.33
## AIC: 14227.33
## (1025 observations deleted due to missingness)Observar el p-value y determinar la significancia de las variables independientes
summary_table <- coef(summary(modelo1))
summary_table## Value Std. Error t value
## educacion -0.12450620 0.01357240 -9.173485
## sexo 0.42914130 0.04462933 9.615679
## zona -0.08566365 0.06126662 -1.398211
## indigena -0.14828724 0.07693982 -1.927315
## disc_grado_adulto 3.10366954 0.06205886 50.011706
## c41_1 -0.46368540 0.19719829 -2.351366
## bajo|medio-bajo -2.31051659 0.43676682 -5.290046
## medio-bajo|medio-alto 2.26109727 0.43812581 5.160840
## medio-alto|alto 9.82087153 0.51515491 19.063919pval <- pnorm(abs(summary_table[, "t value"]),lower.tail = FALSE)* 2
pval## educacion sexo zona
## 4.579669e-20 6.865502e-22 1.620498e-01
## indigena disc_grado_adulto c41_1
## 5.394043e-02 0.000000e+00 1.870461e-02
## bajo|medio-bajo medio-bajo|medio-alto medio-alto|alto
## 1.222856e-07 2.458440e-07 5.035749e-81summary_table <- cbind(summary_table, "p value" = pval)
summary_table## Value Std. Error t value p value
## educacion -0.12450620 0.01357240 -9.173485 4.579669e-20
## sexo 0.42914130 0.04462933 9.615679 6.865502e-22
## zona -0.08566365 0.06126662 -1.398211 1.620498e-01
## indigena -0.14828724 0.07693982 -1.927315 5.394043e-02
## disc_grado_adulto 3.10366954 0.06205886 50.011706 0.000000e+00
## c41_1 -0.46368540 0.19719829 -2.351366 1.870461e-02
## bajo|medio-bajo -2.31051659 0.43676682 -5.290046 1.222856e-07
## medio-bajo|medio-alto 2.26109727 0.43812581 5.160840 2.458440e-07
## medio-alto|alto 9.82087153 0.51515491 19.063919 5.035749e-81Como se observa, todas las variables escogidas tienen un p-value menor a 0.05, por lo tanto, todas son significativas
Cálculo e interpretación de los exponenciales
exp(coef(modelo1))## educacion sexo zona indigena
## 0.8829328 1.5359381 0.9179029 0.8621834
## disc_grado_adulto c41_1
## 22.2795573 0.6289614*Determinamos cuánto explica nuestro modelo con un pseudo R2”
PseudoR2(modelo1, which = c("Nagelkerke"))## Nagelkerke
## 0.502788El modelo explica un 50%
Construimmos la ecuación
primer corte
num_1 = exp(-3.2040681 - ((0.2550297* ##) + (-0.6630363* ##)+ ( 0.1200079* ##) + (-0.2312537* ##) + ()))
denom_1 = 1 + num_1 p_menorigual_muybajo= num_1/denom_1 p_menorigual_muybajo
coef(summary(modelo1))## Value Std. Error t value
## educacion -0.12450620 0.01357240 -9.173485
## sexo 0.42914130 0.04462933 9.615679
## zona -0.08566365 0.06126662 -1.398211
## indigena -0.14828724 0.07693982 -1.927315
## disc_grado_adulto 3.10366954 0.06205886 50.011706
## c41_1 -0.46368540 0.19719829 -2.351366
## bajo|medio-bajo -2.31051659 0.43676682 -5.290046
## medio-bajo|medio-alto 2.26109727 0.43812581 5.160840
## medio-alto|alto 9.82087153 0.51515491 19.063919vemos las probabilidades
head(modelo1$fitted.values)## bajo medio-bajo medio-alto alto
## 1 0.1845723 0.7717369 0.04366706 2.380202e-05
## 5 0.1845723 0.7717369 0.04366706 2.380202e-05
## 8 0.2474697 0.7220423 0.03047158 1.638324e-05
## 9 0.2132337 0.7500127 0.03673379 1.987863e-05
## 14 0.2132337 0.7500127 0.03673379 1.987863e-05
## 15 0.2040509 0.7571746 0.03875340 2.101564e-05¿Cómo se entiende esto?
recodificación de la segunda variable dependiente y segundo modelo
data<-na.omit(data)
summary(data$des_puntaje_adulto)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.2949 30.9944 39.3860 38.3511 47.2970 100.0000modelo <- lm(des_puntaje_adulto~ educacion + sexo +indigena + disc_grado_adulto + c41_1,data)
summary(modelo)##
## Call:
## lm(formula = des_puntaje_adulto ~ educacion + sexo + indigena +
## disc_grado_adulto + c41_1, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -36.546 -5.441 1.583 7.272 37.997
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.67352 1.66132 23.881 < 2e-16 ***
## educacion -0.85879 0.05741 -14.958 < 2e-16 ***
## sexo 2.54003 0.19768 12.849 < 2e-16 ***
## indigena -0.46898 0.33332 -1.407 0.159459
## disc_grado_adulto 11.84285 0.15668 75.588 < 2e-16 ***
## c41_1 -2.74891 0.75418 -3.645 0.000269 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.17 on 11030 degrees of freedom
## Multiple R-squared: 0.41, Adjusted R-squared: 0.4098
## F-statistic: 1533 on 5 and 11030 DF, p-value: < 2.2e-16normalidad
lillie.test(modelo$residuals) ##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: modelo$residuals
## D = 0.070433, p-value < 2.2e-16Dado que tenemos un p-value (2.2e-16) MENOR a (<) 0.05, se rechaza H0 (distribución normal) y podemos acepta H1 que (no hay distribución normal), por lo que nuestro modelo no cumpliaria con los suspuestos de normalidad homocedasticidad
bptest(modelo)##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 493.21, df = 5, p-value < 2.2e-16En este caso con un p-value (2.2e-16) MENOR a(<) que 0.05 se RECHAZA por lo que nuestro modelo tendria un problema de heterocedasticidad.
no-colinealidad
vif(modelo)## educacion sexo indigena disc_grado_adulto
## 1.073845 1.015632 1.004946 1.099283
## c41_1
## 1.017182independencia
durbinWatsonTest (modelo) ## lag Autocorrelation D-W Statistic p-value
## 1 0.008494958 1.982861 0.366
## Alternative hypothesis: rho != 0Interpretación: Dado que la prueba de Durbin-Watson presenta un Pvale (0.83) MAYOR a(>) 0.05, no podemos rechazar la hipotesis nula, por lo que No existe auto-correlación.
linealidad
plot(modelo,1)
autoplot(modelo)
El modelo cumple con el supuesto de linealidad