Modelos de regresión lineal
Econometría1-FCEA-UDELAR
2022-06-19
0.1 Introducción
Veremos las herramientas básicas que necesiarás en el curso de Econometría 1. Veremos como instalar los softwares, como descargar y trabajar con las 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. Irenemos compartiendo códigos, que esperamos te resulten de utilidad, tanto para el curso, como para el posterior desarrollo de tu ejercicio profesional.
0.2 Presentación de R
¿Qué es 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.
0.3 Descarga e instalación de R y Rstudio
Deberás descargar e instalar R y Rstudio.
- Descargar R aquí
- Descargar RStudio aquí
https://www.rstudio.com/products/rstudio/download/
- VIDEO: Tutorial para instalación
https://www.youtube.com/watch?v=ZTzbiAznSjc
- Si queremos trabajar en la nube (sin necesidad de instalar nada en tu PC):
0.4 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 ** > ** , aprece 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 imágen, se puede dividir (vertical u horizontalmente) la interfaz en 4 partes o ventanas:
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 R Markdown o R Notebook.
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: 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.
0.4.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 ggplot()función 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 continue funcionando. Veremos estos mensajes al instalar paquetes.
0.4.2 Librerias 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
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:
Haz click en Packages > Install, y se abre una
ventana emergente.
Escribe el nombre del paquete que quieras instalar, finalmente click
en Install.
Utilizaremos el paquete dplyr.
Algunas funciones del paquete dplyr:
filter(): filtra o elige las observaciones por sus valores.arrange(): reordena 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 resumengroup_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 para aprender R desde la consola que contiene cursos interactivos: install.packages(“swirl”) library(swirl)
0.4.3 Comandos básicos indispensables
- Puedes usar R como una calculadora:
1 / 200 * 30## [1] 0.15(59 + 73 + 2) / 3## [1] 44.66667sin(pi / 2)## [1] 1Puedes crear objetos nuevos usando
<-:nombre_objeto <- valorAlt+ - (signo menos).Los nombres de los objetos deben comenzar con una letra y solo pueden contener letras, números, _ y ..
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.
0.4.4 Clases de objetos
Vectores
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 particularideades 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:
0.5 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.
0.6 ¿Qué es un project?
- Directorio de trabajo (working directory): aquí 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 sesión.
getwd() dice la ruta de donde estoy trabajando.
Para cambiar el directorio: 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.
Hagamos un proyecto:
Haz clic en File > New Project, y después:
0.7 Importar y exportar datos
Utilizamos la libreria
readrdel paquetetidyverse.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()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§ion=10
0.7.1 Para importar datos
- Otra forma de importar datos:
Haz click en Environment > Import Dataset.
Se abre una ventana emergente:
0.7.2 Para exportar datos
write(x, file = “data”) write.table(x, file=paste(“../../..”), “x._”))
x= el nombre del objeto file=paste(“../../..”): ponemos la ruta donde queremos guardar el archivo. x._: es el nombre que le ponemos al archivo para guardarlo.
Guardar en excel: write.csv(x, file=paste(“../../..”), “x.csv”))
0.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.
0.9 Buscando ayuda
Puedes encontrar material de ayuda en:
- Sitio web de la versión en español de “R for Data Science”, de Hadley Wickham y Garrett Grolemund.
- “Machetes” (cheatsheets) de Rstudio.
https://www.rstudio.com/resources/cheatsheets/
- En particular el cheatsheet de importación de datos.
0.9.1 Busquemos ayuda en R
- en la consola escribe
?+ nombre de la función - en la concola escribe
??+ algo parecido al nombre de la función - en la pestaña
Helpde la salida
0.10 Introducción al curso de Econometría 1
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 R studio (software que vemos en el curso). Aquí se encuentran detalles para utilizar R Markdown 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, 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 encontrará 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 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 en 1988.
0.11 Modelo de regresión lineal simple (MRLS)
0.11.1 Cargar y manipular la base
Para cargar los datos se escribe data(“wage1”) # carga la base salarios base <-wage1 # la asigno al objeto 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 (antiguedad)
- nonwhite: =1 si no la persona no es blanca
- female: =1 si la persona es mujer
- married: =1 si la persona es casada
Cambiamos 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
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 4También se puede seleccionar utilizando el nombre de las variables:
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 2Otra forma mas sencilla es utilizar la función select del paquete dplyr, se usa el comando llamado pipe %>%, el mismo enlaza las órdenes de manera secuencial.
Si es la primera vez que se usa recordar instalarlo previamente: #install.packages(“dplyr”)
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 2Se 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.
datos3<-base %>%
select("wage", "educ", "exper", "antiguedad", "married")%>%
filter(married==1)
#View(datos3)0.11.2 Investigar las correlaciones
plot(datos2)
Calculo la matriz de correlaciones de las variables que figuran en datos
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.00000000Existe una correlación positiva entre el salario y la educación (0.4059)
Creamos el diagrama de dispersión entre salario y educación. Utilizaremos 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")
0.11.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}\]
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-16anova(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 ' ' 10.11.4 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.5413593Para consultar la estimación de un coeficiente de regresión se utilizan paréntesis rectos.
ahat <-coef(summary(mod1))[1,1]
bhat <- coef(summary(mod1))[2,1]
ahat## [1] -0.9048516bhat## [1] 0.5413593Miramos los componentes del modelo y los asignamos a nuevos objetos
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$residualsdatos2$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.5914590.11.5 Diagrama de dispersión y recta de regresión
Creamos el gráfico de dispersión entre el salario y la educación e incluímos la estimación de la recta de regresión.
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. Posicionándose encima de cada observación se ven los valores de (x, y) para cada uno de los individuos. Necesita la función plotly.
library(plotly)
plot_ly(data=datos2, x=~educ, y=~wage)0.11.6 Trabajar con los residuos del modelo
Calculamos los residuos del modelo:
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.2914593Comparamos los residuos calculados 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 \]
datos2$resid<- resid
datos2$verif<-datos2$residuos-datos2$residAgrego 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()
Inspección gráfica de los resididuos
plot(datos2$residuos)
Graficamos las funciones de densidad estimadas
library(e1071) # for skewness function
par(mfrow=c(1, 2)) # divide graph area in 2 columns
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")
0.11.7 Predicción e intervalos de confianza
Genero la predicción puntual y el intervalo de confianza para una educación de 10 años, prediction devuelve el valor para la predicción puntual
Calculo el salario esperado para una persona con la educación promedio (\(y/x=\overline{x}\))
mean(datos2$educ)## [1] 12.56274sal_pred<-0.90485+0.054136*12.56
sal_pred## [1] 1.584798Calculo el intervalo de confianza para una persona con la educación promedio (\(x=12.56\))
nuevo <- data.frame(educ=12.56)
future_y<-predict(object=mod1, newdata=nuevo, interval="prediction", level=0.95)Generamos el intervalo de confianza para \(E(y/x)\). Debemos cambiar el argumento de 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$uprAgregamos 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)Generamos el gráfico con los intervalos de confianza para la predicción puntual y para el valor esperado.
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("Precicción de y al 95%")
library(ggplot2)
IC_esp_y_95<-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%")Imprimo los gráficos juntos para poder compararlos mejor. ¿Cuál de los dos tiene mayor amplitud?
#nstall.packages("gridExtra")
library(gridExtra)
grid.arrange(IC_esp_y_95, IC_y, ncol=2,nrow=1)
Veamos el impacto que tiene el nivel de confianza en la amplitud de los intervalos.
library(ggplot2)
#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)
library(gridExtra)
grid.arrange(a, b, ncol=1,nrow=2)
0.12 Modelos de regresión lineal múltiple (MRLM)
Estimaremos modelos de regresión lineal múltiple utilizando la misma base de datos de Wooldridge, wage1 que contiene datos de la Encuesta de población de 1976.
0.12.1 Acceder y conocer la base
Para cargar los datos se escribe data(“wage1”) # carga la base salarios base <-wage1 # la asigno al objeto base con el que trabaremos
Accedemos a la base, la asignamos al objeto base y solicitamos 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 (antiguedad)
- nonwhite: =1 si no la persona no es blanca
- female: =1 si la persona es mujer
- married: =1 si la persona es casada
- lwage: log(salario)
- expersq: exper ^2, el cuadrado de los años de experiencia potencial
- tenursq: ternure^2, el cuadrado de la antiguedad en el trabajo actual
Cambiamos 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. Las asignamos al objeto datos.
datos<-base[,c(1: 7, 22, 23, 24)] Para mirar la dimensión de la base resultante solitamos la dimensión de la misma.
dim(datos)## [1] 526 10Podríamos solicitar el número de columnas (personas en este caso) y de columnas por separado (variables) con las siguientes funciones:
nrow(datos)## [1] 526ncol(datos)## [1] 10Miramos 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 4Solicitamos 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.000.12.2 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.000Para que reconozca el nombre de las variables sin necesidad de citar el nombre del data frame con el signo de $
attach(datos)Exploramos, mediante un gráfico en tres dimensiones, la relación entre el salario, los años de educación y la antiguedad en el cargo que tienen las personas. Se puede cambiar la perspectiva del gráfico con el puntero del mousse.
#install.packages("plotly")
library(plotly)
library(ggplot2)
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)
attach(datos)
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)')
0.12.3 Modelo para el salario
Estimo un primer modelo en el que se explica el nivel del salario en función de los años de educación y de 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-16Ajustamos 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')
0.12.4 Análisis de los residuos
Calculamos la suma de cuadrados de los residuos
SCR <- deviance(mod1) #Suma del cuadrado de los residuos
gl <- df.residual(mod1) # Grados de libertad del modelo
sigma2_gorro <- (SCR/gl)
#otra manera es utilizar la función sigma que devuelve el error estándar de los residuos del modelo.
sigma2_gorro <- (sigma(mod1))^2
SCR## [1] 4998.965gl## [1] 523sigma2_gorro## [1] 9.55825También podemos utilizar la función anova que nos devuelve:
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 ' ' 1Inspeccionamos la normalidad de los residuos mediante análisis gráficos
library(e1071) # for skewness function
plot(density(mod1$residuals), main="Density Plot: residuos", ylab="Frequency", sub=paste("Skewness:", round(e1071::skewness(mod1$residuals), 2))) # density para los 'residuos'
polygon(density(mod1$residuals), col="red")
Solicitamos los siguientes gráficos para inspeccionar la normalidad
de los errores a través de la distribución empírica de los residuos (que
es lo que observamos).
Permite 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-16El 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.
Graficamos los residuos contra cada uno de los regresores
library(ggplot2)
library(gridExtra)
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.
###¿El género 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 resumenes del salario en función del sexo biológico:
library(tidyverse)
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.6Observamos estadísticos de posición mayores para los hombres que para las mujeres.
Asignamos etiquetas a la variable female. La recodificación es en el mismo orden en el que está codificada la variable: 0 es hombre y 1 es mujer.
datos$género<- factor(datos$female,
labels = c("hombre", "mujer"))Observamos el diagrama de dispersión de los salarios por el género.
boxplot(wage ~ género, data = datos)
Estimamos un modelo de regresión múltiple agregando la variable female.
\[\begin{equation} wage_i = \beta _0+\beta _1 educ_i +\beta _2 antiguedad_i+ \beta _3 antiguedad_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-16Utilizamos el comando 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. |
Para concatenar múltiples modelos en una sola tabla. https://jtools.jacob-long.com/reference/export_summs.html
#install.packages("huxtable")
library(huxtable)
library(jtools)
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. | |||
0.12.5 Test de significación global y conjunta
Queremos probar si la experiencia es relevante en el modelo de salarios.
data("wage1") # carga la base salarios
base <-wage1 # la asigno al objeto base con el que trabaremos 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+ \beta _5 female_ +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-16El 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\) parmámetros. En este caso \(q=2\) pues testear si la experiencia no es relevante implica asumir que sus 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+ \beta _3 female_ +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 SCR de ambos modelos. 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.72523Otra manera de calcular el valor del estadístico \(F\) es comparando los coeficientes de determinación \(R2\) de los modelos restringidos (r) y sin restringir (sr).
\(F=\frac{\left(R^2_{nr}-R^2_{r}\right)/q}{1- R^2_{nr}/(n-k-1)}\)
Observación: cuidado con el orden, 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.
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.72523Finalmente, existe una tercera forma de realizar el contraste planteado y es solicitando a R que lo realice por comando con la función linearHypothesis. Requiere el paquete gmodels y car.
Aquí se encuentra la documentaciónd 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- Cuidado! La notación de la función linearHypothesis puede dar lugar a confisión. No estamos testeando si las vaariables 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 comparando 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 utililizar la 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 : los grados de libertad del numerador
- df2 : los grados de libertad del denominador
- lower.tail: Si es VERDADERO, 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 predeterminado es VERDADERO.
alfa=0.05
v.crit.conj <- qf(alfa, 2, ModeloNR$df.residual, lower.tail = FALSE)
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." print(v.crit.conj)## [1] 3.0130570.12.6 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)\) el número de parámetros del modelo no restringido.
\[ \mathrm{H}_{1})\,\mathrm{R}\beta\neq\mathrm{r} \]
\(H_{1})\) implica que al menos una de las igualdades no se cumple. Cuidado!! No es lo mismo que todas las igualdades no se cumplan.
Algunos ejemplos:
Supongamos un modelo con 4 regresores y queremos realizar la siguiente prueba de hipótesis:
\(H_{0}: \beta_{1}=0\) versus \(H_{1})\)\(\beta_{1}\neq 0\)
\[ \mathrm{H}_{0})\,\mathrm{R}\beta=\mathrm{r} \]
\[ \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) \]
Veamos este ejemplo donde 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\), 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 del R para realizar los contrastes.
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 H0 cierta
- R : la matriz de las restricciones lineales.
Para escribir la matriz R, poner los valores de la matriz, leyendo por filas utilizando c(). Cada fila contiene una combinación lineal de los coeficientes de regresión.
Ejemplo para escribir la matriz R en el último ejemplo:
R <- rbind( c(0,0,0,1,0), c(0,0,0,0,1) )
print(R)## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 0 0 1 0
## [2,] 0 0 0 0 1Una manera alternativa de escribir la matriz es especificando todos sus componentes de corrido, pero especificando el número de filas que queremos y que vamos a completar la matriz por filas.
R <- matrix (c(0,0,0,1,0,0,0,0,0,1 ), nrow =2, byrow = T)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-16Utilizamos la función glh para realizar el contraste.
H0 <- glh.test(Modelo, R, 0)Devuelve un objeto tipo lista con los siguientes componentes que los puedo llamar por sus respectivos nombres.
names(H0)## [1] "call" "statistic" "parameter" "p.value" "estimate"
## [6] "null.value" "method" "data.name" "matrix"H0$matrix## (Intercept) female educ exper expersq
## [1,] 0 0 0 1 0
## [2,] 0 0 0 0 1H0$p.value## Estimate
## Estimate 1.332268e-15H0$estimate## Estimate
## [1,] 0.255127617
## [2,] -0.004439636Decidimos como en cualquier contraste con el p-valor.
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."Veamos que el resultado coincide para la prueba de significación inidividual. Escribimos la matriz R indicando los componentes, indicando el 1 en la posición del coeficiente de regresión que queremos testear si es significativo. nrow especifica que quiero una matriz de q filas bynrow le indica que quiero que la matriz se complete por filas
q=1
R1 <- matrix (c(0,0,0,1,0), nrow =q, byrow = T)
View(R1)
H0 <- glh.test(Modelo, R1, 0)Miramos los componentes y los comparamos con los obtenidos en el summary para verificar que es equivalente a la prueba de significación individual para testar 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.5406H0$p.value## Estimate
## Estimate 9.640067e-13H0$estimate## Estimate
## [1,] 0.2551276summary(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-16Finalmente, chequearemos 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))
print(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 1H0 <- glh.test(Modelo, R, 0)H0$statistic## F
## 70.16953H0$p.value## Estimate
## Estimate 0H0$estimate## Estimate
## [1,] -2.114034720
## [2,] 0.556284789
## [3,] 0.255127617
## [4,] -0.004439636summary(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-160.13 Variables ficticias o dummies
0.13.1 Introducció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 (antiguedad)
- nonwhite: =1 si no la persona no es blanca
- female: =1 si la persona es mujer
- married: =1 si la persona es casada
- numdep: =1 número de dependientes que tiene la persona
- lwage: log(salario)
0.13.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-16anova(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 |
0.13.3 Crear variables dummies
Para crear dummies se puede utilizar el paquete fastDummies. Creamos una dummie para cada valor de la variable numdep. Crea una dummie por cada categoría, en este caso son 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 las listo con c() donde listo el nombre de las variables separados por coma.
#dummy_cols(datos, select_columns=c("numdep", "exper"))0.13.4 Crear interacciones entre los regresores
Como ejemplo se van crear las interacciones entre estado civil y género.
Se crean las dummies variables
- mujer casada: =1 si la persona es mujer casada
- hombre casado: =1 si la persona es hombre casado
- mujer soltera: =1 si la persona es mujer soltera
- hombre soltero=1 si la persona es hombre soltero
Para ello es necesario crea antes la variable hombre y la variable soltero que no están en la base como tal. Utilizaremos la librería car:
# 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"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"0.13.5 Estimamos modelos para el salario
El modelo mas básico estima el log del salario en función de la educación de las personas. \[\begin{equation} log(salario_i) = \beta _0+\beta _1 Educ_i + u_i \\ \end{equation}\]
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-16Nos preguntamos si el efecto de la educación es el mismo para hombres y para mujeres, para ello corremos los modelos considerando dos bases separadas por sexo (mujeres y hombres)
library(dplyr)
mujeres<- datos1%>%
select("lwage", "educ", "exper", "tenure") %>%
filter(female==1)
hombres<- datos1%>%
select("lwage", "educ", "exper", "tenure") %>%
filter(female==0)
dim(mujeres)## [1] 252 4dim(hombres)## [1] 274 4Corremos el mismo modelo por un lado con las mujeres y por otro con los hombres:
\[\begin{equation} log(salario_mi) = \beta _0+\beta _1 Educ_mi + u_mi \\ \end{equation}\]
\[\begin{equation} log(salario_hi) = \beta _0+\beta _1 Educ_hi + u_hi \\ \end{equation}\]
sal_mujeres<-lm(lwage~ educ, data = mujeres)
sal_hombres<-lm(lwage~ educ, data = hombres)
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-13summary(sal_hombres)##
## Call:
## lm(formula = lwage ~ educ, data = hombres)
##
## 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-14library(jtools)
library(huxtable)
export_summs(modelo1, sal_mujeres, sal_hombres,
modelo.nombres = c("General","Mujeres","Hombres"))| 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. | |||
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:
datos1$Sexo_biológico<- factor(datos1$female,
labels = c("hombre", "mujer"))Graficamos los diagramas de dispersión del salario con respecto al sexo biológico para ver si existen diferencias:
library(ggplot2)
ggplot(datos1, aes(x=educ, y=lwage, color=Sexo_biológico)) +
geom_point() +
geom_smooth(method='lm', formula=y~x, se=FALSE, col='dodgerblue1') +
theme_light() +
facet_wrap(~Sexo_biológico)+
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_biológico)) +
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")
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-16boxplot(wage ~ Sexo_biológico, data = datos1)
library(tidyverse)
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.42La 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).
0.13.6 Modelos para el salario que consideran interacciones entre regresores
Estimaremos un modelo con interacciones para dar la posibilidad de que el retorno de la educación sea diferente para hombres que para mujeres.
\[\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-16Se 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 hombres.
Estimamos un modelo para el logaritmo del salario, en función del sexo biológico, la educación, la experiencia, la antiguedad 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 hombres? ¿Y para las mujeres?
El rendimiento estimado de la educación para los hombres 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\)), aprox \(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 hombres 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.
names(datos1)
\[\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-16La categoría de referencia en este caso son los hombres 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 antiguedad (tenure), se estima que los hombres 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 hombre soltero Se estima que, a igualdad del resto de las características, una mujer soltera, gana en promedio 11% menos que un hombre soltero.
0.14 Test o contraste de Chow
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 (hombres y mujeres, paíeses, 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 muestra y haciendo regresiones separadas 2. Generando variables ficticias e incluyendo interacciones
Obtengo la base salarios y llamo al paquete stargazer:
#install.packages("wooldridge")
library(wooldridge)
data("wage1")
#install.packages("stargazer")
library(stargazer)
attach(wage1)Implementación del test con la primera forma
Para implementar el contraste estimaremos tres modelos:
- el modelo para la muestra completa (NR)
- para las muestras entre las cuales se desea testear el quiebre estructural
En este caso se estima el modelo para hombres y para mujeres, utilizando 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.01A 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 test, para finalmente compararlo con el valor del estadístico de tabla a un determinado nivel de significación.
\(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.
\(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
Obtenemos las SCR para los tres modelos:
SRCnr<-sum(residuals(modelNR)^2)
SRCn1<-sum(residuals(modeln1)^2)
SRCn2<-sum(residuals(modeln2)^2)
print(SRCnr)## [1] 101.4556print(SRCn1)## [1] 39.0353print(SRCn2)## [1] 49.54719Bajo \(H_0\) cierta \(F\sim F_{(k+1),(n-2k-2)}\)
El número de restricciones \(q=(k+1)\) pues si probar \(H_0\) implica testear que los parámetros asociacos a las interacciones es igual a cero (hay k) y la diferencia en los interceptos (1).
Los grados de libertad del modelo no restringido son \(n-2(k+1)= n-2k-2\). n menos el número de parámetros (2k asociados a los regresores y 2 más de los interceptos de los 2 grupos).
Obtenemos los grados de libertad:
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.
gl <- (nrow(base) - 2*q) #q=(k+1) parámetros
print(q)## [1] 4print(gl)## [1] 518Calculamos el estadístico F en base a los tres modelos:
F<-((SRCnr- (SRCn1+SRCn2))*gl)/((SRCn1+SRCn2)*q)
Fcrit<-qf(0.05, q, gl, lower.tail=FALSE)
print(F)## [1] 18.81934print(Fcrit)## [1] 2.389145Implementación del test con la segunda forma
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 antiguedad. 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: $FF_{(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 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.81934print(Fcrit)## [1] 2.389145Tomamos 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 hombres y mujeres para explicar el salario")} else {
print("No rechazo H_0, no existe evidencia empírica suficiente para especificar modelos diferentes para hombres 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 hombres y mujeres para explicar el salario"0.15 Heteroscedasticidad
0.15.1 Introducción
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 al caso en que la matriz de varianzas y covarianzas de los errores es escalar: \[\begin{align*} \footnotesize{ 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 implica una matriz diagonal, pero varianzas condicionales distintas para diferentes observaciones: \[ \begin{align*} \footnotesize{ 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 implica covarianzas no nulas entre los errores de distintas observaciones: \[\begin{align*} \footnotesize{ 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 se verá 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.
- 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
0.15.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 (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
0.15.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 C.L 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
Cargamos 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 la base “datos” que 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 graph area in 2 columns
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 antiguedad. 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 heteroscedasticiadad.
0.15.3.1 Contraste Breusch-Pagan
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*}\]
Pasos del contraste:
Estimar por MCO el modelo y obtener los residuos \(\widehat{u}_{i}\) y sus cuadrados \(\widehat{u}_{i}^2\)
Alternativa 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}\]
- 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\]
- 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:
- Estimar por MCO la regresión auxiliar: \[\widehat{u}_i^2=\alpha_0+\alpha_1 Z_{1i}+\dots+\alpha_p Z_{pi}+v_i\]
- 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 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.09354417k <- 3
n <- nrow(wage1)
estadistico <- n * R2
estadistico## [1] 49.20423v.critic <- qchisq(0.95,k)
v.critic## [1] 7.814728valorP <- pchisq(q=estadistico, df=k, lower.tail=FALSE)
cbind(estadistico, valorP)## estadistico valorP
## [1,] 49.20423 1.180183e-10Para implementar esta alternativa se puede llamar utilizar directamente la función bptest de la librería lmtest de R.
REVISAR PORQUE NO DA IGUAL!!!!!!!!!!!
#install.packages("lmtest")
library(lmtest)
bptest(modelo)##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 43.096, df = 3, p-value = 2.349e-09El 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.
##Alternativa1** Reescalando.
0.15.3.2 Contraste de White
- Estimar el modelo por MCO y obtener los residuos \(\widehat{u}_i\)
- 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*}\]
- 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 se puede utilizar la misma función, 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-10La 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.96001valorP <- pchisq(q=estadistico, df=9, lower.tail=FALSE)
cbind(estadistico, valorP)## estadistico valorP
## [1,] 63.96001 2.297686e-10Algunos comentarios:
Los contrastes de White y B-P son contrastes generales (no requiere 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
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.
####Contraste de Goldfeld-Quandt Es un contraste válido con muestras pequeñas. 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\)
Pasos del contraste:
Identificar la variable \(X_j\) y ordenar los datos por \(X_j\) (varianza creciente)
Dividir los datos en 3 submuestras, omitiendo \(m\) observaciones centrales
Estimar por MCO el modelo original en las submuestras inicial y final (cada una tiene \((n-m)/2\) observaciones)
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 se significación \(\alpha\) si \(\lambda>c\) donde \(c\) cumple: \[Prob(F_{\frac{n-m}{2}-k-1,\frac{n-m}{2}-k-1}>c)=\alpha\]
La varianza al final debería ser mayor que al principio de la muestra
El cuadrado de los residuos está asociado con \[\sigma^{2}\Rightarrow SCR_{2}>SCR_{1}\]
Observaciónes y comentarios: * Si se sospecha que la varianza depende inversamente de los valores de \(X_j\Rightarrow\) ordenar 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 \end{itemize}
- Excluir `pocas’ observaciones puede provocar que las sub-muestras sean muy parecidas y el test pierde potencia.
Se recomienda 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
Para implementar en R se sugiere utilizar la función gqtest () del paquete lmtest para realizar la prueba Goldfeld-Quandt para determinar si hay heterocedasticidad.
Esta función utiliza la siguiente sintaxis:
gqtest (modelo, orden por, datos, fracción)
dónde:
modelo: El modelo de regresión lineal creado por el comando lm (). order.by: las variables predictoras del modelo. data: el nombre del conjunto de datos. fracción *: el número de observaciones centrales que se eliminarán del conjunto de datos.