#1.- DATOS

library(wooldridge)
## Warning: package 'wooldridge' was built under R version 4.2.3

En esta librería encontraremos todas las bases de datos que contiene el libro. Se deberá utilizar el mismo nombre con el que son referenciadas en el texto.

A modo de ejemplo, trabajaremos a lo largo de todo el docuemnto, con la base wage1 que contiene datos de la Encuesta de población de 1976, recopilados por Henry Farber cuando él y Wooldridge fueron colegas en el MIT (Instituto de Tecnología de Massachusetts) en 1988.

#2.- MANIPULACIÓN DE DATOS

library(wooldridge)
base <- data("wage1")
base <- wage1
names(base)
##  [1] "wage"     "educ"     "exper"    "tenure"   "nonwhite" "female"  
##  [7] "married"  "numdep"   "smsa"     "northcen" "south"    "west"    
## [13] "construc" "ndurman"  "trcommpu" "trade"    "services" "profserv"
## [19] "profocc"  "clerocc"  "servocc"  "lwage"    "expersq"  "tenursq"

Las variables que utilizaremos son las siguientes:

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

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

names(base)[4] <- "antigüedad"
names(base)
##  [1] "wage"       "educ"       "exper"      "antigüedad" "nonwhite"  
##  [6] "female"     "married"    "numdep"     "smsa"       "northcen"  
## [11] "south"      "west"       "construc"   "ndurman"    "trcommpu"  
## [16] "trade"      "services"   "profserv"   "profocc"    "clerocc"   
## [21] "servocc"    "lwage"      "expersq"    "tenursq"

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

base1 <- base[,c(1: 7, 22, 23, 24)]
names(base1)
##  [1] "wage"       "educ"       "exper"      "antigüedad" "nonwhite"  
##  [6] "female"     "married"    "lwage"      "expersq"    "tenursq"
head(base1, n=10)
##     wage educ exper antigüedad 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
## 6   8.75   16     9          8        0      0       1 2.169054      81      64
## 7  11.25   18    15          7        0      0       0 2.420368     225      49
## 8   5.00   12     5          3        0      1       0 1.609438      25       9
## 9   3.60   12    26          4        0      1       0 1.280934     676      16
## 10 18.18   17    22         21        0      0       1 2.900322     484     441
summary(base1)
##       wage             educ           exper         antigüedad    
##  Min.   : 0.530   Min.   : 0.00   Min.   : 1.00   Min.   : 0.000  
##  1st Qu.: 3.330   1st Qu.:12.00   1st Qu.: 5.00   1st Qu.: 0.000  
##  Median : 4.650   Median :12.00   Median :13.50   Median : 2.000  
##  Mean   : 5.896   Mean   :12.56   Mean   :17.02   Mean   : 5.105  
##  3rd Qu.: 6.880   3rd Qu.:14.00   3rd Qu.:26.00   3rd Qu.: 7.000  
##  Max.   :24.980   Max.   :18.00   Max.   :51.00   Max.   :44.000  
##     nonwhite          female          married           lwage        
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :-0.6349  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.: 1.2030  
##  Median :0.0000   Median :0.0000   Median :1.0000   Median : 1.5369  
##  Mean   :0.1027   Mean   :0.4791   Mean   :0.6084   Mean   : 1.6233  
##  3rd Qu.:0.0000   3rd Qu.:1.0000   3rd Qu.:1.0000   3rd Qu.: 1.9286  
##  Max.   :1.0000   Max.   :1.0000   Max.   :1.0000   Max.   : 3.2181  
##     expersq          tenursq       
##  Min.   :   1.0   Min.   :   0.00  
##  1st Qu.:  25.0   1st Qu.:   0.00  
##  Median : 182.5   Median :   4.00  
##  Mean   : 473.4   Mean   :  78.15  
##  3rd Qu.: 676.0   3rd Qu.:  49.00  
##  Max.   :2601.0   Max.   :1936.00
datos1 <- base[,c("wage","educ","exper", "antigüedad" )] 

head(datos1, n = 10)
##     wage educ exper antigüedad
## 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
## 6   8.75   16     9          8
## 7  11.25   18    15          7
## 8   5.00   12     5          3
## 9   3.60   12    26          4
## 10 18.18   17    22         21
summary(datos1)
##       wage             educ           exper         antigüedad    
##  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
library(dplyr) 
## Warning: package 'dplyr' was built under R version 4.2.3
## 
## 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", "antigüedad") 

head(datos2, n=10)
##     wage educ exper antigüedad
## 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
## 6   8.75   16     9          8
## 7  11.25   18    15          7
## 8   5.00   12     5          3
## 9   3.60   12    26          4
## 10 18.18   17    22         21
summary(datos2)
##       wage             educ           exper         antigüedad    
##  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

#3.- CORRELACIONES

plot(datos2)

cor(datos2)
##                 wage        educ      exper  antigüedad
## 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
## antigüedad 0.3468896 -0.05617257  0.4992914  1.00000000

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

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

cor.test(datos2$wage, datos2$educ)
## 
##  Pearson's product-moment correlation
## 
## data:  datos2$wage and datos2$educ
## t = 10.167, df = 524, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3319280 0.4749166
## sample estimates:
##       cor 
## 0.4059033

La correlación entre el salario y la educación es estadísticamente significativa.

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.3
ggplot(datos2, aes(x=educ, y=wage))+
  geom_point(color = "blue", size = 1) + theme_light()+
  ggtitle("Relación entre salaro y educación")

#4.- MODELO DE REGRESIÓN

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 la educación de las personas.

mod1 <- lm(wage ~educ, data=datos2)
summary(mod1)
## 
## 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

El salario independiente de la educación es -0,90 y el coeficiente resulta significativo.

Se estima que un incremento de la educación de un año adicional, con lo demás constante, provoque un incremento del salario de 0.54. La educación es estadísticamente significativa a 5%.

Las variaciones de la educación tan sólo permite explicar el 16.48% de las variaciones del salario (las predicciones no tienen suficiente fiabilidad).

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

Se observa que las pendientes del modelo son conjuntamente significativas (P=0.0000<0.05).

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 objertos del data frame, utilizando el operador $ entre el objeto y el elemento.

mod1$coefficients
## (Intercept)        educ 
##  -0.9048516   0.5413593
ahat <- coef(summary(mod1))[1,1]
ahat
## [1] -0.9048516
bhat <- coef(summary(mod1))[2,1]
bhat
## [1] 0.5413593
predicciones <- predict(mod1)
head(predicciones, n=5)
##        1        2        3        4        5 
## 5.050100 5.591459 5.050100 3.426022 5.591459
resid <- mod1$residuals
head(resid, n = 5)
##          1          2          3          4          5 
## -1.9501003 -2.3514594 -2.0501002  2.5739776 -0.2914593

#5.- DIAGRAMA DE DISPERSIÓN Y RECTA DE REGRESIÓN

ggplot(datos2, aes(x=educ, y=wage))+
  geom_point(color = "blue", size=1)+
  geom_smooth(method = 'lm', formula = y~x, se = TRUE, col = 'red')+
  theme_light()+
  ggtitle("Relación entre salario y educación")

library(plotly)
## Warning: package 'plotly' was built under R version 4.2.3
## 
## 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
ggplotly(data=datos2, x=~educ, y= ~wage)
ggplot(datos2, aes(x=educ, y=wage))+
  geom_point(color = "blue", size=1)+
  geom_smooth(method = 'loess', formula = y~x, se = TRUE, col = 'red')+
  theme_light()+
  ggtitle("Relación entre salario y educación")

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

#6.- RESIDUOS

ggplot(datos2, aes(x = educ, y = wage)) +
  geom_smooth(method = "lm", se = TRUE, 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'

# Para generar el gráfico de los residuos tal cual
plot(resid)

mean(resid)
## [1] -1.19498e-16

La media es próxima a cero (exogeneidad estricta y el estimador es insesgado).

library(e1071)  # para la función skewness
## Warning: package 'e1071' was built under R version 4.2.3
par(mfrow = c(1, 2))  # divide el área de gráficos en 2 columnas

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

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

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

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

#7.- PREDICCIÓN

mean(datos2$educ)
## [1] 12.56274
mod1$coefficients
## (Intercept)        educ 
##  -0.9048516   0.5413593
sal_pred <- 0.90485 + 0.0541359 * 12.56
sal_pred
## [1] 1.584797
nuevo <- data.frame(educ = 12.56)
future_y <- predict(object = mod1, newdata = nuevo, interval = "prediction", level = 0.95)
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
IC_inf_esp_y
## [1] 5.60524
IC_sup_esp_y
## [1] 6.184001
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("Predicción de y al 95%")

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

plot(IC_esp_y)

library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.2.3
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
grid.arrange(IC_esp_y, IC_y, ncol = 2, nrow = 1)

library(jtools)
## Warning: package 'jtools' was built under R version 4.2.3
library(ggstance)
## Warning: package 'ggstance' was built under R version 4.2.3
## 
## Attaching package: 'ggstance'
## The following objects are masked from 'package:ggplot2':
## 
##     geom_errorbarh, GeomErrorbarh
library(broom.mixed)
## Warning: package 'broom.mixed' was built under R version 4.2.3
## 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) 

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

#5.- INFO DE LA SESIÓN

sesion_info <- devtools::session_info()
dplyr::select(
  tibble::as_tibble(sesion_info$packages),
  c(package, loadedversion, source)
)
## # A tibble: 95 × 3
##    package     loadedversion source        
##    <chr>       <chr>         <chr>         
##  1 backports   1.4.1         CRAN (R 4.2.0)
##  2 broom       1.0.5         CRAN (R 4.2.3)
##  3 broom.mixed 0.2.9.4       CRAN (R 4.2.3)
##  4 bslib       0.4.2         CRAN (R 4.2.2)
##  5 cachem      1.0.6         CRAN (R 4.2.2)
##  6 callr       3.7.3         CRAN (R 4.2.3)
##  7 class       7.3-20        CRAN (R 4.2.2)
##  8 cli         3.6.0         CRAN (R 4.2.2)
##  9 codetools   0.2-18        CRAN (R 4.2.2)
## 10 colorspace  2.1-0         CRAN (R 4.2.3)
## # ℹ 85 more rows