Suponga que usted ha decidido investigar la participación en la fuerza laboral femenina de las mujeres casadas mediante el siguiente modelo: \[\begin{equation} infl_i=\beta_0+\beta_1nwifeinc+\beta_2educ+\beta_3exper+\beta_4exper^2+\beta_5age+\beta_6kidslt6+\beta_7kidsge6+u_i \end{equation}\] Donde \(infl\) es una variable dicotómica que toma el valor de 1 si la mujer participó en la fuerza laboral en 1975, fuera del hogar por un salario, y 0 en otro caso; \(nwifeinc\) es el ingreso del esposo; \(educ\) corresponde a los años de educación de la mujer; \(exper\) y \(exper^2\) corresponden a los años de experiencia y experiencia al cuadrado respectivamente, \(age\) es la edad de la mujer; \(kidslt6\) es la cantidad de hijos menores de 6 años y \(kidsge6\) es la cantidad de hijos entre 6 y 18 años.
Para poder obtener los estimadores probit, logit MPL, y testearlas, utilizaremos las siguientes librerias: ## Paquetes a utilizar.
library(wooldridge)#Base de datos de Wooldridge
library(texreg)#Tabla de varias regresiones
library(AER)#Otras herramientas econométricas
library(aod)#Test de Wald
data("mroz")
attach(mroz)
A continuación estimamos el modelo presentado, mediante los 3 tipo de estimación vistos en un comienzo:
M_probit=glm(inlf~nwifeinc+educ+exper+expersq+age+kidslt6+kidsge6,
family = binomial(link = "probit"))
Mlogit=glm(inlf~nwifeinc+educ+exper+expersq+age+kidslt6+kidsge6,
family = binomial(link = "logit"))
MPL=lm(inlf~nwifeinc+educ+exper+expersq+age+kidslt6+kidsge6)
screenreg(l=list(MPL,Mlogit,M_probit),ci.force = F,digits = 4)
##
## ==========================================================
## Model 1 Model 2 Model 3
## ----------------------------------------------------------
## (Intercept) 0.5855 *** 0.4255 0.2701
## (0.1542) (0.8604) (0.5081)
## nwifeinc -0.0034 * -0.0213 * -0.0120 *
## (0.0014) (0.0084) (0.0049)
## educ 0.0380 *** 0.2212 *** 0.1309 ***
## (0.0074) (0.0434) (0.0254)
## exper 0.0395 *** 0.2059 *** 0.1233 ***
## (0.0057) (0.0321) (0.0188)
## expersq -0.0006 ** -0.0032 ** -0.0019 **
## (0.0002) (0.0010) (0.0006)
## age -0.0161 *** -0.0880 *** -0.0529 ***
## (0.0025) (0.0146) (0.0085)
## kidslt6 -0.2618 *** -1.4434 *** -0.8683 ***
## (0.0335) (0.2036) (0.1184)
## kidsge6 0.0130 0.0601 0.0360
## (0.0132) (0.0748) (0.0440)
## ----------------------------------------------------------
## R^2 0.2642
## Adj. R^2 0.2573
## Num. obs. 753 753 753
## AIC 819.5303 818.6044
## BIC 856.5228 855.5969
## Log Likelihood -401.7652 -401.3022
## Deviance 803.5303 802.6044
## ==========================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05
Observamos que si bien los estimadores pueden diferir entre el modelo probit y logit, los retornos marginales serán bastante similares, además, también podemos observar que sus estadísticos son casi idénticos.
A diferencia del MPL, los modelos Probit y Logit no tienen una bondad de ajuste que explique con exactitud la variabilidad de nuestra variable dependiente. Una aproximación bastante utilizada es el Pseudo \(R^2\) o \(R^2\) de McFadden, cuya fórmula es: \[ R^2=1-\frac{ln \hat{L}(M_{libre})}{ln(\hat{L}(M_R))} \] El cual suele ser más bajo que el \(R^2\) que vemos comunmente en MCO. Este estadístico es sencillo de obtener, ya que la regresión nos da la opción de obtener tanto \(ln \hat{L}(M_{libre})\) como \(ln(\hat{L}(M_R))\) mediante las devianzas, cuya fórmula es: \(devianza=-2*log-likehood\). Utilizando este estadístico, podemos armar el \(Pseudo-R^2\)
R2_probit=1-M_probit$deviance/M_probit$null.deviance
R2_logit=1-Mlogit$deviance/Mlogit$null.deviance
paste(round(c(R2_logit,R2_probit),3))
## [1] "0.22" "0.221"
Otra forma de obtener el mismo estadístico, sería mediante la creación del modelo restringido:
probit_restringido=glm(inlf~1,
family = binomial(link = "probit"))
logit_restringido=glm(inlf~1,
family = binomial(link = "logit"))
R2_probit2=1-logLik(M_probit)/logLik(probit_restringido)
R2_logit2=1-logLik(Mlogit)/logLik(logit_restringido)
paste(round(c(R2_probit2,R2_logit2),3))
## [1] "0.221" "0.22"
Como podemos observar, el resultado es el mismo.
Tal como sabemos, test-F ya no es válido para modelos no lineales, y por lo tanto, debemos recurrir a otro tipo de test, entre estos, construiremos la razón de verosimilitud y obtendremos su p-value para la hipótesis de significancia global, luego, mediante un comando estandarizado, obtendremos el test de Wald, cuyo resultado coincidirá con el primer test.
#----Probit----
#Construimos la razón mediante las devianzas.(También, como vimos anteriormente, es similar a )
RV_P=M_probit$null.deviance-M_probit$deviance
#Si lo quisieraamos realizar directamente mediante las funciones y utilizando la fórmula del test, tenemos que:
2*(logLik(M_probit)-logLik(probit_restringido))
## 'log Lik.' 227.142 (df=8)
#Obtenemos los grados de libertad
gdl_probit=M_probit$df.null-M_probit$df.residual
#Y finalmente el p-value
pvalue_probit=pchisq(q=RV_P,df=gdl_probit,lower.tail = F)
#----Logit----
#Ahora realizamos el mismo procedimiento con el modelo logit.
RV_L=Mlogit$null.deviance-Mlogit$deviance
gdl_logit=Mlogit$df.null-Mlogit$df.residual
pvalue_logit=pchisq(q=RV_L, df=gdl_logit, lower.tail=F)
pvalue_logit
## [1] 3.159176e-45
pvalue_probit
## [1] 2.008673e-45
Independiente si realizamos el procedimiento mediante las devianzas, o mediante la fórmula del test de RV, llegaremos al mismo resultado.
El argumento b nos indica los coeficientes, Sigma corresponde a la matriz de varianza y covarianza de los residuos, y finalmente Terms es donde indicamos los coeficientes que quisieramos testear. En este caso, al ser significancia global, nos interesa testear desde las 7 pendientes del modelo.
#probit
wald.test(b=coef(M_probit),Sigma=vcov(M_probit),Terms=1:7)
## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 184.4, df = 7, P(> X2) = 0.0
#logit
wald.test(b=coef(Mlogit),Sigma=vcov(Mlogit),Terms=1:7)
## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 156.6, df = 7, P(> X2) = 0.0