Modelos de regresión lineal

1 Introducción

El objetivo de este material es introducir al estudiante a los conceptos fundamentales de la Econometría, brindando herramientas para la comprensión de trabajos de economía aplicada, así como para el planteo de modelos econométricos sencillos y su estimación y diagnóstico. Se busca cubrir los aspectos básicos que sustentarán el tratamiento más avanzado de algunos tópicos en Econometría II y otras materias optativas. Paralelamente, se aspira a introducir al estudiante en el uso de paquetes estadísticos.

A lo largo de los capítulos se recorrerán los temas del curso de Econometría 1, dictado en la Facultad de Ciencias Económicas (FCEA) de la Universidad de la República (UdelaR). Se verá como instalar y trabajar con el software \(R\), en particular como descargar y trabajar con diferentes bases de datos. Los diferentes capítulos recorrerán los temas del curso implementando ejemplos concretos donde aprenderás a estimar modelos de regresión y a validarlos. Se comienza presentando el Modelo de Regresión Lineal Simple (MRLS), que es utilizado para la introducción de los elementos básicos en la construcción de un modelo econométrico, y los aspectos fundamentales de los distintos métodos de estimación. En base al Modelo de Regresión Lineal Múltiple (MRLM) se presentan en detalle los estimadores mínimo-cuadráticos, enfatizando en sus propiedades estadísticas y los supuestos en que se basan, tanto en muestras finitas como en contexto asintótico. La segunda parte del curso se dedica al análisis de las causas y consecuencias de la violación de cada uno de los supuestos clásicos, así como a las estrategias disponibles para el abordaje de los problemas asociados.

En este material iremos compartiendo códigos, que esperamos te resulten de utilidad, tanto para el curso, como para el posterior desarrollo de tu ejercicio profesional. Contiene ejemplos para aplicar modelos de regresión lineal con los datos del Wooldridge 2010-Introducción a la Econometría (4ta edición), texto de referencia del curso de Econometría 1: https://www.academia.edu/30200962/Introducci%C3%B3n_A_La_Econometr%C3%ADa_4edi_Wooldridge.

2 Presentación de R

2.1 ¿Qué son R y RStudio?

\(R\) es un lenguaje de programación que ofrece una gran variedad de funciones para realizar cálculos estadísticos y generar diversos gráficos a partir de los datos. Sin embargo, el gran potencial está en que, al ser libre y colaborativo, constantemente los usuarios están actualizando y ampliando la cantidad de funciones que presenta. Hoy en día podemos realizar desde operaciones básicas sobre los datos hasta aplicar algoritmos de inteligencia artificial.

Es un lenguaje interpretado, esto significa que los comandos escritos en el teclado son ejecutados directamente sin necesidad de crear ejecutables.

\(R\) es muy flexible, potente y uno de los software mas profesionales que existen actualmente para realizar tareas estadísticas de todo tipo, desde las más elementales, hasta las más avanzadas. Cuenta, además, con la ventaja de permitir el manejo de grandes volúmenes de datos, permite replicar procedimientos anteriores, generar y automatizar informes ahorrando mucho tiempo.

Viene con una librería base, pero existen a libre disposición muchísimas librerías específicas para solucionar problemas puntuales.

\(RStudio\) es una interfaz libre y gratuita que nos permite explotar todo el potencial que tiene el lenguaje de programación \(R\). RStudio es un Ambiente de Desarrollo Integrado (IDE, Integrated Development Environment). Esto quiere decir que expande la funcionalidad y flexibilidad de la interfaz de usuario de R. Se integra todo nuestro proceso de trabajo en un mismo ambiente. Facilita mucho pues es mas amigable e intuitivo.

Las órdenes elementales de \(R\) consisten en expresiones o asignaciones. Una expresión, se evalúa, se imprime el resultado y su valor se pierde.Una asignación en cambio, crea un objeto y no se imprime el resultado.

Cada línea solo puede contener 128 caracteres. Si queremos escribir más, podemos utilizar otra línea que comienza con el símbolo + en lugar de el promt (>).

Mientras que otros programas estadísticos muestran directamente los resultados de un análisis, \(R\) guarda estos resultados como un objeto. \(R\) es un programa orientado a objetos: pueden ser variables, datos, funciones, resultados, etc.

Permite trabajar de manera que nuestro trabajo sea replicable, garantizando la obtención de los mismos resultados en las mismas condiciones. Es una forma de documentar procesos, da garantías y transparencia. En general vamos a cargar una base de datos ya existente, luego vamos a realizar una serie de transformaciones (generar nuevas variables, filtros por determinadas condiciones, cálculos específicos) y vamos a obtener algunos resultados. Recomendamos grabar la base de datos original y la secuencia de órdenes (sintaxis en un archivo tipo Rscript). Nunca vamos a pasar por encima la base original.

A su vez, dentro de \(RStudio\), hay diferentes formatos de archivos (RMarkdown, RScript, RSweave, RNotebook, etc.) y su elección depende del objetivo que tengamos. En nuestro caso, a lo largo del manual trabajaremos con el formato RMarkdown, un tipo de documento de \(RStudio\) que integra texto con código de \(R\) y nos permite generar informes como este a partir de los datos.

2.2 Descarga e instalación de R y Rstudio

Deberás descargar e instalar \(R\) y \(Rstudio\).

• Descarga \(R\) aquí: https://www.r-project.org/

• Descarga \(RStudio\) aquí: https://www.rstudio.com/products/rstudio/download/

• Aquí podrás ver un video tutorial para instalación de https://www.youtube.com/watch?v=ZTzbiAznSjc

• Si queremos trabajar en la nube (sin necesidad de instalar nada en la PC), podremos hacerlo aquí: https://rstudio.cloud

2.3 Características generales

Por defecto, abriendo \(R\) se abre una sola ventana: la consola o ventana de comandos de \(R\). Aquí se introducen los comandos y se visualizan los resultados.

Al final aparece una línea en blanco con el símbolo > (Prompt). A partir de aquí el \(R\) espera que escribamos comandos o instrucciones.

Si en vez de aparece el símbolo > , aparece un + es porque hay una sentencia no finalizada.

Para ejecutar un comando desde la consola basta apretar ENTER.

Se puede trabajar en la consola o crear un programa (script): una secuencia de órdenes podremos guardar, modificar y permitirá replicar y obtener los resultados una y otra vez (exactamente de la misma forma).

Para que \(R\) las ejecute primero debemos seleccionarlas y luego realizar cualquiera de las acciones siguientes: CTRL + R o presionar el ícono Run.

Como se ve en la imagen, se puede dividir (vertical u horizontalmente) la interfaz en 4 partes o paneles:

Panel de Edición: en el cuadrante superior izquierdo vamos a estar creando y modificando nuestro código (script). Aquí también podría haber otro formato de archivo R, como por ejemplo RMarkdown o RNotebook.

Entorno de Variables: en el cuadrante superior derecho veremos todos las bases y los objetos que vayamos cargando y creando. Desde aquí también podremos importar o eliminar base u objetos. Desde la pestaña Historial podremos consultar el historial de comandos y funciones que fuimos utilizando en el proyecto.

Consola: en el cuadrante inferior izquierdo irá apareciendo todo lo que ejecutemos tanto desde el panel de edición como desde el entorno de variables. Aquí también podemos escribir líneas de código que queremos que se ejecuten. No es recomendable pues no se guardan estas órdenes, solo se ejecutan.

Panel de Utilidades: en la ventana inferior derecha se pueden ver varias cosas:

• Files: el directorio donde estamos trabajando.
• Plots: las visualizaciones/gráficos que se van generando.
• Packages: los paquetes de \(R\) disponibles.
• Help: una sección de ayuda donde podemos consultar información de las funciones.
• Viewer: U¿un visor HTML para ver los gráficos interactivos o animados que hayamos hecho.

Cuidado! Un error muy frecuente es pensar que se está escribiendo un comando en la consola, cuando en realidad se está escribiendo en el script.

2.3.1 Errores, mensajes y advertencias

Ver texto rojo en la consola no es necesariamente malo (un error). \(R\) mostrará texto rojo en el panel de la consola en tres situaciones diferentes:

Errores: cuando el texto rojo es un error legítimo, estará precedido por “Error en…” e intentará explicar qué salió mal. Generalmente, cuando hay un error, el código no se ejecutará. Uno muy común es cuando la función que se quiere utilizar necesita instalar previamente un paquete y no se ha hecho. En este caso dirá algo como esto: Error in ggplot(...) : could not find function "ggplot", significa que la función ggplot() no es accesible porque el paquete que contiene la función ggplot2() no se cargó con library(ggplot2).

Advertencias: cuando el texto rojo es una advertencia, estará precedido por “Advertencia:” y \(R\) intentará explicar por qué hay una advertencia. En general, el código seguirá funcionando, pero con algunas advertencias.

Mensajes: cuando el texto rojo no comienza con “Error” o “Advertencia”, es solo un mensaje amistoso. Son mensajes de diagnóstico y no impiden que el código continúe funcionando. Veremos estos mensajes al instalar paquetes por ejemplo.

2.3.2 Librerías y paquetes

• Un paquete es una colección de funciones, datos y documentación que permite extender las capacidades de \(R\) base.
• Las funciones en \(R\) se organizan en paquetes según temas o disciplinas.
• Los paquetes debe instalarse en \(R\) una vez pero es necesario cargarlos en cada sesión (cada vez que abrimos \(R\)).
• Podemos consultar los paquetes disponibles escribiendo la función available.packages()
• También podemos generar nuestras propias funciones e incluso crear un paquete de \(R\).
• En particular Tidyverse es un conjunto de paquetes que comparten una misma filosofía de datos y están diseñados para trabajar conjuntamente de forma natural.

Instalemos el paquete tidyverse:

install.packages("tidyverse")

Una vez instalado el paquete, cargémoslo:

library(tidyverse)

Otra forma de instalar paquetes es la siguiente:

Haz click en Packages > Install, y se abre una ventana emergente. Luego escribe el nombre del paquete que quieras instalar y finalmente click en Install.

Uno de los paquetes que existe dentro de tidyverse es el paquete dplyr. Algunas funciones que contiene este paquete y estaremos usando a lo largo de este documento son:

• filter(): filtra o elige las observaciones por sus valores.

• arrange(): re-ordena las filas.

• select(): selecciona las variables por sus nombres.

• mutate(): crea nuevas variables con transformaciones de variables existentes.

• summarise(): contrae muchos valores en un solo resumen

• group_by(): cambia el alcance de cada función para que actúe ya no sobre todo el conjunto de datos sino de grupo en grupo.

Finalmente te comentamos que existen paquetes que contiene cursos interactivos para aprender \(R\) desde la consola, como lo es el swirl:

install.packages("swirl")
library(swirl)

2.4 Comandos básicos indispensables

• Puedes usar \(R\) como una calculadora:

1 / 200 * 30

## [1] 0.15

(59 + 73 + 2) / 3

## [1] 44.66667

sin(pi / 2)

## [1] 1

• Puedes crear objetos nuevos usando <-:

  • nombre_objeto <- valor

  • Alt + - (signo menos).

  • Los nombres de los objetos deben comenzar con una letra y solo pueden contener letras, números, guiones bajo o puntos.

x <- 3 * 4

• Puedes examinar un objeto escribiendo su nombre:

x 

## [1] 12

• Puedes crear un vector con la función c():

v1 <- c(1, 2, 3, 4.5, "seis")
v1

## [1] "1"    "2"    "3"    "4.5"  "seis"

v2 <- c(3, 7:10)
v2

## [1]  3  7  8  9 10

• Puedes crear matrices.

2.5 Clases de objetos

• Vectores c()

• Matrices y Arrays: Todos los elementos son numéricos matrix()

• Listas: pueden contener elementos de distinto tipo list()

• Factores: se usan para variables categóricas factor()

• Data Frames: se usan para almacenar datos en forma de tablas (filas / columnas). Pueden contener elementos de distinto tipo data.frame()

• Tibbles: dataframes con algunas particularidades por ejemplo no convierten por defecto vectores de texto en factores

• Funciones: argumento, cuerpo, resultado

La clase de un objeto afecta cómo las funciones trabajan con el mismo. Por ejemplo, la función mean() puede generar salidas coherentes con objetos numéricos o lógicos, pero no con caracteres. Lo mismo con la multiplicación.

2.6 ¿Qué es un project?

Primero definamos lo que es un Directorio de trabajo (working directory). El directorio es donde \(R\) busca y guarda los archivos. Al principio de cada sesión podrás especificar en qué directorio trabajarás, pues el \(R\) irá a buscar los archivos allí. A menos que sea cambiado, este directorio se mantendrá hasta el fin de la sesión.

La función getwd() indica la ruta de donde estamos trabajando. Y para cambiar el directorio se utiliza la función setwd(“…/../..”)

Un proyecto permite mantener todos los archivos asociados en un mismo lugar — datos de entrada, scripts, resultados, gráficos. Se puede compartir y evita especificar el directorio de trabajo.

Entonces, hagamos un proyecto:

Haz clic en File > New Project, y después:

2.7 Importar y exportar datos

2.7.1 Para importar datos

Para ello se puede utilizar la librería readr del paquete tidyverse que contiene las funciones explicadas a continuación.

• read_csv() lee archivos delimitados por coma.

• read_csv2() lee archivos separados por punto y coma.

• read_tsv() lee archivos delimitados por tabulaciones.

• read_delim() lee archivos con cualquier delimitador.

El archivo notas.csv, incluye una muestra de los resultados de las personas que rindieron la segunda revisión de Econometría 1 para 2017-2018 y lo pueden encontrar aquí: https://eva.fcea.udelar.edu.uy/course/view.php?id=183&section=10

Otra forma de importar datos es la siguiente: haz click en Environment > Import Dataset.

Se abre una ventana emergente:

2.7.2 Para exportar datos

write(x, file = "data")
write.table(x, file = paste("../../.."), "x._"))

En el código anterior: * x es el nombre del objeto. * en file = paste("../../..") ponemos la ruta donde queremos guardar el archivo. * x._ es el nombre que le ponemos al archivo para guardarlo.

Específicamente, para guardar un archivo en formato .csv, el cual puede ser leido con excel, usamos la función write.csv():

write.csv(x, file=paste("../../.."), "x.csv"))

2.8 Algunas funciones básicas

Una vez importados los datos podemos utilizar estas funciones:

• view(): visualiza los datos como una “tabla”.

• names(): muestra o establece los nombres de un objeto.

• head(): devuelve las primeras 6 partes de un objeto.

• tail(): devuelve las ultimas 6 partes de un objeto.

• dim(): establece la dimensión de un objeto.

• str(): muestra de forma compacta la estructura interna de un objeto.

• summary(): produce resúmenes de los datos o de los resultados de funciones de ajuste de modelos.

2.9 Para buscar ayuda

Puedes encontrar material de ayuda en el sitio web de la versión en español de “R for Data Science”, de Hadley Wickham y Garrett Grolemund. https://es.r4ds.hadley.nz/

Tambien existen los llamados “Machetes” o cheatsheets de Rstudio, los cuales sintetizan en una dos páginas todas las funciones de cada paquete de interés. https://www.rstudio.com/resources/cheatsheets/

2.9.1 Buscar ayuda en R

• En la consola escribe ? + el nombre de la función.
• En la consola escribe ?? + algo parecido al nombre de la función.
• En la pestaña Help como se muestra en la imágen.

2.10 Consejos prácticos para trabajar

• Usá \(R\) todo lo que puedas.

• No desistas, la curva de aprendizaje es empinada, pero vale la pena el esfuerzo.

• No uses otros programas, lo que puedas hacer fuera lo puedes replicar en R.

• Recurrí a los foros y a la ayuda de \(R\) para encontrar las soluciones.

• Prestá atención a los mensajes de error y advertencia.

• Reutilizá sintaxis existente (tuyas o de otros usuarios).

• Aprendé haciendo y equivocándote.

• Recuerda que el código no hace lo que quiere, sino lo que tu le indicas.

• Guardá la base original y el script, con esto sería suficiente para replicar tus resultados.

• Jamás le pases por encima a la base ni a las variables originales.

• Trabaja con \(R\) project, te evitará la definición del directorio de trabajo en \(R\).

• Escribí tus sintaxis en un script y comentalas detalladamente. En \(R\) esto se hace con el uso del símbolo de número o hashtag (#). Todo el texto a la derecha de ese símbolo será ignorado por \(R\). Esto es útil al compartir nuestro código con alguien más, pero también al regresar a código anterior.

3 Introducción al curso de Econometría 1 con R

Este es un documento de trabajo, contiene ejemplos para aplicar modelos de regresión lineal con los datos del Wooldridge 2010-Introducción a la Econometría (4ta edición), texto de referencia del curso de Econometría 1. El mismo fue elaborado en Markdown utilizando \(RStudio\) (software que vemos en el curso). En el siguiente link se encuentran detalles para utilizar \(RMarkdown\) http://rmarkdown.rstudio.com.

El \(R\) permite instalar la librería Wooldridge de la siguiente manera:

if(!require('wooldridge')) {
  install.packages('wooldridge')
  library('wooldridge')
}

Como vimos en la sección anterior, esto es equivalente a lo siguiente:

install.packages("wooldridge") # la instalo si es la primera vez que la llamo.
library(wooldridge) # la llamo una vez que ya la tengo instalada. 

En esta librería encontraremos todas las bases de datos que contiene el libro. Se deberá utilizar el mismo nombre con el que son referenciadas en el texto. A modo de ejemplo, trabajaremos a lo largo de todo el docuemnto, con la base wage1 que contiene datos de la Encuesta de población de 1976, recopilados por Henry Farber cuando él y Wooldridge fueron colegas en el MIT (Instituto de Tecnología de Massachusetts) en 1988.

4 Modelo de regresión lineal simple (MRLS)

4.1 Cargar y manipular la base

Para cargar los datos se escribe

data("wage1") # carga la base salarios
base <- wage1  # la asigno al objeto llamado base con el que trabajaremos

Acceso a la base y solicito los nombres de las variables:

library(wooldridge)
base <- data("wage1")
base <- wage1
#View(base) 
names(base)

##  [1] "wage"     "educ"     "exper"    "tenure"   "nonwhite" "female"  
##  [7] "married"  "numdep"   "smsa"     "northcen" "south"    "west"    
## [13] "construc" "ndurman"  "trcommpu" "trade"    "services" "profserv"
## [19] "profocc"  "clerocc"  "servocc"  "lwage"    "expersq"  "tenursq"

Las variables que utilizaremos son las siguientes:

• wage: salario promedio por hora
• educ: años de educación
• exper: años de experiencia potencial
• ternure: años con el empleador actual (antigüedad)
• nonwhite: es igual a 1 si la persona no es blanca
• female: es igual a 1 si la persona es mujer
• married: es igual a 1 si la persona es casada

En primer lugar queremos cambiar el nombre de la variable que está en la posición 4:

names(base)[4] <- "antiguedad"

Seleccionamos las variables que nos interesan por la posición de las columnas en el data frame. Indicamos que queremos las variables que están desde la posición 1 hasta la 7 y las que están en la posición 22, 23 y 24 de la siguiente forma:

base1 <- base[,c(1: 7, 22, 23, 24)] 

Miramos las observaciones de las primeras 5 personas de la base:

head(base1, n = 5)

##   wage educ exper antiguedad nonwhite female married    lwage expersq tenursq
## 1 3.10   11     2          0        0      1       0 1.131402       4       0
## 2 3.24   12    22          2        0      1       1 1.175573     484       4
## 3 3.00   11     2          0        0      0       0 1.098612       4       0
## 4 6.00    8    44         28        0      0       1 1.791759    1936     784
## 5 5.30   12     7          2        0      0       1 1.667707      49       4

También se pueden seleccionar las variables antes mencionadas utilizando el nombre de las mismas:

datos1 <- base[,c("wage","educ","exper", "antiguedad" )] 

head(datos1, n = 5)

##   wage educ exper antiguedad
## 1 3.10   11     2          0
## 2 3.24   12    22          2
## 3 3.00   11     2          0
## 4 6.00    8    44         28
## 5 5.30   12     7          2

Otra forma mas sencilla aún, es utilizar la función select() del paquete dplyr. Se usa el comando llamado pipe %>%, el cual enlaza las órdenes de manera secuencial.

Si es la primera vez que se usa recordar instalarlo previamente con el comando install.packages("dplyr"). Este paquete es uno de los que conforman el paquete tidyverse, por lo que también funcionaría con el comando install.packages("tidyverse") y luego library(tidyverse).

library(dplyr)

datos2 <- base %>% select("wage", "educ", "exper", "antiguedad")

head(datos2, n=5)

##   wage educ exper antiguedad
## 1 3.10   11     2          0
## 2 3.24   12    22          2
## 3 3.00   11     2          0
## 4 6.00    8    44         28
## 5 5.30   12     7          2

Se puede combinar con un filtro por filas con la función filter() del paquete dplyr. A modo de ejemplo, seleccionamos solo a las personas casadas, es decir aquellas en que la variable married es igual a 1.

datos3 <- base %>% 
  select("wage", "educ", "exper", "antiguedad", "married") %>% 
  filter(married == 1)

4.2 Investigar las correlaciones entre variables

Para investigar si hay correlación entre alguna de las variables se puede realizar un gráfico en el que se presenta la dispersión para cada par de variables.

plot(datos2)

Y también calcular la matriz de correlaciones de las variables que figuran en datos2.

cor(datos2)

##                 wage        educ      exper  antiguedad
## wage       1.0000000  0.40590333  0.1129034  0.34688957
## educ       0.4059033  1.00000000 -0.2995418 -0.05617257
## exper      0.1129034 -0.29954184  1.0000000  0.49929145
## antiguedad 0.3468896 -0.05617257  0.4992914  1.00000000

Existe una correlación positiva entre el salario y la educación (0.4059)

Creamos el diagrama de dispersión entre salario y educación utilizando las funciones de la librería ggplot2.

#install.packages("ggplot2") 
library(ggplot2)

ggplot(datos2, aes(x = educ, y = wage)) + 
  geom_point() + theme_light() +
  ggtitle("Relación entre salario y educación")

4.3 Estimar un modelo de regresión lineal simple (MRLS)

Estimamos un MRLS con el método mínimos cuadrados ordinarios (en adelante MCO) que explique los salarios en función de los años de educación de las personas.

\[\begin{equation} wage_i = \beta _0 + \beta _1 educ_i + u_i \\ \end{equation}\]

Para estimar un MRLS en R se utiliza la función lm() del R base que explicaremos a continuación.

mod1 <- lm(wage ~ educ, data=datos2)
summary(mod1) # para imprimir la salida

## 
## Call:
## lm(formula = wage ~ educ, data = datos2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3396 -2.1501 -0.9674  1.1921 16.6085 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.90485    0.68497  -1.321    0.187    
## educ         0.54136    0.05325  10.167   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.378 on 524 degrees of freedom
## Multiple R-squared:  0.1648, Adjusted R-squared:  0.1632 
## F-statistic: 103.4 on 1 and 524 DF,  p-value: < 2.2e-16

anova(mod1) # para imprimir el análisis de varianzas

## Analysis of Variance Table
## 
## Response: wage
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## educ        1 1179.7 1179.73  103.36 < 2.2e-16 ***
## Residuals 524 5980.7   11.41                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4.3.1 Objeto lm y sus componentes

Los resultados fueron asignados al objeto mod1 de clase lm. Para consultar la clase:

class(mod1)

## [1] "lm"

Dentro del objeto se encuentran 12 elementos que se pueden utilizar llamándolos por los siguientes nombres:

names(mod1)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"

Se puede acceder a ellos como a los objetos del data frame, utilizando el operador $ entre el objeto y el elemento. Por ejemplo, para extraer los coeficientes estimados se escribe lo siguiente:

mod1$coefficients

## (Intercept)        educ 
##  -0.9048516   0.5413593

Para consultar la estimación de un coeficiente de regresión se utilizan paréntesis rectos y se indica la ubicación del mismo dentro de la salida del summary().

ahat <- coef(summary(mod1))[1,1]
ahat

## [1] -0.9048516

bhat <- coef(summary(mod1))[2,1]
bhat

## [1] 0.5413593

Miramos algunos componentes del modelo y los asignamos a nuevos objetos de esta forma:

names(mod1)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"

coeficientes <- mod1$coefficients
ygorro <- mod1$fitted.values
resid <- mod1$residuals

Por último, agregamos al conjunto de datos las predicciones que se obtienen con la función predict().

datos2$predicciones <- predict(mod1) 
head(datos2, 5)

##   wage educ exper antiguedad predicciones
## 1 3.10   11     2          0     5.050100
## 2 3.24   12    22          2     5.591459
## 3 3.00   11     2          0     5.050100
## 4 6.00    8    44         28     3.426022
## 5 5.30   12     7          2     5.591459

4.4 Diagrama de dispersión y recta de regresión

Creamos el gráfico de dispersión entre el salario y la educación, igual al que realizamos de forma exploratoria anteriormente, pero incluimos además la estimación de la recta de regresión con un geom_smooth().

ggplot(datos2, aes(x = educ, y = wage)) + 
  geom_point() +
  geom_smooth(method = 'lm', formula = y ~ x, se = FALSE, col = 'dodgerblue1') +
  theme_light() +
  ggtitle("Relación entre salario y educación")

Es posible obtener un gráfico de dispersión interactivo con el cual posicionándose encima de cada observación se ven los valores de (x, y) para cada uno de los individuos. Para construir dicho gráfico se necesita la función ggplotly() del paquete plotly.

#install.packages("plotly")
library(plotly)

ggplotly(data = datos2, x = ~ educ, y = ~ wage)

4.5 Trabajar con los residuos del modelo

En \(R\) calculamos los residuos del modelo y los agregamos al conjunto de datos (datos2) de la siguiente forma:

datos2$residuos <- datos2$wage - datos2$predicciones
head(datos2, 5)

##   wage educ exper antiguedad predicciones   residuos
## 1 3.10   11     2          0     5.050100 -1.9501003
## 2 3.24   12    22          2     5.591459 -2.3514594
## 3 3.00   11     2          0     5.050100 -2.0501002
## 4 6.00    8    44         28     3.426022  2.5739776
## 5 5.30   12     7          2     5.591459 -0.2914593

Comparamos los residuos calculados manualmente con los que nos dio el modelo. \[\widehat{u}_{i}= y_{i}-(\widehat{\beta}_0 + \widehat{\beta}_1 x_{i}) = y_{i} - \widehat{y}_i \]

resid <- mod1$residuals
datos2$resid <- resid
datos2$verif <- datos2$residuos - datos2$resid

Agrego al gráfico de dispersión los valores estimados de \(y\) (\(\widehat{y}_i\)) en rojo y muestro los residuos (\(\widehat{u}_{i}\))

ggplot(datos2, aes(x = educ, y = wage)) +
  geom_smooth(method = "lm", se = FALSE, color = "lightgrey") +
  geom_segment(aes(xend = educ, yend = predicciones), col = 'red', lty = 'dashed') +
  geom_point() +
  geom_point(aes(y = predicciones), col = 'red') +
  theme_light()

Realizamos un gráfico de dispersión para inspeccionar de forma gráfica los residuos.

plot(datos2$residuos)

Graficamos las funciones de densidad estimadas.

library(e1071)  # para la función skewness
par(mfrow = c(1, 2))  # divide el área de gráficos en 2 columnas

plot(density(datos2$wage), main = "Density Plot: salario", ylab = "Frequency", sub = paste("Skewness:", round(e1071::skewness(datos2$wage), 2)))  # density para 'salario'

polygon(density(datos2$wage), col = "red")

plot(density(datos2$residuos), main = "Density Plot: residuos", ylab = "Frequency", sub = paste("Skewness:", round(e1071::skewness(datos2$residuos), 2)))  # density para los 'residuos'

polygon(density(datos2$residuos), col = "red")

4.6 Predicción e intervalos de confianza

A partir de los datos con los que se viene trabajando, vamos a generar la predicción puntual y el intervalo de confianza para una educación de 10 años.

Primero, calculamos el salario esperado para una persona con la educación promedio: (\(y/x=\overline{x}\)).

mean(datos2$educ)

## [1] 12.56274

sal_pred <- 0.90485 + 0.054136 * 12.56
sal_pred

## [1] 1.584798

En segundo lugar, calculamos el intervalo de confianza para una persona con la educación promedio (\(x = 12.56\)) El argumento interval = prediction devuelve el valor para la predicción puntual.

nuevo <- data.frame(educ = 12.56)
future_y <- predict(object = mod1, newdata = nuevo, interval = "prediction", level = 0.95)

Luego, generamos el intervalo de confianza para \(E(y/x)\). en este caso, debemos cambiar el argumento a interval = "confidence".

future_esp_y <- predict(object = mod1, newdata = nuevo, interval = "confidence", level = 0.95)
future_esp_y <- as.data.frame(future_esp_y)

IC_inf_esp_y <- future_esp_y$lwr
IC_sup_esp_y <- future_esp_y$upr

Agregamos los extremos de los intervalos de confianza a la base con la que estamos trabajando:

nuevos_datos <- cbind(datos2, future_y, IC_inf_esp_y, IC_sup_esp_y)

Finalmente, generamos los gráficos correspondientes con los intervalos de confianza para la predicción puntual (IC_y) y para el valor esperado (IC_esp_y) con el siguiente código:

IC_y <- ggplot(nuevos_datos, aes(x = educ, y = wage)) +
  geom_point() +
  geom_line(aes(y = lwr), color = "red", linetype = "dashed") +
  geom_line(aes(y = upr), color = "red", linetype = "dashed") +
  geom_smooth(method = lm, formula = y ~ x, se = TRUE, level = 0.95, col = 'blue', fill = 'pink2') +
  theme_light() + 
  ggtitle("Predicción de y al 95%")

IC_esp_y <- ggplot(nuevos_datos, aes(x = educ, y = wage)) +
  geom_point() +
  geom_line(aes(y = IC_inf_esp_y), color = "blue", linetype = "dashed") +
  geom_line(aes(y = IC_inf_esp_y), color = "blue", linetype = "dashed") +
  geom_smooth(method = lm, formula = y~x, se = TRUE, level = 0.95, col = 'blue', fill = 'pink2') +
  theme_light() + 
  ggtitle("Predicción de E(y/x) al 95%")

Imprimimos los gráficos uno al lado del otro, para poder compararlos mejor. ¿Cuál de los dos tiene mayor amplitud?

#install.packages("gridExtra")
library(gridExtra)

grid.arrange(IC_esp_y, IC_y, ncol = 2, nrow = 1)

Veamos, a través de la comparación de dos gráficos el impacto que tiene el nivel de confianza en la amplitud de los intervalos. Para ello, tendremos que descargar e instalar algunas librerías. Se presenta primero el código y luego los gráficos obtenidos.

#install.packages("jtools")
#install.packages("ggstance")
#install.packages("broom.mixed")

library(jtools)
library(ggstance)
library(broom.mixed)

a <- plot_summs(mod1, escala = TRUE, plot.distributions = TRUE, inner_ci_level = .90)
b <- plot_summs(mod1, escala = TRUE, plot.distributions = TRUE, inner_ci_level = .7) 

grid.arrange(a, b, ncol = 1, nrow = 2)

5 Modelos de regresión lineal múltiple (MRLM)

Estimaremos modelos de regresión lineal múltiple utilizando la misma base de datos de Wooldridge utilizada como ejemplo en el capítulo anterior, wage1, que contiene datos de la Encuesta de población de 1976.

5.1 Cargar y manipular la base

Accedemos a los datos, los asignamos al objeto llamado base con el que trabajaremos y solicitamos los nombres de las variables:

library(wooldridge)

data("wage1") 
base <- wage1 
names(base) 

##  [1] "wage"     "educ"     "exper"    "tenure"   "nonwhite" "female"  
##  [7] "married"  "numdep"   "smsa"     "northcen" "south"    "west"    
## [13] "construc" "ndurman"  "trcommpu" "trade"    "services" "profserv"
## [19] "profocc"  "clerocc"  "servocc"  "lwage"    "expersq"  "tenursq"

Las variables que utilizaremos en este capítulo son las siguientes:

• wage: salario promedio por hora
• educ: años de educación
• exper: años de experiencia potencial
• ternure: años con el empleador actual (antigüedad)
• nonwhite: es igual a 1 si la persona no es blanca
• female: es igual a 1 si la persona es mujer
• married: es igual a 1 si la persona es casada
• lwage: logaritmo del salario, \(log(wage)\)
• expersq: el cuadrado de los años de experiencia potencial, \(exper^2\)
• tenursq: el cuadrado de la antigüedad en el trabajo actual, \(ternure^2\)

names(base)[4] <- "antiguedad"

Seleccionamos las variables que nos interesan por la posición de las columnas en el data frame. Indicamos que queremos las variables que están desde la posición 1 hasta la 7 y las que están en la posición 22, 23 y 24. Las asignamos al objeto datos.

datos <- base[,c(1: 7, 22, 23, 24)] 

Para mirar la dimensión de la base resultante solicitamos la dimensión de la misma.

dim(datos)

## [1] 526  10

Podríamos solicitar el número de filas (personas en este caso) y de columnas (variables) por separado con las siguientes funciones:

nrow(datos)

## [1] 526

ncol(datos)

## [1] 10

Miramos las observaciones de las primeras 5 personas de la base.

head(datos, n = 5)

##   wage educ exper antiguedad nonwhite female married    lwage expersq tenursq
## 1 3.10   11     2          0        0      1       0 1.131402       4       0
## 2 3.24   12    22          2        0      1       1 1.175573     484       4
## 3 3.00   11     2          0        0      0       0 1.098612       4       0
## 4 6.00    8    44         28        0      0       1 1.791759    1936     784
## 5 5.30   12     7          2        0      0       1 1.667707      49       4

Solicitamos las estadísticas básicas de las variables de la base.

summary(datos)

##       wage             educ           exper         antiguedad    
##  Min.   : 0.530   Min.   : 0.00   Min.   : 1.00   Min.   : 0.000  
##  1st Qu.: 3.330   1st Qu.:12.00   1st Qu.: 5.00   1st Qu.: 0.000  
##  Median : 4.650   Median :12.00   Median :13.50   Median : 2.000  
##  Mean   : 5.896   Mean   :12.56   Mean   :17.02   Mean   : 5.105  
##  3rd Qu.: 6.880   3rd Qu.:14.00   3rd Qu.:26.00   3rd Qu.: 7.000  
##  Max.   :24.980   Max.   :18.00   Max.   :51.00   Max.   :44.000  
##     nonwhite          female          married           lwage        
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :-0.6349  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.: 1.2030  
##  Median :0.0000   Median :0.0000   Median :1.0000   Median : 1.5369  
##  Mean   :0.1027   Mean   :0.4791   Mean   :0.6084   Mean   : 1.6233  
##  3rd Qu.:0.0000   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.: 1.9286  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   : 3.2181  
##     expersq          tenursq       
##  Min.   :   1.0   Min.   :   0.00  
##  1st Qu.:  25.0   1st Qu.:   0.00  
##  Median : 182.5   Median :   4.00  
##  Mean   : 473.4   Mean   :  78.15  
##  3rd Qu.: 676.0   3rd Qu.:  49.00  
##  Max.   :2601.0   Max.   :1936.00

5.2 Investigar las correlaciones entre variables

Calculamos las correlaciones de los pares de variables, indicando que queremos trabajar con 3 decimales:

round(cor(datos, method = "pearson"), 3)

##              wage   educ  exper antiguedad nonwhite female married  lwage
## wage        1.000  0.406  0.113      0.347   -0.039 -0.340   0.229  0.937
## educ        0.406  1.000 -0.300     -0.056   -0.085 -0.085   0.069  0.431
## exper       0.113 -0.300  1.000      0.499    0.014 -0.042   0.317  0.111
## antiguedad  0.347 -0.056  0.499      1.000    0.012 -0.198   0.240  0.326
## nonwhite   -0.039 -0.085  0.014      0.012    1.000 -0.011  -0.062 -0.039
## female     -0.340 -0.085 -0.042     -0.198   -0.011  1.000  -0.166 -0.374
## married     0.229  0.069  0.317      0.240   -0.062 -0.166   1.000  0.271
## lwage       0.937  0.431  0.111      0.326   -0.039 -0.374   0.271  1.000
## expersq     0.030 -0.331  0.961      0.459    0.009 -0.028   0.217  0.023
## tenursq     0.267 -0.069  0.423      0.922   -0.007 -0.176   0.167  0.236
##            expersq tenursq
## wage         0.030   0.267
## educ        -0.331  -0.069
## exper        0.961   0.423
## antiguedad   0.459   0.922
## nonwhite     0.009  -0.007
## female      -0.028  -0.176
## married      0.217   0.167
## lwage        0.023   0.236
## expersq      1.000   0.414
## tenursq      0.414   1.000

Con la función attach(datos) logramos que el R reconozca el nombre de las variables del objeto datos sin necesidad de citar el nombre del data frame y luego el signo de $. A partir de este punto, realizaremos los diferentes comandos teniendo en cuenta que se aplicó esta función.

attach(datos)

Exploramos, mediante un gráfico en tres dimensiones, la relación entre el salario, los años de educación y la antigüedad en el cargo que tienen las personas. Realizamos un gráfico interactivo que nos permite cambiar la perspectiva del gráfico con el puntero del mouse.

library(ggplot2)
library(plotly)

plot_ly(x = educ, y = antiguedad, z = wage, type = "scatter3d", color = wage) %>% 
  layout(scene = list(xaxis = list(title = 'educación (en años)'),
                      yaxis = list(title = 'antiguedad (en años)'),
                      zaxis = list(title = 'Salario (en USD/h)')))

library(scatterplot3d)
graf <- scatterplot3d(x = educ, y = antiguedad, z = wage, pch = 16, 
              cex.lab = 1, highlight.3d = TRUE, type = "h", 
              xlab = 'Años de educación',
              ylab = 'Antiguedad (años)', 
              zlab = 'Salario (USD/h)')

5.3 Estimar un modelo para el salario

Estimamos un primer modelo en el que se explica el nivel del salario en función de los años de educación y de la experiencia. El modelo a estimar es el siguiente:

\[\begin{equation} wage_i = \beta_0 + \beta_1 educ_i + \beta_2 antiguedad_i + u_i \\ \end{equation}\]

mod1 <- lm(wage ~ educ + antiguedad, data = datos)
summary(mod1)

## 
## Call:
## lm(formula = wage ~ educ + antiguedad, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.1438 -1.7288 -0.6372  1.2575 14.7482 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.22162    0.64015   -3.47 0.000563 ***
## educ         0.56914    0.04881   11.66  < 2e-16 ***
## antiguedad   0.18958    0.01871   10.13  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.092 on 523 degrees of freedom
## Multiple R-squared:  0.3019, Adjusted R-squared:  0.2992 
## F-statistic: 113.1 on 2 and 523 DF,  p-value: < 2.2e-16

Ajustamos el plano de regresión estimado.

graf <-  scatterplot3d(x = educ, y = antiguedad, z = wage, pch = 16, 
                     cex.lab = 1, highlight.3d = TRUE, type = "h", 
                     xlab = 'Años de educación',
                     ylab = 'Experiencia (años)', 
                     zlab = 'Salario (USD/h)')

graf$plane3d(mod1, lty.box = "solid", col = 'mediumblue')

5.4 Análisis de los residuos del modelo

Calculamos con la función deviance() la suma de los cuadrados de los residuos

SCR <- deviance(mod1) #Suma del cuadrado de los residuos
SCR

## [1] 4998.965

gl <- df.residual(mod1) #Grados de libertad del modelo
gl

## [1] 523

sigma2_gorro <- (SCR/gl) #Sigma gorro cuadrado
sigma2_gorro

## [1] 9.55825

Otra manera es utilizar la función sigma() que devuelve el error estándar de los residuos del modelo.

sigma2_gorro <- (sigma(mod1))^2 
sigma2_gorro

## [1] 9.55825

También podemos utilizar la función anova() que nos devuelve lo siguiente:

anova(mod1)

## Analysis of Variance Table
## 
## Response: wage
##             Df Sum Sq Mean Sq F value    Pr(>F)    
## educ         1 1179.7 1179.73  123.43 < 2.2e-16 ***
## antiguedad   1  981.7  981.72  102.71 < 2.2e-16 ***
## Residuals  523 4999.0    9.56                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A continuación inspeccionamos la normalidad de los residuos mediante análisis gráficos similar a como lo hicimos con la estimación del modelo lineal simple, con funciones de la librería e1071.

library(e1071)  # librería necesaria para utilizar la función skewness()

plot(density(mod1$residuals), 
     main = "Density Plot: residuos", ylab = "Frequency", 
     sub = paste("Skewness:", round(e1071::skewness(mod1$residuals), 2)))  

polygon(density(mod1$residuals), col = "red")

Solicitamos los gráficos mostrados a continuación para inspeccionar la normalidad de los errores a través de la distribución empírica de los residuos (que es lo que observamos).
Los mismos permiten comparar los cuantiles de los residuos con los cuantiles de una distribución normal. Si el proceso generador de datos fuera normal, los puntos estarían cerca de la recta.

qqnorm(mod1$residuals)
qqline(mod1$residuals)

Realizamos el contraste de normalidad:

shapiro.test(mod1$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  mod1$residuals
## W = 0.88987, p-value < 2.2e-16

El resultado del p-valor indica que debemos rechazar la hipótesis nula de distribución normal de los errores. Esto confirma lo sugerido por el análisis gráfico precedente.

Presentamos dos gráficos en los que se muestran los residuos contra cada uno de los regresores, los cuales se realizan con el siguiente código:

plot1 <- ggplot(data = datos, aes(educ, mod1$residuals)) +
  geom_point() + 
  geom_smooth(color = "firebrick") + 
  geom_hline(yintercept = 0) +
  theme_bw()

plot2 <- ggplot(data = datos, aes(antiguedad, mod1$residuals)) +
  geom_point() + 
  geom_smooth(color = "firebrick") + 
  geom_hline(yintercept = 0) +
  theme_bw()

grid.arrange(plot1, plot2)

El análisis gráfico indica problemas de heteroscedasticidad y de correlación entre los residuos y el nivel de los regresores.

Una manera resumida de verlo es graficando los residuos contra los valores ajustados de \(y\), que es una combinación lineal de los regresores que estamos considerando.

ggplot(data = datos, aes(mod1$fitted.values, mod1$residuals)) +
  geom_point() +
  geom_smooth(color = "firebrick", se = FALSE) +
  geom_hline(yintercept = 0) +
  theme_bw()

Se observan los mismos problemas que mencionamos anteriormente.

5.5 Aplicación en R: ¿El sexo explica los salarios?

Nos preguntamos si el sexo biológico explica el nivel de salarios de las personas. Para responder esta pregunta comenzaremos calculando una tabla donde veremos algunas medidas de resumen del salario en función del sexo biológico:

tabla_salarios <- datos %>%
  group_by(female) %>%
  summarise(cantidad = n(), Mediana = median(wage), Q1 = quantile(wage, 0.25),
            Q3 = quantile(wage, 0.75), Promedio = mean(wage), Desvío = sd(wage),
            Mínimo = min(wage), Máximo = max(wage))

print(tabla_salarios)

## # A tibble: 2 x 9
##   female cantidad Mediana    Q1    Q3 Promedio Desvío Mínimo Máximo
##    <int>    <int>   <dbl> <dbl> <dbl>    <dbl>  <dbl>  <dbl>  <dbl>
## 1      0      274    6     4.14  8.77     7.10   4.16  1.5     25.0
## 2      1      252    3.75  3     5.51     4.59   2.53  0.530   21.6

Observamos estadísticos de posición mayores para los varones que para las mujeres.

Luego asignamos etiquetas a la variable female, con la función factor(). La recodificación es en el mismo orden en el que está codificada la variable: 0 es varon y 1 es mujer.

datos$género <- factor(datos$female, labels = c("varon", "mujer"))

Observamos el diagrama de dispersión de los salarios por el género.

boxplot(wage ~ género, data = datos)

Estimamos nuevamente un modelo de regresión múltiple, pero agregando al modelo anterior, la variable female.

\[\begin{equation} wage_i = \beta _0 + \beta_1 educ_i + \beta_2 antiguedad_i + \beta_3 female_i + u_i \\ \end{equation}\]

mod2 <- lm(wage ~ educ + antiguedad + female, data = datos)
summary(mod2)

## 
## Call:
## lm(formula = wage ~ educ + antiguedad + female, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.5184 -1.8074 -0.4477  1.0270 14.1229 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.84503    0.64774  -1.305    0.193    
## educ         0.53799    0.04709  11.425  < 2e-16 ***
## antiguedad   0.16441    0.01835   8.962  < 2e-16 ***
## female      -1.78839    0.26559  -6.734 4.38e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.968 on 522 degrees of freedom
## Multiple R-squared:  0.3577, Adjusted R-squared:  0.354 
## F-statistic: 96.88 on 3 and 522 DF,  p-value: < 2.2e-16

Utilizamos la función summ() de la librería jtools, que proporciona varias opciones para dar formato a los resúmenes de regresión.

library(jtools)
summ(mod2, scale = TRUE, vifs = TRUE, part.corr = TRUE, confint = TRUE, pvals = TRUE)

Observations	526
Dependent variable	wage
Type	OLS linear regression

F(3,522)	96.88
R²	0.36
Adj. R²	0.35

	Est.	2.5%	97.5%	t val.	p	VIF	partial.r	part.r
(Intercept)	6.75	6.40	7.11	37.21	0.00	NA	NA	NA
educ	1.49	1.23	1.75	11.43	0.00	1.01	0.45	0.40
antiguedad	1.19	0.93	1.45	8.96	0.00	1.05	0.37	0.31
female	-1.79	-2.31	-1.27	-6.73	0.00	1.05	-0.28	-0.24
Standard errors: OLS; Continuous predictors are mean-centered and scaled by 1 s.d.

Si utilizamos la librería jtools en conjunto con el paquete huxtable, podemos concatenar múltiples modelos en una sola tabla con la función export_summs(). En el siguiente link (https://jtools.jacob-long.com/reference/export_summs.html) se detalla como se utiliza y a continuación presentamos el código aplicado a los datos con los que venimos trabajando.

#install.packages("huxtable")
library(huxtable)

mod0 <- lm(wage ~ educ, data = datos)
mod1 <- lm(wage ~ educ + antiguedad, data = datos)
mod2 <- lm(wage ~ educ + antiguedad + female, data = datos)

export_summs(mod0, mod1, mod2)

	Model 1	Model 2	Model 3
(Intercept)	-0.90	-2.22 ***	-0.85
	(0.68)	(0.64)	(0.65)
educ	0.54 ***	0.57 ***	0.54 ***
	(0.05)	(0.05)	(0.05)
antiguedad		0.19 ***	0.16 ***
		(0.02)	(0.02)
female			-1.79 ***
			(0.27)
N	526	526	526
R2	0.16	0.30	0.36
*** p < 0.001; ** p < 0.01; * p < 0.05.

5.6 Test de significación global y conjunta

Queremos probar si la variable experiencia es relevante en el modelo de salarios. Como antes, empezamos cargando los datos:

data("wage1") # carga la base salarios
base <- wage1  # la asigno al objeto base con el que trabajaremos

El estadístico \(F\) relaciona el incremento relativo de la suma de cuadrados de los residuos \(SRC\) al pasar del modelo no restringido (\(nr\)) al modelo restringido (\(r\)).

\(F=\frac{\left(SRC_{R}-SRC_{NR}\right)/q}{SRC_{NR}/(n-k-1)}\)

El estadístico es por definición positivo, pues \(SCR_R\) no puede ser menor que \(SRC{nr}\)

El modelo no restringido (NR) que nos estamos planteando es el siguiente: \[\begin{equation} wage_i = \beta_0 + \beta_1 educ_i + \beta_2 exper_i + \beta_3 expersq_i + \beta_4 tenure_i + \beta_5 female_i + u_i \\ \end{equation}\]

ModeloNR <- lm(wage ~ educ + exper + expersq + tenure + female, data = base)
summary(ModeloNR)

## 
## Call:
## lm(formula = wage ~ educ + exper + expersq + tenure + female, 
##     data = base)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.0152 -1.7066 -0.4399  1.1069 13.5695 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.1200921  0.7120462  -2.977  0.00304 ** 
## educ         0.5300907  0.0485881  10.910  < 2e-16 ***
## exper        0.2048418  0.0344576   5.945 5.08e-09 ***
## expersq     -0.0041266  0.0007489  -5.510 5.65e-08 ***
## tenure       0.1336694  0.0206325   6.479 2.15e-10 ***
## female      -1.7902259  0.2576917  -6.947 1.12e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.878 on 520 degrees of freedom
## Multiple R-squared:  0.3987, Adjusted R-squared:  0.3929 
## F-statistic: 68.94 on 5 and 520 DF,  p-value: < 2.2e-16

El modelo restringido (R) implica suponer que \(H_0\) es cierta. La \(H_0\) incluye \(q\) restricciones de exclusión, involucra restricciones sobre un total de \(k+1\) parámetros. En este caso \(q = 2\) pues testear si la experiencia no es relevante implica asumir que los coeficientes de regresión asociados a \(exper\) y a \(expersq\) son iguales a cero.

El modelo restringido a estimar es el siguiente:

\[\begin{equation} wage_i = \beta_0 + \beta_1 educ_i + \beta_2 tenure_i + \beta_3 female_i + u_i \\ \end{equation}\]

ModeloR <- lm(wage ~ educ + tenure + female, data = base)
summary(ModeloR)

## 
## Call:
## lm(formula = wage ~ educ + tenure + female, data = base)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.5184 -1.8074 -0.4477  1.0270 14.1229 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.84503    0.64774  -1.305    0.193    
## educ         0.53799    0.04709  11.425  < 2e-16 ***
## tenure       0.16441    0.01835   8.962  < 2e-16 ***
## female      -1.78839    0.26559  -6.734 4.38e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.968 on 522 degrees of freedom
## Multiple R-squared:  0.3577, Adjusted R-squared:  0.354 
## F-statistic: 96.88 on 3 and 522 DF,  p-value: < 2.2e-16

Comparamos la suma de los cuadrados de los residuos de ambos modelos y calculamos el estadístico F. Como habíamos visto, la función deviance() devuelve la SCR.

SCR_R <- deviance(ModeloR) 
SCR_NR <- deviance(ModeloNR)

estad.f <- ((SCR_R - SCR_NR)/2) / (SCR_NR/ModeloNR$df.residual)
          
print(estad.f)

## [1] 17.72523

Otra manera de calcular el valor del estadístico \(F\) es comparando los coeficientes de determinación \(R^2\) de los modelos restringidos (r) y sin restringir o no restringido (nr) tal y como se muestra en el código a continuación.

Es estadístico \(F\) se define como: \[F = \frac{\left(R^2_{nr}-R^2_{r}\right)/q}{1- R^2_{nr}/(n-k-1)}\]

R_nr <- summary(ModeloNR)$r.squared
R_r <- summary(ModeloR)$r.squared
estad.f.R <- ((R_nr - R_r)/2) / ((1 - R_nr)/ModeloNR$df.residual)
          
print(estad.f.R)

## [1] 17.72523

Observación: cuidado con el orden de los coeficientes! El \(R^2_{nr}\) va primero en el numerador pues siempre es mayor o igual que el \(R^2_{r}\). Si el estadístico \(F\) da negativo necesariamente hay un error y en general es porque se ha invertido este orden.

Finalmente, existe una tercera forma de realizar el contraste planteado y es solicitando a R que lo realice directamente con la función linearHypothesis(). Requiere los paquete gmodels y car.

Aquí se encuentra la documentación de la función: https://www.rdocumentation.org/packages/car/versions/3.0-13/topics/linearHypothesis

library(gmodels)
library(car)

f.conjunta <- linearHypothesis(ModeloNR,  c("exper=0", "expersq=0"))$F[2]

print(f.conjunta)

## [1] 17.72523

Observación: La notación de la función linearHypothesis() a veces puede dar lugar a confusión. Es importante que quede claro que no estamos testeando si las variables valen cero, no tendría ningún sentido realizar un contraste de hipótesis para ello.

Una vez calculado el valor del estadístico \(F\), por cualquiera de los tres métodos vistos anteriormente, solo resta decidir. Para ello, debemos comparar el valor del estadístico \(F\) con el valor crítico de una distribución \(F\) con los grados de libertad correspondientes.

Bajo \(H_0\) cierta \(F\sim F_{q,n-k-1}\)

La región de rechazo es: \(RR = \left\{\textit{muestras: }F\geq F^{1-\alpha}_{q,n-k-1} \right\}\) siendo \(\alpha\) el nivel de significación de la prueba.

Para encontrar el valor crítico de \(F\) en R, podemos utilizar la función qf(). Esta función devuelve el valor crítico de la distribución \(F\) según el nivel de significancia, los grados de libertad del numerador y los grados de libertad del denominador proporcionados. Utiliza la siguiente sintaxis: qf (p, df1, df2. lower.tail = TRUE) dónde:

• p: es el nivel de significación a utilizar
• df1: son los grados de libertad del numerador
• df2: los grados de libertad del denominador
• lower.tail: Si es verdadero (valor TRUE), se devuelve la probabilidad a la izquierda de p en la distribución \(F\). Si es FALSE, se devuelve la probabilidad a la derecha. El valor por defecto es TRUE.

Mostramos el código aplicado al ejemplo con el que venimos trabajando, que nos mostrará la probabilidad a la izquierda de p.

alfa <-  0.05
v.crit.conj <- qf(alfa, 2, ModeloNR$df.residual, lower.tail = FALSE)

Con el siguiente código, mostramos la interpretación del test realizado, haciendo que se imprima de forma automática (con la función print()), un texto u otro, según el valor crítico obtenido con la función qf() explicada.

if (abs(f.conjunta) >= v.crit.conj){print("F pertenece a la Región crítica o de rechazo. Entonces Rechazo H_0, existe evidencia estadística para afirmar que la experiencia es relevante en el modelo estimado a un 5% de significación.")} else 
  {print("No rechazo H_0, no existe evidencia empírica suficiente para rechazar la hipótesis de que la experiencia no es relevante para explicar el salario de las personas")}

## [1] "F pertenece a la Región crítica o de rechazo. Entonces Rechazo H_0, existe evidencia estadística para afirmar que la experiencia es relevante en el modelo estimado a un 5% de significación."

5.7 Contraste de restricciones lineales utilizando la matriz R

Planteamos el contraste de significación individual (un coeficiente igual a cero), global (todos los coeficientes de la regresión son cero) o conjunta. En este último caso podemos plantear que cualquier combinación lineal de los coeficientes es cero o igual a un determinado número. Los primeros casos se pueden ver como un caso particular de este último. Planteamos \(H_{0})\) de manera general como:

\[ \begin{array}{cccccc} \mathbf{R} & \mathbf{\beta} & = & \mathbf{r} \\ (q \times (k+1)) & ((k+1) \times 1)) & & (q \times 1) \end{array} \]

donde \(q\) es el número de restricciones bajo \(H_{0})\) y \((k+1)\) es el número de parámetros del modelo no restringido.

Por otro lado, \(H_{1})\) implica que al menos una de las igualdades no se cumple. Cuidado! No es lo mismo a que todas las igualdades no se cumplan.

\[\begin{equation} \mathrm{H}_{1})\,\mathrm{R}\beta\neq\mathrm{r} \end{equation}\]

5.7.1 Algunos ejemplos

Ejemplo 1: Supongamos un modelo con 4 regresores con los que queremos realizar la siguiente prueba de hipótesis:

\(H_{0}) \beta_{1}=0\) versus \(H_{1}) \beta_{1}\neq 0\)

\[\begin{equation} \mathrm{H}_{0})\,\mathrm{R}\beta=\mathrm{r} \end{equation}\]

\[ \mathrm{H}_{1})\,\mathrm{R}\beta\neq\mathrm{r} \]

\(\beta' = \left( \begin{array}{ccccc} \beta_0 & \beta_1 & \beta_{2} & \beta_{3} & \beta_{4} \end{array} \right)\)

\[ \begin{array}{ccccc} \mathbf{R} & \mathbf{\beta} &=& \mathbf{r} \\ (1 \times 5)&(5 \times 1) & & (1 \times 1) \end{array} \]

Entonces \[\mathbf{R} = \left( \begin{array}{ccccc} 0 & 1 & 0 & 0 & 0 \end{array} \right) \]

Ejemplo 2: Supongamos ahora que queremos testear si dos coeficientes son iguales a cero.

\(H_{0})\) \(\beta_{3}=\beta_{4}=0\)

\(H_{1})\) algún \(\beta_{j}\neq 0\) con \(j=3,4\)

En este caso \(q=2\), es decir que tenemos 2 restricciones.

\[ \begin{array}{ccccc} \mathbf{R} & \mathbf{\beta} &=& \mathbf{r} \\ (2 \times 5)&(5 \times 1) & & (2 \times 1) \end{array} \]

\(\mathbf{\beta} = \left( \begin{array}{c} \beta_0 \\ \beta_1 \\ \beta_2 \\ \beta_3 \\ \beta_4 \end{array} \right)\)

Entonces \[R = \left( \begin{array}{ccccc} 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{array} \right) \]

Ahora, una vez que comprendimos como escribir las hipótesis de manera matricial. Utilizaremos la función glh.test() de \(R\) para realizar los contrastes. La misma requiere especificar como intputs:

• El nombre del modelo que estoy considerando.
  
• \(q\): el parametro que indica la cantidad de restricciones (numero de filas de la matriz \(R\)).
• \(r\): el vector de los valores que adquieren las restricciones bajo \(H_0\) cierta.
• \(R\): la matriz de las restricciones lineales.

Para escribir la matriz \(R\), se indican los valores de cada fila de la matriz utilizando la función c() y luego se “pegan” con la función rbind(). Cada fila contiene una combinación lineal de los coeficientes de regresión.

Mostramos el código para escribir la matriz \(R\) del ejemplo anterior (ejemplo 2).

R <- rbind( c(0,0,0,1,0), c(0,0,0,0,1) )
R

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    0    0    0    1    0
## [2,]    0    0    0    0    1

Una manera alternativa de escribir dicha matriz, es utilizar la función matrix(), especificando todos sus componentes de corrido dentro de un mismo vector, el número de filas que tendrá (con el comando nrow) y que vamos a completar la matriz por filas (con el comando byrow = T).

R <- matrix(c(0,0,0,1,0,0,0,0,0,1), nrow = 2, byrow = T)
R

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    0    0    0    1    0
## [2,]    0    0    0    0    1

5.7.2 Aplicación en R: estimación del modelo general para los salarios

Estimamos el modelo general para los salarios:

\[\begin{equation} wage_i = \beta_0 + \beta_1 female_i + \beta_2 educ_i + \beta_3 exper_i + \beta_4 expersq_i + u_i \\ \end{equation}\]

Modelo <- lm(wage ~ 1 + female + educ + exper + expersq, data = base)
summary(Modelo)

## 
## Call:
## lm(formula = wage ~ 1 + female + educ + exper + expersq, data = base)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.8380 -1.8998 -0.3573  1.2691 14.2932 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.3192037  0.7388254  -3.139  0.00179 ** 
## female      -2.1140347  0.2625501  -8.052 5.57e-15 ***
## educ         0.5562848  0.0502875  11.062  < 2e-16 ***
## exper        0.2551276  0.0348671   7.317 9.64e-13 ***
## expersq     -0.0044396  0.0007762  -5.720 1.80e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.989 on 521 degrees of freedom
## Multiple R-squared:  0.3501, Adjusted R-squared:  0.3451 
## F-statistic: 70.17 on 4 and 521 DF,  p-value: < 2.2e-16

Utilizamos la función glh() para realizar el contraste. La misma devuelve un objeto tipo lista con los componentes que se devuelven con la función names().

H0 <- glh.test(Modelo, R, 0)
names(H0)

## [1] "call"       "statistic"  "parameter"  "p.value"    "estimate"  
## [6] "null.value" "method"     "data.name"  "matrix"

También podemos acceder a cada uno de los componentes con el símbolo $, como hemos hecho anteriormente en otros ejemplos.

H0$matrix

##      (Intercept) female educ exper expersq
## [1,]           0      0    0     1       0
## [2,]           0      0    0     0       1

H0$p.value

##              Estimate
## Estimate 1.332268e-15

H0$estimate

##          Estimate
## [1,]  0.255127617
## [2,] -0.004439636

Igual que con cualquier contraste decidimos con el p-valor y con el código a continuación vemos que el resultado coincide para la prueba de significación individual.

alpha = 0.05
if (alpha >= H0$p.value){print("F pertenece a la Región crítica o de rechazo. Entonces Rechazo H_0, existe evidencia estadística para afirmar que la variable experiencia es conjuntamente significativa (en nivel y al cuadrado) en el modelo estimado a un 5% de significación.")} else {print("No rechazo H_0, no existe evidencia empírica suficiente para rechazar la hipótesis de  que la experiencia no es relevante para explicar el salario de las personas")}

## [1] "F pertenece a la Región crítica o de rechazo. Entonces Rechazo H_0, existe evidencia estadística para afirmar que la variable experiencia es conjuntamente significativa (en nivel y al cuadrado) en el modelo estimado a un 5% de significación."

Escribimos la matriz \(R\) indicando los componentes, indicando dentro del vector, el \(1\) en la posición del coeficiente de regresión que queremos testear si es significativo. nrow = q especifica que quiero una matriz de \(q\) filas y con byrow = T indicamos que quiero que la matriz se complete por filas.

q <- 1

R1 <- matrix (c(0,0,0,1,0), nrow = q, byrow = T)
R1

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    0    0    0    1    0

H0 <- glh.test(Modelo, R1, 0)

Miramos los componentes y los comparamos con los obtenidos con summary() para verificar que es equivalente a la prueba de significación individual para testear la significación del parámetro de regresión asociado a la experiencia. Recordar que en la significación individual el cuadrado del estadístico \(t\) coincide con el estadístico \(F\).

H0$statistic

##       F 
## 53.5406

H0$p.value

##              Estimate
## Estimate 9.640067e-13

H0$estimate

##       Estimate
## [1,] 0.2551276

summary(Modelo)

## 
## Call:
## lm(formula = wage ~ 1 + female + educ + exper + expersq, data = base)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.8380 -1.8998 -0.3573  1.2691 14.2932 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.3192037  0.7388254  -3.139  0.00179 ** 
## female      -2.1140347  0.2625501  -8.052 5.57e-15 ***
## educ         0.5562848  0.0502875  11.062  < 2e-16 ***
## exper        0.2551276  0.0348671   7.317 9.64e-13 ***
## expersq     -0.0044396  0.0007762  -5.720 1.80e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.989 on 521 degrees of freedom
## Multiple R-squared:  0.3501, Adjusted R-squared:  0.3451 
## F-statistic: 70.17 on 4 and 521 DF,  p-value: < 2.2e-16

Finalmente, chequearemos con el código que sigue, que la prueba de significación global también coincide:

q <- 4

R <- rbind(c(0,1,0,0,0), c(0,0,1,0,0), c(0,0,0,1,0), c(0,0,0,0,1))
R

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    0    1    0    0    0
## [2,]    0    0    1    0    0
## [3,]    0    0    0    1    0
## [4,]    0    0    0    0    1

H0 <- glh.test(Modelo, R, 0)

H0$statistic

##        F 
## 70.16953

H0$p.value

##          Estimate
## Estimate        0

H0$estimate

##          Estimate
## [1,] -2.114034720
## [2,]  0.556284789
## [3,]  0.255127617
## [4,] -0.004439636

summary(Modelo)

## 
## Call:
## lm(formula = wage ~ 1 + female + educ + exper + expersq, data = base)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.8380 -1.8998 -0.3573  1.2691 14.2932 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.3192037  0.7388254  -3.139  0.00179 ** 
## female      -2.1140347  0.2625501  -8.052 5.57e-15 ***
## educ         0.5562848  0.0502875  11.062  < 2e-16 ***
## exper        0.2551276  0.0348671   7.317 9.64e-13 ***
## expersq     -0.0044396  0.0007762  -5.720 1.80e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.989 on 521 degrees of freedom
## Multiple R-squared:  0.3501, Adjusted R-squared:  0.3451 
## F-statistic: 70.17 on 4 and 521 DF,  p-value: < 2.2e-16

6 Variables ficticias o dummies

6.1 Introducción

En esta sección continuaremos trabajando con la base de salarios del texto de Wooldridge. El objeto de este capítulo es aprender a crear variables dummies y a utilizarlas en un MRLM.

Las variables que utilizaremos son las siguientes:

• wage: salario promedio por hora
• educ: años de educación
• exper: años de experiencia potencial
• ternure: años con el empleador actual (antigüedad)
• nonwhite: es igual a 1 si la persona no es blanca
• female: es igual a 1 si la persona es mujer
• married: es igual a 1 si la persona es casada
• numdep: =1 número de dependientes que tiene la persona
• lwage: logaritmo del salario log(salario)

6.2 Modelo para el salario

Continuaremos trabajando con la base de salarios de los temas anteriores (MRLS y MRLM)

library(wooldridge)

data("wage1")
names(wage1)

##  [1] "wage"     "educ"     "exper"    "tenure"   "nonwhite" "female"  
##  [7] "married"  "numdep"   "smsa"     "northcen" "south"    "west"    
## [13] "construc" "ndurman"  "trcommpu" "trade"    "services" "profserv"
## [19] "profocc"  "clerocc"  "servocc"  "lwage"    "expersq"  "tenursq"

\[\begin{equation} los(wage_i) = \beta _0+\beta _1 educ_i +\beta _2 exper_i+ \beta _3 antiguedad_i \beta _4 nro-dependiendientes_i+ +u_i \\ \end{equation}\]

wageModel <- lm(lwage ~ educ + exper + tenure + numdep , data = wage1)
summary(wageModel)

## 
## Call:
## lm(formula = lwage ~ educ + exper + tenure + numdep, data = wage1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.04536 -0.29037 -0.03488  0.29086  1.42843 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.257392   0.112794   2.282   0.0229 *  
## educ        0.093192   0.007566  12.318  < 2e-16 ***
## exper       0.004259   0.001738   2.450   0.0146 *  
## tenure      0.022009   0.003097   7.107 3.92e-12 ***
## numdep      0.009872   0.015763   0.626   0.5314    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4411 on 521 degrees of freedom
## Multiple R-squared:  0.3165, Adjusted R-squared:  0.3113 
## F-statistic: 60.32 on 4 and 521 DF,  p-value: < 2.2e-16

anova(wageModel)

Df	Sum Sq	Mean Sq	F value	Pr(>F)
1	27.6	27.6	142	4.67e-29
1	9.42	9.42	48.4	1.03e-11
1	9.89	9.89	50.8	3.39e-12
1	0.0763	0.0763	0.392	0.531
521	101	0.195

6.3 Crear variables dummies

Para crear dummies se puede utilizar el paquete fastDummies. En este caso crearemos una dummie para cada valor de la variable numdep; se crea una dummie por cada categoría, por lo tanto terminan siendo 7.

#install.packages("fastDummies")
library(fastDummies)

datos1 <- dummy_cols(wage1, select_columns = "numdep")
names(datos1)

##  [1] "wage"     "educ"     "exper"    "tenure"   "nonwhite" "female"  
##  [7] "married"  "numdep"   "smsa"     "northcen" "south"    "west"    
## [13] "construc" "ndurman"  "trcommpu" "trade"    "services" "profserv"
## [19] "profocc"  "clerocc"  "servocc"  "lwage"    "expersq"  "tenursq" 
## [25] "numdep_0" "numdep_1" "numdep_2" "numdep_3" "numdep_4" "numdep_5"
## [31] "numdep_6"

Para crear las dummies para mas de una variable a la vez, lo hacemos con la función c(), listando el nombre de las variables separándolos por comas. Por ejemplo:

dummy_cols(datos, select_columns = c("numdep", "exper"))

6.4 Crear interacciones entre los regresores

Además se pueden realizar interacciones entre variables creando diferentes dummies para cada una de las categorías obtenidas al cruzar las variables. Como ejemplo, mostraremos como crear las interacciones entre las variables que indican el estado civil y el género.

Para ello es necesario crear antes la variable que indica que la persona es varon y la variable que indica que la persona es soltera, que no están en la base como tal. Utilizaremos la librería car, específicamente la función recode():

# install.packages("car")
library(car)

datos1$male <- recode(datos1$female, "1 = 0; 0 = 1")
datos1$single <- recode(datos1$married, "1=0; 0=1")

names(datos1)

##  [1] "wage"     "educ"     "exper"    "tenure"   "nonwhite" "female"  
##  [7] "married"  "numdep"   "smsa"     "northcen" "south"    "west"    
## [13] "construc" "ndurman"  "trcommpu" "trade"    "services" "profserv"
## [19] "profocc"  "clerocc"  "servocc"  "lwage"    "expersq"  "tenursq" 
## [25] "numdep_0" "numdep_1" "numdep_2" "numdep_3" "numdep_4" "numdep_5"
## [31] "numdep_6" "male"     "single"

Crearemos entonces las variables dummies:

• mujer casada: es igual a 1 si la persona es mujer casada.
• varon casado: es igual a 1 si la persona es varon casado.
• mujer soltera: es igual a 1 si la persona es mujer soltera.
• varon soltero: es igual a 1 si la persona es varon soltero.

attach(datos1)

marrmale <- married*male
marrfem <- married*female
singfem <- single*female

datos1 <- cbind(datos1, marrmale, marrfem, singfem)
names(datos1)

##  [1] "wage"     "educ"     "exper"    "tenure"   "nonwhite" "female"  
##  [7] "married"  "numdep"   "smsa"     "northcen" "south"    "west"    
## [13] "construc" "ndurman"  "trcommpu" "trade"    "services" "profserv"
## [19] "profocc"  "clerocc"  "servocc"  "lwage"    "expersq"  "tenursq" 
## [25] "numdep_0" "numdep_1" "numdep_2" "numdep_3" "numdep_4" "numdep_5"
## [31] "numdep_6" "male"     "single"   "marrmale" "marrfem"  "singfem"

Observación: Recordemos que con la función attach() se puede acceder a los objetos del data frame que indicamos como argumento, solamente con el nombre de los mismos. Luego, con el signo \(*\) logramos las interacciones antes explicadas. Finalmente, la función cbind() combina por columnas varios vectores, matrices o data frames.

6.5 Estimamos modelos para el salario

El modelo mas básico es aquel en el que se estima el logaritmo del salario en función de la educación de las personas. Se define como: \[\begin{equation} log(salario_i) = \beta_0 + \beta_1 Educ_i + u_i \\ \end{equation}\]

En R, se estima con el siguiente codigo:

modelo1 <- lm(lwage ~ educ, data = datos1)
summary(modelo1)

## 
## Call:
## lm(formula = lwage ~ educ, data = datos1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.21158 -0.36393 -0.07263  0.29712  1.52339 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.583773   0.097336   5.998 3.74e-09 ***
## educ        0.082744   0.007567  10.935  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4801 on 524 degrees of freedom
## Multiple R-squared:  0.1858, Adjusted R-squared:  0.1843 
## F-statistic: 119.6 on 1 and 524 DF,  p-value: < 2.2e-16

Nos preguntamos si el efecto de la educación es el mismo para varones y para mujeres. Para ello debemos comparar dos modelos considerando dos bases separadas por sexo (una para mujeres y otra para varones).

El modelo para las mujeres se define como: \[\begin{equation} log(salario_mi) = \beta_0 + \beta_1 Educ_mi + u_mi \\ \end{equation}\]

Análogamente, el modelo para los varones se define como:

\[\begin{equation} log(salario_hi) = \beta _0+\beta _1 Educ_hi + u_hi \\ \end{equation}\]

En primer lugar creamos las dos nuevas bases (mujeres y varones):

mujeres <- datos1 %>% 
  select("lwage", "educ", "exper", "tenure") %>% 
  filter(female == 1)

varones <- datos1 %>% 
  select("lwage", "educ", "exper", "tenure") %>% 
  filter(female == 0)

dim(mujeres)

## [1] 252   4

dim(varones)

## [1] 274   4

A continuación corremos el mismo modelo pero por un lado con la base mujeres y por otro con la base varones:

sal_mujeres <- lm(lwage ~ educ, data = mujeres)
sal_varones <- lm(lwage ~ educ, data = varones)

summary(sal_mujeres)

## 
## Call:
## lm(formula = lwage ~ educ, data = mujeres)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.02673 -0.24397 -0.06163  0.21415  1.21924 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.46589    0.12890   3.614 0.000364 ***
## educ         0.07716    0.01026   7.520 9.82e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.402 on 250 degrees of freedom
## Multiple R-squared:  0.1845, Adjusted R-squared:  0.1812 
## F-statistic: 56.55 on 1 and 250 DF,  p-value: 9.824e-13

summary(sal_varones)

## 
## Call:
## lm(formula = lwage ~ educ, data = varones)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.11585 -0.34240 -0.01708  0.32659  1.34740 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.825955   0.127813   6.462 4.75e-10 ***
## educ        0.077228   0.009731   7.936 5.47e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4828 on 272 degrees of freedom
## Multiple R-squared:  0.188,  Adjusted R-squared:  0.185 
## F-statistic: 62.99 on 1 and 272 DF,  p-value: 5.471e-14

library(jtools)
library(huxtable)

export_summs(modelo1, sal_mujeres, sal_varones, 
              modelo.nombres = c("General", "Mujeres", "Varones"))

	Model 1	Model 2	Model 3
(Intercept)	0.58 ***	0.47 ***	0.83 ***
	(0.10)	(0.13)	(0.13)
educ	0.08 ***	0.08 ***	0.08 ***
	(0.01)	(0.01)	(0.01)
N	526	252	274
R2	0.19	0.18	0.19
*** p < 0.001; ** p < 0.01; * p < 0.05.

Observando los resultados, el retorno de la educación parece ser el mismo para ambos géneros. Eso no significa que el salario sea el mismo.

Etiquetamos las categorías para los gráficos y luego graficamos los diagramas de dispersión del salario con respecto al sexo biológico para ver si existen diferencias:

datos1$Sexo_biologico <- factor(datos1$female, labels = c("varón", "mujer"))

ggplot(datos1, aes(x = educ, y = lwage, color = Sexo_biologico)) + 
  geom_point() +
  geom_smooth(method = 'lm', formula = y ~ x, se = FALSE, col = 'dodgerblue1') +
  theme_light() +
  facet_wrap(~ Sexo_biologico)+
  ggtitle("Relación entre salario y educación según el sexo biológico")

Otra manera, estima las rectas de regresión para cada submuestra y las muestra en un mismo gráfico:

#Plot interaction model
ggplot(datos1, aes(x = educ, y = lwage, col = Sexo_biologico)) +
  geom_point(alpha = 0.05) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(y = "lwage", x = "educ", title = "Salario en función de la educación")

Finalmente, para contrastar si el retorno de la educación es el mismo según el sexo biológico estimamos el siguiente modelo: \[\begin{equation} log(salario_i) = \beta_0 + \beta_1 Mujer_i + u_i \\ \end{equation}\]

model_sexo <- lm(lwage ~ female, data = datos1)
summary(model_sexo)

## 
## Call:
## lm(formula = lwage ~ female, data = datos1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.05123 -0.31774 -0.04889  0.35548  1.65773 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.81357    0.02981  60.830   <2e-16 ***
## female      -0.39722    0.04307  -9.222   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4935 on 524 degrees of freedom
## Multiple R-squared:  0.1396, Adjusted R-squared:  0.138 
## F-statistic: 85.04 on 1 and 524 DF,  p-value: < 2.2e-16

boxplot(wage ~ Sexo_biologico, data = datos1)

tabla_salarios <- datos1 %>%
  group_by(female) %>%
  summarise(Cantidad = n(), Promedio = mean(lwage))

print(tabla_salarios)

## # A tibble: 2 x 3
##   female Cantidad Promedio
##    <int>    <int>    <dbl>
## 1      0      274     1.81
## 2      1      252     1.42

La diferencia en los valores promedio es igual a \(1.416353 - 1.813570 = -0.397217\), que coincide con la estimación del coeficiente de regresión simple del logaritmo del salario contra el sexo biológico (variable female*).

6.6 Modelos para el salario que consideran interacciones entre regresores

Estimaremos ahora un modelo con interacciones para dar la posibilidad de que el retorno de la educación sea diferente para varones que para mujeres. El modelo es el siguiente:

\[\begin{equation} log(salario_i) = \beta_0 + \beta_1 Female_i + \beta_2 Educ_i + \beta_3 Female_i*Educ_i + u_i \\ \end{equation}\]

modelo4 <- lm(lwage ~ female + educ + female*educ, data = datos1)
summary(modelo4)

## 
## Call:
## lm(formula = lwage ~ female + educ + female * educ, data = datos1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.02673 -0.27468 -0.03721  0.26221  1.34740 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.260e-01  1.181e-01   6.997 8.08e-12 ***
## female      -3.601e-01  1.854e-01  -1.942   0.0527 .  
## educ         7.723e-02  8.988e-03   8.593  < 2e-16 ***
## female:educ -6.408e-05  1.450e-02  -0.004   0.9965    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4459 on 522 degrees of freedom
## Multiple R-squared:  0.3002, Adjusted R-squared:  0.2962 
## F-statistic: 74.65 on 3 and 522 DF,  p-value: < 2.2e-16

Se observa que la interacción no es significativa, esto es equivalente al resultado visto antes que mostraba que el rendimiento de la educación era el mismo para mujeres que para varones.

Otra alternativa: estimamos un modelo para el logaritmo del salario, en función del sexo biológico, la educación, la experiencia, la antigüedad y sus efectos cuadráticos.

\[\begin{equation} log(salario_i) = \beta_0 + \beta_1 Female_i + \beta_2 Educ_i + \beta_3 exper_i + \beta_4 expersq_i + \beta_5 tenure_i + \beta_6 tenuresq_i + \beta_7 Female_i*Educ_i + u_i \\ \end{equation}\]

modelo5 <- lm(lwage ~ female + educ + exper + expersq + tenure + tenursq + educ*female, data = datos1)
summary(modelo5)

## 
## Call:
## lm(formula = lwage ~ female + educ + exper + expersq + tenure + 
##     tenursq + educ * female, data = datos1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.83265 -0.25261 -0.02374  0.25396  1.13584 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.3888060  0.1186871   3.276  0.00112 ** 
## female      -0.2267886  0.1675394  -1.354  0.17644    
## educ         0.0823692  0.0084699   9.725  < 2e-16 ***
## exper        0.0293366  0.0049842   5.886 7.11e-09 ***
## expersq     -0.0005804  0.0001075  -5.398 1.03e-07 ***
## tenure       0.0318967  0.0068640   4.647 4.28e-06 ***
## tenursq     -0.0005900  0.0002352  -2.509  0.01242 *  
## female:educ -0.0055645  0.0130618  -0.426  0.67028    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4001 on 518 degrees of freedom
## Multiple R-squared:  0.441,  Adjusted R-squared:  0.4334 
## F-statistic: 58.37 on 7 and 518 DF,  p-value: < 2.2e-16

¿Cuál es el rendimiento estimado de la educación para los varones? ¿Y para las mujeres?

El rendimiento estimado de la educación para los varones es \(\widehat{\beta}_2=0.082\), es decir, \(8.2\%\). Para las mujeres es (\(\widehat{\beta}_2-\widehat{\beta}_7=0.082-0.0056=0.0764\)), aproximadamente \(7.6\%\).

La diferencia en el rendimiento de la educación para ambos sexos se estima en poco mas de medio punto porcentual (\(-0.56\%\)). Como vimos, esta diferencia no es grande desde el punto de vista económico ni tampoco es estadísticamente significativa (observar el p-valor). En consecuencia, no se ha encontrado suficiente evidencia empírica para rechazar la hipótesis de que el rendimiento en la educación sea el mismo para varones y para mujeres.

Estimamos un modelo para un último modelo para el salario incluyendo las interacciones entre el sexo biológico y el estado civil.

\[\begin{equation} log(salario_i) = \beta_0 + \beta_1 marrmale_i + \beta_2 marrfem_i + \beta_3 singfem_i + \beta_4 educ_i + \beta_5 exper_i + \beta_6 expersq_i + \beta_7 tenure_i + \beta_8 tenuresq_i + u_i \\ \end{equation}\]

modelo6 <- lm(lwage ~ marrmale + marrfem + singfem + educ + exper + expersq + tenure + tenursq, data = datos1)
summary(modelo6)

## 
## Call:
## lm(formula = lwage ~ marrmale + marrfem + singfem + educ + exper + 
##     expersq + tenure + tenursq, data = datos1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.89697 -0.24060 -0.02689  0.23144  1.09197 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.3213781  0.1000090   3.213 0.001393 ** 
## marrmale     0.2126757  0.0553572   3.842 0.000137 ***
## marrfem     -0.1982676  0.0578355  -3.428 0.000656 ***
## singfem     -0.1103502  0.0557421  -1.980 0.048272 *  
## educ         0.0789103  0.0066945  11.787  < 2e-16 ***
## exper        0.0268006  0.0052428   5.112 4.50e-07 ***
## expersq     -0.0005352  0.0001104  -4.847 1.66e-06 ***
## tenure       0.0290875  0.0067620   4.302 2.03e-05 ***
## tenursq     -0.0005331  0.0002312  -2.306 0.021531 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3933 on 517 degrees of freedom
## Multiple R-squared:  0.4609, Adjusted R-squared:  0.4525 
## F-statistic: 55.25 on 8 and 517 DF,  p-value: < 2.2e-16

La categoría de referencia en este caso son los varones solteros, por lo tanto la estimación de los coeficientes de las dummies debe interpretarse como diferencial respecto a esta categoría.

Manteniendo constantes los niveles de educación, experiencia y antigüedad (tenure), se estima que los varones casados, en promedio, ganan aproximadamente \(21,3\%\) más que los solteros. Se estima que, a igualdad del resto de las características, una mujer casada, gana en promedio \(19,8\%\) menos que un varón soltero. Se estima que, a igualdad del resto de las características, una mujer soltera, gana en promedio \(11\%\) menos que un varón soltero.

7 Test o contraste de Chow

7.1 Introducción

Se conoce como contraste de Chow a la prueba que evalúa si existen diferencias en los parámetros entre subgrupos definidos en la muestra (por ejemplo: varones y mujeres, países, etc.). En principio se supone que la muestra se divide en dos, aunque se puede generalizar a más grupos como en el caso de los países. Es especialmente utilizado en series temporales, donde se analiza si a partir de cierto período existió un cambio en los parámetros .

La prueba se puede realizar de dos formas:

1. Dividiendo el espacio en dos muestras y haciendo regresiones separadas.

2. Generando variables ficticias e incluyendo interacciones.

Con el código a continuación, nuevamente obtengo el conjunto de datos salarios y llamo al paquete stargazer:

#install.packages("wooldridge")
library(wooldridge)
data("wage1")

#install.packages("stargazer")
library(stargazer)

attach(wage1)

7.2 Implementación del test con la forma 1

Para implementar el contraste estimaremos tres modelos: el modelo para la muestra completa (\(NR\)), y los modelos para las dos muestras entre las cuales se desea testear el quiebre estructural

En este caso se estima el modelo para varones y para mujeres, utilizando la indexación de nuestra base de datos, wage1[female == 1,]. Los resultados de estos modelos se muestran en una tabla de regresión utilizando la librería stargazer.

modelNR <- lm(lwage ~ 1 + educ + exper + tenure)
modeln1 <- lm(lwage ~ 1 + educ + exper + tenure, data = wage1[female == 1,])
modeln2 <- lm(lwage ~ 1 + educ + exper + tenure, data = wage1[female == 0,] )

Imprimimos los resultados:

stargazer(modelNR, modeln1, modeln2, type = "text", title = "Results", align = TRUE)

## 
## Results
## ===========================================================================================
##                                               Dependent variable:                          
##                     -----------------------------------------------------------------------
##                                                      lwage                                 
##                               (1)                     (2)                     (3)          
## -------------------------------------------------------------------------------------------
## educ                       0.092***                0.080***                0.096***        
##                             (0.007)                 (0.010)                 (0.009)        
##                                                                                            
## exper                       0.004**                  0.002                 0.008***        
##                             (0.002)                 (0.002)                 (0.002)        
##                                                                                            
## tenure                     0.022***                 0.010*                 0.018***        
##                             (0.003)                 (0.005)                 (0.004)        
##                                                                                            
## Constant                   0.284***                 0.356**                 0.322**        
##                             (0.104)                 (0.141)                 (0.139)        
##                                                                                            
## -------------------------------------------------------------------------------------------
## Observations                  526                     252                     274          
## R2                           0.316                   0.212                   0.365         
## Adjusted R2                  0.312                   0.202                   0.358         
## Residual Std. Error    0.441 (df = 522)        0.397 (df = 248)        0.428 (df = 270)    
## F Statistic         80.391*** (df = 3; 522) 22.233*** (df = 3; 248) 51.836*** (df = 3; 270)
## ===========================================================================================
## Note:                                                           *p<0.1; **p<0.05; ***p<0.01

A partir de las suma del cuadrado de los residuos obtenidos en los modelos anteriores (\(SCR\)) y la obtención de los grados de libertad del test, se puede obtener el estadístico calculado asociado al mismo, para finalmente compararlo con el valor del estadístico de tabla a un determinado nivel de significación.

El estadístico \(F\) se define como:

\(F = \frac{\left(SCR_{r} - SCR_{nr}\right)/q}{SCR_{nr}/(n-k-1)}\)

La \(SCR\) del modelo no restringido (\(SCR_{nr}\)) puede obtenerse de dos regresiones separadas: una para cada grupo (\(SCR_1\) y \(SCR_2\)), obtenidas con la cantidad de observaciones del grupo 1 (\(n_1\)) y del grupo 2 (\(n_2\)) respectivamente. Esto es: \(SCR_{nr}=SCR_1+SCR_2\)

Entonces, tenemos que el estadístico de contraste es:

\(F=\frac{\left(SCR_{r}-(SRC_{1}+SRC_{2})\right)/(k+1)}{(SRC_{1}+SRC_{2})/(n-2(k+1))} = \frac{\left(SCR_{r}-(SRC_{1}+SRC_{2})\right)}{(SRC_{1}+SRC_{2})}*\frac{n-2(k+1))}{k+1}\)

A este estadístico se le conoce como estadístico de Chow.

Obtenemos las \(SCR\) para los tres modelos con el código a continuación:

SRCnr <- sum(residuals(modelNR)^2)
SRCn1 <- sum(residuals(modeln1)^2)
SRCn2 <- sum(residuals(modeln2)^2)

print(SRCnr)

## [1] 101.4556

print(SRCn1)

## [1] 39.0353

print(SRCn2)

## [1] 49.54719

Bajo \(H_0\) cierta \(F\sim F_{(k+1),(n-2k-2)}\).

El número de restricciones es \(q=(k+1)\), pues probar \(H_0\) implica testear que los parámetros asociados a las interacciones es igual a cero (hay \(k\)) y la diferencia en los interceptos (hay 1 más).

Los grados de libertad del modelo no restringido son \(n\) menos el número de parámetros; \(2k\) asociados a los regresores y \(2\) más de los interceptos de los 2 grupos (\(n-2(k+1) = n-2k-2\)).

En R obtenemos los grados de libertad de la siguiente manera:

q <- length(coefficients(modeln1)) #podría ser cualquiera de los tres modelos en realidad, en este caso q=4 pues son 4 coeficientes de regresión que tienen los modelos.
q

## [1] 4

gl <- (nrow(base) - 2*q)  #q = (k+1) parámetros
gl

## [1] 518

Luego calculamos el estadístico \(F\) en base a los tres modelos definidos anteriormente:

F <- ((SRCnr - (SRCn1 + SRCn2))*gl) / ((SRCn1 + SRCn2)*q)
F

## [1] 18.81934

Fcrit <- qf(0.05, q, gl, lower.tail = FALSE)
Fcrit

## [1] 2.389145

7.3 Implementación del test con la forma 2

Una versión alternativa del contraste utiliza el mismo modelo restringido, pero escribe el modelo no restringido incorporando la variable que supongo discrimina y sus interacciones con cada uno de los regresores. En este caso el sexo biológico (female) y las interacciones con la educación, con la experiencia y con la antigüedad.

El modelo no restringido a estimar será el siguiente: \(y_i= \beta_0 + \beta^*_0 fem_i +\beta_1 educ_{i} +\beta^*_1 educ_{i} fem_i+\beta_2 exper_{i} +\beta^*_2 exper_{i} fem_i+\beta_3 tenure_{i} +\beta^*_3 tenure_{ki} fem_i+u_i\)

En este caso la hipótesis a contrastar es: \(H_0)\ \beta^*_0 =\beta^*_1 = \beta^*_2=\beta^*_3=0\)

Para realizar el test se puede utilizar cualquiera de las versiones del contraste \(F\), por ejemplo: \(F=\frac{\left(R^2_{nr}-R^2_{r}\right)/(k+1)}{(1- R^2_{nr})/(n-2k-2)}\)

Bajo \(H_0\) cierta y suponiendo normalidad de los residuos \(F\sim F_{(k+1),(n-2k-2)}\).

ModeloR <- lm(lwage ~ 1 + educ + exper + tenure)
ModeloNR <- lm(lwage ~ 1 + educ + educ*female + exper + exper*female + tenure + tenure*female)

q <- 4 #nro de restricciones (4 parámetros)
gl <- 526 - 6 - 2 #grados libertad en el denominador (n-2k-2)

R_nr <- summary(ModeloNR)$r.squared
R_r <- summary(ModeloR)$r.squared
estad.f.R <- ((R_nr - R_r)/q)/((1 - R_nr)/(gl))

Calculamos el valor crítico del estadístico \(F\):

Fcrit <- qf(0.05, q, gl, lower.tail = FALSE)
print(estad.f.R)

## [1] 18.81934

print(Fcrit)

## [1] 2.389145

Tomamos la decisión:

if (abs(F) >= Fcrit){print("F pertenece a la región crítica o de rechazo. Entonces rechazo H_0, existe evidencia estadística para afirmar que es necesario trabajar con modelos diferentes para varones y mujeres para explicar el salario")} else {
            print("No rechazo H_0, no existe evidencia empírica suficiente para especificar modelos diferentes para varones y para mujeres")}

## [1] "F pertenece a la región crítica o de rechazo. Entonces rechazo H_0, existe evidencia estadística para afirmar que es necesario trabajar con modelos diferentes para varones y mujeres para explicar el salario"

8 Especificación

8.1 Introducción

En este capítulo tratamos lo principales problemas a considerar en la especificación de la parte sistemática del modelo. En particular haremos énfasis en la elección de las variables independientes a incluir y en la elección de la forma funcional.

Los problemas de especificación a considerar son:

• Omisión de variables relevantes.
• Inclusión de variables irrelevantes.
• Adopción de una forma funcional incorrecta.
• Errores de medida en las variables.

La elección de una forma funcional incorrecta, en general, produce sesgo en los estimadores. Es usual estimar inicialmente un modelo lineal en los parámetros. Aunque este puede:

• ser un modelo no lineal que fue linealizado (transformación logarítmica por ejemplo), o
• incluir no linealidades en las variables como potencias e interacciones entre variables.

En la práctica son transformaciones que entran linealmente como una nueva variable y el modelo continúa siendo lineal en los parámetros.

En algunos casos el modelo teórico indica que el modelo a ser estimado es no lineal en los parámetros y no linealizable y en ese caso será necesario trabajar con métodos diferentes que veremos mas adelante.

8.2 Evaluación de problemas de especificación.

En las siguientes secciones veremos algunos métodos para detectar problemas.

8.2.1 Análisis de los residuos del modelo

Como siempre, partimos el análisis en R de este capitulo, llamando al conjunto de datos de salarios con el que venimos trabajando:

library(wooldridge)
data("wage1")

Estimamos el modelo de regresión múltiple:

modelo <- lm(wage ~ educ + exper + tenure)

Analizamos los residuos del modelo estimado anteriormente:

residuos <- resid(modelo)
hist(residuos, freq = FALSE, main = "Distribución de los residuos", ylab = "Densidad")
lines(density(residuos), lwd = 2, col = 'red')

Creamos la función histDenNorm() para superponer la distribución empírica de los residuos con la teórica normal:

histDenNorm <- function (x, ...) {
   hist(x, ...) # Histograma
   lines(density(x), col = "blue", lwd = 2) # Densidad
   x2 <- seq(min(x), max(x), length = 40)
   f <- dnorm(x2, mean(x), sd(x))
   lines(x2, f, col = "red", lwd = 2) # Normal
   legend("topright", c("Histograma", "Densidad", "Normal"), box.lty = 0,
          lty = 1, col = c("black", "blue", "red"), lwd = c(1, 2, 2))
}

Agregamos la curva de distribución normal para comparar visualmente utilizando la función creada anteriormente:

histDenNorm(residuos, prob = TRUE, main = "Distribución de los residuos")

La no normalidad de los residuos puede revelar problemas, y puede evaluarse con el histograma de residuos o contrastarse con la prueba de Jarque-Bera (1987).

El contraste de Jarque-Bera es un contraste de validez asintótica basado en los residuos de la regresión MCO. Se contrasta la hipótesis nula de normalidad:

\[\begin{equation} H_0)\ u\sim Normal \\ H_1)\ u\nsim Normal \end{equation}\]

Se computa el coeficiente de asimetría (\(S\)) y el coeficiente de kurtosis (\(K\)) de la distribución de los errores, para obtener el siguiente estadístico:

\[\begin{equation} JB=\frac{n}{6}\left[S^{2}+\frac{(K-3)^{2}}{4}\right] \end{equation}\]

Dado que una distribución normal tiene \(S=0\) y \(K=3\), el estadístico toma valores altos (en valor absoluto) cuando la distribución se aleja de la Normal. Entonces, bajo \(H_0\) se tiene que: \[\begin{equation} JB\overset{a}{\sim}\chi^{2}_2 \end{equation}\]

#install.packages("tseries")
library(tseries)

jarque.bera.test(residuos)

## 
##  Jarque Bera Test
## 
## data:  residuos
## X-squared = 650.81, df = 2, p-value < 2.2e-16

Se rechaza la hipótesis de normalidad de los errores.

8.2.2 Identificación de datos atípicos e influyentes

Una forma es visualizar el gráfico de los residuos:

par(mfrow=c(1,2))
plot(residuos) # dispersión
boxplot(residuos) # gáfico de cajas

También se puede utilizar la función outlierTest() del paquete car como se muestra a continuación, como resultado se obtienen las filas con valores mas extremos o atípicos.

#install.packages(car)
library(car)

outlierTest(modelo)

##     rstudent unadjusted p-value Bonferroni p
## 15  4.866985         1.5041e-06   0.00079117
## 229 4.849068         1.6401e-06   0.00086269
## 440 4.811042         1.9690e-06   0.00103570
## 186 4.746360         2.6795e-06   0.00140940
## 112 4.061856         5.6179e-05   0.02955000
## 59  4.029082         6.4339e-05   0.03384300

También se puede calcular la distancia de cook. Este indicador se calcula para una de las observaciones (individuos), mide el cambio en el valor estimado de \(y\) para todas las observaciones con y sin la presencia de la observación \(i\). Se busca medir cuánto impacta la observación o individuo \(i\) en los valores ajustados del MRLM.

Mas información de la distancia de cook se encuentra aquí: http://r-statistics.co/Outlier-Treatment-With-R.html

cooksd <- cooks.distance(modelo)
plot(cooksd, pch = "*", cex = 2)
abline(h = 4*mean(cooksd, na.rm = T), col = "red") 
text(x = 1:length(cooksd) + 1, y = cooksd, labels = ifelse(cooksd > 4*mean(cooksd, na.rm = T), names(cooksd), ""), col = "red")

Se pueden ver y analizar las filas influyentes de los datos originales. Se extrae y examina cada fila influyente solicitando ver todos los valores de los regresores con la segunda orden.

influential <- as.numeric(names(cooksd)[(cooksd > 4*mean(cooksd, na.rm = T))])  # influential row numbers
head(wage1[influential, ])  # influential observations

wage	educ	exper	tenure	nonwhite	female	married	numdep	smsa	northcen	south	west	construc	ndurman	trcommpu	trade	services	profserv	profocc	clerocc	servocc	lwage	expersq	tenursq
18.2	17	22	21	0	0	1	0	1	0	0	1	0	0	0	0	0	0	1	0	0	2.9	484	441
22.2	12	31	15	0	0	1	1	1	0	0	1	0	0	0	0	0	0	1	0	0	3.1	961	225
17.3	16	14	0	0	0	1	1	1	0	0	1	0	0	0	0	0	1	1	0	0	2.85	196	0
10	8	13	0	1	0	0	0	1	0	0	0	0	0	0	1	0	0	0	0	0	2.3	169	0
21.6	18	8	8	0	1	0	0	1	0	0	0	0	0	0	0	0	1	1	0	0	3.07	64	64
20	14	26	23	0	0	1	2	1	0	0	0	0	0	0	0	0	0	1	0	0	2.99	676	529

8.2.3 Contraste Reset

Este contraste tiene potencia para detectar formas funcionales incorrectas, pero no tiene potencia para detectar omisión de variables que se relacionen linealmente con los regresores incluidos (problema que veremos mas adelante como endogeneidad).

Se considera un modelo extendido que incluye además de los regresores con los que estamos trabajando las potencias de los valores ajustados de \(y\): \[\begin{equation} y=\beta_{0}+\beta_{1}x_{1}+\ldots+\beta_{k}x_{k}+\delta_{2}\hat{y}^{2}+\ldots+\delta_{m}\hat{y}^{m}+v \end{equation}\]

Las potencias de \(\hat{y}\) son funciones no lineales de las \(x_{j}\) Si el modelo inicial era incorrecto, agregar estos términos debe mejorar el ajuste significativamente. Esta es la escencia del contraste.

Observación: debería determinarse el orden del polinomio en \(\hat{y}\), aunque es sencillo ensayar con distintos órdenes.

Las hipótesis del contraste son las siguientes:

\[\begin{align*} H_{0}&)\ \delta_{2}=\ldots=\delta_{m}=0\\ H_{1}&)\ \text{Algún} \ \delta_{i}\neq0, \ i=2,\ldots,m \end{align*}\]

El estadístico de contraste es: \[\begin{equation} E=\frac{\left(SCR_{R}-SCR_{SR}\right)/(m-1)}{SCR_{SR}/(n-k-m)}\sim F_{m-1,n-k-m} \end{equation}\]

Es muy sencillo de aplicar pues no requiere especificar un modelo alternativo. Sin embargo, es poco constructivo, en el sentido que no brinda detalles de cómo proceder si el modelo original es rechazado.

CUIDADO!: Si se omitieron variables relevantes que están correlacionadas con las \(x_{j}\) es probable que lo estén con sus productos y potencias.

Para solicitarlo a R utilizaremos el comando resettest(formula, potencias, type, ... ),
del paqete lmtest, dónde:

• formula es la formula del modelo o el modelo que se testea.
• potencias es a que potencias se va a incluir.
• type indica si se incluyen las potencias del valor ajustado (fitted) o de los regresores (regresors)

#install.packages("lmtest")
library(lmtest)

resettest(modelo, potencia = 2:3, tipo = "ajustado")

## 
##  RESET test
## 
## data:  modelo
## RESET = 11.566, df1 = 2, df2 = 520, p-value = 1.217e-05

Si rechazo la hipótesis nula debo considerar un modelo alternativo.

9 Heteroscedasticidad

9.1 Introducción

En esta sección se analiza el supuesto de varianza condicional constante del término de error: \[Var(u|X_i)=\sigma^2 \ \forall \ i=1,\dots,n\]

Llamamos perturbaciones esféricas al caso en que la matriz de varianzas y covarianzas de los errores es escalar: \[\begin{align*} Var(u|X)=E(uu'|X)=\sigma^2I_n= \begin{bmatrix} \sigma^2 & 0 & \dots & 0 \\ 0 & \sigma^2 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \sigma^2 \end{bmatrix} \end{align*} \]

La heteroscedasticidad implica una matriz diagonal pero varianzas condicionales distintas para diferentes observaciones: \[ \begin{align*} E(uu'|X)=\sigma^2\Omega_n= \begin{bmatrix} \sigma^2_{1} & 0 & \dots & 0 \\ 0 & \sigma^2_{2} & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \sigma^2_{n} \end{bmatrix},\text{ con } \sigma^2_{i} \neq \sigma^2_{j} \text{ para algún }i\neq j \end{align*} \] La autocorrelación implica covarianzas no nulas entre los errores de distintas observaciones: \[\begin{align*} E(uu'|X)=\sigma^2\Omega= \begin{bmatrix} \sigma^2_{1} & \sigma_{12} & \dots & \sigma_{1n} \\ \sigma_{21} & \sigma^2_{2} & \dots & \sigma_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ \sigma_{n1} & \sigma_{n2} & \dots & \sigma^2_{n} \end{bmatrix},\text{ con } \sigma_{ij} \neq 0 \text{ para algún }i\neq j \end{align*} \]

En esta sección entonces, veremos como realizar contrastes para testear heteroscedasticidad y como estimar un MRLM en presencia de la misma. La heteroscedasticidad puede ser vista como un problema de parámetros no constantes, donde \(\sigma^2\) cambia para distintas observaciones o grupos.

Algunas aclaraciones:

• La varianza condicional es conocida como función scedástica.
• La heteroscedasticidad es el cambio de la varianza en la distribución del término de error condicional a los valores de los regresores, lo que implica que se trata de una variabilidad del error asociada a los regresores.
• Hay casos de varianza variable no asociada a los regresores, y eso NO es heteroscedasticidad.

9.2 Causas y consecuencias

Algunas causas comunes de heteroscedasticidad son:

• Varianza asociada al nivel de la variable dependiente.
• Agregación de datos.
• Outliers.
• Omisión de variables relevantes correlacionadas con los regresores.
• Forma funcional incorrecta (o transformación incorrecta de los datos).
• Errores de medida asociados a un regresor.

¿Cuáles son las consecuencias de la heteroscedasticidad?

• No produce sesgo ni inconsistencia de los estimadores \(MCO\).
• Los estimadores \(MCO\) dejan de ser \(MELI\) (Mejores Estimadores Lineales Insesgados) ya que no se cumple Gauss-Markov.
• La estimación de \(Var(\widehat{\beta}_{MCO}|X)\) es sesgada, y esto invalida la inferencia: las pruebas \(t\) y \(F\), los distintos contrastes asintóticos, los intervalos de confianza.

9.3 Métodos de detección

Si no hay información a priori sobre la heteroscedasticidad, una estrategia es estimar el modelo por \(MCO\) y analizar sus residuos.

Una manera es mediante la inspección gráfica. En general se realizan gráficos de dispersión de los residuos o sus cuadrados:

• En el tiempo si son series temporales.
• Contra algún regresor o combinación lineal de regresores (como \(\widehat{y}\)) en sección cruzada.
• Se trata de ver si \(\widehat{u}^{2}\) está sistemáticamente relacionado con \(\widehat{y}\) o con algún regresor del modelo.

9.3.1 Aplicación en R

Comenzamos cargando la base de salarios de Wooldridge:

library(wooldridge)
data("wage1")
attach(wage1)

Estimamos el siguiente modelo de regresión lineal múltiple: \[\begin{equation} wage_i = \beta_0 + \beta_1 educ_i + \beta_2 exper_i + \beta_3 antiguedad_i + u_i \\ \end{equation}\]

modelo <- lm(wage ~ educ + exper + tenure)

Calculamos los residuos y las predicciones del modelo:

uhat <- resid(modelo)
uhat2 <- (uhat)^2
yhat <- predict(modelo)

Generamos una nueva base datos, la cual contiene los residuos, el cuadrado de los residuos y las predicciones y los regresores utilizados en el modelo estimado.

datos <- cbind(wage1, uhat, uhat2, yhat)

Análisis gráfico de la heteroscedasticidad a través de la relación entre los cuadrados de los residuos y los regresores del modelo:

par(mfrow=c(1, 2))  # divide el área del gráfico en 2 columnas

plot(educ, uhat2, main = "Residuos vs educación")
plot(exper, uhat2, main = "Residuos vs experiencia")

plot(tenure, uhat2, main = "Residuos vs antiguedad")
plot(yhat, uhat2, main = "Residuos en función de las predicciones")

Se observa una varianza creciente de los residuos a mayores años de educación. A simple vista no se observan cambios en la varianza de los residuos que dependan de la experiencia y la antigüedad. Las predicciones del modelo de regresión lineal son una combinación lineal de los regresores, entonces al graficar los residuos contra las predicciones también se observan cambios en la varianza con respecto a los niveles de las predicciones. Parece existir heteroscedasticidad.

9.4 Contraste Breusch-Pagan (BPG)

El contraste puede incluir cierto “modelo” para la heteroscedasticidad.

Breusch y Pagan (1979) consideran heteroscedasticidad de la forma: \[\sigma_{i}^{2}=h(z_{i}'\alpha)=h(\alpha_0+\alpha_1Z_{1i}+\dots+\alpha_{p}Z_{pi})\] donde:

• Las \(p\) variables \(Z\) suelen coincidir con los \(k\) regresores \(X\).
• Las \(Z\) pueden incluir nuevas variables o excluir algunas \(X\).
• Se recomienda parsimonia (\(p\) bajo) y que \(Z\) incluya pocas variables nuevas.
• \(h\) es una función no especificada.

Godfrey (1978) propone una forma específica para la función \(h\): \[\sigma_{i}^{2}\text{=}e^{(\alpha_{0}+\alpha_{1}Z_{1i}+\dots+\alpha_{p}Z_{pi})}\]

Se busca contrastar (con cualquier función \(h\)): \[\begin{align*} H_{0}) & \ \alpha_1=\dots=\alpha_{p}=0 \\ H_{1}) & \ \text{Algún }\alpha_j\neq0,\,\ \text{con } j=1,\dots,p \end{align*}\]

9.4.1 Pasos del contraste

Estimar por \(MCO\) el modelo y obtener los residuos \(\widehat{u}_{i}\) y sus cuadrados \(\widehat{u}_{i}^2\)

Alternativa 1:

1. Calcular la estimación \(MV\) de la varianza del error bajo \(H_{0}\) cierta: \[\widetilde{\sigma}^{2}=\frac{1}{n}\sum_{i=1}^{n}\widehat{u}_{i}^{2}\]
2. Estimar por MCO la siguiente regresión auxiliar, y obtener la \(SCE_{aux}\): \[\left(\frac{\widehat{u}_{i}}{\widetilde{\sigma}}\right)^{2}=\alpha_0+\alpha_1 Z_{1i}+\dots+\alpha_p Z_{pi}+v_i\]
3. Bajo \(H_{0}\) y suponiendo \(v_i\sim N(0;\sigma^{2})\) para la regresión auxiliar: \[\frac{SCE_{aux}}{2}\overset{a}{\sim}\chi_{p}^{2}\]

Alternativa 2, sin reescalar:

1. Estimar por MCO la regresión auxiliar: \[\widehat{u}_i^2=\alpha_0+\alpha_1 Z_{1i}+\dots+\alpha_p Z_{pi}+v_i\]
2. Usar el estadístico \(ML\), que bajo \(H_0\) sigue una distribución \(\chi^2\): \[nR_{aux}^{2}\overset{a}{\sim}\chi_{p}^{2}\]

Comenzaremos por la segunda alternativa pues es un poco mas sencilla de implementar.

Estimamos la regresión auxiliar:

aux <- lm((uhat)^2 ~ educ + exper, tenure, data = wage1) 

R2 <- summary(aux)$r.squared
R2

## [1] 0.09354417

k <- 3
n <- nrow(wage1)
estadistico <- n * R2
estadistico

## [1] 49.20423

v.critic <- qchisq(0.95, k)
v.critic

## [1] 7.814728

valorP <- pchisq(q = estadistico, df = k, lower.tail = FALSE)
cbind(estadistico, valorP)

##      estadistico       valorP
## [1,]    49.20423 1.180183e-10

Para implementar esta alternativa se puede utilizar directamente la función bptest() de la librería lmtest de R.

#install.packages("lmtest")
library(lmtest)

bptest(modelo)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 43.096, df = 3, p-value = 2.349e-09

El \(p\)-valor es muy pequeño, por lo tanto rechazo la hipótesis nula y en coherencia con lo sugerido en el análisis gráfico existe heteroscedasticidad.

Ahora realicemos la alternativa 1, reescalando.

9.5 Contraste de White

Los pasos a seguir para realizar el contraste de White son los siguientes:

1. Estimar el modelo por \(MCO\) y obtener los residuos \(\widehat{u}_i\).

2. Estimar por \(MCO\) una regresión auxiliar de dichos residuos al cuadrado sobre todos los regresores, sus cuadrados y sus productos cruzados: \[\begin{align*} \widehat{u}_i^2=&\alpha_0+\alpha_1X_{1i}+\dots+\alpha_kX_{ki} +\alpha_{k+1}X_{1i}^2+\dots+\alpha_{2k}X_{ki}^2 +\alpha_{2k+1}X_{1i}X_{2i}+\dots+\alpha_{3k}X_{1i}X_{ki} \\ &+\alpha_{3k+1}X_{2i}X_{3i}+\dots+\alpha_{4k-1}X_{2i}X_{ki} +\dots+\alpha_{2k+k(k-1)/2}X_{(k-1)i}X_{ki}+v_{i} \end{align*}\]

3. Usar el \(R_{aux}^{2}\) para contrastar: \[\begin{align*} H_{0}) & \ E(u_i^2)=\sigma^{2} & H_0) \ &\alpha_1=\dots=\alpha_{2k+k(k-1)/2}=0 \\ H_{1}) & \ E(u_i^2)=\sigma_{i}^{2}\ &\ H_1) \ &\text{Algún }\alpha_j\neq0, \text{ con } j=1,\dots,2k+k(k-1)/2 \end{align*}\]

Para utilizar el test de White en R se puede utilizar la misma función que para el test de Breusch-Pagan, pero con el argumento varformula. Se debe escribir la fórmula con los términos no lineales, pues los lineales los considera por defecto.

bptest(modelo, varformula = ~ educ * exper + educ * tenure + exper * tenure + I(educ^2) + I(exper^2) + I(tenure^2), data = wage1)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 63.96, df = 9, p-value = 2.298e-10

La conclusión es la misma, se rechaza la hipótesis que asume homoscedasticidad, pues se encuentra evidencia de heteroscedasticidad.

#white_lm(ModelNR, interactions=T)

Verificamos que da lo mismo que hacerlo en forma manual:

aux2 <- lm(resid(modelo)^2 ~ educ + exper + tenure + educ * exper + educ * tenure + exper * tenure + I(educ^2) + I(exper^2) + I(tenure^2), data = wage1) 

R2 <- summary(aux2)$r.squared

estadistico <- n * R2
estadistico

## [1] 63.96001

valorP <- pchisq(q = estadistico, df = 9, lower.tail = FALSE)
cbind(estadistico, valorP)

##      estadistico       valorP
## [1,]    63.96001 2.297686e-10

Algunos comentarios:

• Los contrastes de White y BPG son contrastes generales (no requieren definir una forma de la heteroscedasticidad a priori), aunque asumen una forma implícita. Sus resultados son condicionales a que la heteroscedasticidad pueda ser ajustada de la forma propuesta.

• Los contrastes de White y de BPG son constructivos. Si se rechaza \(H_0\) el propio contraste sugiere qué variables intervienen en la forma de la heteroscedasticidad y se podría estimar utilizando mínimos cuadrados generalizados factibles (\(MCGF\)).

• Los contrastes de White y de BPG tienen validez asintótica. Se requiere una muestra grande, y si no se tiene deberán utilizarse otros contrastes.

• Con varias variables en el modelo original, se produce un aumento cuadrático en el número de variables en la regresión auxiliar. Ello hace que se pierden muchos grados de libertad y la prueba pierde potencia.

Una alternativa es aplicar la versión modificada (Wooldridge, 2004). Se utilizan potencias de los valores ajustados del modelo en lugar de los cuadrados de los regresores y sus productos cruzados. Como el valor ajustado es una combinación lineal de los regresores, las potencias de \(\widehat{y}\) son una combinación particular de esos cuadrados y productos cruzados. Se muestra que la regresión auxiliar de este contraste es una versión restringida de la del contraste original de White. Se conserva el espíritu del test de White, sin perder los grados de libertad.

9.6 Contraste de Goldfeld-Quandt

El contraste de Goldfeld-Quandt es válido con muestras pequeñas y funciona eliminando algunas observaciones ubicadas en el centro del conjunto de datos y luego probando para ver si la dispersión de los residuos es diferente de los dos conjuntos de datos resultantes que están a ambos lados de las observaciones centrales.

Se supone que la varianza del término de error depende de una variable \(X_j\)

Se busca contrastar un hipótesis nula de homoscedasticidad frente a una alternativa de heteroscedasticidad asociada a la variable \(X_j\)

9.6.1 Pasos del contraste

1. Identificar la variable \(X_j\) y ordenar los datos por \(X_j\) (varianza creciente).

2. Dividir los datos en 3 submuestras, omitiendo \(m\) observaciones centrales.

3. Estimar por \(MCO\) el modelo original en las submuestras inicial y final (cada una tiene \((n-m)/2\) observaciones)

4. Calcular la \(SCR\) en cada regresión y construir el estadístico de prueba, que bajo \(H_0\) sigue una distribución \(F\): \[\lambda=\frac{SCR_2}{SCR_1}=\frac{\widehat{\sigma}_2^2}{\widehat{\sigma}_1^2}\sim F_{\frac{n-m}{2}-k-1,\frac{n-m}{2}-k-1}\] Rechazamos \(H_0\) al nivel de significación \(\alpha\) si \(\lambda>c\) donde \(c\) cumple: \[Prob\bigg(F_{\frac{n-m}{2}-k-1,\frac{n-m}{2}-k-1}>c\bigg)=\alpha\]

5. La varianza al final debería ser mayor que al principio de la muestra.

6. El cuadrado de los residuos está asociado con \[\sigma^{2}\Rightarrow SCR_{2}>SCR_{1}\].

9.6.2 Observaciones y comentarios

Algunas observaciones y comentarios respecto al contraste de Goldfeld-Quandt son los que siguen:

• Si se sospecha que la varianza depende inversamente de los valores de \(X_j\Rightarrow\) se deberán ordenar los datos en forma decreciente y proceder de igual modo.

• Excluir “demasiadas” observaciones aumenta la potencia del test (más diferencia entre submuestras), pero se pierden muchos grados de libertad.

• Excluir “pocas” observaciones puede provocar que las submuestras sean muy parecidas y el test pierde potencia.

• Se recomienda entonces considerar \(m=n/3\).

• Los resultados son condicionales a elección de la variable \(X_j\).

• El procedimiento se puede repetir con todas las variables candidatas.

9.6.3 Aplicación en R

Para implementar la prueba Goldfeld-Quandt para determinar si hay heterocedasticidad en R, se sugiere utilizar la función gqtest() del paquete lmtest para realizar .

Esta función utiliza la sintaxis: gqtest (modelo, orden por, datos, fracción), dónde:

• modelo: es el modelo de regresión lineal creado por el comando lm().
• order.by: son las variables predictoras del modelo.
• data: es el nombre del conjunto de datos.
• fracción: es el número de observaciones centrales que se eliminarán del conjunto de datos.

gqtest(modelo, order.by = ~ educ, fraction = 2, data = wage1) 

## 
##  Goldfeld-Quandt test
## 
## data:  modelo
## GQ = 1.9534, df1 = 258, df2 = 258, p-value = 5.091e-08
## alternative hypothesis: variance increases from segment 1 to 2

Para estimar corrigiendo la heteroscedasticidad de los errores, se puede utilizar la función lm_robust() del paquete estimatr:

#install.packages("estimatr")
library(estimatr)

modelo <- lm(wage ~ educ + exper + tenure, data = wage1)
modelo_robusto <- lm_robust(wage ~ educ + exper + tenure, data = wage1)

Obtendremos las estimaciones con \(MCO\) y \(MCO\) robusto utilizando la función modelsummary() del paquete modelsummary:

#install.packages("modelsummary")
library(modelsummary) 

modelsummary(list("MCO" = modelo, "MCO robusto" = modelo_robusto), gof_omit = 'IC|Log|Adj|p\\.value|statistic|se_type', stars = TRUE)

	MCO	MCO robusto
(Intercept)	-2.873***	-2.873***
	(0.729)	(0.812)
educ	0.599***	0.599***
	(0.051)	(0.061)
exper	0.022+	0.022*
	(0.012)	(0.011)
tenure	0.169***	0.169***
	(0.022)	(0.030)
Num.Obs.	526	526
R2	0.306	0.306
F	76.873
RMSE	3.08
Std.Errors		HC2
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

10 Multicolinealidad

Con multicolinealidad exacta no se puede obtener la estimación \(MCO\), sin embargo es frecuente tener multicolinealidad aproximada. En este caso es posible obtener el estimador \(MCO\) utilizando la expresión habitual: \(\widehat{\beta}_{MCO}=(X'X)^{-1}X'y\).

Con esto, no se viola ninguno de los supuestos clásicos, por lo tanto los estimadores \(MCO\) siguen siendo \(MELI\) (Mejores Estimadores Lineales Insesgados) por el teorema de Gauss-Markov. Sin embargo, las varianzas y covarianzas de los estimadores se incrementan mucho en presencia de multicolinealidad aproximada. Esto implica que los estimadores se vuelven menos precisos e implica:

• Intervalos de confianza grandes.
• En los contrastes de significación individual tiende a no rechazarse \(H_0\).
• Cambios mínimos en la matriz \(X'X\) afectan fuertemente las estimaciones.

El problema es equivalente al caso en que se cuenta con pocas observaciones (que Goldberger llama irónicamente micronumerosidad).

10.1 Efectos de la multicolinealidad

10.2 Multicolinealidad exacta: aplicación en R

library(wooldridge)
data("wage1")
attach(wage1)

Generamos una variable llamada educ1, como combinación lineal de educación, que no considera los primeros 6 años de educación primaria:

educ1 <- educ - 6
wage2 <- cbind(wage1, educ1)

mod <- lm(wage ~ educ+ educ1 + exper + tenure, data = wage2)
summary(mod)

## 
## Call:
## lm(formula = wage ~ educ + educ1 + exper + tenure, data = wage2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.6068 -1.7747 -0.6279  1.1969 14.6536 
## 
## Coefficients: (1 not defined because of singularities)
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.87273    0.72896  -3.941 9.22e-05 ***
## educ         0.59897    0.05128  11.679  < 2e-16 ***
## educ1             NA         NA      NA       NA    
## exper        0.02234    0.01206   1.853   0.0645 .  
## tenure       0.16927    0.02164   7.820 2.93e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.084 on 522 degrees of freedom
## Multiple R-squared:  0.3064, Adjusted R-squared:  0.3024 
## F-statistic: 76.87 on 3 and 522 DF,  p-value: < 2.2e-16

Observar que el R quita automáticamente la variable educ1 para poder estimar por MCO.

Para mostrar los efectos de la multicolinealidad aproximada se utilizará la simulación del siguiente libro: https://fhernanb.github.io/libro_regresion/multicoli.html

Se simulan dos variables, por un lado \(X_1\) con resultados de una distribución uniforme con recorrido \((0, 10)\) y una variable \(X_2\) con multicolinealidad aproximada con \(X_1\).

\[X_2 = 2*X_1 + ruido\]

Si no estuviera el término del ruido, la multicolinealidad sería exacta.

Se simula \(y = -3+2*X_1 -4*X_2\). Esto hace que se conozca el verdadero proceso generador de datos de \(y\) y los verdaderos valores de los parámetros serían:

\(\beta_0 = -3\), \(\beta_1 = 2\), \(\beta_2 = -4\) y \(\sigma = 2\)

gen_dat <- function(n) {
  x1 <- runif(n=n, min=0, max=10)
  x2 <- x1 * 2 + rnorm(n=n, sd=0.01)  #x2 es el doble de x1 + ruido
  y <- rnorm(n=n, mean= - 3 + 2 * x1 - 4 * x2, sd=2)
  data.frame(y, x1, x2)
}

En total se simulan 40 observaciones y se guardan en dos data frames diferentes, con 20 observaciones cada uno. Las primeras 20 van al data frame datos1 y las restantes 20 al data frame datos2. Se fija la semilla en el valor 12345 para que el lector pueda replicar el ejemplo y obtener los mismos resultados.

set.seed(12345)
datos <- gen_dat(n = 40)
datos1 <- datos[1:20,]
datos2 <- datos[21:40,]

Se estiman los modelos con cada uno de los data frame:

mod1 <- lm(y ~ x1 + x2, data = datos1)
mod2 <- lm(y ~ x1 + x2, data = datos2)

modelsummary(list("Datos1" = mod1, "Datos2" = mod2), gof_omit = 'IC|Log|Adj|p\\.value|statistic|se_type', stars = TRUE)

	Datos1	Datos2
(Intercept)	-1.560	-2.402+
	(1.135)	(1.140)
x1	100.994	-36.608
	(94.085)	(81.561)
x2	-53.574	15.323
	(47.073)	(40.776)
Num.Obs.	20	20
R2	0.986	0.982
F	617.474	461.293
RMSE	2.34	2.30
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

Se pueden observar las diferencias entre las estimaciones en ambos modelos, sin embargo, son dos muestras del mismo proceso generador de datos. A su vez, observar la diferencia con los verdaderos valores de los parámetros: \(\beta_0 = -3\), \(\beta_1 = 2\), \(\beta_2 = -4\) y \(\sigma = 2\)

Tenemos el típico de multicolinealidad; modelo globalmente significativo, pero ninguno de los regresores es relevante (son no significativos).

summary(mod1)

## 
## Call:
## lm(formula = y ~ x1 + x2, data = datos1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.9549 -1.7802  0.5549  1.5669  4.0319 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   -1.560      1.135  -1.374    0.187
## x1           100.994     94.085   1.073    0.298
## x2           -53.574     47.073  -1.138    0.271
## 
## Residual standard error: 2.341 on 17 degrees of freedom
## Multiple R-squared:  0.9864, Adjusted R-squared:  0.9848 
## F-statistic: 617.5 on 2 and 17 DF,  p-value: < 2.2e-16

summary(mod2)

## 
## Call:
## lm(formula = y ~ x1 + x2, data = datos2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7246 -1.8012  0.2883  1.2949  3.9046 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   -2.402      1.140  -2.107   0.0503 .
## x1           -36.608     81.561  -0.449   0.6592  
## x2            15.323     40.776   0.376   0.7117  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.303 on 17 degrees of freedom
## Multiple R-squared:  0.9819, Adjusted R-squared:  0.9798 
## F-statistic: 461.3 on 2 and 17 DF,  p-value: 1.545e-15

Ahora bien, podemos tener multicolinealidad sin observar este síntoma. Veamos esta simulación donde la variable \(X_4\) no depende del resto en el proceso generador de datos:

gen_dat <- function(n) {
  x1 <- sample(5:25, size = n, replace = TRUE)
  x2 <- rpois(n, lambda = 5)
  x3 <- rbinom(n, size = 10, prob = 0.4)
  x4 <- rbeta(n = n, shape1 = 0.5, shape2 = 0.7)
  ruido1 <- runif(n = n, min = -0.5, max = 0.5)
  ruido2 <- runif(n = n, min = -0.5, max = 0.5)
  x5 <- x1 - 3 * x2 + ruido1
  x6 <- x2 - 4 * x3 + ruido2
  y <- rnorm(n = n, mean = - 3 + 2 * x1 - 4 * x2, sd = 2)
  data.frame(y, x1, x2, x3, x4, x5, x6)
}

datos3 <- gen_dat(n = 30)

No tenemos el síntoma típico que se da en la inconsistencia entre el test \(F\) y los test \(t\) de significación individual, sin embargo el modelo tiene problemas de multicolinealidad aproximada.

10.2.1 Detección de la multicolinealidad

Calcularemos los factores de inflación de varianzas (\(FIV\))

Se define el Factor de Incremento de la Varianza (\(FIV\)) como: \[FIV_{j}=\frac{Var(\widehat{\beta}_{j})}{Var(\widehat{b}_{j})}=\frac{1}{1-R_{j}^{2}}\] Si no existe multicolinealidad el \(R_{j}^{2}\) va a ser cercano a cero y en este caso: \[Var(\widehat{\beta}_{j})=\frac{\sigma^{2}}{SCT_j}= Var(\widehat{b}_{j})\]

El \(FIV\) puede ser interpretado como el aumento de la varianza del estimador \(MCO\) con respecto al caso en que no existe relación entre el regresor \(j\) y el resto de los regresores del modelo.

El \(FIV\) muestra en qué medida se agranda la varianza del estimador como consecuencia de la no ortogonalidad de los regresores.

\[ FIV_{j}=\frac{V(\hat{\beta}_{j})}{V(\hat{b}_{j})}=\frac{1}{1-R_{j}^{2}} \]

Otro concepto a definir es el Factor de Tolerancia de la Varianza (\(TOL\)):

\[ TOL_{j}=1-R_{j}^{2} \]

\(FIV_{j}>10\) ó \(TOL_{j}<0.1\) son indicadores de multicolinealidad grave.

Para calcular los factores de inflación de varianzas en R, se utiliza la función vif() del paquete car. Esta función devuelve el \(FIV\) de cada uno de los regresores del modelo. Obviamente son los mismos ambos conjuntos de datos y por ende ambos modelos.

library(car)

vif(mod1)

##       x1       x2 
## 294997.4 294997.4

vif(mod2)

##       x1       x2 
## 172810.7 172810.7

Se puede observar que los mismos son muy superiores al nivel crítico 10 pues existe multicolinealidad aproximada.

Ajustamos un nuevo modelo, con el otro subconjunto de datos:

mod3 <- lm(y ~ ., data = datos3)
vif(mod3)

##         x1         x2         x3         x4         x5         x6 
##  540.50557 1408.06675  298.35108    1.05934 1979.04151  292.87883

summary(mod3)

## 
## Call:
## lm(formula = y ~ ., data = datos3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4429 -1.2756 -0.2266  0.9174  3.9270 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -2.73053    1.46522  -1.864  0.07520 . 
## x1            4.72654    1.27255   3.714  0.00114 **
## x2          -12.46779    3.88771  -3.207  0.00391 **
## x3            1.05875    4.29016   0.247  0.80726   
## x4            0.06121    1.02178   0.060  0.95275   
## x5           -2.73795    1.28350  -2.133  0.04380 * 
## x6            0.29167    1.07228   0.272  0.78804   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.711 on 23 degrees of freedom
## Multiple R-squared:  0.9922, Adjusted R-squared:  0.9901 
## F-statistic: 486.2 on 6 and 23 DF,  p-value: < 2.2e-16

Ya no tenemos la inconsistencia entre la prueba \(F\) con las pruebas de significación individuales, sin embargo el modelo tiene multicolinealidad aproximada y eso se ve claramente en los FIV.

11 Endogeneidad