Laboratorio 6
Federico Daverio
rm(list=ls())
setwd("C:/DAVE2/CIDE/Doc/Econometria II/Lab_6")
options(scipen=999)
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.0.5## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.4 v dplyr 1.0.7
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 2.0.1 v forcats 0.5.1## Warning: package 'ggplot2' was built under R version 4.0.5## Warning: package 'tibble' was built under R version 4.0.5## Warning: package 'tidyr' was built under R version 4.0.5## Warning: package 'readr' was built under R version 4.0.5## Warning: package 'dplyr' was built under R version 4.0.5## Warning: package 'forcats' was built under R version 4.0.5## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()## Warning: package 'janitor' was built under R version 4.0.5##
## Attaching package: 'janitor'## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.test## Warning: package 'modelsummary' was built under R version 4.0.5## Loading required package: dfidx## Warning: package 'dfidx' was built under R version 4.0.5##
## Attaching package: 'dfidx'## The following object is masked from 'package:stats':
##
## filter## Warning: package 'sandwich' was built under R version 4.0.5## Loading required package: msmModelos con datos truncados
La regresión truncada se utiliza para modelar variables dependientes para las cuales algunas de las observaciones no se incluyen en el análisis debido al valor de la variable dependiente.
## Loading required package: foreign## Loading required package: truncreg## Warning: package 'truncreg' was built under R version 4.0.3## Loading required package: maxLik## Warning: package 'maxLik' was built under R version 4.0.5## Loading required package: miscTools##
## Please cite the 'maxLik' package as:
## Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
##
## If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
## https://r-forge.r-project.org/projects/maxlik/## Loading required package: boot##
## Attaching package: 'boot'## The following object is masked from 'package:msm':
##
## cavEjemplos
Ejemplo 1. Un estudio de estudiantes en un programa especial GATE (educación para dotados y talentosos) desea modelar el rendimiento en función de las habilidades del lenguaje y el tipo de programa en el que el estudiante está actualmente inscrito. Una de las principales preocupaciones es que los estudiantes deben tener un puntaje mínimo de 40 para ingresar al programa especial. Por lo tanto, la muestra se trunca en una puntuación de logro de 40.
Ejemplo 2. Un investigador tiene datos para una muestra de estadounidenses cuyos ingresos están por encima de la línea de pobreza. Por tanto, se trunca la parte inferior de la distribución del ingreso. Si el investigador tuviera una muestra de estadounidenses cuyos ingresos estuvieran en la línea de pobreza o por debajo de ella, la parte superior de la distribución de ingresos se truncaría. En otras palabras, el truncamiento es el resultado de muestrear solo una parte de la distribución de la variable de resultado.
Datos
Sigamos con el ejemplo 1 de arriba. Tenemos un archivo de datos hipotéticos, truncreg.dta, con 178 observaciones. La variable de resultado se llama achiv y la variable de puntuación de la prueba de idioma se llama langscore. La variable prog es una variable predictiva categórica con tres niveles que indican el tipo de programa en el que se matricularon los estudiantes.
Veamos los datos. Siempre es una buena idea comenzar con estadísticas descriptivas.
## id achiv langscore prog
## Min. : 3.00 Min. :41.00 Min. :31.00 general : 40
## 1st Qu.: 55.25 1st Qu.:47.00 1st Qu.:47.50 academic:101
## Median :102.50 Median :52.00 Median :56.00 vocation: 37
## Mean :103.62 Mean :54.24 Mean :54.01
## 3rd Qu.:151.75 3rd Qu.:63.00 3rd Qu.:61.75
## Max. :200.00 Max. :76.00 Max. :67.00Hacemos un analisis grafica:
Y unos boxplot:
Aproximamos con un polinomio:
ggplot(dat, aes(x = langscore, y = achiv)) +
geom_point() +
stat_smooth(method = "loess") +
facet_grid(. ~ prog, margins=TRUE)## `geom_smooth()` using formula 'y ~ x'
Posibles metodos de analisis
Regresión OLS: puede analizar estos datos mediante la regresión OLS. La regresión MCO no ajustará las estimaciones de los coeficientes para tener en cuenta el efecto de truncar la muestra en 40, y los coeficientes pueden estar gravemente sesgados. Esto puede conceptualizarse como un error de especificación del modelo (Heckman, 1979).
Regresión truncada: la regresión truncada aborda el sesgo introducido cuando se usa la regresión OLS con datos truncados. Tenga en cuenta que con la regresión truncada, la varianza de la variable de resultado se reduce en comparación con la distribución que no está truncada. Además, si la parte inferior de la distribución está truncada, entonces la media de la variable truncada será mayor que la media de la variable no truncada; si el truncamiento es de arriba, la media de la variable truncada será menor que la de la variable no truncada. Estos tipos de modelos también pueden conceptualizarse como modelos de selección de Heckman, que se utilizan para corregir el sesgo de selección muestral.
Regresión censurada: a veces, los conceptos de truncamiento y censura se confunden. Con los datos censurados tenemos todas las observaciones, pero no conocemos los valores “verdaderos” de algunas de ellas. Con el truncamiento, algunas de las observaciones no se incluyen en el análisis debido al valor de la variable de resultado. Sería inapropiado analizar los datos de nuestro ejemplo utilizando un modelo de regresión censurado.
Analisis
A continuación, usamos la función truncreg en el paquete truncreg para estimar un modelo de regresión truncado. El argumento de punto indica dónde se truncan los datos y la dirección indica si están truncados a la izquierda o a la derecha.
##
## Call:
## truncreg(formula = achiv ~ langscore + prog, data = dat, point = 40,
## direction = "left")
##
## BFGS maximization method
## 57 iterations, 0h:0m:0s
## g'(-H)^-1g = 2.5E-05
##
##
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 11.29942 6.77173 1.6686 0.09519 .
## langscore 0.71267 0.11446 6.2264 0.0000000004773 ***
## progacademic 4.06267 2.05432 1.9776 0.04797 *
## progvocation -1.14422 2.66958 -0.4286 0.66821
## sigma 8.75368 0.66647 13.1343 < 0.00000000000000022 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: -591.31 on 5 DfEn la tabla de coeficientes, tenemos los coeficientes de regresión truncados, el error estándar de los coeficientes, las pruebas z de Wald (coeficiente / se) y el valor p asociado con cada prueba z (mostrado como valores t).
La estadística auxiliar / sigma es equivalente al error estándar de estimación en la regresión OLS. El valor de 8,76 se puede comparar con la desviación estándar de logro que fue de 8,96. Esto muestra una modesta reducción. La salida también contiene una estimación del error estándar de sigma.
La variable langscore es estadísticamente significativa. Un aumento de unidad en la puntuación de lenguaje conduce a un aumento de .71 en el rendimiento previsto. Una de las variables indicadoras de prog también es estadísticamente significativa. En comparación con los programas generales, los programas académicos son aproximadamente 4.07 más altos.
Para determinar si el prog en sí es estadísticamente significativo, podemos probar modelos considerandolo y no conisderandolo para la prueba de dos grados de libertad de esta variable.
## 'log Lik.' 0.02516651 (df=3)La prueba de chi-cuadrado de dos grados de libertad indica que el prog es un predictor estadísticamente significativo del logro. Podemos obtener las medias esperadas para cada programa en la media de langscore reparametrizando el modelo.
Podemos centrar nuestra regresiòn:
dat <- within(dat, {mlangscore <- langscore - mean(langscore)})
malt <- truncreg(achiv ~ 0 + mlangscore + prog, data = dat, point = 40)
summary(malt)##
## Call:
## truncreg(formula = achiv ~ 0 + mlangscore + prog, data = dat,
## point = 40)
##
## BFGS maximization method
## 21 iterations, 0h:0m:0s
## g'(-H)^-1g = 4.3E-07
##
##
##
## Coefficients :
## Estimate Std. Error t-value Pr(>|t|)
## mlangscore 0.71259 0.11448 6.2248 0.000000000482 ***
## proggeneral 49.78926 1.89714 26.2443 < 0.00000000000000022 ***
## progacademic 53.85340 1.15012 46.8242 < 0.00000000000000022 ***
## progvocation 48.65315 2.14049 22.7299 < 0.00000000000000022 ***
## sigma 8.75545 0.66684 13.1299 < 0.00000000000000022 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: -591.31 on 5 DfObservamos que todo lo que ha cambiado es que el intercepto se ha ido y las puntuaciones del programa ahora son los valores esperados cuando langscore está en su media para cada tipo de programa.
También podríamos calcular los intervalos de confianza de arranque si quisiéramos. Primero, definimos una función que devuelve los parámetros de interés y luego usamos la función de arranque para ejecutar el bootstrap.
f <- function(data, i) {
require(truncreg)
m <- truncreg(formula = achiv ~ langscore + prog, data = data[i, ], point = 40)
as.vector(t(summary(m)$coefficients[, 1:2]))
}
set.seed(10)
(res <- boot(dat, f, R = 1200, parallel = "snow", ncpus = 4))##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = dat, statistic = f, R = 1200, parallel = "snow",
## ncpus = 4)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 11.2994158 -0.0750432302 5.91023663
## t2* 6.7717257 0.0099996401 0.86430871
## t3* 0.7126732 0.0019041192 0.09730007
## t4* 0.1144602 0.0002004845 0.01366713
## t5* 4.0626698 -0.0485288252 1.94566248
## t6* 2.0543191 0.0076934779 0.24494415
## t7* -1.1442162 0.0425573767 2.81757637
## t8* 2.6695799 0.0219657954 0.29789986
## t9* 8.7536778 -0.0900719907 0.54763875
## t10* 0.6664744 -0.0083607460 0.07505858Podríamos usar el error estándar bootstrap para obtener una aproximación normal para una prueba de significancia e intervalos de confianza para cada parámetro. Sin embargo, en su lugar obtendremos el percentil y los intervalos de confianza del 95 por ciento ajustados al sesgo, utilizando la función boot.ci.
# basic parameter estimates with percentile and bias adjusted CIs
parms <- t(sapply(c(1, 3, 5, 7, 9), function(i) {
out <- boot.ci(res, index = c(i, i + 1), type = c("perc", "bca"))
with(out, c(Est = t0, pLL = percent[4], pUL = percent[5], bcaLL = bca[4],
bcaLL = bca[5]))
}))
# add row names
row.names(parms) <- names(coef(m))
# print results
parms## Est pLL pUL bcaLL bcaLL
## (Intercept) 11.2994158 -0.47667113 22.8722200 -0.55662632 22.8711476
## langscore 0.7126732 0.52197893 0.9055977 0.52286226 0.9065408
## progacademic 4.0626698 -0.06878459 7.8897378 0.04873108 7.9296538
## progvocation -1.1442162 -6.79210189 4.3581860 -6.92519441 4.0693050
## sigma 8.7536778 7.62374802 9.7726403 7.86169145 9.9453352Las conclusiones son las mismas que las de las pruebas del modelo predeterminado. Podemos calcular una estimación aproximada del grado de asociación del modelo general, correlacionando achiv con el valor predicho y elevando al cuadrado el resultado.
## [1] 0.5524392El r cuadro serà:
## [1] 0.305189El valor calculado de .31 es una estimación aproximada del R2 que encontraría en una regresión de MCO. La correlación al cuadrado entre los valores de aptitud académica observados y pronosticados es de aproximadamente 0,31, lo que indica que estos predictores explicaron más del 30% de la variabilidad en la variable de resultado.
Cuidado
La función truncreg está diseñada para funcionar cuando el truncado está en la variable de resultado en el modelo. Es posible tener muestras truncadas en función de uno o más predictores. Por ejemplo, modelar el GPA de la universidad en función del GPA de la escuela secundaria (HSGPA) y los puntajes del SAT implica una muestra que se trunca en función de los predictores, es decir, solo los estudiantes con puntajes HSGPA y SAT más altos son admitidos en la universidad.
Debemos tener cuidado con el valor que se usa como valor de truncamiento, ya que afecta la estimación de los coeficientes y los errores estándar. En el ejemplo anterior, si hubiéramos usado point = 39 en lugar de point = 40, los resultados habrían sido ligeramente diferentes. No importa que no haya valores de 40 en nuestra muestra.
Modelos Censurados
Hasta ahora nos enfocamos sobre todo a limitaciones en el tipo de variables (binaria, ordenada, truncada, etc).
Ahora regresamos a una variable continua, pero vemos lo que pasa si la muestra no es representativa para la poblacion.
Definicion (Sesgo de seleccion): Hay un sesgo de seleccion si la probabilidad de una observacion en particular de estar en la muestra depende del fenomeno que estamos analizando.
Ejemplos de sesgo de seleccion:
Los problemas de sesgo de seleccion estan casi siempre presentes en Econometria aplicada
Modelos de seleccion Heckit o Correccion de Heckman permite corregir el sesgo de seleccion.
La correccion de Heckman, un metodo estadistico en dos pasos que permite corregir las muestras no seleccionadas al azar.
#La seleccion: formalmente Para formalizar este problema de seleccion, consideramos dos ecuaciones. Primero tenemos la ecuacion de seleccion:
\[z_i^\ast=w_i^{T}\gamma+v_i\] donde \(z_i^\ast\) es una variable no observada. Y tenemos la regresion de interes: \[y_i= x_i^T\beta+u_i\] \(y_i\) es observado solo si \(z_i^\ast >0\)
El punto crucial es que se observa \(y_i\) solo si \(z_i^\ast >0\)
Si solo tomamos en cuenta los datos observados,date cuenta que \[E[y_i|y_i \quad observado]=x_i^T\beta+E[u_i|v_i>-w_i^{T}\gamma]\]
Por lo tanto la estimacion \(y_i=x_i^T\beta+u_i\) genera un sesgo en \(\beta\) debido a la exclusion de la lambda de Heckman(similar a un problema de variable omitida).
#Idea de la correccion de Heckman - Suponga que un investigador desea estimar los determinantes de la oferta salarial, pero tiene acceso a observaciones solo para los que trabajan. - Dado que las personas que trabajan son seleccionadas no aleatoriamente de la poblacion, la estimacion de los salarios de la subpoblacion que tiene trabajo puede introducir un sesgo.
La correccion de Heackman se lleva a cabo en dos etapas.
En la primera etapa, el investigador formula un modelo para la probabilidad de tener trabajo. La especificacion canonica es una regresion probit (Ecuacion de seleccion o Probit selection equation). La estimacion del modelo produce resultados que se pueden utilizar para predecir la probabilidad de empleo para cada individuo.
En la segunda etapa, el investigador corrige para el error de seleccion mediante la incorporacion de una transformacion de estas probabilidades individuales predichas como una variables explicativa adicional (Modelo de interes o Outcome equation)
El modelo Heckit se puede estimar por MV o por el Heckman’s two-step procedure.
-Heckman’s two-step procedure en ecuaciones
Defines \(D_i=1_{\text{observacion i es parte de la muestra}}\)
El modelo consiste en dos ecuaciones (similar al two-part model), primero estimamos la ecuacion de seleccion con un modelo probit
\[P(D=1)=\Phi(X_1\gamma)\] donde \(X_1\) son un set de regresores.
Segundo estimamos la ecuacion principal (de interes) por MCO e incluimos un termino de correccion:
\[y=X_2\beta+\sigma_{12}\lambda(X_1\hat{\lambda})+v\]
donde \[\lambda(X_1\hat{\lambda}) \quad \text{es el inverse Mills ratio.}\]
si \(\sigma_{12}\neq0\) es necesario hacer la correcion por sesgo de seleccion.
Al asumir el supuesto de normalidad bivariada de los errores tenemos una expresion de la forma \[E[y_i|y_i \quad observado]=x_i^T\beta+h(\rho_{uv},\sigma_u,\sigma_v)\] en concreto:
\[E[y_i|y_i\quad\text{es observado}]=x_i^T\beta+\rho_{ uv}\sigma_u\lambda_i\left(\dfrac{w_i^T\hat{\lambda}}{\sigma_v}\right)\]
equivalente a: \[E[y_i|y_i\quad\text{es observado}]=x_i^T\beta+\beta_{\lambda}\lambda_i\left(\dfrac{-w_i^T\hat{\lambda}}{\sigma_v}\right)\] y el inverse Mills ratio asume la forma \[\lambda(X_1\hat{\lambda})=\dfrac{\phi(X_1\hat{\lambda})}{\Phi(X_1\hat{\lambda})}\]
Para que \(\beta_{\lambda}\) sea \(0\) necesitas que \(\rho_{uv}\) se cero. Por ese motiva se plantea la prueba de hipotesis \(H_0:\rho=0\) vs \(H_a: \rho\neq 0\).
El problema de seleccion puede verse como un problema de variables omitida \[y|z_i^\ast>0= E[y_i|y_i \text{observado}]+ u_i\] \[y|z_i^\ast>0= x_i^T\beta+\beta_{\lambda}\lambda_i\left(\dfrac{w_i^T\hat{\lambda}}{\sigma_v}\right)+ u_i\]
recuerda que \(\beta_{\lambda}=\rho_{uv}\sigma_u\)
- Teoricamente hay dos condiciones en las cuales, el estimador MCO puede dar resultados sin sesgo:
- Si \(\beta_{\lambda} = 0\), lo cual es equivalente a decir que \(\rho_{uv} = 0\), no tememos un sesgo, porque el termino desaparece en la ecuacion. Cuando usamos el estimador de Heckman, se puede hacer una simple prueba estadistica para ver si \(\beta_{\lambda}= 0\) que es equivalente a que \(\rho_{uv}=0\).
- La otra situacion en la cual no tenemos un sesgo con MCO es cuando \(\lambda_i\) no tiene correlacion con las xs. Si consideramos el ejemplo del mercado laboral donde queremos estimar el salario de mujeres. La condicion implicaria que la probabilidad de trabajar (margen extensivo) es independiente de las caracteristicas que determinan el salario. Tambien se puede decir que los determinantes de la decision de trabajar son independientes de los factores que explican el salario! Cosa que no es muy probable o si?!
#Ejemplo: Ecuacion de la oferta salarial para mujeres casadas Se aplica la correccion de seleccion muestral a los datos sobre las mujeres casadas en MROZ.RAW.
###Paquete \(\texttt{sampleSelection}\)
#install.packages("sampleSelection") #Para Heckit
#install.packages("mvtnorm")
#install.packages("wooldridge")
library("wooldridge")## Warning: package 'wooldridge' was built under R version 4.0.5## Warning: package 'sampleSelection' was built under R version 4.0.5Vease documentacion en R: https://www.rdocumentation.org/packages/sampleSelection/versions/1.2-6
#Datos
## [1] 753 22## 'data.frame': 753 obs. of 22 variables:
## $ inlf : int 1 1 1 1 1 1 1 1 1 1 ...
## $ hours : int 1610 1656 1980 456 1568 2032 1440 1020 1458 1600 ...
## $ kidslt6 : int 1 0 1 0 1 0 0 0 0 0 ...
## $ kidsge6 : int 0 2 3 3 2 0 2 0 2 2 ...
## $ age : int 32 30 35 34 31 54 37 54 48 39 ...
## $ educ : int 12 12 12 12 14 12 16 12 12 12 ...
## $ wage : num 3.35 1.39 4.55 1.1 4.59 ...
## $ repwage : num 2.65 2.65 4.04 3.25 3.6 ...
## $ hushrs : int 2708 2310 3072 1920 2000 1040 2670 4120 1995 2100 ...
## $ husage : int 34 30 40 53 32 57 37 53 52 43 ...
## $ huseduc : int 12 9 12 10 12 11 12 8 4 12 ...
## $ huswage : num 4.03 8.44 3.58 3.54 10 ...
## $ faminc : num 16310 21800 21040 7300 27300 ...
## $ mtr : num 0.721 0.661 0.692 0.781 0.622 ...
## $ motheduc: int 12 7 12 7 12 14 14 3 7 7 ...
## $ fatheduc: int 7 7 7 7 14 7 7 3 7 7 ...
## $ unem : num 5 11 5 5 9.5 7.5 5 5 3 5 ...
## $ city : int 0 1 0 0 1 1 0 0 0 0 ...
## $ exper : int 14 5 15 6 7 33 11 35 24 21 ...
## $ nwifeinc: num 10.9 19.5 12 6.8 20.1 ...
## $ lwage : num 1.2102 0.3285 1.5141 0.0921 1.5243 ...
## $ expersq : int 196 25 225 36 49 1089 121 1225 576 441 ...
## - attr(*, "time.stamp")= chr "25 Jun 2011 23:03"##
## 0 1
## 325 428Un data.frame con 753 observaciones de 22 variables:
inlf: =1 if in lab frce, 1975
hours: hours worked, 1975
kidslt6: # kids < 6 years
kidsge6: # kids 6-18
age: woman’s age in yrs
educ: years of schooling
wage: est. wage from earn, hrs
repwage: rep. wage at interview in 1976
hushrs: hours worked by husband, 1975
husage: husband’s age
huseduc: husband’s years of schooling
huswage: husband’s hourly wage, 1975
faminc: family income, 1975
mtr: fed. marg. tax rte facing woman
motheduc: mother’s years of schooling
fatheduc: father’s years of schooling
unem: unem. rate in county of resid.
city: =1 if live in SMSA
exper: actual labor mkt exper
nwifeinc: (faminc - wage*hours)/1000
lwage: log(wage)
expersq: exper^2
Ecuacion de interes y Ecuacion de seleccion
En este ejemplo la ecuacion de oferta salarial es la estandar, con log(wage) como la variable dependiente y educ, exper y \(exper^2\) como las variables explicativas.
Con el fin de probar y corregir el sesgo de seleccion muestral, debido a la inobservabilidad de la oferta salarial para mujeres que no estan trabajando, se necesita estimar un modelo probit para la participacion en la fuerza laboral. Ademas de las variables de educacion y experiencia, se incluyen los factores otros ingresos, edad,numero de hijos pequenos y numero de hijos mayores. El hecho de que estas cuatro variables esten excluidas de la ecuacion de oferta salarial es un supuesto: se supone que nwifeinc, age, kidslt6 y kidsge6 no tienen efecto sobre la oferta salarial.
#Estimacion ## Forma 1: Metodo de dos pasos o “Heckman’s two-step procedure”
HeckmanTS<-selection(inlf ~ educ+ exper + I(exper^2)+age + nwifeinc + kidslt6 + kidsge6, lwage ~ educ+ exper + I(exper^2),
data = mroz, method = "2step")
summary(HeckmanTS)## --------------------------------------------
## Tobit 2 model (sample selection model)
## 2-step Heckman / heckit estimation
## 753 observations (325 censored and 428 observed)
## 15 free parameters (df = 739)
## Probit selection equation:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.270077 0.508593 0.531 0.59556
## educ 0.130905 0.025254 5.183 0.000000281331104 ***
## exper 0.123348 0.018716 6.590 0.000000000083424 ***
## I(exper^2) -0.001887 0.000600 -3.145 0.00173 **
## age -0.052853 0.008477 -6.235 0.000000000761007 ***
## nwifeinc -0.012024 0.004840 -2.484 0.01320 *
## kidslt6 -0.868328 0.118522 -7.326 0.000000000000621 ***
## kidsge6 0.036005 0.043477 0.828 0.40786
## Outcome equation:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5781032 0.3050062 -1.895 0.05843 .
## educ 0.1090655 0.0155230 7.026 0.00000000000483 ***
## exper 0.0438873 0.0162611 2.699 0.00712 **
## I(exper^2) -0.0008591 0.0004389 -1.957 0.05068 .
## Multiple R-Squared:0.1569, Adjusted R-Squared:0.149
## Error terms:
## Estimate Std. Error t value Pr(>|t|)
## invMillsRatio 0.03226 0.13362 0.241 0.809
## sigma 0.66363 NA NA NA
## rho 0.04861 NA NA NA
## --------------------------------------------###Comentario: No hay evidencia de un problema de seleccion muestral en la estimacion de la ecuacion de oferta salarial. El coeficiente de \(\hat{\lambda}\) tiene un estadistico \(t\) muy pequeno \((0.241)\), asi que no se puede rechazar \(H_0: \rho=0\).
Forma 2: Metodo de Maxima Verosimilitud
###################################################
HeckmanML<-selection(inlf ~ educ+ exper + I(exper^2)+age + nwifeinc + kidslt6 + kidsge6,lwage ~ educ+ exper + I(exper^2),data = mroz, method = "ml")
summary(HeckmanML)## --------------------------------------------
## Tobit 2 model (sample selection model)
## Maximum Likelihood estimation
## Newton-Raphson maximisation, 3 iterations
## Return code 1: gradient close to zero (gradtol)
## Log-Likelihood: -832.8851
## 753 observations (325 censored and 428 observed)
## 14 free parameters (df = 739)
## Probit selection equation:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2664491 0.5089578 0.524 0.60077
## educ 0.1313414 0.0253823 5.175 0.000000294669116 ***
## exper 0.1232818 0.0187242 6.584 0.000000000086808 ***
## I(exper^2) -0.0018863 0.0006004 -3.142 0.00175 **
## age -0.0528287 0.0084792 -6.230 0.000000000780903 ***
## nwifeinc -0.0121321 0.0048767 -2.488 0.01307 *
## kidslt6 -0.8673987 0.1186509 -7.311 0.000000000000693 ***
## kidsge6 0.0358724 0.0434753 0.825 0.40957
## Outcome equation:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5526963 0.2603785 -2.123 0.0341 *
## educ 0.1083502 0.0148607 7.291 0.000000000000793 ***
## exper 0.0428368 0.0148785 2.879 0.0041 **
## I(exper^2) -0.0008374 0.0004175 -2.006 0.0452 *
## Error terms:
## Estimate Std. Error t value Pr(>|t|)
## sigma 0.66340 0.02271 29.215 <0.0000000000000002 ***
## rho 0.02661 0.14708 0.181 0.856
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------###Comentario Como una alternativa al metodo anterior de estimacion de dos pasos esta la estimacion completa de maxima verosimilitud. Es mas complicado puesto que requiere obtener la distribucion conjunta de y y s.
Usando OLS
##
## Call:
## lm(formula = lwage ~ educ + exper + I(exper^2), data = mroz)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.08404 -0.30627 0.04952 0.37498 2.37115
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5220406 0.1986321 -2.628 0.00890 **
## educ 0.1074896 0.0141465 7.598 0.000000000000194 ***
## exper 0.0415665 0.0131752 3.155 0.00172 **
## I(exper^2) -0.0008112 0.0003932 -2.063 0.03974 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6664 on 424 degrees of freedom
## (325 observations deleted due to missingness)
## Multiple R-squared: 0.1568, Adjusted R-squared: 0.1509
## F-statistic: 26.29 on 3 and 424 DF, p-value: 0.000000000000001302##Comparamos MCO, Heckit y MV
## (Intercept) educ exper I(exper^2)
## -0.5220405615 0.1074896401 0.0415665091 -0.0008111931## (Intercept) educ exper I(exper^2)
## -0.5781031866 0.1090655213 0.0438873379 -0.0008591142## (Intercept) educ exper I(exper^2)
## -0.5526962913 0.1083501918 0.0428368191 -0.0008374258###Comentario Practicamente no hay grandes diferencias en los coeficientes asociados a las pendientes estimadas con los diversa tecnicas estimaciòn. Los rendimientos estimados de la educaciòn difieren por sòlo una decima parte de un punto porcentual.
## educ
## 0.001575881Confrontamos con stargazer:
##
## Please cite as:## Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.2. https://CRAN.R-project.org/package=stargazerstargazer(OLSw, HeckmanTS, HeckmanML,
title="Regresiòn salario mujeres casadas", type="text",
df=FALSE, digits=4)##
## Regresiòn salario mujeres casadas
## ==============================================================
## Dependent variable:
## ------------------------------------------
## lwage
## OLS selection
## (1) (2) (3)
## --------------------------------------------------------------
## educ 0.1075*** 0.1091*** 0.1084***
## (0.0141) (0.0155) (0.0149)
##
## exper 0.0416*** 0.0439*** 0.0428***
## (0.0132) (0.0163) (0.0149)
##
## I(exper2) -0.0008** -0.0009* -0.0008**
## (0.0004) (0.0004) (0.0004)
##
## Constant -0.5220*** -0.5781* -0.5527**
## (0.1986) (0.3050) (0.2604)
##
## --------------------------------------------------------------
## Observations 428 753 753
## R2 0.1568
## Adjusted R2 0.1509
## Log Likelihood -832.8851
## rho 0.0486 0.0266 (0.1471)
## Inverse Mills Ratio 0.0323 (0.1336)
## Residual Std. Error 0.6664
## F Statistic 26.2862***
## ==============================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Ejemplo Green 2008
Estimaciòn en dos pasos
summary( heckit( lfp ~ age + I( age^2 ) + faminc + kids + educ,
wage ~ exper + I( exper^2 ) + educ + city, Mroz87 ) )## --------------------------------------------
## Tobit 2 model (sample selection model)
## 2-step Heckman / heckit estimation
## 753 observations (325 censored and 428 observed)
## 14 free parameters (df = 740)
## Probit selection equation:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.156806923 1.402085958 -2.965 0.003127 **
## age 0.185395096 0.065966659 2.810 0.005078 **
## I(age^2) -0.002425897 0.000773540 -3.136 0.001780 **
## faminc 0.000004580 0.000004206 1.089 0.276544
## kidsTRUE -0.448986740 0.130911496 -3.430 0.000638 ***
## educ 0.098182281 0.022984120 4.272 0.0000219 ***
## Outcome equation:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.9712003 2.0593505 -0.472 0.637
## exper 0.0210610 0.0624646 0.337 0.736
## I(exper^2) 0.0001371 0.0018782 0.073 0.942
## educ 0.4170174 0.1002497 4.160 0.0000356 ***
## city 0.4438379 0.3158984 1.405 0.160
## Multiple R-Squared:0.1264, Adjusted R-Squared:0.116
## Error terms:
## Estimate Std. Error t value Pr(>|t|)
## invMillsRatio -1.098 1.266 -0.867 0.386
## sigma 3.200 NA NA NA
## rho -0.343 NA NA NA
## --------------------------------------------Estimaciòn por maxima verosimilitud
summary( selection( lfp ~ age + I( age^2 ) + faminc + kids + educ,
wage ~ exper + I( exper^2 ) + educ + city, Mroz87 ) )## --------------------------------------------
## Tobit 2 model (sample selection model)
## Maximum Likelihood estimation
## Newton-Raphson maximisation, 5 iterations
## Return code 8: successive function values within relative tolerance limit (reltol)
## Log-Likelihood: -1581.258
## 753 observations (325 censored and 428 observed)
## 13 free parameters (df = 740)
## Probit selection equation:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.119691981 1.400516371 -2.942 0.003368 **
## age 0.184015424 0.065867312 2.794 0.005345 **
## I(age^2) -0.002408697 0.000772297 -3.119 0.001886 **
## faminc 0.000005680 0.000004416 1.286 0.198782
## kidsTRUE -0.450614870 0.130185426 -3.461 0.000568 ***
## educ 0.095280799 0.023153419 4.115 0.000043 ***
## Outcome equation:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.9630242 1.1982209 -1.638 0.102
## exper 0.0278683 0.0615514 0.453 0.651
## I(exper^2) -0.0001039 0.0018388 -0.056 0.955
## educ 0.4570051 0.0732299 6.241 0.000000000733 ***
## city 0.4465290 0.3159209 1.413 0.158
## Error terms:
## Estimate Std. Error t value Pr(>|t|)
## sigma 3.1084 0.1138 27.307 <0.0000000000000002 ***
## rho -0.1320 0.1651 -0.799 0.424
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------Ejemplo con una variable binaria como Y para un modelo de selecciòn
Estimamos la probabilidad de educación de las mujeres sobre sus posibilidades de obtener un salario alto (> $ 5 / h en 1975 USD), utilizando datos de PSID. utilizaremos educaciòn como una variable explicativa y añadiremos la edad, numero de niñes, y entradas no de trabajo como restricciones de exclusion.
data(Mroz87)
m <- selection(lfp ~ educ + age + kids5 + kids618 + nwifeinc,
wage >= 5 ~ educ, data = Mroz87 )
summary(m)## --------------------------------------------
## Tobit 2 model with binary outcome (sample selection model)
## Maximum Likelihood estimation
## BHHH maximisation, 8 iterations
## Return code 8: successive function values within relative tolerance limit (reltol)
## Log-Likelihood: -653.2037
## 753 observations (325 censored and 428 observed)
## 9 free parameters (df = 744)
## Probit selection equation:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.430362 0.475966 0.904 0.366
## educ 0.156223 0.023811 6.561 0.00000000010018220 ***
## age -0.034713 0.007649 -4.538 0.00000661066866332 ***
## kids5 -0.890560 0.112663 -7.905 0.00000000000000969 ***
## kids618 -0.038167 0.039320 -0.971 0.332
## nwifeinc -0.020948 0.004318 -4.851 0.00000149356272461 ***
## Outcome equation:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.5213 0.5611 -8.058 0.00000000000000308 ***
## educ 0.2879 0.0369 7.800 0.00000000000002090 ***
## Error terms:
## Estimate Std. Error t value Pr(>|t|)
## rho 0.1164 0.2706 0.43 0.667
## --------------------------------------------Ejemplo utilizando numeros aleatorio
## Warning: package 'mvtnorm' was built under R version 4.0.5nObs <- 1000
sigma <- matrix( c( 1, -0.7, -0.7, 1 ), ncol = 2 )
errorTerms <- rmvnorm( nObs, c( 0, 0 ), sigma )
myData <- data.frame( no = c( 1:nObs ), x1 = rnorm( nObs ), x2 = rnorm( nObs ),
u1 = errorTerms[ , 1 ], u2 = errorTerms[ , 2 ] )
myData$y <- 2 + myData$x1 + myData$u1
myData$s <- ( 2 * myData$x1 + myData$x2 + myData$u2 - 0.2 ) > 0
myData$y[ !myData$s ] <- NA
myOls <- lm( y ~ x1, data = myData)
summary( myOls )##
## Call:
## lm(formula = y ~ x1, data = myData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.72495 -0.70243 0.03919 0.67678 2.71326
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.53595 0.05973 25.72 <0.0000000000000002 ***
## x1 1.35165 0.05700 23.71 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9664 on 471 degrees of freedom
## (527 observations deleted due to missingness)
## Multiple R-squared: 0.5442, Adjusted R-squared: 0.5432
## F-statistic: 562.3 on 1 and 471 DF, p-value: < 0.00000000000000022## Tobit 2 model## --------------------------------------------
## Tobit 2 model (sample selection model)
## 2-step Heckman / heckit estimation
## 1000 observations (527 censored and 473 observed)
## 8 free parameters (df = 993)
## Probit selection equation:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2170 0.0584 -3.715 0.000214 ***
## x1 1.8337 0.1103 16.625 < 0.0000000000000002 ***
## x2 0.9555 0.0752 12.706 < 0.0000000000000002 ***
## Outcome equation:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.11248 0.10806 19.55 <0.0000000000000002 ***
## x1 1.00408 0.08015 12.53 <0.0000000000000002 ***
## Multiple R-Squared:0.5858, Adjusted R-Squared:0.584
## Error terms:
## Estimate Std. Error t value Pr(>|t|)
## invMillsRatio -0.9109 0.1294 -7.037 0.00000000000366 ***
## sigma 1.0537 NA NA NA
## rho -0.8645 NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------Ejemplo utilizand numeros aleatorios con una estimaciòn 2SLS
library( "mvtnorm" )
nObs <- 1000
sigma <- matrix( c( 1, 0.5, 0.1, 0.5, 1, -0.3, 0.1, -0.3, 1 ), ncol = 3 )
errorTerms <- rmvnorm( nObs, c( 0, 0, 0 ), sigma )
myData <- data.frame( no = c( 1:nObs ), x1 = rnorm( nObs ), x2 = rnorm( nObs ),
u1 = errorTerms[ , 1 ], u2 = errorTerms[ , 2 ], u3 = errorTerms[ , 3 ] )
myData$w <- 1 + myData$x1 + myData$u1
myData$y <- 2 + myData$w + myData$u2
myData$s <- ( 2 * myData$x1 + myData$x2 + myData$u3 - 0.2 ) > 0
myData$y[ !myData$s ] <- NA
myHeckit <- heckit( s ~ x1 + x2, y ~ w, data = myData )
summary( myHeckit ) ## --------------------------------------------
## Tobit 2 model (sample selection model)
## 2-step Heckman / heckit estimation
## 1000 observations (525 censored and 475 observed)
## 8 free parameters (df = 993)
## Probit selection equation:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.10310 0.05690 -1.812 0.0703 .
## x1 1.80891 0.10945 16.528 <0.0000000000000002 ***
## x2 0.91502 0.07534 12.146 <0.0000000000000002 ***
## Outcome equation:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.19327 0.09431 12.65 <0.0000000000000002 ***
## w 1.35630 0.03652 37.14 <0.0000000000000002 ***
## Multiple R-Squared:0.775, Adjusted R-Squared:0.774
## Error terms:
## Estimate Std. Error t value Pr(>|t|)
## invMillsRatio 0.03004 0.10256 0.293 0.77
## sigma 0.87328 NA NA NA
## rho 0.03440 NA NA NA
## --------------------------------------------Es sesgadpo, para corregir
## --------------------------------------------
## Tobit 2 model (sample selection model)
## 2-step Heckman / heckit estimation
## 1000 observations (525 censored and 475 observed)
## 8 free parameters (df = 993)
## Probit selection equation:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.10310 0.05690 -1.812 0.0703 .
## x1 1.80891 0.10945 16.528 <0.0000000000000002 ***
## x2 0.91502 0.07534 12.146 <0.0000000000000002 ***
## Outcome equation:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.21705 0.11113 19.95 <0.0000000000000002 ***
## w 0.87572 0.04286 20.43 <0.0000000000000002 ***
## Multiple R-Squared:0.6929, Adjusted R-Squared:0.6916
## Error terms:
## Estimate Std. Error t value Pr(>|t|)
## invMillsRatio -0.5133 0.1185 -4.33 0.0000164 ***
## sigma 1.0616 NA NA NA
## rho -0.4835 NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------Hemos corregido el sesgo
Ejemplo tobit
N <- 500
library(mvtnorm)
vc <- diag(3)
vc[lower.tri(vc)] <- c(0.9, 0.5, 0.6)
vc[upper.tri(vc)] <- vc[lower.tri(vc)]
eps <- rmvnorm(N, rep(0, 3), vc)
xs <- runif(N)
ys <- xs + eps[,1] > 0
xo1 <- runif(N)
yo1 <- xo1 + eps[,2]
xo2 <- runif(N)
yo2 <- xo2 + eps[,3]
a <- selection(ys~xs, list(yo1 ~ xo1, yo2 ~ xo2))
summary(a)## --------------------------------------------
## Tobit 5 model (switching regression model)
## Maximum Likelihood estimation
## Newton-Raphson maximisation, 13 iterations
## Return code 8: successive function values within relative tolerance limit (reltol)
## Log-Likelihood: -875.5618
## 500 observations: 146 selection 1 (FALSE) and 354 selection 2 (TRUE)
## 10 free parameters (df = 490)
## Probit selection equation:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.01261 0.09841 0.128 0.898
## xs 1.10186 0.15617 7.055 0.00000000000591 ***
## Outcome equation 1:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.00186 0.16275 -0.011 0.991
## xo1 1.21147 0.15952 7.595 0.000000000000158 ***
## Outcome equation 2:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1005 0.1681 0.598 0.55
## xo2 1.0614 0.1656 6.410 0.000000000342 ***
## Error terms:
## Estimate Std. Error t value Pr(>|t|)
## sigma1 1.03365 0.09765 10.585 <0.0000000000000002 ***
## sigma2 0.90042 0.05252 17.145 <0.0000000000000002 ***
## rho1 0.94996 0.02531 37.539 <0.0000000000000002 ***
## rho2 0.28230 0.32700 0.863 0.388
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------Ejemplo por una regresiòn con tratamiento
Estimamos el efecto del tratamiento sobre los ingresos
- resutlados de la seleccion: la participacion al tratamiento Ñ variables explicativas de seleccion: edad, educacion,, desempleo en el 74 y el 75, raza
- \(Y\): Logaritmo de los ingresos reales en 1978
- Variables explicativas por la \(Y\):tratamiento, edad, educacion, raza. Las variablesde desempleo son utilizadas como restriccion de explucion.
data(Treatment, package="Ecdat")
a <- heckit(treat~poly(age,2) + educ + u74 + u75 + ethn,
log(re78)~treat + poly(age,2) + educ + ethn,
data=Treatment)
print(summary(a))
b<-lm(log(re78)~treat + poly(age,2) + educ + ethn,
data=Treatment)
stargazer::stargazer(a,b,type = text)Tobit
El modelo tobit, también llamado modelo de regresión censurado, está diseñado para estimar relaciones lineales entre variables cuando hay censura por la izquierda o por la derecha en la variable dependiente (también conocida como censura desde abajo y arriba, respectivamente). La censura desde arriba tiene lugar cuando los casos con un valor en o por encima de algún umbral, todos toman el valor de ese umbral, de modo que el valor real puede ser igual al umbral, pero también puede ser mayor. En el caso de censurar desde abajo, se censuran valores aquellos que caen en o por debajo de algún umbral.
## Loading required package: GGally## Warning: package 'GGally' was built under R version 4.0.5## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2## Loading required package: VGAM## Warning: package 'VGAM' was built under R version 4.0.5## Loading required package: stats4## Loading required package: splines##
## Attaching package: 'VGAM'## The following object is masked from 'package:sampleSelection':
##
## probit## The following object is masked from 'package:wooldridge':
##
## wine## The following objects are masked from 'package:boot':
##
## logit, simplex## The following objects are masked from 'package:mlogit':
##
## has.intercept, lrtest## The following object is masked from 'package:modelsummary':
##
## Max## The following object is masked from 'package:tidyr':
##
## fillEjemplos
Ejemplo 1. En la década de 1980, existía una ley federal que restringía las lecturas del velocímetro a no más de 85 mph. Entonces, si quisiera intentar predecir la velocidad máxima de un vehículo a partir de una combinación de caballos de fuerza y tamaño del motor, obtendría una lectura no superior a 85, independientemente de qué tan rápido se desplazara realmente el vehículo. Este es un caso clásico de censura por la derecha (censura desde arriba) de los datos. De lo único que estamos seguros es de que esos vehículos viajaban al menos a 85 mph.
Ejemplo 2. Un proyecto de investigación está estudiando el nivel de plomo en el agua potable en el hogar en función de la edad de la casa y los ingresos familiares. El kit de análisis de agua no puede detectar concentraciones de plomo por debajo de 5 partes por mil millones (ppb). La EPA considera peligrosos los niveles superiores a 15 ppb. Estos datos son un ejemplo de censura por la izquierda (censura desde abajo).
Ejemplo 3. Considere la situación en la que tenemos una medida de aptitud académica (en una escala de 200-800) que queremos modelar utilizando los puntajes de las pruebas de lectura y matemáticas, así como el tipo de programa en el que está inscrito el estudiante (académico, general , o vocacional). El problema aquí es que los estudiantes que responden correctamente a todas las preguntas de la prueba de aptitud académica reciben una puntuación de 800, aunque es probable que estos estudiantes no sean “verdaderamente” iguales en aptitud. Lo mismo ocurre con los estudiantes que responden incorrectamente a todas las preguntas. Todos estos estudiantes tendrían una puntuación de 200, aunque es posible que no todos tengan las mismas aptitudes.
datos
El conjunto de datos contiene 200 observaciones. La variable de aptitud académica es apta, los resultados de las pruebas de lectura y matemáticas son lectura y matemáticas respectivamente. La variable prog es el tipo de programa en el que se encuentra el estudiante, es una variable categórica (nominal) que toma tres valores, académico (prog = 1), general (prog = 2) y vocacional (prog = 3). La variable id es una variable de identificación.
Ahora veamos los datos de forma descriptiva. Tenga en cuenta que en este conjunto de datos, el valor más bajo de apt es 352. Es decir, ningún estudiante recibió una puntuación de 200 (la puntuación más baja posible), lo que significa que aunque fue posible censurar desde abajo, no ocurre en el conjunto de datos.
## id read math prog
## Min. : 1.00 Min. :28.00 Min. :33.00 Length:200
## 1st Qu.: 50.75 1st Qu.:44.00 1st Qu.:45.00 Class :character
## Median :100.50 Median :50.00 Median :52.00 Mode :character
## Mean :100.50 Mean :52.23 Mean :52.65
## 3rd Qu.:150.25 3rd Qu.:60.00 3rd Qu.:59.00
## Max. :200.00 Max. :76.00 Max. :75.00
## apt
## Min. :352.0
## 1st Qu.:575.5
## Median :633.0
## Mean :640.0
## 3rd Qu.:705.2
## Max. :800.0analisis grafica:
f <- function(x, var, bw = 15) {
dnorm(x, mean = mean(var), sd(var)) * length(var) * bw
}
p <- ggplot(dat, aes(x = apt, fill=prog))
p + stat_bin(binwidth=15) +
stat_function(fun = f, size = 1,
args = list(var = dat$apt))## Warning: Multiple drawing groups in `geom_function()`. Did you use the correct
## `group`, `colour`, or `fill` aesthetics?
Al observar el histograma anterior, podemos ver la censura en los valores de apt, es decir, hay muchos más casos con puntuaciones de 750 a 800 de los que cabría esperar al observar el resto de la distribución. A continuación se muestra un histograma alternativo que resalta aún más el exceso de casos donde apt = 800. En el histograma a continuación, la opción de rupturas produce un histograma donde cada valor único de apt tiene su propia barra (estableciendo rupturas iguales a un vector que contiene valores desde el mínimo de apt hasta el máximo de apt). Debido a que apt es continuo, la mayoría de los valores de apt son únicos en el conjunto de datos, aunque cerca del centro de la distribución hay algunos valores de apt que tienen dos o tres casos. El pico en el extremo derecho del histograma es la barra para los casos en los que apt = 800, la altura de esta barra en relación con todas las demás muestra claramente el exceso de casos con este valor.
## Warning: Multiple drawing groups in `geom_function()`. Did you use the correct
## `group`, `colour`, or `fill` aesthetics?
Analizamos las relaciones bivariadas en nuestro dataset:
## read math apt
## read 1.0000000 0.6622801 0.6451215
## math 0.6622801 1.0000000 0.7332702
## apt 0.6451215 0.7332702 1.0000000Con ggplot:
En la primera fila de la matriz de diagramas de dispersión que se muestra arriba, vemos los diagramas de dispersión que muestran la relación entre read y apt, así como también entre math y apt. Tenga en cuenta la colección de casos en la parte superior de estos dos diagramas de dispersión, esto se debe a la censura en la distribución de apt.
Utilizamos el Tobit:
##
## Call:
## vglm(formula = apt ~ read + math + prog, family = tobit(Upper = 800),
## data = dat)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 209.55956 32.54590 6.439 0.00000000012 ***
## (Intercept):2 4.18476 0.05235 79.944 < 0.0000000000000002 ***
## read 2.69796 0.61928 4.357 0.00001320835 ***
## math 5.91460 0.70539 8.385 < 0.0000000000000002 ***
## proggeneral -12.71458 12.40857 -1.025 0.305523
## progvocational -46.14327 13.70667 -3.366 0.000761 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: mu, loglink(sd)
##
## Log-likelihood: -1041.063 on 394 degrees of freedom
##
## Number of Fisher scoring iterations: 5
##
## No Hauck-Donner effect found in any of the estimatesLa tabla denominada coeficientes proporciona los coeficientes, sus errores estándar y el estadístico z.
Los coeficientes de regresión de Tobit se interpretan de manera similar a los coeficientes de regresión de OLS; sin embargo, el efecto lineal está en la variable latente sin censura, no en el resultado observado.
Para un aumento de una unidad en la lectura, hay un aumento de 2.6981 puntos en el valor predicho de apto.
Un aumento de una unidad en matemáticas se asocia con un aumento de 5,9146 unidades en el valor predicho de apto.
Los términos para prog tienen una interpretación ligeramente diferente. El valor previsto de apt es -46,1419 puntos más bajo para los estudiantes en un programa vocacional que para los estudiantes en un programa académico.
El coeficiente etiquetado como “(Intercepción): 1” es la intersección o constante del modelo.
El coeficiente etiquetado como “(Intercepción): 2” es una estadística auxiliar. Si exponenciamos este valor, obtenemos una estadística que es análoga a la raíz cuadrada de la varianza residual en la regresión MCO. El valor de 65.6773 puede comparado con la desviación estándar de aptitud académica que fue 99.21, una reducción sustancial.
La probabilidad logarítmica final, -1041.0629, se muestra hacia la parte inferior de la salida, se puede usar en comparaciones de modelos anidados.
A continuación, calculamos los valores p para cada uno de los coeficientes del modelo. Calculamos el valor p para cada coeficiente usando los valores z y luego lo mostramos en una tabla con los coeficientes. Los coeficientes para lectura, matemáticas y prog = 3 (vocacional) son estadísticamente significativos.
ctable <- coef(summary(m))
pvals <- 2 * pt(abs(ctable[, "z value"]), df.residual(m), lower.tail = FALSE)
cbind(ctable, pvals)## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 209.559557 32.54589921 6.438893 0.00000000012034812930792
## (Intercept):2 4.184759 0.05234618 79.943922 0.00000000000000000000000
## read 2.697959 0.61927743 4.356625 0.00001320834607070725204
## math 5.914596 0.70538721 8.384892 0.00000000000000005077232
## proggeneral -12.714581 12.40856959 -1.024661 0.30552301628167571889705
## progvocational -46.143271 13.70667208 -3.366482 0.00076133432325494214063
## pvals
## (Intercept):1 0.00000000035058389605921694833975954530558283295249566435813903808593750000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## (Intercept):2 0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001299833
## read 0.00001686815319946651512469232292446008614206220954656600952148437500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## math 0.00000000000000091224343995323472862254976512019766232697293162345886230468750000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## proggeneral 0.30615170215254988717035189438320230692625045776367187500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
## progvocational 0.00083619122671922717707215788607300055446103215217590332031250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000Podemos probar la significancia del tipo de programa en general ajustando un modelo sin programa en él y usando una prueba de razón de verosimilitud.
m2 <- vglm(apt ~ read + math, tobit(Upper = 800), data = dat)
(p <- pchisq(2 * (logLik(m) - logLik(m2)), df = 2, lower.tail = FALSE))## [1] 0.003155176El LRT con dos grados de libertad está asociado con un valor p de 0,0032, lo que indica que el efecto general del prog es estadísticamente significativo.
Podemos caclular intervalo de confianzas al 95%:
b <- coef(m)
se <- sqrt(diag(vcov(m)))
cbind(LL = b - qnorm(0.975) * se, UL = b + qnorm(0.975) * se)## LL UL
## (Intercept):1 145.770767 273.348348
## (Intercept):2 4.082163 4.287356
## read 1.484198 3.911721
## math 4.532062 7.297129
## proggeneral -37.034931 11.605768
## progvocational -73.007854 -19.278687Y hacer una anlisis grafica:
dat$yhat <- fitted(m)[,1]
dat$rr <- resid(m, type = "response")
dat$rp <- resid(m, type = "pearson")[,1]
par(mfcol = c(2, 3))
with(dat, {
plot(yhat, rr, main = "Fitted vs Residuals")
qqnorm(rr)
plot(yhat, rp, main = "Fitted vs Pearson Residuals")
qqnorm(rp)
plot(apt, rp, main = "Actual vs Pearson Residuals")
plot(apt, yhat, main = "Actual vs Fitted")
})
El gráfico de la parte inferior derecha muestra los valores predichos, o ajustados, graficados contra los reales. Esto puede resultar particularmente útil al comparar modelos. Podemos calcular la correlación entre estos dos, así como la correlación al cuadrado, para tener una idea de cuán preciso nuestro modelo predice los datos y qué parte de la varianza en el resultado se explica por el modelo.
#Referencias: 1. https://stats.idre.ucla.edu/r/dae/truncated-regression/