Índice
1. Introducción 2. Acceder y conocer la base 3. Estimar el modelo 4. Predicción puntual e intervalos de confianza

Introducción

Este es un documento de trabajo, contiene ejemplos para aplicar modelos de regresión lineal simple con los datos del libro de texto del curso de Econometría 1, Wooldridge 2010-Introducción a la Econometría (4ta edución). 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’) }

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 ella encontrará todas las bases de datos que contiene el libro, utilizando 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 eran colegas en el MIT en 1988.

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

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:

Cambio el nombre de la variable que está en la posición 4

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

Selecciono las variables que me interesan por la posición de las columnas en el data frame. Indico que quiero 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)]

Miro 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 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          2

Otra forma mas sencilla de seleccionar la base 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: #install.packages(“dplyr”)

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
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

Lo puedo combinar con un filtro por filas con la función filter del paquete dplyr. Selecciono solo a las personas casadas.

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

#View(datos3)
  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.00000000

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

Estimar el modelo

Por el momento vamos a trabajar solo con el salario (wage) y educación (educ) en un MRLS.

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")

Estimo un MRLS con el método Mínimos cuadrados ordinarios (en adelante MCO)

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

Interpretar la estimación de los parámetros.

Incluyo la recta de regresión al gráfico de dispersió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, me paro en cada observación y veo los valores. Necesita la función plotly

library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
plot_ly(data=datos2, x=~educ, y=~wage)
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Miro los componentes del modelo y los asigno 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$residuals

Para consultar la estimación de un coeficiente de regresión se utilizan paréntesis rectos.

#View(coeficientes)
beta.estimate <- coeficientes["educ"]
print(beta.estimate)
##      educ 
## 0.5413593

Agrego a mis datos una nueva variable con los valores estimados según el MRLS.

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

Calculo los residuos

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

Comparo los residuos calculados con los que me dio el modelo

datos2$resid<- resid
datos2$verif<-datos2$residuos-datos2$resid

Agrego al gráfico de dispersión los ygorro en rojo y muestro los residuos

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()
## `geom_smooth()` using formula 'y ~ x'

Inspección gráfica de los resididuos

plot(datos2$residuos)

Inspecciono 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")

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

mean(datos2$educ)
## [1] 12.56274
sal_prom<-0.90485+0.054136*12.56
nuevo <- data.frame(educ=10)
future_y<-predict(object=mod1, newdata=nuevo, interval="prediction", level=0.95)

Genera el intervalo de confianza para el valor esperado de \(y\). 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$upr

Genera el data frame completo, agregando los extremos de los intervalos de confianza

nuevos_datos <- cbind(datos2, future_y, IC_inf_esp_y, IC_sup_esp_y)
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)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
grid.arrange(IC_esp_y_95, IC_y, ncol=2,nrow=1)

Veo el efecto de los niveles de confianza.

library(ggplot2)

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

library(jtools)
library(ggstance)
## 
## Attaching package: 'ggstance'
## The following objects are masked from 'package:ggplot2':
## 
##     geom_errorbarh, GeomErrorbarh
library(broom.mixed)
## Registered S3 methods overwritten by 'broom':
##   method            from  
##   tidy.glht         jtools
##   tidy.summary.glht jtools
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)