Taller Propiedades estimador MCO
Universidad del Norte
Luisa F. Llorente de la Ossa
21/8/2021
Usted quiere estimar el retorno a la educación. Para ello plantea una ecuación de Mincer de la forma
\[ln(w_i)=\alpha+\beta_1{si}+\beta_2{expi}+\beta_3{exp^2}+γMujer_i+e_i\]
se carga la base de datos y se añade la variable de ln(w)
load("C:/Users/jebet/Downloads/wagew.rda")
colnames(wagew)<-c("sexo","edad","nivel educativo","ultimo_grado","ingreso_mensual","años_escolaridad")
ingreso=log(wagew$ingreso_mensual)
#insertamos la columna de ln(w) a la dataframe1. Análisis gráfico
1.1 Grafico de dispersión
library(ggplot2)
wagew<-cbind(wagew,ln_w=c(ingreso))#Añadir columna de ln(w) a la dataframe
Datos <- wagew[apply(wagew, 1, function(row) all(row !=0 )),]#Sacar datos de 0
#Crear graficolibrary(ggplot2)
ggplot(Datos,aes(x=Datos$años_escolaridad,y=Datos$ln_w))+geom_point(shape=21)+
labs(
title = paste("Grafico de dispersion"),
subtitle = "Relacion años de escolaridad e ingreso",
caption = "Autor: Luisa Llorente",
x = "Años de escolaridad",
y = "Logaritmo natural del ingreso")+theme_linedraw()+stat_smooth(method=lm,
col="#E90F0F",
se=FALSE,
size=1)## `geom_smooth()` using formula 'y ~ x'
Del anterior grafico se denota que hay una relación proporcional entre ambas variables, es decir, una relación positiva entre el logaritmo del ingreso laboral mensual y los años de escolaridad,por tanto, se puede intuir que a mayor educación se devengan salarios más altos. Asimismo,los datos muestra una tendencia de distribución uniforme.
1.2 Gráfico- dsitribución empírica, densidad, de ln(w) para hombres y para mujeres (rango de valores del eje X entre 12 y 16)
library(ggplot2)
Sexo<- Datos$sexo
Sexo<-factor(Datos$sexo, labels = c("Masculino","Femenino"))
ggplot(Datos,aes(x=Datos$ln_w,fill=Sexo))+
geom_density()+scale_x_continuous(limits=c(12,16))+
labs(
title=paste("Grafico de densidades superpuesto"),
x = "Logaritmo del ingreso laboral mensual",
y = "Densidad")
En el anterior gráfico de densidad se puede visualizar la distribución de los datos para la muestra de hombres y mujeres. Se observa que la distribución de los ingresos para ambos sexos es similar.A excepción de una leve diferencia en la concetración de los datos de los hombres en el centro de la disribución, no se evidencian diferencias importantes.
2. Estimación del modelo y reporte de los resultados
library(fastDummies)## Warning: package 'fastDummies' was built under R version 4.0.5Datos<-dummy_cols(Datos,select_columns=c("sexo"))
exp=(Datos$edad-Datos$años_escolaridad-6)
Datos<-cbind(Datos,exp=c(exp))
Ingreso<-Datos$ln_w
añosEscolaridad<-Datos$años_escolaridad
Experiencia<-Datos$exp
Mujer<-Datos$sexo_2
Rg1<-lm(Ingreso~añosEscolaridad+Experiencia+Mujer,data=Datos)
summary(Rg1)##
## Call:
## lm(formula = Ingreso ~ añosEscolaridad + Experiencia + Mujer,
## data = Datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.8015 -0.2507 0.0296 0.3089 3.5420
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.3047249 0.0064639 1903.60 <2e-16 ***
## añosEscolaridad 0.1214209 0.0004283 283.46 <2e-16 ***
## Experiencia 0.0161493 0.0001231 131.16 <2e-16 ***
## Mujer -0.2173294 0.0030715 -70.76 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5596 on 140096 degrees of freedom
## Multiple R-squared: 0.3687, Adjusted R-squared: 0.3687
## F-statistic: 2.727e+04 on 3 and 140096 DF, p-value: < 2.2e-16library(sjPlot)## Warning: package 'sjPlot' was built under R version 4.0.5## Registered S3 methods overwritten by 'parameters':
## method from
## as.double.parameters_kurtosis datawizard
## as.double.parameters_skewness datawizard
## as.double.parameters_smoothness datawizard
## as.numeric.parameters_kurtosis datawizard
## as.numeric.parameters_skewness datawizard
## as.numeric.parameters_smoothness datawizard
## print.parameters_distribution datawizard
## print.parameters_kurtosis datawizard
## print.parameters_skewness datawizard
## summary.parameters_kurtosis datawizard
## summary.parameters_skewness datawizardlibrary (Rcpp)## Warning: package 'Rcpp' was built under R version 4.0.5tab_model(Rg1,show.est=TRUE,show.se=TRUE,show.std=TRUE,show.df=TRUE,show.p=TRUE,show.fstat=TRUE,digits=2,title="Resumen del modelo",p.style="numeric_stars")| Ingreso | ||||||||
|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | std. Error | std. Beta | standardized std. Error | CI | standardized CI | p | df |
| (Intercept) | 12.30 *** | 0.01 | 0.00 | 0.00 | 12.29 – 12.32 | -0.00 – 0.00 | <0.001 | 140096.00 |
| añosEscolaridad | 0.12 *** | 0.00 | 0.66 | 0.00 | 0.12 – 0.12 | 0.65 – 0.66 | <0.001 | 140096.00 |
| Experiencia | 0.02 *** | 0.00 | 0.30 | 0.00 | 0.02 – 0.02 | 0.29 – 0.30 | <0.001 | 140096.00 |
| Mujer | -0.22 *** | 0.00 | -0.15 | 0.00 | -0.22 – -0.21 | -0.16 – -0.15 | <0.001 | 140096.00 |
| Observations | 140100 | |||||||
| R2 / R2 adjusted | 0.369 / 0.369 | |||||||
|
||||||||
A la vista del resultado analítico podemos afirmar que el ajuste del modelo es aceptable, ya que el valor de R2 es 0,369, por tanto, las variables escogidas explican los ingresos aproximadamente en un 37% o que el 36,9% de la variabilidad del logaritmo natural de los ingresos mensuales a su promedio es explicado por el modelo de regresión ajustado. Asimismo, Se observa que todas las variables son significativas al 99% y como P-Valor<0,05 se rechaza la hipotesis nula.
Cuando no se tienen en cuenta los años de eduación,la experiencia y la categoría de sexo femenino se obtiene un intercepto estimado de 12.30. También, hay una relación positiva entre las variables de los años de educación y experencia con el logaritmo natural de los ingresos mensuales, en este sentido, en promedio, por cada año adicional de educación los ingresos mensuales aumentan en 12%, mientras que por cada año adicional de experencia los ingresos aumentan en un 2%. Sin embargo, se observa una relación negativa entre el hecho de ser mujer y los ingresos, situación que resulta alarmanete,pues, lo que que se observa en el modelo es que las trabajdaoras registran una penalidad salarial de 22%, obteniendo así menos de ingresos.
3. ¿El coeficiente estimado \(\beta_1\) tiene interpretación causal?
La relación entre el salario y la formación academica ha sido sido tema de investigación de diversos autores a nivel nacional, esto con el fin de entender la tasa de retornos de la educación y los salarios. Por ejemplo, Fukusaki, G. Y. (2007), Nordin M.; Persson, I. (2010) y Kasakis, P.; Faggian, A (2017) argumentan que la experiencia y los años de preparación es determinante en el salario en Colombia, sin embargo, también se debe tener en cuenta otras variables como la zona geografica y el sector económico. Por su parte, hala-Nieto, Y. & Carrera-Canizales, X. (2018) confirman que hay una relación evidente entre educación y salario, al encontrar en un modelo sobre determinantes del salario en la ciudad de Bogotá y argumentar lo siguiente “por un aumento en un año de educación, el salario mensual aumentaría en 12.5675%”. Por ende, se comprueba que hay una relación causal positiva entre los años de escolaridad y el nivel salarial.
De esta forma, teniendo en cuenta lo anterior, los años de escolaridad no implican una correlación con los salarios para el caso de Colombia resultando ser una variable apropiada en el modelo estimado, así,una mayor preparación academica es relevante para obtener mayores ingresos y es un aspecto valorado en el mercado nacional. No obtante, esto varía dependiendo de aspectos como los países. Por ejemplo, en Alemania se privilegia la experencia sobre la eduación e incluso esto puede ser un factor que a nivel general varía dependiendo de la zona geografica, el tipo de actividad y si se desarrolla en una zona rural o urbana.
4. Insesgamiento e inferencia: ejercicio de simulación
1. Estime el modelo con la muestra de datos. ¿Como se compara este coeficiente estimado con el de la estimacion inicial?
library(mvtnorm)## Warning: package 'mvtnorm' was built under R version 4.0.5N<-10000
coefs<-cbind("hat_beta_1" = numeric(1000), "hat_beta_2" = numeric(1000))
set.seed(1)
X <- rmvnorm(N, c(50, 100), sigma = cbind(c(10, 2.5), c(2.5, 10)))
u <- rnorm(N, sd = 5)
Y <- 5 + 0.5 * X[, 1] + 3 * X[, 2] + u
xdf<-data.frame(X,Y)
dfs<-sample(c(TRUE,FALSE),nrow(xdf),replace=TRUE,prob=c(0.01,0.99))
dfs<-xdf[dfs,]
## Posteriormente, se hace la regresion:
Mod5<- lm(Y~X1, data = dfs)
summary(Mod5)##
## Call:
## lm(formula = Y ~ X1, data = dfs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.7058 -5.4948 0.6255 5.8224 31.4948
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 263.0701 18.0856 14.55 < 2e-16 ***
## X1 1.2899 0.3644 3.54 0.000645 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.55 on 87 degrees of freedom
## Multiple R-squared: 0.1259, Adjusted R-squared: 0.1159
## F-statistic: 12.53 on 1 and 87 DF, p-value: 0.000645En este modelo se observa que el error estándar es mayor al de la estimacion inicial. sin embargo, las variables explicativas siguen conservando los niveles de significancia al 99%.
2. Estime el modelo anterior para 1000 muestras diferentes. Guarde el coeficiente estimado en cada regresion, calcule el promedio de los coeficientes estimados, y grafique la distribucion empirica. ¿Podria decir que el estimador es insesgado? ¿por qué?
## Primer paso; ciclo que genere y guarde lass estimaciones.
nrsim<-1000
for (i in 1:nrsim){
dfs<-sample(c(TRUE,FALSE),nrow(xdf),replace=TRUE,prob=c(0.01,0.99))
dfs<-xdf[dfs,]
ols<-lm(Y~X1,data=dfs)
coefs[i,]<-coef(ols)[-1]
}
coefs.df<-data.frame(coefs)
## Segundo paso: Calculo del promedio.
library(dplyr)## Warning: package 'dplyr' was built under R version 4.0.5##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionmean(coefs.df$hat_beta_1)## [1] 1.255209msd1<-coefs.df%>%summarise(b1m=mean(coefs.df$hat_beta_1),
b1sd=sd(coefs.df$hat_beta_1))
## Tercer paso: Grafico de la distribucion empirica.
ggplot(coefs.df,aes(x=hat_beta_1))+geom_density(fill="lightblue",alpha=0.5)+
geom_vline(xintercept=mean(coefs.df$hat_beta_1),linetype=4)+theme_minimal()+
labs(title="Distribución empirica",subtitle="1000 muestras",x=expression(beta_hat[1])) 
Como el valor promedio estimado difiere del parametro poblacional conocido se concluye que el estimador no es insesgado.
3. Vuelva a simular los datos, pero esta vez cambie el valor 2.5 que define la covarianza entre x1 y x2 por 0. Repita el punto 2. En este caso, ¿podria decir que el estimador es insesgado? ¿por que?
# Primer paso
N<-10000
coefs<-cbind("hat_beta_1" = numeric(1000), "hat_beta_2" = numeric(1000))
set.seed(1)
X <- rmvnorm(N, c(50, 100), sigma = cbind(c(10, 0), c(0, 10)))
u <- rnorm(N, sd = 5)
Y <- 5 + 0.5 * X[, 1] + 3 * X[, 2] + u
xdf<-data.frame(X,Y)
dfs<-sample(c(TRUE,FALSE),nrow(xdf),replace=TRUE,prob=c(0.01,0.99))
dfs<-xdf[dfs,]
# Segundo paso
Mod4<- lm(Y~X1, data = dfs)
summary(Mod4)##
## Call:
## lm(formula = Y ~ X1, data = dfs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.4290 -6.2288 0.8539 5.8100 31.5650
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 300.4344 18.2341 16.477 <2e-16 ***
## X1 0.5385 0.3668 1.468 0.146
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.788 on 87 degrees of freedom
## Multiple R-squared: 0.02417, Adjusted R-squared: 0.01295
## F-statistic: 2.155 on 1 and 87 DF, p-value: 0.1457nrsim<-1000
for (i in 1:nrsim){
dfs<-sample(c(TRUE,FALSE),nrow(xdf),replace=TRUE,prob=c(0.01,0.99))
dfs<-xdf[dfs,]
ols<-lm(Y~X1,data=dfs)
coefs[i,]<-coef(ols)[-1]
}
coefs.df<-data.frame(coefs)
mean(coefs.df$hat_beta_1)## [1] 0.518473msd2<-coefs.df%>%summarise(b1m=mean(coefs.df$hat_beta_1),
b1sd=sd(coefs.df$hat_beta_1))
# Grafico
ggplot(coefs.df,aes(x=hat_beta_1))+geom_density(fill="lightblue",alpha=0.4)+
geom_vline(xintercept=mean(coefs.df$hat_beta_1),linetype=4)+theme_minimal()+
labs(title="Distribución empirica",subtitle="1000 muestras",x=expression(beta_hat[1])) 
Se observa que el valor estimado no se aleja mucho del parámetro poblacional conocido, por ende, el estimador es insesgado.
Varianza del estimador
4. Tome una muestra de n = 20 y estime nuevamente la varianza del estimador y el error estandar. Compare lo obtenido en el ejercicio de demostracion. Explique la diferencia.
library(mvtnorm)
N<-10000
coefs<-cbind("hat_beta_1" = numeric(1000), "hat_beta_2" = numeric(1000))
set.seed(1)
X <- rmvnorm(N, c(50, 100), sigma = cbind(c(10, 2.5), c(2.5, 10)))
u <- rnorm(N, sd = 5)
Y <- 5 + 0.5 * X[, 1] + 3 * X[, 2] + u
xdf<-data.frame(X,Y)
dfs<-sample(c(TRUE,FALSE),nrow(xdf),replace=TRUE,prob=c(0.002,0.998))
dfs<-xdf[dfs,]
Mod1<- lm(Y~X1+X2, data = dfs)
summary(Mod1)##
## Call:
## lm(formula = Y ~ X1 + X2, data = dfs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.302 -3.582 1.010 3.594 7.756
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.4143 56.3604 0.096 0.925196
## X1 0.2463 0.7140 0.345 0.736644
## X2 3.1040 0.6149 5.048 0.000373 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.804 on 11 degrees of freedom
## Multiple R-squared: 0.7448, Adjusted R-squared: 0.6984
## F-statistic: 16.05 on 2 and 11 DF, p-value: 0.0005468Cuando se reduce el tamaño de la muestra tanto la varianza del estimador como el error estandar tienden a aumentar, lo que conlleva a una pérdida de eficiencia del modelo, puesto que, su precisión para explicar la variable dependiente disminuye, de esta forma, se observa que un menor número de observaciones conlleva a que los intervalos de confianza y el p-valor sean menos significativos Asimismo, \(x_1\) obtiene una reducción de su estimación a diferencia de \(x_2\) que presenta un aumento.
5. En el ejercicio de demostracion, triplique la covarianza de x1 y x2. Estime nuevamente la varianza del estimador y el error estandar. Compare y explique.
N<-10000
coefs<-cbind("hat_beta_1" = numeric(1000), "hat_beta_2" = numeric(1000))
set.seed(1)
X <- rmvnorm(N, c(50, 100), sigma = cbind(c(10, 7.5), c(7.5, 10)))
u <- rnorm(N, sd = 5)
Y <- 5 + 0.5 * X[, 1] + 3 * X[, 2] + u
xdf<-data.frame(X,Y)
dfs<-sample(c(TRUE,FALSE),nrow(xdf),replace=TRUE,prob=c(0.01,0.99))
dfs<-xdf[dfs,]
Mod2<- lm(Y~X1+X2, data = dfs)
summary(Mod2)##
## Call:
## lm(formula = Y ~ X1 + X2, data = dfs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.8316 -3.9653 -0.0361 3.9530 15.3369
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.7816 21.4429 0.083 0.933975
## X1 1.1247 0.2813 3.998 0.000135 ***
## X2 2.7125 0.2854 9.505 4.6e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.184 on 86 degrees of freedom
## Multiple R-squared: 0.7744, Adjusted R-squared: 0.7691
## F-statistic: 147.6 on 2 and 86 DF, p-value: < 2.2e-16Tomando en cuenta el modelo original, se puede observar que triplicar la covarianza no influyó de forma importante en la significancia, sin embargo, esto conllevó a un aumento de error estándar, y, por ende, repercutió de forma negativa en la precisión del modelo.
6. Realice el ejercicio anterior, pero con una muestra de n=20. Compare y explique.
N<-10000
coefs<-cbind("hat_beta_1" = numeric(1000), "hat_beta_2" = numeric(1000))
set.seed(1)
X <- rmvnorm(N, c(50, 100), sigma = cbind(c(10, 7.5), c(7.5, 10)))
u <- rnorm(N, sd = 5)
Y <- 5 + 0.5 * X[, 1] + 3 * X[, 2] + u
xdf<-data.frame(X,Y)
dfs<-sample(c(TRUE,FALSE),nrow(xdf),replace=TRUE,prob=c(0.002,0.998))
dfs<-xdf[dfs,]
Mod3<- lm(Y~X1+X2, data = dfs)
summary(Mod3)##
## Call:
## lm(formula = Y ~ X1 + X2, data = dfs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.302 -3.582 1.010 3.594 7.756
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.8713 60.3651 -0.048 0.96292
## X1 0.1269 1.0588 0.120 0.90674
## X2 3.2465 0.9636 3.369 0.00626 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.804 on 11 degrees of freedom
## Multiple R-squared: 0.7662, Adjusted R-squared: 0.7237
## F-statistic: 18.02 on 2 and 11 DF, p-value: 0.0003378Con una muestra de n=20 y la nueva covarianza se observa una diminución del R2,un aumento del P-valor en las variables seguido de un inremento del error estándar (pérdida de precisión) y la varianza.
El ejercicio propuesto genera las siguientes conclusiones.Cumplir el suspuesto de la insesgadez es relevante para que los estimadores sean más consistentes y no diverjan significativamente del valor real; tener muestras grandes y que estén distribuidas de forma normal facilita que haya un menor error estándar y una menor varianza,y, por tanto, el modelo sea más preciso y eficiente.
Referencias
Fukusaki, G. Y. (2007). Retornos a la educación superior en el mercado laboral: ¿vale la pena el esfuerzo
Chala-Nieto, Y. & Carrera-Canizales, X. (2018). Relación de los años de educación, experiencia y género, con el salario en Bogotá para el año 2014. Trabajo de Grado. Universidad Católica de Colombia. Facultad de Ciencias Económicas y Administrativas. Programa de Economía. Bogotá, Colombia