1. Repaso básico del modelo lineal


Una regresión lineal con una variable explicativa es un método para analizar cómo una variable X (independiente) genera un cambio en una vartiable Y (variable dependiente). El modelo general de regresiones se explica a continuación:



La notación puede variar en algunos libros: por ejemplo la constante a veces es escrita como \(\alpha\), u otras veces como \(\beta_1\); En el caso del error este también es escrito como \(\varepsilon_i\). Será deber del estudiante estar atento para identificar qué elemento es cada uno, pero la constante es fácil de identificar pues no va acompañada de una variable explicativa \(X_i\).


…pero falta aclarar algo importante en la imagen explicada


\(\beta_0\) y \(\beta_1\) son desconocidos. Son medidas poblacionales que debemos estimar a partir de la muestra que tenemos disponibles. Por lo tanto la notación de la pendiente estimada por la muuestra es nuestro amigo beta gorro: \(\hat{\beta}\).


2. Estimación de mínimos cuadrados ordinarios (MCO)


La característica principal los MCO es que elige los coeficientes de regresión de tal forma que la recta de regresión estimada se encuentre lo más cercana posible a los datos observados, y la cercanía está medida por la suma de los errores al cuadrado que se cometen con la predicción de Y dado X.

Sean \(b_0\) y \(b_1\) algunos de los estimadores de \(\beta_0\) y \(\beta_1\). La recta de regresión basada en esos estimadores es \(b_0 + b_1X\) por lo que el valor de \(Y_i\) previsto mediante esta esta es \(b_0 + b_1X_i\). Por lo tanto el error se puede calcular como \(u_i = Y_i -(b_0 + B-1X_i)\).

La suma de todos los errores al cuadrado para n predicciones es:

\[\sum_{i=1}^{n}(Y_i - b_0 - b_1X_i)^2\]

MCO tiene su propia notación. El estimador MCO de \(\beta_0\) es \(\hat\beta_0\), y el estrimador \(\beta_1\) es \(\hat\beta_1\). lOS RESIDUOS \(\hat u_i = Y_i - \hat Y_i\). En este sentido los gorros muestran que son homogos muestrales de \(\beta_0\) y \(\beta_1\).

Para obtener \(\hat \beta_1\) se computa:

\[\hat \beta_1 = \frac{\sum_{i=1}^n(x_i - \bar x)(y_i - \bar y)}{\sum_{i=1}^n(x_i - \bar x)^2} = \frac{S_{xy}}{S^2_x}\] Y para la constante:

\[\hat \beta_0 = \bar y - \hat \beta_1X_i\]


2.1 Medidas de ajuste (\(R^2\) y Error estándar)


Las medidas de ajuste las podemos entender como mecanismos para observar en qué medida esta regresión lineal descibre correctamente. El \(R^2\) y el error estándar miden la bondad de ajuste de la recta de regresión MCO a los datos.

2.1.1 \(R^2\)

Se popuede sefinir como la proporción de varianza muestral de \(Y_i\) explicada por (o predicha por) \(X_i\). Matemáticamente el \(R^2\) puede escribirse como el cociente entre la suma explicada (de cuadrados) y la suma total de cuadrados. La suma explicada es la suma de las desviaciones al cuasdrado de los valores prerdichos de \(Y_i\) (\(\bar Y_i\)), respecto de su media y la suma total es la suma de los cuadrados de las desviaciones de \(Y_i\) respecto de su media:

\[ SE = \sum_{i=1}^n(\hat Y_i - \bar Y)^2\]

\[ST = \sum_{i=1}^n(Y_i - \bar Y)^2\]

Por lo tanto:

\[R^2 = \frac{SE}{ST}\]


De manera alternativa, \(R^2\) puede ser expresado como el cociente entre la varianza \(Y_i\) no explicada por \(X_i\). Para eso podemos tomar la suma de cuadrados de residuos, o suma resicual, que es:

\[SR = \sum_{i=1}^n \hat u_i^2\]

Por lo tanto el \(R^2\) puede se reescrito como:

\[R² = 1 - \frac{SR}{ST}\]

El \(R^2\) toma valores entre 0 y 1, donde 0 es que el modelo no explica nada de la variación \(Y_i\) sobre \(X_i\) y 1 es que explica completamente dicha variación. Habitualmente los valores nunca son muy cercanos a 0 o 1, sino que están en un intermedio. Si un \(R^2\) de 0.8 implica que el modelo tiene un muy buen rendimiento de la variable explicativa sobre \(Y_i\).


2.1.2 Error estandar


El error estandar es un estimador de la desviación típica del error de regresión \(u_i\). Como \(u_i\) e \(Y_i\) tienen las mismas medidas, el ESR es una medida de dispersión de las observacionesen torno a una recta de regresión. Ya que los \(U_i\) no están disponibles por ser medidas poblacionales, se usan sus medidas muestrales \(\hat u_i\):

\[ESR = S_{\hat u}, donde s_{\hat u}^2 = \frac{1}{n - 2} \sum_{i=1}^n = \frac{SR}{n - 2}\] La razón de usar \(n - 2\) es que corrige una desvación a la baja introducido al estimar dos coeficientes de regresión. Esto se denomina Corección por grados de libertad.


Ejemplo: Ingreso y edad


Para ejemplificar todo lo que hemos expuesto hasta ahora, analizaremos un pequeño caso: La relación entre ingreso percapita (variable dependiente) en relación a los años de escolaridad (variable independiente). Para esto usaremos los datos de la encuesta CASEN 2017, y para que la base de datos sea manejable sólo haremos el ejemplo de la comuna de La Cisterna, Santiago de Chile.

Computamos la regresión y obtenemos los siguientes resultados:


ejemplo <-lm(ytotcor ~ Exp, data = casen)

stargazer::stargazer(ejemplo, type = "html")
Dependent variable:
ytotcor
Exp -8,494.868**
(3,666.977)
Constant 864,623.600***
(89,282.810)
Observations 119
R2 0.044
Adjusted R2 0.036
Residual Std. Error 511,859.100 (df = 117)
F Statistic 5.367** (df = 1; 117)
Note: p<0.1; p<0.05; p<0.01


Lo que nos está diciendo esta regresión es que el cambio marginal de la experiencia en una unidad hace que el ingreso total baje 8.494 pesos chilenos. Esto es … raro. La experiencia nos debería llevar a tener mejor sueldo, pero esto demuestra lo contrario. Si observamos la significancia, esta es significativa en p < 0.05, y el error estandar (que es el bajos que sale entre paréntesis debajo del coeficiente Exp es de 3.666.977), que es grande pero tampoco es gigante y major al mismo coeficiente.

De manera gráfica la regresión muestra lo siguiente:


library(ggplot2)
ggplot(data = casen, aes(x = Exp, y = ytotcor)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE, color = "firebrick") +
  theme_bw() + labs(x = "", y = "")


Cada punto del gráfico es el valor y en su respectiva x (puntos de experiencia). La línea roja marca la pendiente que el estimador MCO calcula con los datos y el sombreado muestra el error estandar de la regresión. Los puntos muy disparados hacia arriba se pueden considerar como outliers, Y el estimador de Mínimos cuadrados ordinarios es sensible a ellos.

¿Qué nos muestra esta regresión? que en las dos puntas la regresión aumenta sus errores estandar, implicando la presencia de valores que hacen variar mucho el coeficiente obtenido, y que en particular hay una gran cantidad de outliers en los valores x de menos experiencia (sobre todo entre 15 a 25 años de experiencia) que tienen ingresos muy altos. Es posible pensar que esos valores son los que hacen que la pendiente se vuelva negativa, pero aún queda mucho por explorar. ya seguiremos analizando este ejemplo.

2.2 ¿Por qué usar Mínimos cuadrados ordinarios?


Podemos resumir en las siguientes razones:

  1. Son el lenguaje común entre economistas, y es el método dominante en todas las ciencias sociales.

  2. Bajo los supuestos que veremos ahora, MCO es insesgado y consistente.

  3. En el marco de otros estimadores insesgados, MCO es uno de los más eficiente. Siempre y cuando se mantengan algunas condiciones.


2.3 Supuestos del modelo lineal


Cuando se hace una estimación por Mínimos cuadrados ordinarios, no sólo se buscan los valores de \(\hat \beta_0\) y \(\hat \beta_1\), sino que hacer inferencia a ver si \(\hat Y_i\) se acerca al valor real de Y, es decir \(E(Y | X_i)\). Esta inferencia se realiza bao ciertos supuestos de \(X_i\) y \(u_i\). Los supuestos son críticos para una interpretación válida de los resultados.

Según el libro de Gujaragit (2007) el modelo lineal tiene 10 supuestos:

  1. El modelo de regresión es lineal en los parámetros

\[Y_i = \beta_o + \beta_1X_i + u_i\]

  1. Los valores de X son fijos en un muestreo repetido

Esto quiere decir que los valores de X no son estocásticos. Recordar que es un análisis de regresión condicional a cierto valor x.

  1. El valor medio del error \(u_i = 0\)

Cada Y poblacional correspondiente a determinado X está distribuído alrededor de la media. Los valores por encima y debajo son los errores, y el promedio del vaor medio es cero (ya que los valores positivos y negativos se deben cancelar hasta llegar a cero)

  1. Las varianzas de \(u_i\) son iguales (Homocedasticidad)

Homocedástico se refiere a “Igual disperción”. La varianza de \(u_i\) para cada \(X_i\) es algún número positivo constante igual a \(\sigma^2\). En otras palabras:

\[var(u_i | X) = E[u_i - E(u_i) | X]^2\] \[ = E(u_i^2 | X )\] \[ = \sigma^2 \]

  1. No Existe correlación entre errores

Este es un supuesto de no correlación lineal. Si existen patrones en la distribución de errores, están autocorrelacionados. Si eso pasa, \(Y\) ya no sólo depende de X, sino también de \(u_i\).

  1. La varianza entre \(u_i\) y \(X_i\) es cero (\(E[U_iX_i] = 0\))

Para que la estimación por MCO esté correcta, los errores \(u_i\) y \(x_i\) no están correlacionados. Si lo está, no es posible determinar el efecto diferenciador entre \(u_i\) y \(x_i\).

De todas formas, buenas noticias: Este supuesto se cumple automáticamente si la variable X no es estocástica y se mantiene el supuesto 3.

  1. El número de observaciones n debe ser mayor al número de parámetros por estimar

Se necesitan al menos dos pares de observaciones por parámetro.

  1. Variables en valores de X

No todos los valores de X en una muesta dada deben ser iguales. Si todos on iguales implica que \(X = \bar X\), y la estimación de \(\beta_0\) y \(\beta_1\) sería imposible.

  1. El modelo de regresión está correctamente especificado

Con este supuesto no habría sesgo de especificación por error. Y la teoría econométrica asume que el modelo está especificado correctamente.

Sin embargo se puede testear a partir de algunas preguntas: ¿Cuáles variables están incluídas en el modelo? ¿Cuál es la forma funcional del modelo? ¿Es lineal envariables, en parámetros o ambos? ¿Cuáles son los supuestos probabilísticos considerados para \(Y_i\), \(X_i\) y \(u_i\) en el modelo?.

El sesgo consiste en escoger la forma funcional equivocada. En este sentido, la construcción de modelos econométricos es más un arte que una ciencia.

  1. No hay multicolinealidad perfecta

Quiere decir que no hay relaciones perfectamente lineales entre variables explicativas. En algunos casos también se busca la colinealidad imperfecta.

Si tienen relación pero no es lineal no importa para efectos de MCO.

Por el contrario, si hay colinealidad perfecta, los coeficientes son indeterminados y los errores son infinitos.

Si hay colinealidad pero no es perfecta, los coeficientes se pueden determinar pero los errores estandar serán altos.

¿Por qué puede haber multicolinealidad?

Puede ser por método de recolección de muestra (muestras con intervalos limitados), restricciones del modelo o en la población de muestra, errores de especificación del modelo o porque está sobredeterminado.


3. Testeo de supuestos de modelo lineal


3.1 Prueba de heterocedasticidad (White Test)

El test de White es una forma intuitiva y rápida de analizar si tenemos heterocedasticidad en nuestra regresión. Con el paquete lmtest y la función bptest lo analizaremos rápidamente:


bptest(ejemplo, ~ Exp + I(Exp^2), data = casen)
## 
##  studentized Breusch-Pagan test
## 
## data:  ejemplo
## BP = 3.2018, df = 2, p-value = 0.2017


Vemos que el test estadístico \(X^2\) da 3.2018, los grados de libertad son 2 y el p-valor es de 0.2017. La forma de interpretar esto es la siguiente:

El test de White usa la siguiente forma: La hipótesis nula (\(H_0\)) es que la homocedasticidad está presente (lo que buscamos!), y la alternativa (\(H_A\)), es que la heterocedasticidad está presente. Como el p-valor no es menos a 0.05, fallamos al rechazar la hipótesis nula. Es decir: al tener valores mayores de 0.05 implica que no tenemos suficiente evidencia para decir que la heterocedasticidad está presente en el modelo de regresión, por lo tanto está correcta.

Si el modelo hubiese tenido heterocedasticidad ¿Qué podríamos hacer? podemos hacer la regresión lineal robusta con la función rlm.


3.2 Test de Ramsey para errores de especificación de Modelos


Comos señala la página oficial del paquete lmtest de R, el test Reset de Ramsey es un diagnóstico para las corrección de las formas funcionales de modelos de estimación. Su principal característica es observar si hay variables no estipuladas que afectan a mi variable dependiente, de modo de poder analizar la pertinencia del modelo.

El test de ramsey tiene como hipótesis nula que no hay errores de especificación en su modelo. La hipótesis alternativa es que sí tiene errores de específicación.


resettest(formula = ejemplo, power = 2:3, type = "regressor")
## 
##  RESET test
## 
## data:  ejemplo
## RESET = 0.52495, df1 = 2, df2 = 115, p-value = 0.593


Al igual que el tetst de White, la Hipótesis nula es que no hay errores de especificación del modelo. Al realizarlo vemos que tenemos un p-valor de 0.513, lo que indica que no tenemos evidencia para creer que hay errores de especificación, por lo que el modelo sí está bien especificado.


3.3 Plot de diagnóstico para la regresión


par(mfrow = c(2,2))

plot(ejemplo)


¿Cómo interpretar esto?

La eventual presencia de patrones de movimiento de la línea indicarían algún problema del modelo lineal.

En nuestro caso se ve relativamente bien pero tiene un paulatino levantamiento en la zona derecha que muestra que no es perfecta. Eso puede ser una señal de heterocedascitidad a pesar del test anterior.

Habitualmente una forma de arreglar esto es adecuar la escala de la variable dependiente, como por ejemplo transformarla a logaritmo raíz cuadrada.

Sólo por curiosidad, probamos esta técnica en nuestr modelo:


ejemplo2 <- lm(log(ytotcor) ~ Exp, data = casen)

plot(ejemplo2, 3)


En nuestro caso pasarlo a logaritmo ha incrementado la diferencia en la medida que aumentaba la experiencia, dando mayores señales de heterocedasticidad.

Los outliers pueden ser identificados examinando los residuos estandarizados, los cuales son los residuos dividido por su error estandar estimado. Los residuos estandarizados pueden ser interpretados como el número de errores estandar fuera de la linea de regresión.

Por otra parte, los valores influyentes son valores que, en su inclusión o exclusión, pueden alterar los resultados del análisis de regresión. Ese valor está asociado a un residuo grande. En este sentido, es importante señalar que no todos los outliers son influyentes.

Para medir esto se ha desarrollado la distancia de Cook para determinar la influencia de un valor. Esta métrica define influencia como la combinación de apalancamiento (leverage) y tamaño de residuo.

Una forma de ver la influencia con la distancia de cook es que exceda \(\frac{4}{n - p - 1}\), donde \(n\) es el número de observaciones y \(p\) es el número de variables predictoras.

plot(ejemplo, 4)

plot(ejemplo, 5)


En el gráfico residuals vs leverage, los valores outliers que son influyentes se ubican en la esquina superior derecha y la esquina inferior derecha. En nuestro caso, casi no tenemos valores de ese tipo, salvo tres claramente marcados y uno en particular.


3.4 Regresión lineal multivariada y testeo de multicolinealidad

En este caso agregaremos otras variables que pueden estar explicando la variabilidad del ingreso sumando la edad y el sexo de la persona:


ingreso_multi <-lm(ytotcor ~ Exp + edad + SexDum , data = casen)

stargazer::stargazer(ingreso_multi, type = "html")
Dependent variable:
ytotcor
Exp -68,349.590***
(18,571.080)
edad 65,486.560***
(19,946.570)
SexDum -38,404.650
(90,452.380)
Constant -528,968.400
(440,910.200)
Observations 119
R2 0.127
Adjusted R2 0.104
Residual Std. Error 493,338.800 (df = 115)
F Statistic 5.576*** (df = 3; 115)
Note: p<0.1; p<0.05; p<0.01


Para terminar con los análisis de supuestos, analizaremos la multicolinealidad perfecta. Su aparición es compleja pues implicaría que no hay ortogonalidad en los regresores, lo que haría engañosa la interpretación de esta regresión.

Para testearla haremos una matriz de correlaciones entre las variables explicativas:

casen_matrix <-casen%>%
  dplyr::select(Exp, edad, SexDum)

casen_matrix%>%
  cor(method = "pearson")%>%
  round(digits = 2) ->matcor

matcor
##          Exp  edad SexDum
## Exp     1.00  0.98  -0.01
## edad    0.98  1.00  -0.01
## SexDum -0.01 -0.01   1.00


Como esta regresión es pequeña y sólo tiene 3 regresores es fácil de observar los problemas de colinealidad casi perfecta qe hay entre edad y experiencia (se entiende, pues experiencia fue creada a partir de esta variable). Pero cuando tenemos 5, 6 o 10 regresores vale la pena tener herramientas gráficas para observarla. las mostramos a continuación:


corrplot(matcor, type = "upper", order = "hclust", tl.col = "black")


De esta forma tenemos herramientas detalladas para analizar la multicolinealidad.

4. Regresiones con una variable dependiente Binaria (Logit y Probit)


Cuando se trata de una variable dependiente binaria, hay que interpretar la función de regresión como una predicción de probabilidad.

Si Y es una variable binaria 0 o 1, su valor esperado (o media) es la probabilidad de \(Y = 1\), osea:

\[E(Y) = 0 * Pr(Y = 0) + 1 * Pr(Y = 1) = Pr (Y = 1)\]

n el contexto de regresión, el valor esperado está condicionado al valores de los regresores, por lo tanto la probabilidad está condicionada a \(X\).

Por ejemplo, el modelo de probabilidad lineal (que no vimos en clases y no profundizaré más que por temas pedagógicos) está dado por:

\[Pr(Y = 1 | X_1, X_2,..., X_kn)= \beta_0 + \beta_1X_1 + \beta_2X_2 + ... \beta_kX_k\]

Donde el coeficiente \(\beta_1\) es la probabilidad de \(Y = 1\) cuando \(X_1\), ceteris paribus. Esos coeficientes se estiman por MCO ya que es una probabilidad lineal.

¿Cuál es el problema de esto? que raras veces las probabilidades son lineales, además que hace que valores muy bajos de \(X\) den números negativos y se supone que son entre 0 y 1. No tiene sentido.

Es por eso que se usan otros modelos de regresión, como el Logit y Probit.

Para esta sección usaremos la encuesta CEP de Diciembre de 2022, para analizar la probabilidad de aprobar el gobierno de Gabriel Boric.


4.1 Regresión Probit


Este es un modelo de regresión no-lineal para variables dependientes binarias, que se expresa en la siguiente fórmula:

\[Pr(Y = 1 | X) = \Phi (\beta_0 + \beta_1X)\] Done \(\Phi\) es una función de distribución de probabilidades acumulada normal estandar.

EN el caso de tener múltiples regresores siemplemente se agregan a la función de distribución de probabilidad:

\[Pr(Y = 1 | X_1, X_2) = \Phi(\beta_0, \beta_1X_1 + \beta_2X_2)\]

Puntualmente ¿Qué provca la varación de \(X\) en \(Y\)? Lo que provoca es un cambio en la esperanza condicional; es decir, en la probabilidad condicional de que \(Y = 1\).

Para ejemplificar este modelo, vamos a analizar la probabilidad de obtener el Ingreso ético familiar a pasar de la proporción de hijos dentro del número familiar. Esto lo haremos en base a la encuesta CASEN 2017:

el código resumido (con las variables ya construídas) es el siguiente:


CasenFiltra <-readRDS(file = "casenfiltrada1.RDS")

Asigna_probit <-glm(AsignaFam ~ partninos, family = binomial(link = "probit"), data = CasenFiltra)

stargazer::stargazer(Asigna_probit, type = "html")
Dependent variable:
AsignaFam
partninos 0.595***
(0.116)
Constant -2.038***
(0.024)
Observations 17,893
Log Likelihood -1,986.193
Akaike Inf. Crit. 3,976.386
Note: p<0.1; p<0.05; p<0.01


Vemos que la regresión muestra que la variable elegida es significativa al 99%, y que el cambio de proporción de niños en el hogar en una unidad aumenta un 59% de unidades probit para obtener el Ingreso ético familiar. Suena raro eso de “Unidades probit” no? ya lo resolveremos cuando veamos efectos marginales.


4.2 Regresión Logit


Es una función similar a la probit, sólo que la función de distribución de probabilidad es logística estándar:

\[Pr(Y = 1 | X_1, X_2,..., X_k) = F(\beta_0 + \beta_1 + \beta_2X_2 + ... + \beta_kX_k) = \frac{1}{1 + e^{(\beta_0 + \beta_1X_1 + \beta_2X_2 + .. + \beta_kX_k)}}\]


Donde \(F\) Es la función de distribución de probabilidad logística. También se puede estimar por método de máxima verosimilitud.

Al igual que Probit, el Logit es más fácil de interpretar calculando las probabilidades estimadas y diferencias entre probabilidades estimadas.

Utilizando el mismo ejemplo del modelo Probit, para ver el cambio de probabilidad de obtener el Ingreso Ético familiar en base a la proporción de niños en el hogar.

AsignaLogit<- glm(AsignaFam ~ partninos, data = CasenFiltra, family = binomial(link = "logit"))

stargazer::stargazer(AsignaLogit, type = "html")
Dependent variable:
AsignaFam
partninos 1.364***
(0.266)
Constant -3.850***
(0.058)
Observations 17,893
Log Likelihood -1,986.521
Akaike Inf. Crit. 3,977.042
Note: p<0.1; p<0.05; p<0.01


En esta regresion vemos que el cambio es mucho más acentuado que el probit. por cada unidad de cambio en la proporción de niños en el núcleo familiar aumenta 1.4 veces unidades logit la probabilidad de obtener el ingreso ético familiar, al 99% de confianza.


4.3 Regresión Probit y Logit multivariada y efectos marginales

Al igual que las regresiones lineales multivariadas, podemos agregas tantas variables explicativas deseamos considerando la cantidad de observaciones que tenemos.

Para tener resultados que sean interpretables es importante transformar estas regresiones en predicciones para valores específicos (como lo hace la función predict en R), pero también del cambio marginal promedio a partir del cálculo de efectos marginales.

dejamos el código de la regresión multivariada Probit y Logit donde buscamos observar si el Ingreso Percápita total, la escolaridad del jefe de familia, la condición de ruralidad, cantidad de adultos mayores en la familia, si el jefe del hogar es mujer, la cantidad de integrantes en el núcleo y la presencia de discapacitados afecta a la probabilidad de obtener el Ingreso ético familiar:


Asigna_probit_multi <-glm(AsignaFam ~ PercapTotal + escojefe + partninos + rural + partmujer + partadultmayor + jefemujermonoparental + jefemujer + intergrantesnucl + Disca, family = binomial(link = "probit"), data = CasenFiltra)


Asigna_logit_multi <- glm(AsignaFam ~ PercapTotal + escojefe + partninos + rural + partmujer + partadultmayor + jefemujermonoparental + jefemujer + intergrantesnucl + Disca, family = binomial(link = "logit"), data = CasenFiltra)

stargazer::stargazer(Asigna_probit_multi, type = "html")
Dependent variable:
AsignaFam
PercapTotal 0.00000**
(0.00000)
escojefe -0.003
(0.009)
partninos 0.500**
(0.217)
rural 0.047
(0.065)
partmujer -0.089
(0.133)
partadultmayor 0.040
(0.099)
jefemujermonoparental 0.155
(0.116)
jefemujer 0.166
(0.105)
intergrantesnucl 0.220***
(0.031)
Disca 0.069
(0.068)
Constant -2.797***
(0.176)
Observations 8,256
Log Likelihood -950.762
Akaike Inf. Crit. 1,923.524
Note: p<0.1; p<0.05; p<0.01
stargazer::stargazer(Asigna_logit_multi, type = "html")
Dependent variable:
AsignaFam
PercapTotal 0.00000**
(0.00000)
escojefe -0.007
(0.020)
partninos 1.172**
(0.479)
rural 0.088
(0.152)
partmujer -0.143
(0.311)
partadultmayor 0.018
(0.238)
jefemujermonoparental 0.277
(0.257)
jefemujer 0.422*
(0.235)
intergrantesnucl 0.475***
(0.067)
Disca 0.153
(0.163)
Constant -5.521***
(0.412)
Observations 8,256
Log Likelihood -951.220
Akaike Inf. Crit. 1,924.441
Note: p<0.1; p<0.05; p<0.01


En estos resultados vemos algunas diferencias que son relevantes. En términos de significancia estadística, el modelo Probit destaca el Percapita total, la proporción de niños en núcleo familiar y la cantida de integrantes por núcleo familiar. En la magnitud del cambio el Ingreso percapita total tiende a cero, por lo que no hay mayor cambio, mientras que la concentración de niños lllega hasta 50% unidades probit más, y el aumento en un integrante más por núcleo lleva a 22% unidades Probit más de probabilidad de tener la asignación. Por otra parte el modelo Logit agrega que si el jefe de hogar es mujer aumenta considerablemente la probabilidad de obtener este Ingreso especial, así como los valores de magnitud tanto en la proporción de niños y el número de integrantes del núcleo familiar es mayor.

¿Cómo se puede justificar esta diferencia? Es imposrtante recordar que los modelos Probit, al distribuir con funciones de probabilidad acumulada estándar, son algo màs sensibles a los outliers, y quizás esa sea la razón. Sin embargo siempre es importante testear ambos modelos para observar la consistencia de la estimación y del modelo.

Ahora pasaremos a cómo interpretar esta información. Es muy importante poder traspasar estas regresiones probabilidades para poder interpretar los resultados.

A continuación dejamos los efectos marginales promedio de la regresión Probit:


test1 <-margins(Asigna_probit_multi)

test2 <-as_tibble(summary(test1))
test2 %>% dplyr::select(factor, AME, lower, upper,p) 
factorAMElowerupperp
Disca0.004   -0.0038  0.0118  0.315  
escojefe-0.000185-0.00117 0.0007980.712  
intergrantesnucl0.0128  0.00906 0.0165  1.4e-11
jefemujer0.00966 -0.00231 0.0216  0.114  
jefemujermonoparental0.00897 -0.00427 0.0222  0.184  
partadultmayor0.00231 -0.00898 0.0136  0.689  
partmujer-0.00519 -0.0203  0.00993 0.501  
partninos0.029   0.00414 0.0539  0.0223 
PercapTotal7.48e-088.25e-091.41e-070.0276 
rural0.00271 -0.00472 0.0101  0.474  


Esto nos entrega una tabla regumen con los efectos marginales promedios (AME) por cada variable explicativa, así como el límite inferior (lower) y el máximo (upper) de cada una de ellas.

Luego de eso creamos un gráfico con los efectos marginales de cada uno:


p1 <- ggplot(data = test2, aes(x = reorder(factor, AME),
                              y = AME, ymin = lower, ymax = upper))

p1 + geom_hline(yintercept = 0, color = "blue") +
  geom_pointrange() + coord_flip() +
  labs(x = NULL, y = "Average Marginal Effect", title = "Efectos Marginales Promedio para Asignación de Ingreso Ético Familiar  \n en base a modelo Probit") 


Así, tenemos una forma gráfica de ver qué tan distantes del 0 y en qué sentido se encuentra el cálculo de probabilidad promedio en cada una. Confirmamos entonces que por cada cambio unitario en el porcentaje de niños en el grupo familiar aumenta casi un 3% la probabilidad de asignación familiar, seguido por la cantida dde integrante del número, y por otro lado la cantidad de mujeres tiene un efecto negativo en la probabilidad de obtener el beneficio.

Veamos con el modelo Logit:


test3 <-margins(Asigna_logit_multi)

test4 <-as_tibble(summary(test3))
test4 %>% dplyr::select(factor, AME, lower, upper,p) 
factorAMElowerupperp
Disca0.00383 -0.00414 0.0118  0.346   
escojefe-0.000165-0.00116 0.0008330.746   
intergrantesnucl0.0119  0.00836 0.0154  3.79e-11
jefemujer0.0105  -0.00105 0.0221  0.0747  
jefemujermonoparental0.00693 -0.00569 0.0196  0.282   
partadultmayor0.000456-0.0112  0.0121  0.939   
partmujer-0.00358 -0.0188  0.0116  0.644   
partninos0.0293  0.00559 0.053   0.0154  
PercapTotal7.3e-08 9.99e-091.36e-070.0232  
rural0.00221 -0.00526 0.00968 0.562   


p1 <- ggplot(data = test4, aes(x = reorder(factor, AME),
                              y = AME, ymin = lower, ymax = upper))

p1 + geom_hline(yintercept = 0, color = "blue") +
  geom_pointrange() + coord_flip() +
  labs(x = NULL, y = "Average Marginal Effect", title = "Efectos Marginales Promedio para Asignación de Ingreso Ético Familiar  \n en base a modelo Logit") 

Finalmente, vemos que las grandes diferencias que habíamos visto anteriormente entre los dos modelos no eran tal, y de hecho se comportan de forma muy fimilar una del otro.

4.4 Diferencias entre Logit, Probit y probabilidad lineal

**Las regresiones logit y probit son casi iguales. Se usa más la logit porque era más rápida de calcular a mano que la Probit. Luego del avance de capacidades computacionales su diferencia se volvió irrelevante.


Entonces: ¿Cuáles usar?


  1. Los tres son aproximaciones de funciones de regresión poblacional desconocidas. es decir \(E(Y | X) = Pr(Y = 1 | X)\).
  2. Modelos de probabilidad lineal es + fácil pero es malo, pues no capta la naturaleza no lineal de la verdadera función poblacional.
  3. Modelos Probit y Logit dan resultados muy similares.


4.4 Estimación por Máxima verosimilitud


La estimación por máxima verosimilitud es una de las varias formas de estimar modelos Logit y Probit; otro ejemplo es la estimación por Mínimos Cuadrados no lineales pero que no hablaremos acá. Lo fundamental es entender es que la estimación por máxima verosimilitud hace posible estimar una función Logit o Probit.

Primero se parte con una función de verosimilitud, que es una distribución de probabilidad conjunta de los datos considerada como una función de los coeficientes desconocifod. Por otro lado el estimador de máxima verosimilitud (EMV) de los coeficientes desconocidos está compuesto por los valore sde los coeficientes que maximizan la función de verosimilitud. En este sentido, el EMV son los valores de los parámetros que “más probablemente” hayan generado los datos.


4.5 Medidas de ajuste: Pseudo \(R^2\)


el Pseudo \(R^2\) mide el ajuste del modelo mediante la función de verosimilitud. Entonces lo que hace esta medida es medir la calidad del ajuste mediante la comparación del valor de la función de verosimilitud maximizada con todas las variables explicativas con el valor de la función de verosimilitud sin regresores.

Si quiero obtener el pseudo \(R^2\) de la regesión se realiza con el siguiente código:


# Carga de librerías y datos de ejemplo
library(pscl)

# Cálculo del Pseudo R-cuadrado y adición a la tabla
tabla <- summary(Asigna_probit_multi)$coef
tabla$Pseudo_R2 <- pR2(Asigna_probit_multi)
## fitting null model for pseudo-r2
# Mostrar la tabla con el Pseudo R-cuadrado
tabla
## [[1]]
## [1] -2.796668
## 
## [[2]]
## [1] 1.288443e-06
## 
## [[3]]
## [1] -0.003186777
## 
## [[4]]
## [1] 0.5000359
## 
## [[5]]
## [1] 0.0467108
## 
## [[6]]
## [1] -0.08945153
## 
## [[7]]
## [1] 0.03977568
## 
## [[8]]
## [1] 0.1545014
## 
## [[9]]
## [1] 0.1663019
## 
## [[10]]
## [1] 0.2196523
## 
## [[11]]
## [1] 0.06886231
## 
## [[12]]
## [1] 0.1763539
## 
## [[13]]
## [1] 5.687503e-07
## 
## [[14]]
## [1] 0.008632594
## 
## [[15]]
## [1] 0.2173457
## 
## [[16]]
## [1] 0.0652261
## 
## [[17]]
## [1] 0.1327903
## 
## [[18]]
## [1] 0.09922296
## 
## [[19]]
## [1] 0.1160955
## 
## [[20]]
## [1] 0.1048455
## 
## [[21]]
## [1] 0.03073727
## 
## [[22]]
## [1] 0.06846321
## 
## [[23]]
## [1] -15.85828
## 
## [[24]]
## [1] 2.265394
## 
## [[25]]
## [1] -0.3691563
## 
## [[26]]
## [1] 2.300648
## 
## [[27]]
## [1] 0.7161366
## 
## [[28]]
## [1] -0.6736299
## 
## [[29]]
## [1] 0.4008717
## 
## [[30]]
## [1] 1.330812
## 
## [[31]]
## [1] 1.586161
## 
## [[32]]
## [1] 7.146125
## 
## [[33]]
## [1] 1.005829
## 
## [[34]]
## [1] 1.232317e-56
## 
## [[35]]
## [1] 0.02348853
## 
## [[36]]
## [1] 0.7120112
## 
## [[37]]
## [1] 0.02141156
## 
## [[38]]
## [1] 0.473907
## 
## [[39]]
## [1] 0.5005466
## 
## [[40]]
## [1] 0.6885146
## 
## [[41]]
## [1] 0.1832507
## 
## [[42]]
## [1] 0.1127028
## 
## [[43]]
## [1] 8.926178e-13
## 
## [[44]]
## [1] 0.3144977
## 
## $Pseudo_R2
##           llh       llhNull            G2      McFadden          r2ML 
## -950.76207146 -996.50825216   91.49236140    0.04590647    0.01102075 
##          r2CU 
##    0.05138491


Acá podemos ver seis tipos de cálculos de pseudo \(R^2\) y nos indican que el modelo se ajusta terriblemente mal. No explica nada.


5. Regresión Tobit para modelos censurados


En algunas ocasiones ocurre que, por alguna razón de la naturaleza de los datos, los métodos de Mínimos Cuadrados Ordinarios puede no ser correcto. En este caso hablamos de Datos truncados y Datos censurados.

-Los Datos censurados Existen cuando un determinado límite en la variable respuesta (superior, inferior o ambos) a partir del cual todas las observaciones se les asigna el mismo valor. Ejemplo: un instrumento de medida como una balanza tiene un límita de deteccion de 0; eso se llama censura inferior. Por otra parte, cuando se pregunta por nivel de ingresos de las personas, la escala se divide en múltples intervalos y la última es “más de 5 millones” o lo que sea, hace que el tope llegue a esa categoría, sabiendo que existen datos sobre ella; eso se llama censura superior.

  • Los datos truncados Son datos que a partir de cierto límite no existen observaciones. La diferencia fundamental con los censurados es que, en estos últimos, las observaciones sí existen en la población latente pero no se puede captar en el muestreo.

¿Por qué importan tanto estas consideraciones? porque como se usa método de mínimos cuadrados, la pendiente de la recta varía drásticamente si no se tienen en consideración estas características.

Por ejemplo este modelo OLS clásico, y a la derecha el modelo censurado donde, por alguna razón, no se puede mejor bajo el y=30:



Cuando existe esta situación, los investigadores tienen que tomar una decisión: o todos los valores \(<30 = 30\) o \(<30 = 0\). Esto cambia profundamente la pendiente como vemos a continuación:



En el caso del punto de corte = 0 la pendiente aumenta su inclinación por el peso de los 0. En el segundo es el efecto contrario. Esto puede implicar problemas de heterocedascticidad y e independencia de observaciones.

Es por eso que existe el modelo Tobit

En la regresión de Tobit (1958) se considera que existe una variable latente \(Y∗\) no observable y una variable Y observable, formada por la parte no censurada de \(Y*\). El objetivo es ser capaz de estimar parámetros de \(Y∗\) empleando solo una muestra de la parte observable. Para ello, asumiendo censura por el límite inferior (misma idea para la superior), considera que el valor esperado de la variable censurada Y puede definirse como:

\[E(y) = [P(No censurado)* E(y | y > \tau)] + [P(censurado)* E(y | y = \tau_y)]\]

\[y_i=\left\{\begin{array}{ll}y_i^* & \text { if } y_i^*>\tau \\ \tau_y & \text { if } y_i^* \leq \tau\end{array}\right.\] Donde \(P\) es la proobabilidad, \(\tau\) límite de censura, \(\tau_y\) el valor que se le asigna a la variable latente \(Y*\) cuano se le aplica la censura \(y | y > \tau\) y \(y | y > \tau_y\) son el condicional de que Y sea mayor que el límite de censura y el valor asignado a la censura respectivamente.

En definitiva, el modelo Tobit permite hacer inferencias de \(Y*\) teniendo sólamente \(Y\).

Hay dos formas de hacer regresiones Tobit. Usaremos el ejemplo depuntajes en la universidad que van desde 200 a 800. Se entiende que 200 es un corte arbitrario, porque probablemente hayan personas que linealmente puedan sacar menos de 200 en base a sus conocimientos. Por otro lado 800 puntos es el máximo, pero muchas de estas personas que sacan ese puntaje probablemente saben más de ello. Estos son los límites superiores e inferiores.

datos <- read.csv("https://stats.idre.ucla.edu/stat/data/tobit.csv")

ggplot(data = datos, aes(x = apt)) +
  geom_histogram(binwidth = 10) +
  theme_bw() +
  labs(title = "Distribución de la aptitud académica\nde los estudiantes")


La exploración de datos muestra que nadie llega a los 200 puntos, por lo que la censura inferior es irrelevante, sin embargo la censura superior sí lo es, ya que mucha gente está en los 800 puntos, que es mucho mayor a lo esperado en una distribución normal.

modelo_tobit <- vglm(apt ~ read + math + prog, tobit(Upper = 800), data = datos)

summary(modelo_tobit)
## 
## Call:
## vglm(formula = apt ~ read + math + prog, family = tobit(Upper = 800), 
##     data = datos)
## 
## Coefficients: 
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1  209.55956   32.54590   6.439 1.20e-10 ***
## (Intercept):2    4.18476    0.05235  79.944  < 2e-16 ***
## read             2.69796    0.61928   4.357 1.32e-05 ***
## math             5.91460    0.70539   8.385  < 2e-16 ***
## 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 estimates


La forma de interpretar es la misma que un MCO, pero referidos a la variable latente no censurada.

También se puede hacer con la función censReg:

modelo_tobit2 <- censReg(apt ~ read + math + prog, right = 800, data = datos)
summary(modelo_tobit2)
## 
## Call:
## censReg(formula = apt ~ read + math + prog, right = 800, data = datos)
## 
## Observations:
##          Total  Left-censored     Uncensored Right-censored 
##            200              0            183             17 
## 
## Coefficients:
##                 Estimate Std. error t value  Pr(> t)    
## (Intercept)    209.56597   32.77196   6.395 1.61e-10 ***
## read             2.69794    0.61881   4.360 1.30e-05 ***
## math             5.91448    0.70982   8.332  < 2e-16 ***
## proggeneral    -12.71476   12.40646  -1.025 0.305434    
## progvocational -46.14390   13.72419  -3.362 0.000773 ***
## logSigma         4.18474    0.05301  78.946  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Newton-Raphson maximisation, 8 iterations
## Return code 1: gradient close to zero (gradtol)
## Log-likelihood: -1041.063 on 6 Df

Nótese que hay unas pequeñas diferencias pero probablemente se deban a sus maneras particulares de maximizar la función de verosimilitd.