Rincón de Práctica R#
Probit/Logit usando R#
Ejemplo 1 (Nivel usuario).#
Nota: El enfoque del siguiente ejemplo estará solo en ilustrar el computo de efectos marginales manualmente (no en la parte de estimación vía MLE).
Primero simular datos y luego
comparar resultados entre los modelos LPM, Logit y Probit. Notar que, como vimos en la teoría, no es posible realizar una comparación directa entre los coeficientes estimados de los tres modelos.
# 1. Simular datos
set.seed(123)
nobs <- 100
x <- rnorm(nobs)
y <- rbinom(n=nobs, size=1, prob=pnorm(x))
# 1. Comparar LMP / Logit / Probit
lpm <- lm(y~x)
logit <- glm(y~x, family=binomial(link="logit"))
probit<- glm(y~x, family=binomial(link="probit"))
library(stargazer)
stargazer(lpm, logit, probit, type="text")
Please cite as:
Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
============================================================
Dependent variable:
----------------------------------------
y
OLS logistic probit
(1) (2) (3)
------------------------------------------------------------
x 0.303*** 2.083*** 1.186***
(0.044) (0.447) (0.234)
Constant 0.603*** 0.683** 0.376**
(0.040) (0.269) (0.152)
------------------------------------------------------------
Observations 100 100 100
R2 0.326
Adjusted R2 0.319
Log Likelihood -45.679 -45.875
Akaike Inf. Crit. 95.358 95.750
Residual Std. Error 0.400 (df = 98)
F Statistic 47.369*** (df = 1; 98)
============================================================
Note: *p<0.1; **p<0.05; ***p<0.01
Empleando la media de \(x\) y los resultados del logit, ilustrar el computo manual del efecto marginal:
\[f_i*\hat{\beta}_1=f(\hat{\beta}_0+\hat{\beta}_1 x_i)*\hat{\beta}_1=\frac{exp(-\hat{\beta}_0-\hat{\beta}_1 x_i)}{(1+exp(-\hat{\beta}_0-\hat{\beta}_1 x_i))^2}\hat{\beta}_1\]
# Efectos marginales
b0 <- logit$coefficients[1]
b1 <- logit$coefficients[2]
xbar <- mean(x)
(exp(-(b0+(b1*xbar)))/(1 +exp(-(b0+(b1*xbar))))^2)*b1 # Efecto marginal
# Comparación usando mfx
library(mfx)
logitmfx(y~x, data=data.frame(cbind(y,x)))
(Intercept): 0.433155566019768
Call:
logitmfx(formula = y ~ x, data = data.frame(cbind(y, x)))
Marginal Effects:
dF/dx Std. Err. z P>|z|
x 0.433156 0.085583 5.0612 4.166e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Es decir, como se esperaba, ambos resultados concuerdan.