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}\).
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\]
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.
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\).
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.
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.
Podemos resumir en las siguientes razones:
Son el lenguaje común entre economistas, y es el método dominante en todas las ciencias sociales.
Bajo los supuestos que veremos ahora, MCO es insesgado y consistente.
En el marco de otros estimadores insesgados, MCO es uno de los más eficiente. Siempre y cuando se mantengan algunas condiciones.
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:
- El modelo de regresión es lineal en los parámetros
\[Y_i = \beta_o + \beta_1X_i + u_i\]
- 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.
- 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)
- 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 \]
- 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\).
- 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.
- 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.
- 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.
- 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.
- 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.
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.
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.
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.
Normal Q-Q: Usado para examinar si los residuos están normalmente distribuídos. Es bueno su los puntos de los residuos siguen la linea recta dibujada. En nuestro caso se comienza a levantar así que hay indicios de que los residuos no están normalmente distribuídos.
Scale-location:: Usado para chequear la homogeneidad de la varianza de los residuos (u homocedasticidad). La línea horizontal con dispersión de puntos equivalente es un buen indicio de homogeneidad de varianza.
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.
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.
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.
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.
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.
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)
| factor | AME | lower | upper | p |
|---|---|---|---|---|
| Disca | 0.004 | -0.0038 | 0.0118 | 0.315 |
| escojefe | -0.000185 | -0.00117 | 0.000798 | 0.712 |
| intergrantesnucl | 0.0128 | 0.00906 | 0.0165 | 1.4e-11 |
| jefemujer | 0.00966 | -0.00231 | 0.0216 | 0.114 |
| jefemujermonoparental | 0.00897 | -0.00427 | 0.0222 | 0.184 |
| partadultmayor | 0.00231 | -0.00898 | 0.0136 | 0.689 |
| partmujer | -0.00519 | -0.0203 | 0.00993 | 0.501 |
| partninos | 0.029 | 0.00414 | 0.0539 | 0.0223 |
| PercapTotal | 7.48e-08 | 8.25e-09 | 1.41e-07 | 0.0276 |
| rural | 0.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)
| factor | AME | lower | upper | p |
|---|---|---|---|---|
| Disca | 0.00383 | -0.00414 | 0.0118 | 0.346 |
| escojefe | -0.000165 | -0.00116 | 0.000833 | 0.746 |
| intergrantesnucl | 0.0119 | 0.00836 | 0.0154 | 3.79e-11 |
| jefemujer | 0.0105 | -0.00105 | 0.0221 | 0.0747 |
| jefemujermonoparental | 0.00693 | -0.00569 | 0.0196 | 0.282 |
| partadultmayor | 0.000456 | -0.0112 | 0.0121 | 0.939 |
| partmujer | -0.00358 | -0.0188 | 0.0116 | 0.644 |
| partninos | 0.0293 | 0.00559 | 0.053 | 0.0154 |
| PercapTotal | 7.3e-08 | 9.99e-09 | 1.36e-07 | 0.0232 |
| rural | 0.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.
**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?
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.
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.
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.
¿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.