Laboratorio 8
Federico Daverio
rm(list=ls())
setwd("C:/DAVE2/CIDE/Doc/Econometria II/Lab_7")
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: msmRepaso
OLS
datos<-as_tibble(wooldridge::wage1) %>%
rename(SalarioHora=wage,
Educación=educ,
Experiencia=exper,
Antiguedad=tenure,
Género=female,
EdoCivil=married) %>%
mutate(Género = recode_factor(Género, `1` = "Mujer", `0` = "Hombre"),
EdoCivil = recode_factor(EdoCivil, `1` = "Casad@", `0` = "Solter@"))Estimacion MCO
Añadimos màs regresores:
##
## Call:
## lm(formula = SalarioHora ~ Género + Educación, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.9890 -1.8702 -0.6651 1.0447 15.4998
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.65055 0.65232 -2.530 0.0117 *
## GéneroHombre 2.27336 0.27904 8.147 0.00000000000000276 ***
## Educación 0.50645 0.05039 10.051 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.186 on 523 degrees of freedom
## Multiple R-squared: 0.2588, Adjusted R-squared: 0.256
## F-statistic: 91.32 on 2 and 523 DF, p-value: < 0.00000000000000022Comparamos modelos:
##
## ===================================================================
## Dependent variable:
## -----------------------------------------------
## SalarioHora
## (1) (2)
## -------------------------------------------------------------------
## GéneroHombre 2.512*** 2.273***
## (0.303) (0.279)
##
## Educación 0.506***
## (0.050)
##
## Constant 4.588*** -1.651**
## (0.219) (0.652)
##
## -------------------------------------------------------------------
## Observations 526 526
## R2 0.116 0.259
## Adjusted R2 0.114 0.256
## Residual Std. Error 3.476 (df = 524) 3.186 (df = 523)
## F Statistic 68.537*** (df = 1; 524) 91.315*** (df = 2; 523)
## ===================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Cambiamos el nivel de referencia:
datos <- datos %>%
mutate(Género = relevel(Género, ref = "Hombre"),EdoCivil=relevel(EdoCivil,ref="Solter@"))
modelo1.1<-lm(SalarioHora~Género+Educación,data = datos)
summary(modelo1.1)##
## Call:
## lm(formula = SalarioHora ~ Género + Educación, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.9890 -1.8702 -0.6651 1.0447 15.4998
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.62282 0.67253 0.926 0.355
## GéneroMujer -2.27336 0.27904 -8.147 0.00000000000000276 ***
## Educación 0.50645 0.05039 10.051 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.186 on 523 degrees of freedom
## Multiple R-squared: 0.2588, Adjusted R-squared: 0.256
## F-statistic: 91.32 on 2 and 523 DF, p-value: < 0.00000000000000022Obtenemos el grafico:
interceptos <- c(coef(modelo1.1)["(Intercept)"], #intercepto de hombres
coef(modelo1.1)["(Intercept)"] + coef(modelo1.1)["GéneroMujer"]) #intercepto de mujeres
rectas.df<- data.frame(interceptos = interceptos,
slopes = rep(coef(modelo1.1)["Educación"], 2),
Género = levels(datos$Género))Y lo ploteamos:
qplot(x = Educación, y = SalarioHora, color = Género, data = datos) +
geom_abline(aes(intercept = interceptos,
slope = slopes,
color = Género), data = rectas.df)
Erores std robustos y clusterizados
Calculamos los errores robustos:
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.622817 0.728684 0.8547 0.3931
## GéneroMujer -2.273362 0.270203 -8.4135 0.0000000000000003820 ***
## Educación 0.506452 0.059896 8.4556 0.0000000000000002784 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Comparamos:
##
## =======================================================
## Dependent variable:
## -----------------------------------
## SalarioHora
## OLS coefficient
## test
## (1) (2)
## -------------------------------------------------------
## GéneroMujer -2.273*** -2.273***
## (0.279) (0.270)
##
## Educación 0.506*** 0.506***
## (0.050) (0.060)
##
## Constant 0.623 0.623
## (0.673) (0.729)
##
## -------------------------------------------------------
## Observations 526
## R2 0.259
## Adjusted R2 0.256
## Residual Std. Error 3.186 (df = 523)
## F Statistic 91.315*** (df = 2; 523)
## =======================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Si queremos errores clusterizados:
## Warning: package 'clubSandwich' was built under R version 4.0.5## Registered S3 method overwritten by 'clubSandwich':
## method from
## bread.mlm sandwich## Coef. Estimate SE t-stat p-val (naive-t) Sig.
## 1 (Intercept) 0.623 0.5213 1.19 0.44365
## 2 GéneroMujer -2.273 0.0192 -118.44 0.00537 **
## 3 Educación 0.506 0.0408 12.42 0.05113 .Dependiente Binaria
El modelo Probit
mydata <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")
## convert rank to a factor (categorical variable)
mydata$rank <- factor(mydata$rank)
## view first few rows
head(mydata)## admit gre gpa rank
## 1 0 380 3.61 3
## 2 1 660 3.67 3
## 3 1 800 4.00 1
## 4 1 640 3.19 4
## 5 0 520 2.93 4
## 6 1 760 3.00 2Utilizamos el probit:
myprobit <- glm(admit ~ gre + gpa + rank, family = binomial(link = "probit"),
data = mydata)
## model summary
summary(myprobit)##
## Call:
## glm(formula = admit ~ gre + gpa + rank, family = binomial(link = "probit"),
## data = mydata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6163 -0.8710 -0.6389 1.1560 2.1035
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.386836 0.673946 -3.542 0.000398 ***
## gre 0.001376 0.000650 2.116 0.034329 *
## gpa 0.477730 0.197197 2.423 0.015410 *
## rank2 -0.415399 0.194977 -2.131 0.033130 *
## rank3 -0.812138 0.208358 -3.898 0.0000971 ***
## rank4 -0.935899 0.245272 -3.816 0.000136 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 499.98 on 399 degrees of freedom
## Residual deviance: 458.41 on 394 degrees of freedom
## AIC: 470.41
##
## Number of Fisher Scoring iterations: 4Tanto gre, gpa como los tres términos para rango son estadísticamente significativos. Los coeficientes de regresión probit dan el cambio en la puntuación z o índice probit para un cambio de una unidad en el predictor.
Para un aumento de una unidad en gre, la puntuación z aumenta en 0,001.
Por cada aumento de una unidad en gpa, la puntuación z aumenta en 0,478.
Las variables indicadoras de rango tienen una interpretación ligeramente diferente. Por ejemplo, haber asistido a una institución de pregrado de rango 2, en comparación con una institución con rango 1 (el grupo de referencia), disminuye el puntaje z en 0.415.
El mdoelo logit
mydata <- read.csv("https://stats.idre.ucla.edu/stat/data/binary.csv")
mydata$rank <- factor(mydata$rank)
mylogit <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial")
summary(mylogit)##
## Call:
## glm(formula = admit ~ gre + gpa + rank, family = "binomial",
## data = mydata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6268 -0.8662 -0.6388 1.1490 2.0790
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.989979 1.139951 -3.500 0.000465 ***
## gre 0.002264 0.001094 2.070 0.038465 *
## gpa 0.804038 0.331819 2.423 0.015388 *
## rank2 -0.675443 0.316490 -2.134 0.032829 *
## rank3 -1.340204 0.345306 -3.881 0.000104 ***
## rank4 -1.551464 0.417832 -3.713 0.000205 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 499.98 on 399 degrees of freedom
## Residual deviance: 458.52 on 394 degrees of freedom
## AIC: 470.52
##
## Number of Fisher Scoring iterations: 4A continuación, vemos los residuos de desviación, que son una medida del ajuste del modelo. Esta parte de la salida muestra la distribución de los residuos de desviación para casos individuales usados en el modelo. A continuación, analizamos cómo utilizar resúmenes de la estadística de desviación para evaluar el ajuste del modelo.
La siguiente parte del resultado muestra los coeficientes, sus errores estándar, el estadístico z (a veces llamado estadístico z de Wald) y los valores p asociados. Tanto gre como gpa son estadísticamente significativos, al igual que los tres términos para rango. Los coeficientes de regresión logística dan el cambio en las probabilidades logarítmicas del resultado para un aumento de una unidad en la variable predictora.
Por cada cambio de una unidad en gre, las probabilidades logarítmicas de admisión (frente a no admisión) aumentan en 0,002.
Para un aumento de una unidad en el gpa, las probabilidades logarítmicas de ser admitido en la escuela de posgrado aumentan en 0,804.
Las variables indicadoras de rango tienen una interpretación ligeramente diferente. Por ejemplo, haber asistido a una institución de pregrado con rango 2, en comparación con una institución con rango 1, cambia las probabilidades logarítmicas de admisión en -0,675.
Debajo de la tabla de coeficientes se encuentran los índices de ajuste, incluidos los residuales nulos y de desviación y el AIC. Más adelante mostramos un ejemplo de cómo puede usar estos valores para ayudar a evaluar el ajuste del modelo.
Caclualmos ahora los irr
## (Intercept) gre gpa rank2 rank3 rank4
## 0.0185001 1.0022670 2.2345448 0.5089310 0.2617923 0.2119375Ahora podemos decir que para un aumento de una unidad en el gpa, las probabilidades de ser admitido en la escuela de posgrado (en lugar de no ser admitido) aumentan en un factor de 2,23.
Multinomiales
Multinomial logit
No varian los regresores
library(foreign)
ml <- read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta")
test1<- multinom(prog ~ ses + write, data = ml)## # weights: 15 (8 variable)
## initial value 219.722458
## iter 10 value 179.985215
## final value 179.981726
## converged## Call:
## multinom(formula = prog ~ ses + write, data = ml)
##
## Coefficients:
## (Intercept) sesmiddle seshigh write
## academic -2.851973 0.5332914 1.1628257 0.05792480
## vocation 2.366097 0.8246384 0.1802176 -0.05567514
##
## Std. Errors:
## (Intercept) sesmiddle seshigh write
## academic 1.166437 0.4437319 0.5142215 0.02141092
## vocation 1.174251 0.4901237 0.6484508 0.02333135
##
## Residual Deviance: 359.9635
## AIC: 375.9635Vemos diferencias:
##
## ==============================================
## Dependent variable:
## ----------------------------
## academic vocation
## (1) (2)
## ----------------------------------------------
## sesmiddle 0.533 0.825*
## (0.444) (0.490)
##
## seshigh 1.163** 0.180
## (0.514) (0.648)
##
## write 0.058*** -0.056**
## (0.021) (0.023)
##
## Constant -2.852** 2.366**
## (1.166) (1.174)
##
## ----------------------------------------------
## Akaike Inf. Crit. 375.963 375.963
## ==============================================
## Note: *p<0.1; **p<0.05; ***p<0.01\(\beta_4^1\) mide el efecto de cambios en el score de escritura sobre el log de la razon de probabilidades de elegir el programa academico respecto a el programa general (si aumenta el score en write es mas probable que un alumno elija el programa academico en vez de el programa general)
\(\beta_4^2\) mide el efectode cambios en el score de escritura sobre el log de la razon de probabilidades de elegir el programa vocacional respecto a el programa general(si aumenta el score en write es menos probable que un alumno elija el programa vocacional en lugar de el programa general)
## (Intercept) sesmiddle seshigh write
## academic 0.05773033 1.704533 3.198960 1.0596353
## vocation 10.65572511 2.281056 1.197478 0.9458464Ante un cambio unitario en score de write,es un 5% mas probable que se elija un programa academico respecto aun programa general.
Logit condicioncinal
data("Electricity", package = "mlogit")
Electr <- mlogit.data(Electricity, id.var = "id", choice = "choice",varying = 3:26, shape = "wide", sep = "")
head(Electr,10)## ~~~~~~~
## first 10 observations out of 17232
## ~~~~~~~
## choice id alt pf cl loc wk tod seas chid idx
## 1 FALSE 1 1 7 5 0 1 0 0 1 1:1
## 2 FALSE 1 2 9 1 1 0 0 0 1 1:2
## 3 FALSE 1 3 0 0 0 0 0 1 1 1:3
## 4 TRUE 1 4 0 5 0 1 1 0 1 1:4
## 5 FALSE 1 1 7 0 0 1 0 0 2 2:1
## 6 FALSE 1 2 9 5 0 1 0 0 2 2:2
## 7 TRUE 1 3 0 1 1 0 1 0 2 2:3
## 8 FALSE 1 4 0 5 0 0 0 1 2 2:4
## 9 FALSE 1 1 9 5 0 0 0 0 3 3:1
## 10 FALSE 1 2 7 1 0 1 0 0 3 3:2
##
## ~~~ indexes ~~~~
## chid id alt
## 1 1 1 1
## 2 1 1 2
## 3 1 1 3
## 4 1 1 4
## 5 2 1 1
## 6 2 1 2
## 7 2 1 3
## 8 2 1 4
## 9 3 1 1
## 10 3 1 2
## indexes: 1, 1, 2##
## Call:
## mlogit(formula = choice ~ pf + cl + loc + wk + tod + seas | 0,
## data = Electr, method = "nr")
##
## Frequencies of alternatives:choice
## 1 2 3 4
## 0.22702 0.26393 0.23816 0.27089
##
## nr method
## 4 iterations, 0h:0m:2s
## g'(-H)^-1g = 9.34E-07
## gradient close to zero
##
## Coefficients :
## Estimate Std. Error z-value Pr(>|z|)
## pf -0.6252278 0.0232223 -26.924 < 0.00000000000000022 ***
## cl -0.1082991 0.0082442 -13.136 < 0.00000000000000022 ***
## loc 1.4422429 0.0505571 28.527 < 0.00000000000000022 ***
## wk 0.9955040 0.0447801 22.231 < 0.00000000000000022 ***
## tod -5.4627587 0.1837125 -29.735 < 0.00000000000000022 ***
## seas -5.8400308 0.1866779 -31.284 < 0.00000000000000022 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-Likelihood: -4958.6Modelos de conteo
poisson
p <- read.csv("https://stats.idre.ucla.edu/stat/data/poisson_sim.csv")
p <- within(p, {
prog <- factor(prog, levels=1:3, labels=c("General", "Academic",
"Vocational"))
id <- factor(id)
})
with(p, tapply(num_awards, prog, function(x) {
sprintf("M (SD) = %1.2f (%1.2f)", mean(x), sd(x))
}))## General Academic Vocational
## "M (SD) = 0.20 (0.40)" "M (SD) = 1.00 (1.28)" "M (SD) = 0.24 (0.52)"
HAcemos la regresiòn:
##
## Call:
## glm(formula = num_awards ~ prog + math, family = "poisson", data = p)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2043 -0.8436 -0.5106 0.2558 2.6796
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.24712 0.65845 -7.969 0.0000000000000016 ***
## progAcademic 1.08386 0.35825 3.025 0.00248 **
## progVocational 0.36981 0.44107 0.838 0.40179
## math 0.07015 0.01060 6.619 0.0000000000362501 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 287.67 on 199 degrees of freedom
## Residual deviance: 189.45 on 196 degrees of freedom
## AIC: 373.5
##
## Number of Fisher Scoring iterations: 6El coeficiente de matemáticas es .07. Esto significa que el recuento de registros esperado para un aumento de una unidad en matemáticas es .07.
La variable indicadora progAcademic se compara entre prog = “Académico” y prog = “General”, el recuento de registros esperado para prog = “Académico” aumenta en aproximadamente 1,1. La variable indicadora prog.Vocational es la diferencia esperada en el recuento de registros ((aproximadamente .37)) entre prog = “Vocational” y el grupo de referencia (prog = “General”).
Bn1 y 2
##
## Attaching package: 'MASS'## The following object is masked from 'package:dfidx':
##
## select## The following object is masked from 'package:dplyr':
##
## selectdat <- read.dta("https://stats.idre.ucla.edu/stat/stata/dae/nb_data.dta")
dat <- within(dat, {
prog <- factor(prog, levels = 1:3, labels = c("General", "Academic", "Vocational"))
id <- factor(id)
})
summary(m1 <- glm.nb(daysabs ~ math + prog, data = dat))##
## Call:
## glm.nb(formula = daysabs ~ math + prog, data = dat, init.theta = 1.032713156,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1547 -1.0192 -0.3694 0.2285 2.5273
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.615265 0.197460 13.245 < 0.0000000000000002 ***
## math -0.005993 0.002505 -2.392 0.0167 *
## progAcademic -0.440760 0.182610 -2.414 0.0158 *
## progVocational -1.278651 0.200720 -6.370 0.000000000189 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.0327) family taken to be 1)
##
## Null deviance: 427.54 on 313 degrees of freedom
## Residual deviance: 358.52 on 310 degrees of freedom
## AIC: 1741.3
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.033
## Std. Err.: 0.106
##
## 2 x log-likelihood: -1731.258La variable matemática tiene un coeficiente de -0,006, que es estadísticamente significativo. Esto significa que por cada aumento de una unidad en matemáticas, el recuento de registros esperado del número de días ausentes disminuye en 0,006.
La variable indicadora que se muestra como progAcademic es la diferencia esperada en el recuento de registros entre el grupo 2 y el grupo de referencia (prog = 1). El recuento de registros esperado para el nivel 2 de prog es 0.44 menor que el recuento de registros esperado para el nivel 1. La variable indicadora de progVocational es la diferencia esperada en el recuento de registros entre el grupo 3 y el grupo de referencia. es 1,28 menor que el recuento de registros esperado para el nivel 1.
Modelos truncados
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.
## 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/Ploteamos:
ggplot(dat, aes(x = langscore, y = achiv)) +
geom_point() +
stat_smooth(method = "loess") +
facet_grid(. ~ prog, margins=TRUE)## `geom_smooth()` using formula 'y ~ x'
Regresion truncada.
##
## 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 DfLa 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.
modelos censurados
Tobit
## 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:AER':
##
## tobit## The following object is masked from 'package:lmtest':
##
## lrtest## The following object is masked from 'package:car':
##
## logit## 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':
##
## filldat <- read.csv("https://stats.idre.ucla.edu/stat/data/tobit.csv")
ggpairs(dat[, c("read", "math", "apt")])
Aplicamos el modelo:
##
## 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 estimatesLos 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.
Heckit
## Warning: package 'sampleSelection' was built under R version 4.0.5##
## Attaching package: 'sampleSelection'## The following object is masked from 'package:VGAM':
##
## probitEstimacion 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
## --------------------------------------------Estimacion por MV:
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
## --------------------------------------------